text:
% Scrip: mcgrid.m
% Calculates Magnitude shift map for two specific time periods
% Uses view_mcgrid to plot the results
%
% J. Woessner
% last update: 22.01.04

report_this_filefun(mfilename('fullpath'));

global no1 bo1 inb1 inb2 valeg valeg2 CO valm1

valeg = 1;
valm1 = min(a.Magnitude);
prf = NaN;
if sel == 'in'
    % Set the grid parameter
    %Initial values
    dx = 1;
    dy = 1;
    ni = 150;
    Nmin = 150;
    ra = 50;
    fMaxRadius = 5;
    fSplitTime = 2000.4;

    % cut catalog at mainshock time:
%     l = a.Date > maepi(1,3);
%     a = a.subset(l);

    % cat at selecte magnitude threshold
    l = a.Magnitude < valm1;
    a(l,:) = [];
    newt2 = a;

    ho2=true;
    timeplot
    ho2=false;

    %The definitions in the following line were present in the initial bvalgrid.m file.
    %stan2 = NaN; stan = NaN; prf = NaN; av = NaN;

    % make the interface
    % creates a dialog box to input grid parameters
    %
    figure_w_normalized_uicontrolunits(...
        'Name','Grid Input Parameter',...
        'NumberTitle','off', ...
        'MenuBar','none', ...
        'NextPlot','new', ...
        'units','points',...
        'Visible','off', ...
        'Position',[ wex+200 wey-200 650 250]);
    axis off
%     labelList2=[' Automatic Mcomp (max curvature) | Fixed Mc (Mc = Mmin) | Automatic Mcomp (90% probability) | Automatic Mcomp (95% probability) | Best (?) combination (Mc95 - Mc90 - max curvature) | Constant Mc'];
%     labelPos = [0.2 0.8  0.6  0.08];
%     hndl2=uicontrol(...
%         'Style','popup',...
%         'Position',labelPos,...
%         'Units','normalized',...
%         'String',labelList2,...
%         'Callback','inb2 =get(hndl2,''Value''); ');
%
%     set(hndl2,'value',5);

    % creates a dialog box to input grid parameters
    %

    oldfig_button = uicontrol('BackGroundColor',[.60 .92 .84], ...
        'Style','checkbox','string','Plot in Current Figure',...
        'Position',[.78 .7 .20 .08],...
        'Units','normalized');

    set(oldfig_button,'value',1);

    freq_field=uicontrol('Style','edit',...
        'Position',[.30 .60 .12 .08],...
        'Units','normalized','String',num2str(ni),...
        'Callback','ni=str2double(get(freq_field,''String'')); set(freq_field,''String'',num2str(ni));set(tgl2,''value'',0); set(tgl1,''value'',1)');

    freq_field0=uicontrol('Style','edit',...
        'Position',[.30 .50 .12 .08],...
        'Units','normalized','String',num2str(ra),...
        'Callback','ra=str2double(get(freq_field0,''String'')); set(freq_field0,''String'',num2str(ra)) ; set(tgl2,''value'',1); set(tgl1,''value'',0)');

    freq_field2=uicontrol('Style','edit',...
        'Position',[.30 .40 .12 .08],...
        'Units','normalized','String',num2str(dx),...
        'Callback','dx=str2double(get(freq_field2,''String'')); set(freq_field2,''String'',num2str(dx));');

    freq_field3=uicontrol('Style','edit',...
        'Position',[.30 .30 .12 .080],...
        'Units','normalized','String',num2str(dy),...
        'Callback','dy=str2double(get(freq_field3,''String'')); set(freq_field3,''String'',num2str(dy));');

    freq_field4=uicontrol('Style','edit',...
        'Position',[.6 .30 .12 .080],...
        'Units','normalized','String',num2str(fSplitTime),...
        'Callback','fSplitTime=str2double(get(freq_field4,''String'')); set(freq_field4,''String'',num2str(fSplitTime));');

    freq_field7=uicontrol('Style','edit',...
        'Position',[.30 .20 .12 .080],...
        'Units','normalized','String',num2str(Nmin),...
        'Callback','Nmin=str2double(get(freq_field7,''String'')); set(freq_field7,''String'',num2str(Nmin));');

    freq_field8=uicontrol('Style','edit',...
        'Position',[.6 .60 .12 .080],...
        'Units','normalized','String',num2str(fMaxRadius),...
        'Callback','fMaxRadius=str2double(get(freq_field8,''String'')); set(freq_field8,''String'',num2str(fMaxRadius));');

    tgl1 = uicontrol('Style','radiobutton',...
        'string','Number of Events:',...
        'Position',[.05 .60 .2 .0800], 'Callback','set(tgl2,''value'',0)',...
        'Units','normalized');

    set(tgl1,'value',0);

    tgl2 =  uicontrol('Style','radiobutton',...
        'string','OR: Constant Radius',...
        'Position',[.05 .50 .2 .080], 'Callback','set(tgl1,''value'',0)',...
        'Units','normalized');
    set(tgl2,'value',1);

    create_grid =  uicontrol('Style','radiobutton',...
        'string','Calculate a new grid', 'Callback','set(load_grid,''value'',0), set(prev_grid,''value'',0)','Position',[.78 .55 .2 .080],...
        'Units','normalized');

    set(create_grid,'value',1);

    prev_grid =  uicontrol('Style','radiobutton',...
        'string','Reuse the previous grid', 'Callback','set(load_grid,''value'',0),set(create_grid,''value'',0)','Position',[.78 .45 .2 .080],...
        'Units','normalized');

    load_grid =  uicontrol('Style','radiobutton',...
        'string','Load a previously saved grid', 'Callback','set(prev_grid,''value'',0),set(create_grid,''value'',0)','Position',[.78 .35 .2 .080],...
        'Units','normalized');

    save_grid =  uicontrol('Style','checkbox',...
        'string','Save selected grid to file',...
        'Position',[.78 .22 .2 .080],...
        'Units','normalized');

    close_button=uicontrol('Style','Pushbutton',...
        'Position',[.60 .05 .15 .12 ],...
        'Units','normalized','Callback','close;done','String','Cancel');

    go_button1=uicontrol('Style','Pushbutton',...
        'Position',[.20 .05 .15 .12 ],...
        'Units','normalized',...
        'Callback','tgl1 =get(tgl1,''Value'');tgl2 =get(tgl2,''Value'');prev_grid = get(prev_grid,''Value'');create_grid = get(create_grid,''Value''); load_grid = get(load_grid,''Value''); save_grid = get(save_grid,''Value''); oldfig_button = get(oldfig_button,''Value''); close,sel =''ca'', mc_grid',...
        'String','Go');

    text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0.10 0.98 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.l ,...
        'FontWeight','bold',...
        'String','Please choose an Mc estimation option   ');
    txt3 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0.30 0.75 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.l ,...
        'FontWeight','bold',...
        'String',' Grid Parameter');
    txt5 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[-0.1 0.4 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Spacing in x (dx) in deg:');

    txt6 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[-0.1 0.3 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Spacing in y (dy) in deg:');

    txt7 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[-0.1 0.18 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Min. No. of events > Mc:');

    txt9 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0.42 0.28 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.s ,...
        'FontWeight','bold',...
        'String','Time window:');

    txt11 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0.42 0.62 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.s ,...
        'FontWeight','bold',...
        'String','Max. Radius /[km]:');

    set(gcf,'visible','on');
    watchoff

end   % if nargin ==0

% get the grid-size interactively and
% calculate the b-value in the grid by sorting
% thge seimicity and selectiong the ni neighbors
% to each grid point

if sel == 'ca'

    %In the following line, the program .m is called, which creates a rectangular grid from which then selects,
    %on the basis of the vector ll, the points within the selected poligon.

    % get new grid if needed
    if load_grid == 1
        [file1,path1] = uigetfile(['*.mat'],'previously saved grid');
        if length(path1) > 1
            think
            load([path1 file1])
        end
        plot(newgri(:,1),newgri(:,2),'k+')
    elseif load_grid ==0  &&  prev_grid == 0
        selgp
        if length(gx) < 2  ||  length(gy) < 2
            errordlg('Selection too small! (Dx and Dy are in degreees! ');
            return
        end
    elseif prev_grid == 1
        plot(newgri(:,1),newgri(:,2),'k+')
    end
    %   end

    gll = ll;

    if save_grid == 1

        sFile = ['*.mat'];
        sPath = pwd;
        [file1,path1] = uiputfile([sFile], 'Grid File Name?');
        sSaveFile = [path1 file1];
        save(sSaveFile, 'newgri', 'dx', 'dy', 'gx', 'gy', 'xvect', 'yvect', 'tmpgri', 'll', 'fSplitTime', 'Nmin', 'ra', 'a');
    end

    %   selgp
    itotal = length(newgri(:,1));
%     if length(gx) < 4 | length(gy) < 4
%         errordlg('Selection too small! (Dx and Dy are in degreees! ');
%         return
%     end

    zmap_message_center.set_info(' ','Running... ');think
    %  make grid, calculate start- endtime etc.  ...
    %
    t0b = min(a.Date)  ;
    n = a.Count;
    teb = a(n,3) ;
    tdiff = round((teb - t0b)*365/par1);
    loc = zeros(3, length(gx)*length(gy));

    % loop over  all points
    %
    i1 = 0.;
    mRes =[];
    mResult = [];
    allcount = 0.;
    wai = waitbar(0,' Please Wait ...  ');
    set(wai,'NumberTitle','off','Name','Rate change grid - percent done');
    drawnow

    % loop over all points
    for i= 1:length(newgri(:,1))
        i/length(newgri(:,1));
        % Grid node point
        x = newgri(i,1);y = newgri(i,2);
        allcount = allcount + 1.;

        % calculate distance from center point and sort with distance
        l = sqrt(((a.Longitude-x)*cos(pi/180*y)*111).^2 + ((a.Latitude-y)*111).^2) ;
        [s,is] = sort(l);
        mCat = a(is(:,1),:) ;       % re-orders matrix to agree row-wise
        % Use Radius to determine grid node catalogs
        l3 = l <= ra;
        mCat = a.subset(l3);      % new data per grid point (mCat) is sorted in distance

        %         % Select earthquakes in non-overlapping rectangles
        %         vSel = (a.Longitude >= (newgri(i,1)-dx/2)) & (a.Longitude < (newgri(i,1)+dx/2)) &...
        %             (a.Latitude >= (newgri(i,2)-dy/2)) & (a.Latitude < (newgri(i,2)+dy/2));
        %         % Select earthquakes in overlapping rectangles
        %         vSel = (a.Longitude >= (newgri(i,1)-dx)) & (a.Longitude < (newgri(i,1)+dx)) &...
        %             (a.Latitude >= (newgri(i,2)-dy)) & (a.Latitude < (newgri(i,2)+dy));
        %         mCat = a.subset(vSel);

        % Initialize
        fMinTime = min(mCat(:,3));
        fMaxTime = max(mCat(:,3));

        % Select data from 2 time periods
        vSelT = mCat(:,3) < fSplitTime;
        mCat1 = mCat(vSelT,:);
        mCat2 = mCat(~vSelT,:);
        % Length of catalog (resolution)
        nNumevents1 = length(mCat1(:,1));
        nNumevents2 = length(mCat2(:,1));

        % Minimum bin: essential for computing difference in FMD
        fMinBin = roundn(min(mCat(:,6)),-1);
        fMaxBin = roundn(max(mCat(:,6)),-1);

        if (length(mCat1(:,1)) >=Nmin & length(mCat2(:,1)) >= Nmin)
            % Compute change in FMD normalized by time period
            % Time periods
            fPeriod1 = max(mCat1(:,3))-min(mCat1(:,3));
            fPeriod2 = max(mCat2(:,3))-min(mCat2(:,3));
            [vFMD1, vBin1] = hist(mCat1(:,6),fMinBin:0.1:fMaxBin);
            [vFMD2, vBin2] = hist(mCat2(:,6),fMinBin:0.1:fMaxBin);
            fChFMD = max(cumsum(abs(vFMD2./fPeriod2-vFMD1./fPeriod1)));
            % Calculate shift
            [fMshift, fProbability, fAICc, mProblikelihood, bH] = calc_loglikelihood_dM2(mCat1, mCat2);
            % Check for validity of model using KS-Test result to produce validated result of magnitude shift
            if (bH == 1 | fMshift == 0)
                fMshift_valid = NaN;
            else
                fMshift_valid = fMshift;
            end
            % Calculate Mc
            [mResult1, fMls1, fMc1, fMu1, fSigma1, mDatPredBest1,...
                    vPredBest1, fBvalue1, fAvalue1, bH1] = calc_McEMR_kstest(mCat1, 0.1);
            [mResult2, fMls2, fMc2, fMu2, fSigma2, mDatPredBest2,...
                    vPredBest2, fBvalue2, fAvalue2, bH2] = calc_McEMR_kstest(mCat2, 0.1);
            % Mc change
            fdMc = fMc2-fMc1;
            if (bH1 ==0  &&  bH2 == 0)
                fdMc_val = fdMc;
            else
                fdMc_val = NaN;
            end
            % Calculate Utsu-Test
            [dA, fProbEqual, fProbDifferent] = calc_Utsu(fBvalue1, fBvalue2, nNumevents1, nNumevents2);
            % Create result matrix
            mResult = [mResult; i fMshift fProbability fAICc bH nNumevents1 nNumevents2 fMshift_valid fChFMD fMc1 fMc2 fdMc fdMc_val dA fProbEqual fProbDifferent bH1 fBvalue1 fAvalue1 bH2 fBvalue2 fAvalue2];
        else
            % Create result matrix
            mResult = [mResult; i NaN NaN NaN NaN nNumevents1 nNumevents2 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN];
        end
        waitbar(allcount/itotal)
    end  % for newgr

    % Save the data to rcval_grid.mat
    save  mcgrid.mat mResult gx gy dx dy par1 tdiff t0b teb a main faults mainfault coastline yvect xvect tmpgri ll bo1 newgri gll ra maepi fSplitTime
    disp('Saving data to mcgrid.mat in current directory')

    close(wai)
    watchoff

    normlap2=ones(length(tmpgri(:,1)),1)*nan;
    % Magnitude shift
    normlap2(ll)= mResult(:,2);
    mMagShift = reshape(normlap2,length(yvect),length(xvect));

    % KS-Test-Value
    normlap2(ll)= mResult(:,5);
    mHkstest = reshape(normlap2,length(yvect),length(xvect));

    %%% Resolution parameters
    % Number of events per grid node in first time period
    normlap2(ll)= mResult(:,6);
    mNumevents1 = reshape(normlap2,length(yvect),length(xvect));

    % Number of events per grid node in first time period
    normlap2(ll)= mResult(:,7);
    mNumevents2 = reshape(normlap2,length(yvect),length(xvect));

    % Validated magnitude shift
    normlap2(ll)= mResult(:,8);
    mMagShift_valid = reshape(normlap2,length(yvect),length(xvect));

    % Absolute FMD difference
    normlap2(ll)= mResult(:,9);
    mChFMD = reshape(normlap2,length(yvect),length(xvect));

    % Mc Periode 1
    normlap2(ll)= mResult(:,10);
    mMc1 = reshape(normlap2,length(yvect),length(xvect));

    % Mc Periode 2
    normlap2(ll)= mResult(:,11);
    mMc2 = reshape(normlap2,length(yvect),length(xvect));

    % Mc change
    normlap2(ll)= mResult(:,12);
    mdMc = reshape(normlap2,length(yvect),length(xvect));

    % Mc change validated by KS-Test
    normlap2(ll)= mResult(:,13);
    mdMc_val = reshape(normlap2,length(yvect),length(xvect));

    % Difference of AIC for Utsu-Test
    normlap2(ll)= mResult(:,14);
    mdAIC_Utsu = reshape(normlap2,length(yvect),length(xvect));

    % Probability of stationarity (favoring stationarity)
    normlap2(ll)= mResult(:,15);
    mStationary1_Utsu = reshape(normlap2,length(yvect),length(xvect));

    % Probability of stationarity (favoring non-stationarity)
    normlap2(ll)= mResult(:,16);
    mStationary2_Utsu = reshape(normlap2,length(yvect),length(xvect));

    % Data to plot first map
    re3 = mMagShift;
    lab1 = 'Magnitude shift';

    % View the map
    view_mcgrid

end   % if sel = na

% Load exist b-grid
if sel == 'lo'
    [file1,path1] = uigetfile(['*.mat'],'b-value gridfile');
    if length(path1) > 1
        think
        load([path1 file1])

        normlap2=ones(length(tmpgri(:,1)),1)*nan;

        % Magnitude shift
        normlap2(ll)= mResult(:,2);
        mMagShift = reshape(normlap2,length(yvect),length(xvect));

        % KS-Test-Value
        normlap2(ll)= mResult(:,5);
        mHkstest = reshape(normlap2,length(yvect),length(xvect));

        %%% Resolution parameters
        % Number of events per grid node in first time period
        normlap2(ll)= mResult(:,6);
        mNumevents1 = reshape(normlap2,length(yvect),length(xvect));

        % Number of events per grid node in first time period
        normlap2(ll)= mResult(:,7);
        mNumevents2 = reshape(normlap2,length(yvect),length(xvect));

        try
            % Validated magnitude shift
            normlap2(ll)= mResult(:,8);
            mMagShift_valid = reshape(normlap2,length(yvect),length(xvect));

            % Absolute FMD difference
            normlap2(ll)= mResult(:,9);
            mChFMD = reshape(normlap2,length(yvect),length(xvect));
            % Mc Periode 1
            normlap2(ll)= mResult(:,10);
            mMc1 = reshape(normlap2,length(yvect),length(xvect));

            % Mc Periode 2
            normlap2(ll)= mResult(:,11);
            mMc2 = reshape(normlap2,length(yvect),length(xvect));

            % Mc change
            normlap2(ll)= mResult(:,12);
            mdMc = reshape(normlap2,length(yvect),length(xvect));

            % Mc change
            normlap2(ll)= mResult(:,13);
            mdMc_val = reshape(normlap2,length(yvect),length(xvect));

            % Difference of AIC for Utsu-Test
            normlap2(ll)= mResult(:,14);
            mdAIC_Utsu = reshape(normlap2,length(yvect),length(xvect));

            % Probability of stationarity (favoring stationarity)
            normlap2(ll)= mResult(:,15);
            mStationary1_Utsu = reshape(normlap2,length(yvect),length(xvect));

            % Probability of stationarity (favoring non-stationarity)
            normlap2(ll)= mResult(:,16);
            mStationary2_Utsu = reshape(normlap2,length(yvect),length(xvect));

        catch
            disp('Validated magnitude shift map not available');
        end
        % Initial map set to relative rate change
        re3 = mMagShift;
        lab1 = 'Magnitude shift';

        old = re3;
        % Plot
        view_mcgrid;
    else
        return
    end
end

author: CelsoReyes


text:
function em = errormatrix(fet,clu)

cluIx = unique(clu);
nClu = length(cluIx);
nDim = size(fet,2);

mu = zeros(nClu,nDim);
sigma = zeros(nDim,nDim,nClu);
p = zeros(1,nClu);
for ii=1:nClu
    f = fet(clu==cluIx(ii),:);
    % Mean fet vector
    mu(ii,:) = mean(f);
    % Covariance. We add 1e-5 to the diagonal to ensure that it will be
    % positive definite
    sigma(:,:,ii) = cov(f)+1e-5*eye(nDim,nDim);
    p(ii) = sum(clu==cluIx(ii));
end
obj = gmdistribution(mu,sigma,p);
%keyboard

em = zeros(nClu,nClu);

for ii=1:nClu
    p = posterior(obj,fet(clu==cluIx(ii),:));
    %keyboard
    em(ii,:) = mean(p);
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% iirfilt() -  (high|low|band)-pass filter data using an elliptic IIR filter
%              Multiple data channels and epochs are supported. Forward and
%              reverse filtering are used to avoid phase distortions.
% Usage:
%  >> [smoothdata] = iirfilt(data,srate,locutoff,hicutoff);
%  >> [smoothdata,filtwts] = iirfilt(data,srate,locutoff,hicutoff, ...
%                                    epochframes, trans_bw, revfilt,rp,rs)
% Inputs:
%   data        = (channels,frames*epochs) data to filter
%   srate       = data sampling rate (Hz)
%   locutoff    = low-edge frequency in pass band (Hz).  If 0, lowpass only.
%   hicutoff    = high-edge frequency in pass band (Hz). If 0, highpass only.
%   epochframes = frames per epoch (filter each epoch separately)
%                 {default|0: data is 1 epoch}
%   trans_bw    = width (in Hz) of transition interval between stop and
%                 pass bands {default: 1 Hz, or if pass band f < 5 Hz, f/3 Hz}.
%                 Enter 0 to use default.
%   revfilt     = [0|1] reverse filter (i.e. bandpass filter to notch
%                 or band-reject filter). {0}
%   rp          = ripple amplitude in dB in the pass band {default: 0.0025 dB,
%                 or if pass band f < 5 Hz, 0.01 dB}
%   rs          = ripple amplitude in dB in the stop band {default: 40 dB,
%                 or if pass band f < 5 Hz, 30 dB}
%   causal      = ['on'|'off'] use causal filter. Default 'off'.
%
% Note: Requires the MATLAB Signal Processing Toolbox.
%
% Authors: Maksym Pozdin (mpozdin.ece04@gtalumni.org, IOL/ONRC,2004),
%          with Arnaud Delorme and Scott Makeig (SCCN/INC/UCSD, La Jolla CA)
%
% See also: eegfilt(), eegfiltfft()

% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% Notes:
% Currently, with cutoff less then 5 Hz, both HPF and LPF use more 'relaxed'
% filter parameters to achieve best results.

function [smoothdata,filtwts] = iirfilt(data,srate,locutoff,hicutoff,epochframes, trans_bw, revfilt, rp, rs, causal)

if nargin<4
    fprintf('');
    help iirfilt
    return
end

data = double(data);

if exist('ellipord') ~= 2 | exist('ellip') ~= 2
   error('*** ellip() requires the Matlab Signal Processing toolbox. ***');
end

[chans frames] = size(data);
if chans > 1 & frames == 1,
    help iirfilt
    error('input data should be a row vector.');
end
nyq            = srate*0.5;  % Nyquist frequency
MINFREQ = 5;

if locutoff>0 & hicutoff > 0 & locutoff > hicutoff,
    error('locutoff > hicutoff ???');
end
if locutoff < 0 | hicutoff < 0,
   error('locutoff | hicutoff < 0 ???');
end

if locutoff>nyq,
    error('locutoff cannot be > srate/2');
end

if hicutoff>=nyq
   hicutoff = 0;
end
if nargin < 6
    trans_bw = 0;
end;

if trans_bw==0

    if locutoff < MINFREQ & locutoff > 0
        trans_bw=locutoff/3;
    elseif hicutoff < MINFREQ & hicutoff > 0
        trans_bw=hicutoff/3;
        %rp=0.01;
        %rs=30;
    else
        % Default transition bandwidth is 1 Hz
        trans_bw=1;
    end;
end

if nargin<7
   revfilt = 0;
end
if nargin<8 || isempty(rp)
   % Ripple in the passband
   rp=0.0025;
end
if nargin<9 || isempty(rs)
   % Ripple in the stopband
   rs=40;
end
if nargin<10
   causal = 'off';
end

if nargin<5
        epochframes = 0;
end
if epochframes ==0,
    epochframes = frames;    % default
end
epochs = fix(frames/epochframes);
if epochs*epochframes ~= frames,
    error('epochframes does not divide frames.\n');
end

if locutoff > 0 & hicutoff ==0 & revfilt==0              % highpass filter

    ws=(locutoff-trans_bw)/nyq;
    wp=(locutoff)/nyq;
    [N,wn] = ellipord(wp,ws,rp,rs);
    fprintf('HPF has been designed.\n');
    fprintf('HPF has cutoff of %1.1f Hz, transition bandwidth of %1.1f Hz and its order is %1.1f\n',locutoff, trans_bw,N);
    [b,a]=ellip(N,rp,rs,wn,'high');

elseif hicutoff > 0 & locutoff==0 & revfilt==0                %  lowpass filter

    ws=(hicutoff+trans_bw)/nyq;
    wp=(hicutoff)/nyq;
    [N,wn] = ellipord(wp,ws,rp,rs);
    fprintf('LPF has been designed.\n');
    fprintf('LPF has cutoff of %1.1f Hz, transition bandwidth of %1.1f Hz and its order is %1.1f\n',hicutoff, trans_bw,N);
    [b,a]=ellip(N,rp,rs,wn);

elseif hicutoff > 0 & locutoff> 0 & revfilt==0               %  bandpass filter

    %first part LPF
    ws=(hicutoff+trans_bw)/nyq;
    wp=(hicutoff)/nyq;
    [N,wn] = ellipord(wp,ws,rp,rs);
    fprintf('BPF has been designed.\n');
    fprintf('LPF has cutoff of %1.1f Hz, transition bandwidth of %1.1f Hz and its order is %1.1f\n',hicutoff, trans_bw,N);
    %[zl, pl, kl]=ellip(N,rp,rs,wn);
    %[lpf_sos,lpf_g] = zp2sos(zl,pl, kl);
    [bl,al]=ellip(N,rp,rs,wn);

    %second part HPF
    ws=(locutoff-trans_bw)/nyq;
    wp=(locutoff)/nyq;
    [N,wn] = ellipord(wp,ws,rp,rs);
    fprintf('HPF has cutoff of %1.1f Hz, transition bandwidth of %1.1f Hz and its order is %1.1f\n',locutoff, trans_bw,N);
    %[zh, ph, kh]=ellip(N,rp,rs,wn);
    %[hpf_sos,hpf_g] = zp2sos(zh,ph, kh);
    [bh,ah]=ellip(N,rp,rs,wn,'high');
    %help fvtool
    %b=conv(bh,bl);a=conv(ah,al);
    b.bl=bl;b.bh=bh; a.al=al;a.ah=ah;

elseif hicutoff > 0 & locutoff> 0 & revfilt==1               %  bandreject filter

    %first part LPF
    ws=(locutoff+trans_bw)/nyq;
    wp=(locutoff)/nyq;
    [N,wn] = ellipord(wp,ws,rp,rs);
    fprintf('BRF has been designed.\n');
    fprintf('LPF has cutoff of %1.1f Hz, transition bandwidth of %1.1f Hz and its order is %1.1f\n',locutoff, trans_bw,N);
    [bl,al]=ellip(N,rp,rs,wn);

    %second part HPF
    ws=(hicutoff-trans_bw)/nyq;
    wp=(hicutoff)/nyq;
    [N,wn] = ellipord(wp,ws,rp,rs);
    fprintf('HPF has cutoff of %1.1f Hz, transition bandwidth of %1.1f Hz and its order is %1.1f\n',hicutoff, trans_bw,N);
    [bh,ah]=ellip(N,rp,rs,wn,'high');
    b.bl=bl;b.bh=bh; a.al=al;a.ah=ah;

else
    error('You must provide a non-0 low or high cut-off frequency');
end

smoothdata = zeros(chans,frames);
for e = 1:epochs                % filter each epoch, channel
    for c=1:chans
        if isstruct(a) & isstruct(b) & revfilt==0            %BPF - filter with LPF and HPF in series

            if strcmpi(causal, 'on')
                smoothdata1(c,(e-1)*epochframes+1:e*epochframes) = filter(b.bl,a.al,data(c,(e-1)*epochframes+1:e*epochframes));
                smoothdata( c,(e-1)*epochframes+1:e*epochframes) = filter(b.bh,a.ah,smoothdata1(c,(e-1)*epochframes+1:e*epochframes));
            else
                smoothdata1(c,(e-1)*epochframes+1:e*epochframes) = filtfilt(b.bl,a.al,data(c,(e-1)*epochframes+1:e*epochframes));
                smoothdata( c,(e-1)*epochframes+1:e*epochframes) = filtfilt(b.bh,a.ah,smoothdata1(c,(e-1)*epochframes+1:e*epochframes));
            end;

        elseif isstruct(a) & isstruct(b) & revfilt==1         %BRF - filter with LPF and HPF in parallel
            if strcmpi(causal, 'on')
                smoothdata1(c,(e-1)*epochframes+1:e*epochframes) = filter(b.bl,a.al,data(c,(e-1)*epochframes+1:e*epochframes));
                smoothdata2(c,(e-1)*epochframes+1:e*epochframes) = filter(b.bh,a.ah,data(c,(e-1)*epochframes+1:e*epochframes));
            else
                smoothdata1(c,(e-1)*epochframes+1:e*epochframes) = filtfilt(b.bl,a.al,data(c,(e-1)*epochframes+1:e*epochframes));
                smoothdata2(c,(e-1)*epochframes+1:e*epochframes) = filtfilt(b.bh,a.ah,data(c,(e-1)*epochframes+1:e*epochframes));
            end;
            smoothdata=smoothdata1+smoothdata2;               %combine final results
        else
            if strcmpi(causal, 'on')
                smoothdata(c,(e-1)*epochframes+1:e*epochframes) = filter(b,a,data(c,(e-1)*epochframes+1:e*epochframes));
            else
                smoothdata(c,(e-1)*epochframes+1:e*epochframes) = filtfilt(b,a,data(c,(e-1)*epochframes+1:e*epochframes));
            end;
        end;
        if epochs == 1
            if rem(c,20) ~= 0
                fprintf('.');
            else
                fprintf('%d',c);
            end
        end
    end
end
fprintf('\n');

%% To test my filter: draws Frequency response
% N=10001; n=1:N;
% xx=zeros(size(n));
% xx(ceil(length(xx)/2))=1;
% if isstruct(a) & isstruct(b)
%     yy1=filtfilt(b.bl,a.al,xx);
%     yy=filtfilt(b.bh, a.ah,yy1);
% else
%     yy=filtfilt(b,a,xx);
% end;
% YY=fft(yy,100000);
% ff=-nyq:2*nyq/length(YY):nyq-2*nyq/length(YY);
% figure;plot(ff,abs(fftshift(YY)));grid on;title('Frequency response');

author: PatternRecognition

repo: OpenBMI


text:
function [XXs,YYs] = Delequalsamples(XXs,YYs)

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    zerosindex        = YYs == 0;
    XXs(zerosindex,:) = [];
    YYs(zerosindex)   = [];
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function tests = test_spm_ncTpdf
% Unit Tests for spm_ncTpdf
%__________________________________________________________________________
% Copyright (C) 2018 Wellcome Trust Centre for Neuroimaging

% $Id: test_spm_ncTpdf.m 7260 2018-02-19 10:55:53Z guillaume $

tests = functiontests(localfunctions);

function test_spm_ncTpdf_1(testCase)
exp = spm_Tpdf(0,1);
act = spm_ncTpdf(0,1,0);
tol = 1e-12;
testCase.verifyEqual(act, exp,'AbsTol',tol);

exp = spm_Tpdf(-2:0.5:2,2);
act = spm_ncTpdf(-2:0.5:2,2,0);
tol = 1e-12;
testCase.verifyEqual(act, exp,'AbsTol',tol);

exp = spm_Tpdf(0,1:4);
act = spm_ncTpdf(0,1:4,0);
tol = 1e-12;
testCase.verifyEqual(act, exp,'AbsTol',tol);

exp = spm_Tpdf(-2:2,1:5);
act = spm_ncTpdf(-2:2,1:5,0);
tol = 1e-12;
testCase.verifyEqual(act, exp,'AbsTol',tol);

function test_spm_ncTpdf_2(testCase)
exp = 0.193064705260108; % nctpdf(0,1,1)
act = spm_ncTpdf(0,1,1);
tol = 1e-12;
testCase.verifyEqual(act, exp,'AbsTol',tol);

exp = 2.072225240523640e-04; % nctpdf(1,10,5);
act = spm_ncTpdf(1,10,5);
tol = 1e-12;
testCase.verifyEqual(act, exp,'AbsTol',tol);

exp = [... % nctpdf(-1:3,6:10,-2:2);
   0.240403340723731
   0.233509118657019
   0.227607580145303
   0.223756811547761
   0.219731001980809]';
act = spm_ncTpdf(-1:3,6:10,-2:2);
tol = 1e-12;
testCase.verifyEqual(act, exp,'AbsTol',tol);

author: spm


text:
%%***********************************************************************
%%  Read in a problem in SDPA sparse format.
%%
%%  [blk,At,C,b] = read_sdpa(fname)
%%
%%  Input: fname = name of the file containing SDP data in
%%                 SDPA foramt.
%%  Important: the data is assumed to contain only
%%             semidefinite and linear blocks.
%%
%%***********************************************************************
%% SDPNAL+
%% Copyright (c) 2014 by
%% Liuqin Yang, Defeng Sun, and Kim-Chuan Toh
%%***********************************************************************
   function [blk,At,C,b] = read_sdpa(fname)

%%
%%  Open the file for input
%%
   compressed = 0;
   if exist(fname)
      fid = fopen(fname,'r');
   elseif exist([fname,'.Z']);
      compressed = 1;
      unix(['uncompress ',fname,'.Z']);
      fid = fopen(fname,'r');
   elseif exist([fname,'.gz']);
      compressed = 2;
      unix(['gunzip ',fname,'.gz']);
      fid = fopen(fname,'r');
   else
      fprintf(['*** Problem "',fname,'" not found, please specify the correct path or problem name. \n']);
      blk = []; At = []; C = []; b = [];
      return;
   end
%%
%%  Clean up special characters and comments from the file
%%
   [datavec,count] = fscanf(fid,'%c');
   linefeeds = findstr(datavec,char(10));
   comment_chars = '*"=';
   cumidx = [];
   for i=1:length(comment_chars)
      idx = findstr(datavec,comment_chars(i));
      cumidx = [cumidx,idx];
   end
   for j=length(cumidx):-1:1
      if (cumidx(j)==1) || (strcmp(datavec(cumidx(j)-1),char(10)))
         datavec(cumidx(j):linefeeds(min(find(cumidx(j)<linefeeds))))='';
      else
         datavec(cumidx(j):linefeeds(min(find(cumidx(j)<linefeeds)))-1)='';
      end
   end
   special_chars = ',{}()';
   cumidx=[];
   for i=1:length(special_chars)
      idx = findstr(datavec,special_chars(i));
      cumidx = [cumidx,idx];
   end
   datavec(cumidx) = blanks(length(cumidx));
   clear linefeeds;
%%
%% Close the file
%%
   fclose('all');
   if compressed==1; unix(['compress ',fname]); end;
   if compressed==2; unix(['gzip ',fname]); end;
%%
%%  Next, read in basic problem size parameters.
%%
   datavec = sscanf(datavec,'%f');
   if size(datavec,1) < size(datavec,2); datavec = datavec'; end;
   m = datavec(1);
   numblk  = datavec(2);
   blksize = datavec(2+[1:numblk]);
   if size(blksize,1) > size(blksize,2); blksize = blksize'; end
%%
%% Get input  b.
%%
   idxstrt = 2+numblk;
   b = datavec(idxstrt+[1:m]);
   idxstrt = idxstrt+m;
   b = -b;
%%
%% Construct blk
%%
   deblksize = 80;
   spblkidxtmp = find( (blksize>1) & (blksize < deblksize) );
   spblkidxtmp = sort(spblkidxtmp);
   deblkidx = find( (blksize<=1) | (blksize >= deblksize) );
   denumblk = length(deblkidx);
   linblkidx = zeros(1,denumblk);
   for p = 1:denumblk
      n = blksize(deblkidx(p));
      if (n > 1);
         blk{p,1} = 's'; blk{p,2} = n;
         n2 = n*(n+1)/2;
         At{p,1} = sparse(n2,m);
         C{p,1} = sparse(n,n);
      else
         linblkidx(p) = p;
         blk{p,1} = 'l'; blk{p,2} = abs(n);
         At{p,1} = sparse(abs(n),m);
         C{p,1} = sparse(abs(n),1);
      end
   end
   if ~isempty(spblkidxtmp)
      maxnumblk = 200;
      spnumblk = ceil(length(spblkidxtmp)/maxnumblk);
      for q = 1:spnumblk
         if (q < spnumblk)
            spblkidxall{q} = spblkidxtmp([(q-1)*maxnumblk+1: q*maxnumblk]);
         else
            spblkidxall{q} = spblkidxtmp([(q-1)*maxnumblk+1: length(spblkidxtmp)]);
         end
         tmp = blksize(spblkidxall{q});
         blk{denumblk+q,1} = 's';
         blk{denumblk+q,2} = tmp;
         n2 = sum(tmp.*(tmp+1))/2;
         At{denumblk+q,1} = sparse(n2,m);
         C{denumblk+q,1} = sparse(sum(tmp),sum(tmp));
      end
   else
      spnumblk = 0;
   end
   linblkidx(denumblk+[1:spnumblk]) = zeros(1,spnumblk);
%%
%% Construct single blocks of A,C
%%
   len = length(datavec);
   Y = reshape(datavec(idxstrt+1:len),5,(len-idxstrt)/5)';
   clear datavec;
   Y = sortrows(Y,[1 2]);
   matidx = [0; find(diff(Y(:,1)) ~= 0); size(Y,1)];
%%
   for k = 1:length(matidx)-1
      idx = [matidx(k)+1 : matidx(k+1)];
      Ytmp  = Y(idx,1:5);
      matno = Ytmp(1,1);
      Ytmp2 = Ytmp(:,2);
      for p = 1:denumblk
         n  = blksize(deblkidx(p));
         idx = find(Ytmp2 == deblkidx(p));
         ii = Ytmp(idx,3); jj = Ytmp(idx,4); vv =Ytmp(idx,5);
         len = length(idx);
         if (n > 1)
            idxtmp = find(ii > jj);
            if ~isempty(idxtmp);
               tmp = jj(idxtmp);
               jj(idxtmp) = ii(idxtmp); ii(idxtmp) = tmp;
            end
            tmp = -sparse(ii,jj,vv,n,n);
            tmp = tmp + triu(tmp,1)';
         else
            tmp = -sparse(ii,ones(len,1),vv,abs(n),1);
         end
         if (matno == 0)
            C{p,1} = tmp;
         else
            if (n > 1)
               At{p,1}(:,matno) = svec(blk(p,:),tmp,1);
            else
               At{p,1}(:,matno) = tmp;
            end
         end
      end
   end
%%
%% Construct big sparse block of A,C
%%
if (spnumblk > 0)
   Y1 = Y(:,1);
   diffY1 = find(diff([-1; Y1; inf]));
   for kk = 1:length(diffY1)-1
      idx = [diffY1(kk) : diffY1(kk+1)-1];
      matno = Y1(diffY1(kk));
      Ytmp  = Y(idx,1:5);
      Ytmp2 = Ytmp(:,2);
      maxYtmp2  = Ytmp2(length(Ytmp2));
      minYtmp2  = Ytmp2(1);
      diffYtmp2 = Ytmp2(find(diff([-1; Ytmp2])));
      for q = 1:spnumblk
         spblkidx    = spblkidxall{q};
         maxspblkidx = spblkidx(length(spblkidx));
         minspblkidx = spblkidx(1);
         count = 0;
         if (minYtmp2 <= maxspblkidx) & (maxYtmp2 >= minspblkidx)
            tmpblksize = blksize(spblkidx);
            n = sum(tmpblksize);
            cumspblksize = [0 cumsum(tmpblksize)];
            n2 = sum(tmpblksize.*(tmpblksize+1))/2;
            idx = zeros(n2,1); offset = zeros(n2,1);
            for t = [1:length(diffYtmp2)]
  	       p = find(spblkidx == diffYtmp2(t));
               if ~isempty(p)
	          idxtmp = find(Ytmp2 == spblkidx(p));
                  len = length(idxtmp);
     	          idx(count+[1:len]) = idxtmp;
                  offset(count+[1:len]) = cumspblksize(p);
                  count = count + len;
               end
            end
            idx = idx(1:count); offset = offset(1:count);
            ii = Ytmp(idx,3)+offset; jj = Ytmp(idx,4)+offset; vv = Ytmp(idx,5);
            idxtmp = find(ii > jj);
            if ~isempty(idxtmp);
               tmp = jj(idxtmp);
               jj(idxtmp) = ii(idxtmp); ii(idxtmp) = tmp;
            end
	    idxeq = find(ii==jj);
            tmp = spconvert([ii jj -vv; jj ii -vv; n n 0]) ...
                  + spconvert([ii(idxeq) jj(idxeq) vv(idxeq); n n 0]);
            if (matno == 0)
               C{denumblk+q,1} = tmp;
            else
               At{denumblk+q,1}(:,matno) = svec(blk(denumblk+q,:),tmp,1);
            end
         end
      end
   end
end
%%
%% put all linear blocks together as a single linear block
%%
   idx = find(linblkidx);
   if (length(idx) > 1)
      sdpidx = find(linblkidx==0);
      blktmp = 0; Atmp = []; Ctmp = [];
      for k = 1:length(idx)
          tmp = linblkidx(idx(k));
          blktmp = blktmp+blk{tmp,2};
          Atmp = [Atmp; At{tmp}];
          Ctmp = [Ctmp; C{tmp}];
      end
      At = At(sdpidx); C = C(sdpidx); blk = blk(sdpidx,:);
      len = length(sdpidx);
      blk(2:len+1,:) = blk;
      blk{1,1} = 'l'; blk{1,2} = blktmp;
      At(2:len+1,1) = At; C(2:len+1,1) = C;
      At{1,1} = Atmp; C{1,1} = Ctmp;
   end
%%******************************************************************

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function cost = computecost(f,c,Q,x,p)

if isempty(f)
    f = 0;
end
cost = f+c'*x+x'*Q*x;
if ~isequal(p.K.m,0)
    top = size(p.F_struc,1)-sum(p.K.m.^2)+1;
    for i = 1:length(p.K.m)
        X = p.F_struc(top:top + p.K.m(i)^2-1,:)*[1;x];
        X = reshape(X,p.K.m(i),p.K.m(i));
        cost = cost + p.K.maxdetgain(i)*sum(real(log(real(eig(X)))));
        top = top + p.K.m(i)^2;
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function options = rbfOptions

% RBFOPTIONS Default options for RBF network.

% MLTOOLS

options.outFunc = 'linear';
options.activeFunc = 'gaussian';
options.hiddenDim = 20;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [vert,conn,tria,tnum] = smooth2(varargin)
%SMOOTH2 "hill-climbing" mesh-smoothing for two-dimensional,
%2-simplex triangulations.
%   [VERT,EDGE,TRIA,TNUM] = SMOOTH2(VERT,EDGE,TRIA,TNUM) re-
%   turns a "smoothed" triangulation {VERT,TRIA}, incorpora-
%   ting "optimised" vertex coordinates and mesh topology.
%
%   VERT is a V-by-2 array of XY coordinates in the triangu-
%   lation, EDGE is an array of constrained edges, TRIA is a
%   T-by-3 array of triangles, and TNUM is a T-by-1 array of
%   part indices. Each row of TRIA and EDGE define an eleme-
%   nt. VERT(TRIA(II,1),:), VERT(TRIA(II,2),:) and VERT(TRIA
%   (II,3),:) are the coordinates of the II-TH triangle. The
%   edges in EDGE are defined in a similar manner. NUM is an
%   array of part indexing, such that TNUM(II) is the index
%   of the part in which the II-TH triangle resides.
%
%   [VERT,EDGE,TRIA,TNUM] = SMOOTH2(... ,OPTS) passes an ad-
%   ditional options structure OPTS, containing user-defined
%   parameters, including:
%
% - OPTS.VTOL = {+1.0E-02} -- relative vertex movement tole-
%   rance, smoothing is converged when (VNEW-VERT) <= VTOL *
%   VLEN, where VLEN is a local effective length-scale.
%
% - OPTS.ITER = {+32} -- max. number of smoothing iterations
%
% - OPTS.DISP = {+ 4} -- smoothing verbosity. Set to INF for
%   quiet execution.
%
%   See also REFINE2, TRICOST, TRIDEMO

%   This routine is loosely based on the DISTMESH algorithm,
%   employing a "spring-based" analogy to redistribute mesh
%   vertices. Such an approach is described in: P.O. Persson
%   and Gilbert Strang. "A simple mesh generator in MATLAB."
%   SIAM review 46(2) 2004, pp: 329--345. Details of the al-
%   gorithm used here are somewhat different, with an alter-
%   ative spring-based update employed, in addition to hill-
%   climbing element quality guarantees, and vertex density
%   controls.

%-----------------------------------------------------------
%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 21/07/2017
%-----------------------------------------------------------

    vert = []; conn = []; tria = [] ;
    tnum = [];
    opts = []; hfun = []; harg = {} ;

%---------------------------------------------- extract args
    if (nargin>=+1), vert = varargin{1}; end
    if (nargin>=+2), conn = varargin{2}; end
    if (nargin>=+3), tria = varargin{3}; end
    if (nargin>=+4), tnum = varargin{4}; end
    if (nargin>=+5), opts = varargin{5}; end

   [opts] = makeopt(opts) ;

%---------------------------------------------- default CONN
    if (isempty(conn))

       [edge] = tricon2(tria);

        ebnd = edge(:,4) < +1;              %-- use bnd edge
        conn = edge(ebnd,1:2);

    end

%---------------------------------------------- default TNUM
    if (isempty(tnum))
        tnum = ones(size(tria, 1), 1) ;
    end

%---------------------------------------------- basic checks
    if ( ~isnumeric(vert) || ...
         ~isnumeric(conn) || ...
         ~isnumeric(tria) || ...
         ~isnumeric(tnum) || ...
         ~isstruct (opts) )
        error('smooth2:incorrectInputClass' , ...
            'Incorrect input class.') ;
    end

%---------------------------------------------- basic checks
    if (ndims(vert) ~= +2 || ...
        ndims(conn) ~= +2 || ...
        ndims(tria) ~= +2 || ...
        ndims(tnum) ~= +2 )
        error('smooth2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

    if (size(vert,2)~= +2 || ...
        size(conn,2)~= +2 || ...
        size(tria,2)~= +3 || ...
        size(tnum,2)~= +1 || ...
        size(tria,1)~= size(tnum,1) )
        error('smooth2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

    nvrt = size(vert,1) ;

%---------------------------------------------- basic checks
    if (min(min(conn(:,1:2))) < +1 || ...
            max(max(conn(:,1:2))) > nvrt )
        error('smooth2:invalidInputs', ...
            'Invalid EDGE input array.') ;
    end

    if (min(min(tria(:,1:3))) < +1 || ...
            max(max(tria(:,1:3))) > nvrt )
        error('smooth2:invalidInputs', ...
            'Invalid TRIA input array.') ;
    end

%---------------------------------------------- output title
    if (~isinf(opts.disp))
        fprintf(1,'\n') ;
        fprintf(1,' Smooth triangulation...\n') ;
        fprintf(1,'\n') ;
        fprintf(1,[...
' -------------------------------------------------------\n', ...
'      |ITER.|          |MOVE(X)|          |DTRI(X)|     \n', ...
' -------------------------------------------------------\n', ...
             ] ) ;
    end

%---------------------------------------------- polygon bnds
    node = vert; PSLG = conn; part = {};

    pmax = max(tnum(:));
    for ppos = +1 : pmax

        tsel = tnum == ppos ;
        tcur = tria(tsel,:) ;

       [ecur,tcur] ...
           = tricon2 (tcur) ;

        ebnd = ecur(:,4)==0 ;

        same = setset2( ...
            PSLG,ecur(ebnd,1:2));

        part{ppos} = find(same) ;

    end

%---------------------------------------------- inflate bbox
    vmin = min(vert,[],1);
    vmax = max(vert,[],1);

    vdel = vmax - 1.*vmin;
    vmin = vmin - .5*vdel;
    vmax = vmax + .5*vdel;

    vbox = [
        vmin(1), vmin(2)
        vmax(1), vmin(2)
        vmax(1), vmax(2)
        vmin(1), vmax(2)
           ] ;
    vert = [vert ; vbox] ;

%---------------------------------------------- DO MESH ITER
    tnow =  tic ;

    tcpu = struct('full',0.,'dtri',0., ...
        'tcon',0.,'iter',0.,'undo',0., ...
        'keep',0.) ;

    for iter = +1 : opts.iter

    %------------------------------------------ inflate adj.
        ttic = tic ;

       [edge,tria] = tricon2(tria,conn) ;

        tcpu.tcon = ...
            tcpu.tcon + toc(ttic) ;

    %------------------------------------------ compute scr.
        oscr = triscr2(vert,tria) ;

    %------------------------------------------ vert. iter's
        ttic =  tic ;

        nvrt = size(vert,1);
        nedg = size(edge,1);

        IMAT = sparse( ...
           edge(:,1),(1:nedg)',+1,nvrt,nedg) ;
        JMAT = sparse( ...
           edge(:,2),(1:nedg)',+1,nvrt,nedg) ;

        EMAT = IMAT + JMAT ;

        vdeg = sum(EMAT,2) ;                %-- vertex |deg|
        free = (vdeg == 0) ;

        vold = vert ;
        for isub = +1 : max(+2,min(+8,iter))

        %-- compute HFUN at vert/midpoints
            hvrt = evalhfn(vert, ...
                edge,EMAT,hfun,harg) ;

            hmid = hvrt(edge(:,1),:) ...
                 + hvrt(edge(:,2),:) ;
            hmid = hmid * +.5 ;

        %-- calc. relative edge extensions
            evec = vert(edge(:,2),:) ...
                 - vert(edge(:,1),:) ;
            elen = ...
                sqrt(sum(evec.^2,2)) ;

            scal = +1.0-elen./hmid ;
            scal = min (+1.0, scal);
            scal = max (-1.0, scal);

        %-- projected points from each end
            ipos = vert(edge(:,1),:) ...
                 -.67*[scal,scal].*evec;
            jpos = vert(edge(:,2),:) ...
                 +.67*[scal,scal].*evec;

           %scal = ...                      %-- nlin. weight
           %   max(abs(scal).^.5,eps^.75);
            scal = ...
               max(abs(scal).^ 1,eps^.75);

        %-- sum contributions edge-to-vert
            vnew = ...
            IMAT*([scal,scal] .* ipos) ...
          + JMAT*([scal,scal] .* jpos) ;

            vsum = max(EMAT*scal,eps^.75);

            vnew = vnew ./ [vsum,vsum] ;

        %-- fixed points. edge projection?
            vnew(conn(:),1:2) = ...
            vert(conn(:),1:2) ;

            vnew(vdeg==0,1:2) = ...
            vert(vdeg==0,1:2) ;

        %-- reset for the next local iter.
            vert = vnew ;

        end

        tcpu.iter = ...
            tcpu.iter + toc(ttic) ;

    %------------------------------------------ hill-climber
        ttic = tic ;

    %-- unwind vert. upadte if score lower
        nscr = ones(size(tria,1),1);
        btri = true(size(tria,1),1);

        umax = + 8 ;

        for undo = +1 : umax

            nscr(btri) = triscr2( ...
                vert,tria(btri,:)) ;

        %-- TRUE if tria needs "unwinding"
            smin = +.70 ;
            smax = +.90 ;
            sdel = .025 ;

            stol = smin+iter*sdel;
            stol = min (smax,stol) ;

            btri = nscr <= stol ...
                 & nscr <  oscr ;

            if (~any(btri)), break; end

        %-- relax toward old vert. coord's
            ivrt = ...
               unique(tria(btri,1:3));

            bvrt = ...
               false(size(vert,1),1) ;
            bvrt(ivrt) = true;

            if (undo ~= umax)
                bnew =  +.75 ^ undo ;
                bold =  +1.0 - bnew ;
            else
                bnew =  +0.0 ;
                bold =  +1.0 - bnew ;
            end

            vert(bvrt,:) = ...
                bold * vold(bvrt,:) ...
              + bnew * vert(bvrt,:) ;

            btri = any( ...
               bvrt(tria(:,1:3)),2) ;

        end

        oscr = nscr ;

        tcpu.undo = ...
            tcpu.undo + toc(ttic) ;

    %------------------------------------- test convergence!
        ttic = tic ;

        vdel = ...
            sum((vert-vold).^2,2) ;

        evec = vert(edge(:,2),:) ...
             - vert(edge(:,1),:) ;
        elen = ...
            sqrt(sum(evec.^2,2)) ;

        hvrt = evalhfn(vert, ...
            edge,EMAT,hfun,harg) ;

        hmid = hvrt(edge(:,1),:) ...
             + hvrt(edge(:,2),:) ;
        hmid = hmid * 0.5 ;
        scal = elen./hmid ;

        emid = vert(edge(:,1),:) ...
             + vert(edge(:,2),:) ;
        emid = emid * 0.5 ;

    %------------------------------------- |deg|-based prune
        keep = false(size(vert,1),1);
        keep(vdeg>+4) = true ;

        keep(conn(:)) = true ;
        keep(free(:)) = true ;

    %------------------------------------- 'density' control
        lmax = +5. / +4.  ;
        lmin = +1. / lmax ;

        less = scal<=lmin ;
    	more = scal>=lmax ;

        vbnd = false(size(vert,1),1);
        vbnd(conn(:,1)) = true ;
        vbnd(conn(:,2)) = true ;

        ebad = vbnd(edge(:,1)) ...     %-- not at boundaries
             | vbnd(edge(:,2)) ;

        less(ebad(:)) = false;
        more(ebad(:)) = false;

    %------------------------------------- force as disjoint
        lidx = find (less) ;

        for lpos = 1 : length(lidx)
            epos = lidx(lpos,1);
            inod = edge(epos,1);
            jnod = edge(epos,2);
        %--------------------------------- if still disjoint
            if (keep(inod) && ...
                keep(jnod) )

            keep(inod) = false ;
            keep(jnod) = false ;

            else

            less(epos) = false ;

            end
        end

        ebad = ...
           keep(edge(less,1)) ...
         & keep(edge(less,2)) ;

        more(ebad(:)) = false ;

    %------------------------------------- reindex vert/tria
        redo = ...
            zeros(size(vert,1),1);
        itop = ...
            length(find(keep));
        iend = ...
            length(find(less));

        redo(keep) = (1:itop)';

        redo(edge(less,1)) = ...        %-- to new midpoints
            (itop+1 : itop+iend)';
        redo(edge(less,2)) = ...
            (itop+1 : itop+iend)';

        vnew =[vert(keep,:) ;
               emid(less,:) ;
              ] ;
        tnew = redo(tria(:,1:3)) ;

        ttmp = sort(tnew,2) ;           %-- filter collapsed
        okay = all( ...
            diff(ttmp,1,2)~=0,2) ;
        okay = ...
            okay & ttmp(:,1) > 0 ;
        tnew = tnew(okay,:) ;

    %------------------------------------- quality preserver
        nscr = ...
            triscr2  (vnew,tnew) ;

        stol = +0.80 ;

        tbad = nscr < stol ...
             & nscr < oscr(okay) ;

        vbad = ...
            false(size(vnew,1),1);
        vbad(tnew(tbad,:)) = true;

    %------------------------------------- filter edge merge
        lidx = find (less) ;

        ebad = ...
        vbad(redo(edge(lidx,1))) | ...
        vbad(redo(edge(lidx,2))) ;

        less(lidx(ebad)) = false ;

        keep(edge(...
        lidx(ebad),1:2)) =  true ;

    %------------------------------------- reindex vert/conn
        redo = ...
            zeros(size(vert,1),1);
        itop = ...
            length(find(keep));
        iend = ...
            length(find(less));

        redo(keep) = (1:itop)';

        redo(edge(less,1)) = ...
            (itop+1 : itop+iend)';
        redo(edge(less,2)) = ...
            (itop+1 : itop+iend)';

        vert =[vert(keep,:);
               emid(less,:);
               emid(more,:);
              ] ;
        conn = redo(conn(:,1:2)) ;

        tcpu.keep = ...
            tcpu.keep + toc(ttic) ;

    %------------------------------------- build current CDT
        ttic = tic ;

       [vert,conn,tria,tnum] = ...
            deltri2 (vert, ...
                conn,node,PSLG,part) ;

        tcpu.dtri = ...
            tcpu.dtri + toc(ttic) ;

    %------------------------------------- dump-out progess!
        vdel = vdel./(hvrt.*hvrt) ;
        move = vdel > opts.vtol^2 ;
        nmov = ...
            length(find(move));

        ntri = size(tria,1) ;

        if (mod(iter,opts.disp)==+0)
            fprintf(+1, ...
            '%11i %18i %18i\n', ...
            [iter,nmov,ntri]) ;
        end

    %------------------------------------- loop convergence!
        if (nmov == +0), break; end

    end

    tria = tria( :,1:3);

%----------------------------------------- prune unused vert
    keep = false(size(vert,1),1);
    keep(tria(:)) = true ;
    keep(conn(:)) = true ;

    redo = zeros(size(vert,1),1);
    redo(keep) = ...
        (+1:length(find(keep)))';

    conn = redo(conn);
    tria = redo(tria);

    vert = vert(keep,:);

    tcpu.full = ...
        tcpu.full + toc(tnow) ;

    if (opts.dbug)
%----------------------------------------- print debug timer
        fprintf(1,'\n') ;
        fprintf(1,' Mesh smoothing timer...\n');
        fprintf(1,'\n') ;
        fprintf(1, ...
        ' FULL: %f \n', tcpu.full);
        fprintf(1, ...
        ' DTRI: %f \n', tcpu.dtri);
        fprintf(1, ...
        ' TCON: %f \n', tcpu.tcon);
        fprintf(1, ...
        ' ITER: %f \n', tcpu.iter);
        fprintf(1, ...
        ' UNDO: %f \n', tcpu.undo);
        fprintf(1, ...
        ' KEEP: %f \n', tcpu.keep);
        fprintf(1,'\n') ;
    end

    if (~isinf(opts.disp)), fprintf(1,'\n'); end

end

function [hvrt] = evalhfn(vert,edge,EMAT,hfun,harg)
%EVALHFN eval. the spacing-fun. at mesh vertices.

    if (~isempty (hfun))
        if (isnumeric(hfun))
            hvrt = hfun * ...
              ones(size(vert,1),1) ;
        else
            hvrt = feval( ...
                hfun,vert,harg{:}) ;
        end
    else

%-- no HFUN - HVRT is mean edge-len. at vertices!
        evec = vert(edge(:,2),:) - ...
               vert(edge(:,1),:) ;
        elen = sqrt(sum(evec.^2,2)) ;

        hvrt = (EMAT*elen) ...
            ./ max(sum(EMAT,2),eps) ;

        free = true(size(vert,1),1) ;
        free(edge(:,1)) = false;
        free(edge(:,2)) = false;

        hvrt(free)      =  +inf;

    end

end

function [opts] = makeopt(opts)
%MAKEOPT setup the options structure for SMOOTH2.

    if (~isfield(opts,'iter'))
        opts.iter = +32;
    else
    if (~isnumeric(opts.iter))
        error('smooth2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.iter)~= +1)
        error('smooth2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    if (opts.iter <= +0)
        error('smooth2:invalidOptionValues', ...
            'Invalid OPT.ITER selection.') ;
    end
    end

    if (~isfield(opts,'disp'))
        opts.disp = + 4;
    else
    if (~isnumeric(opts.disp))
        error('smooth2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.disp)~= +1)
        error('smooth2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    if (opts.disp <= +0)
        error('smooth2:invalidOptionValues', ...
            'Invalid OPT.DISP selection.') ;
    end
    end

    if (~isfield(opts,'vtol'))
        opts.vtol = +1.0E-02;
    else
    if (~isnumeric(opts.vtol))
        error('smooth2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.vtol)~= +1)
        error('smooth2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    if (opts.vtol <= 0.)
        error('smooth2:invalidOptionValues', ...
            'Invalid OPT.VTOL selection.') ;
    end
    end

    if (~isfield(opts,'dbug'))
        opts.dbug = false;
    else
    if (~islogical(opts.dbug))
        error('refine2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.dbug)~= +1)
        error('refine2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    end

end

author: CHLNDDEV


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%	SurfBox-MATLAB (c)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%	Yue Lu and Minh N. Do
%%
%%	Department of Electrical and Computer Engineering
%%	Coordinated Science Laboratory
%%	University of Illinois at Urbana-Champaign
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%	NXdiamondmapping.m
%%
%%	First created: 04-20-05
%%	Last modified: 04-13-06
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function m = NXdiamondmapping(N, beta)

%% Get the mapping kernel for the fan-shaped frequency response
%
%   See also ......
%

if rem(N, 2) ~= 1
    error('N must be an odd integer so that we have a zero-phase mapping function.');
end

%% Use the kaiser window to truncate and smooth the ideal sinc sequences

%% w = window(@kaiser, N, beta);
w = calculateKaiserwindow(N, beta);

ind = [1 : N]' - (N+1)/2;
w = w .* sinc(ind / 2);
w(N+1) = 0;

%% Get the mapping kernel
m = zeros(2*N-1);
[k1, k2] = meshgrid(-N+1:N-1, -N+1:N-1);
a = k1 + k2 + (N+1)/ 2;
a(a <= 0) = N + 1;
a(a > N) = N + 1;
a = w(a);
b = k1 - k2 + (N+1)/ 2;
b(b <= 0) = N + 1;
b(b > N) = N + 1;
b = w(b);

m = a .* b;
%% If you want to check the frquency response of this mapping kernel.
% mapp = m;
% dispfreqd2(CRISPfilter(mapp, [1;1]));

%%	This software is provided "as-is", without any express or implied
%%	warranty. In no event will the authors be held liable for any
%%	damages arising from the use of this software.

author: thomaskuestner


text:
% Interface for querying a database
%
% Author: Jonathan Karr
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 1/21/2011
classdef Database < handle
    properties (Abstract = true, SetAccess = protected)
        hostName
        schema
        userName
        password
    end

    methods
        function this = Database(varargin)
            switch nargin
                case 1
                    this.hostName = varargin{1}.hostName;
                    this.schema = varargin{1}.schema;
                    this.userName = varargin{1}.userName;
                    this.password = varargin{1}.password;
                case 4
                    this.hostName = varargin{1};
                    this.schema = varargin{2};
                    this.userName = varargin{3};
                    this.password = varargin{4};
                otherwise
                    throw(MException('Database:error', 'invalid options'));
            end
        end
    end

    methods (Abstract = true)
        result = query(this)
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function Ascale = cmultiscale(scale, A)
% Multiplies the 2D slices in a 3D array by individual scalars.
%
% function Ascale = cmultiscale(scale, A)
%
% Basically the same as multiscale but is compatible with structs with
% fields real and imag, which means this can be used for automatic
% differentiation with complex arrays in older Matlab versions.
%
% See also: manoptADhelp multiscale

% This file is part of Manopt: www.manopt.org.
% Original author: Xiaowen Jiang, Aug. 31, 2021.
% Contributors: Nicolas Boumal
% Change log:

    if isnumeric(A) && isnumeric(scale)

        Ascale = multiscale(scale, A);

    else

        A = tocstruct(A);
        scale = tocstruct(scale);

        assert(ndims(A.real) == 3, ...
           'cmultiscale is only well defined for 3D arrays.');

        [~, ~, N] = size(A.real);

        assert(numel(scale.real) == N, ...
           ['scale must be a vector whose length equals the third ' ...
            'dimension of A, that is, the number of 2D matrix slices ' ...
            'in the 3D array A. It can also be a struct with fields ' ...
            'real and imag with size as stated above.']);

        scale.real = reshape(scale.real, 1, 1, N);
        scale.imag = reshape(scale.imag, 1, 1, N);
        Ascale = cdottimes(A, scale);

    end

end

% Test code
% n = 3; m = 5; N = 17;
% A = randn(n, m, N);
% scale = randn(N, 1);
% Z = multiscale(scale, A) - cmultiscale(scale, A);
% norm(Z(:))
% A = randn(n, m, N) + 1i*randn(n, m, N);
% Z = multiscale(scale, A) - cmultiscale(scale, A);
% norm(Z(:))
% scale = randn(N, 1) + 1i*randn(N, 1);
% Z = multiscale(scale, A) - cmultiscale(scale, A);
% norm(Z(:))
% B.real = real(A); B.imag = imag(A);
% Z = cmultiscale(scale, B);
% Zr = real(multiscale(scale, A)) - Z.real;
% Zi = imag(multiscale(scale, A)) - Z.imag;
% norm(Zr(:))
% norm(Zi(:))
% scalebis = tocstruct(scale);
% Z = cmultiscale(scalebis, B);
% Zr = real(multiscale(scale, A)) - Z.real;
% Zi = imag(multiscale(scale, A)) - Z.imag;
% norm(Zr(:))
% norm(Zi(:))

author: NicolasBoumal

repo: manopt


text:
dbcnn.name = 'dbcnn' ;
dbcnn.opts = {...
  'type', 'bcnn', ...
  'modela', 'D:\dbcnn\data\models\imagenet-vgg-verydeep-16.mat', ...
  'layera', 30,...
  'modelb', 'D:\dbcnn\data\models\scnn.mat', ...
  'layerb', 18,...
  'shareWeight', false,...
  };

opts.setupNameList = {'dbcnn'};
opts.encoderList = {{dbcnn}};
opts.datasetList = {{'live', 1}};
%opts.learningRate = 1e-6;
opts.momentum = 0.9;
opts.batchSize = 8;
opts.numEpoch = 30;
opts.dataset = opts.datasetList{1,1}{1};

srcc = zeros(1,10);
plcc = zeros(1,10);

switch opts.dataset
    case 'live'
        subdirec = 'data\checkgpu\live-seed-01';
        modelpath = 'models\LIVE_models';
        datapath = 'data\checkgpu\live-seed-01';
        opts.learningRate = 1e-6;
    case 'csiq'
        subdirec = 'data\checkgpu\csiq-seed-01';
        modelpath = 'models\CSIQ_models';
        datapath = 'data\checkgpu\csiq-seed-01';
        opts.learningRate = 1e-6;
    case 'tid'
        subdirec = 'data\checkgpu\tid-seed-01';
        modelpath = 'models\TID_models';
        datapath = 'data\checkgpu\tid-seed-01';
        opts.learningRate = 1e-5;
    case 'mlive'
        subdirec = 'data\checkgpu\mlive-seed-01';
        modelpath = 'models\MLIVE_models';
        datapath = 'data\checkgpu\mlive-seed-01';
        opts.learningRate = 1e-5;
    case 'clive'
        subdirec = 'data\checkgpu\clive-seed-01';
        modelpath = 'models\Challen_models';
        datapath = 'data\checkgpu\clive-seed-01';
        opts.learningRate = 1e-5;
end

for split = 1:10
    mkdir(subdirec);
    rmdir(subdirec,'s');
    imdbpath = fullfile(subdirec,'imdb');
    mkdir(imdbpath);
    copyfile(fullfile(modelpath,num2str(split),'imdb-seed-1.mat'),...
        fullfile(imdbpath,'imdb-seed-1.mat'));
    [options, imdb] = run_experiments_bcnn_train(opts);
    deploy_each_dagnet(opts.numEpoch, datapath);
    options.ftpath = fullfile(subdirec,'fine-tuned-model');
    [srcc(1,split),plcc(1,split),index] = find_bestmodel(opts.numEpoch, imdb, options);
end

repo: DBCNN


function [ min_cost  ] = get_min_cost( compressResultsDir, qFile, lamda, epsilon_soft )
%GET_MIN_COST Summary of this function goes here
%   Detailed explanation goes here

    q_file_t = fopen(qFile);
    observe_value = fscanf(q_file_t, '%d');
    fclose(q_file_t);
    observe_value(find(observe_value > 100)) = 100;  % remove some unusual counts
    q_value = (1 - mapminmax(observe_value', 0, 1))';  % actually same as in CVPR_W 2013

    fileExt = '*.txt';
    files = dir(fullfile(compressResultsDir,fileExt));

    % get the union cost
    saved_ID = [];
    for i = 0:length(files)-1

        fileCnt = num2str(i);
        fileName = [compressResultsDir, fileCnt, '.txt'];

        file_t = fopen(fileName);
        point_ID = fscanf(file_t, '%d');
        fclose(file_t);

        saved_ID = union(saved_ID, point_ID);
    end

    % sum the q_value * if_compress
    min_cost = 0;
    for i = 1:1:length(saved_ID)
        min_cost = min_cost + q_value(saved_ID(i));
    end

%     % sum the soft part
%     for i = 1:length(files)
%         min_cost = min_cost + lamda*epsilon_soft(i,:);
%     end

end


text:
function angout = npi2pi(angin,units,approach)

    %NPI2PI Truncates angles into the -180 deg to 180 deg range
    %
    %
    %  ang = NPI2PI(angin) transforms input angles into the
    %  -180 to 180 degree range.
    %
    %  ang = NPI2PI(angin,'units') uses the units defined by the
    %  input string 'units'.  If omitted, default units of 'degrees'
    %  are assumed.
    %
    %  ang = NPI2PI(angin,'units','method') uses the method
    %  defined by the corresponding input string.  Valid methods are:
    %  'exact' for the exact transformation; 'inward' where all angles
    %  are shifted epsilon towards the origin before the -180 to 180 degree
    %  transformation; 'outward' where all angles are shifted epsilon
    %  away from the origin before the -180 to 180 degree transformation.
    %  If omitted, default method of 'exact' is assumed.
    %
    %  See also ZERO22PI, ANGLEDIM

    %  Copyright 1996-2000 Systems Planning and Analysis, Inc. and The MathWorks, Inc.
    %  Written by:  E. Byrns, E. Brown
    %   $Revision: 1399 $    $Date: 2006-08-11 11:19:27 +0200 (Fr, 11 Aug 2006) $

    %report_this_filefun(mfilename('fullpath'));

    if nargin == 0
        error('Incorrect number of arguments')
    elseif nargin == 1
        units = [];     approach = [];
    elseif nargin == 2
        approach = [];
    end

    %  Empty argument tests

    if isempty(approach);   approach = 'exact';    end
    if isempty(units);      units    = 'degrees';  end

    %  Convert inputs to radians

    angin = angledim(angin,units,'radians');

    %  Exact approach -- Eliminates the use of ATAN2 function.

    %  Approximate approach -- Some inconsistencies with ATAN2 function
    %  across platforms when an exact multiple of -pi is used.  Some platforms
    %  atan2(-1,0) = -pi, and for some platforms atan2(-1,0) = pi;

    switch lower(approach)
        case 'exact'

            %  Exact approach is not straightforward because of this mapping behavior:
            %       -3pi maps to -pi;   -2pi maps to -pi;  -pi maps to -pi;
            %        pi  maps to pi;     2pi maps to pi;   3pi maps to pi
            %  Note that the mapping point changes when the sign on pi changes.

            angout = pi*((abs(angin)/pi) - ...
                2*ceil(((abs(angin)/pi)-1)/2)) .* sign(angin);

        case 'inward'

            %  Move data epsilon towards (inward) the origin.  Eliminates any
            %  points which start identically on a multiple of pi.  Then
            %  use the atan2 function.

            epsilon = epsm('radians');
            angin   = angin*(1 - epsilon);
            angout  = atan2(sin(angin),cos(angin));

        case 'outward'

            %  Move data epsilon towards (away from) the origin.  Eliminates any
            %  points which start identically on a multiple of pi.  Then
            %  use the atan2 function.

            epsilon = epsm('radians');
            angin   = angin * (1 + epsilon);
            angout  = atan2(sin(angin),cos(angin));

        otherwise
            error('Unrecognized approach string')
    end

    angout = angledim(angout,'radians',units);  %  Convert to the original units


text:
classdef CSMaster < handle
    %CS_LAB superclass
    %
    % (c) Thomas Kuestner
    % ---------------------------------------------------------------------

    properties
        type % reconstruction type
        kSpace % save kSpace in class to reduce memory overhead
        iNOUTER % outer loop iterations
        iNINNER % inner loop iterations
        p % L-p norm for reweighted minimum norm
        lambda % lagrange multiplier for noise stability
        lambdaCalib % lagrange multiplier for calibration
        lambdaTV % lagrange multiplier for TV
        lambdaESPReSSo % lagrange multiplier for ESPReSSo constraint
        lambdaMC % lagrange multiplier for motion correction
%         oversampling % oversampling along space directions
%         dimension % 2D (x-y or x-y-t), 3D (x,y,z) or 4D (x,y,z,t) problem
%         dim % nPha x nFreq x nZ x nTime x nCha
        measPara % struct containing all measurement parameters (dim, dimension, ...)
        kernelSize % kernel window dimension for convolution with calibration data
        calibSize % k-space center needed for initial estimate of W or GRAPPA calibration
        fullMask % mask of acquired data points
        currSlice % current slice, just needed for 2D with several slices
        currCha % current channel, just needed for 2D/3D with several channels
        trafo % struct of (tight frame) transformation in sparse basis
        espresso % espresso parameters
        window % cut out window options
        lSense % apply sense mask for multicoil reconstruction
        dRecontime % reconstruction time
    end

    methods
        function obj = CSMaster(iNOUTER, iNINNER, p, lambda, measPara)
            % constructor
            obj.iNOUTER = iNOUTER;
            obj.iNINNER = iNINNER;
            obj.p       = p;
            obj.lambda  = lambda;
            obj.measPara = measPara;
        end
    end

end


text:
function H= stimutil_drawArrow(direction, varargin)
%STIMUTIL_DRAWARROW - Draw a thick arrow
%
%Synopsis:
% H= stimutil_drawArrow(DIRECTION, 'Prop1', VALUE1, ...)
%
%Arguments:
% DIRECTION: Derection of the arrow. Format is either a scalar
%   corresponding to a clock (3, 6, 9, 12), or a string ('right',
%   'down', 'left', 'up').
% Optional Properties:
%  'arrow_size': Size of the arrow in the units of the axis
%  'arrow_width': Width of the arrow's line, relative to arrow_size.
%  'arrow_blunt': Width of the arrow's blunt edges, relative to arrow_size.
%     Use 0 for sharp edges. Default: 0 if 'cross'=0, else 0.1.
%  'arrow_color': Default: 0.8*[1 1 1]
%  'arrow_edgecolor': Default 'none'.
%  'cross': Switch to superimpose a fiaxtion cross; Default: 0
%  'cross_size': Size of the cross, relative to arrow_size
%  'cross_width': Width of the cross' line, relative to arrow_size
%  'cross_color': Default: Color of the figure.
%  'cross_edgecolor': Default: 'none'.
%
%Returns
% H: Struct of handles of the graphical objects

% 11-2007 Benjamin Blankertz

props= {'arrow_size'        0.3           '!DOUBLE[1]'
        'arrow_width'       0.2           '!DOUBLE[1]'
        'arrow_blunt'       0             '!DOUBLE[1]'
        'arrow_tip_blunt'   0             '!DOUBLE[1]'
        'arrow_color'       0.8*[1 1 1]   '!DOUBLE[1 3]'
        'arrow_edgecolor'   'none'        'CHAR'
        'cross'             0             'BOOL'
        'cross_size'        0.2           '!DOUBLE[1]'
        'cross_width'       0.03          '!DOUBLE[1]'
        'cross_color'       0*[1 1 1]     '!DOUBLE[1 3]'
        'cross_edgecolor'   'none'        'CHAR'
       };

if nargin==0,
  H= props; return;
end

opt= opt_proplistToStruct(varargin{:});
opt= opt_setDefaults(opt, props, 1);

x= opt.arrow_size;
w= x*opt.arrow_width;
b= x*opt.arrow_blunt;
t= x*opt.arrow_tip_blunt;

%xpr= [0 x 0 0 -x+w -x+w 0];
%ypr= [x 0 -x -w -w w w];
xpr= [b x x b -b -b -x+w -x+w -b -b];
ypr= [x t -t -x -x -w -w w w x];
switch(direction),
 case {'up','tongue',0,12},
  xp= ypr;
  yp= xpr;
 case {'right',3},
  xp= xpr;
  yp= ypr;
 case {'down','foot',6},
  xp= -ypr;
  yp= -xpr;
 case {'left',9},
  xp= -xpr;
  yp= -ypr;
 otherwise,
  error('invalid direction: use 3, 6, 9, 12; or ''right'', ''down'', ...');
end

H.arrow= patch(xp, yp, opt.arrow_color);
set(H.arrow, 'EdgeColor', opt.arrow_edgecolor);

if opt.cross,
  c= x*opt.cross_size;
  v= x*opt.cross_width;
  H.cross= patch([v v c c v v -v -v -c -c -v -v], ...
                 [c v v -v -v -c -c -v -v v v c], opt.cross_color);
  set(H.cross, 'EdgeColor', opt.cross_edgecolor);
end

author: bbci


text:
% nlsS = getNLSStruct( extra, dispOn, zoom)
%
% extra.tVec    : defining TIs
%               (not called TIVec because it looks too much like T1Vec)
% extra.T1Vec   : defining T1s
% dispOn        : 1 - display the struct at the end
%                 0 (or omitted) - no display
% zoom          : 1 (or omitted) - do a non-zoomed search
%                 x>1 - do an iterative zoomed search (x-1) times (slower search)
%                 NOTE: When zooming in, convergence is not guaranteed.
%
% Data Model    : a + b*exp(-TI/T1)
%
% written by J. Barral, M. Etezadi-Amoli, E. Gudmundson, and N. Stikov, 2009
%  (c) Board of Trustees, Leland Stanford Junior University

function nlsS = getNLSStruct(extra, dispOn, zoom)

nlsS.tVec = extra.tVec(:);
nlsS.N = length(nlsS.tVec);
nlsS.T1Vec = extra.T1Vec(:);
nlsS.T1Start = nlsS.T1Vec(1);
nlsS.T1Stop = nlsS.T1Vec(end);
nlsS.T1Len = length(nlsS.T1Vec);

% The search algorithm to be used
nlsS.nlsAlg = 'grid'; % Grid search

% Display the struct so that the user can see it went ok
if nargin < 2
  dispOn = 0;
end

% Set the number of times you zoom the grid search in, 1 = no zoom
% Setting this greater than 1 will reduce the step size in the grid search
% (and slow down the fit significantly)

if nargin < 3
	nlsS.nbrOfZoom = 2;
else
	nlsS.nbrOfZoom = zoom;
end

if nlsS.nbrOfZoom > 1
    nlsS.T1LenZ = 21; % Length of the zoomed search
end

if dispOn
  % Display the structure for inspection
  nlsS
end

% Set the help variables that can be precomputed:
% alpha is 1/T1,
% theExp is a matrix of exp(-TI/T1) for different TI and T1,
% rhoNormVec is a vector containing the norm-squared of rho over TI,
% where rho = exp(-TI/T1), for different T1's.
switch(nlsS.nlsAlg)
  case{'grid'}
    alphaVec = 1./nlsS.T1Vec;
    nlsS.theExp = exp( -nlsS.tVec*alphaVec' );
    nlsS.rhoNormVec = ...
        sum( nlsS.theExp.^2, 1)' - ...
        1/nlsS.N*(sum(nlsS.theExp,1)').^2;
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [behav] = processConvertAuxilToAccel(fbasename,varargin)
% USAGE
%
%   processConvertAuxilToAccel(fbasename,varargin)
%
%
%  INPUTS
%    fbasename   basename of the dat file
%
%
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'syncDatFile' name of binary file where sync signal is stored
%                   (default = filebasename_digitalin.dat)
%     'accelFile'   basename of .dat with accelerometry data
%     'recSampFq'  sampling freqeuncy of the sync signal (default= 20kHz)
%     'syncChan'    sync channel to read (default = 1)
%     'numChans'    number of channels in the sync file (default = 1)
%     'accelSampFq'   sampling frequency of the final pos file (after
%                   interpolation)
%    =========================================================================
%
% OUTPUT
%    behav a behavior struct with the following fields added
%    behav.accel.x - acceleration in the x direction
%    behav.accel.y - acceleration in the x direction
%    behav.accel.z - acceleration in the x direction
%    behav.accel.ts - timestamps of downsamples acceleration signal

accelFile = 'auxiliary.dat'; % default intan file
recSampFq  = 20000;
numChans = 3;
accelSampFq = 1250;

for i = 1:2:length(varargin)
  if ~isa(varargin{i},'char')
    error(['Parameter ' num2str(i+3) ' is not a property .']);
  end
  switch(lower(varargin{i}))
    case 'accelfile'
      syncDatFile = varargin{i+1};
      if ~isa(duration,'char')
        error('Incorrect value for property ''syncDatFile'' .');
      end
    case 'recsampFq'
      recSampFq = varargin{i+1};
      if ~isa(frequency,'numeric')
        error('Incorrect value for property ''recSampFq'' .');
      end
    case 'numchans'
      numChans = varargin{i+1};
      if ~isa(numChans,'numeric')
        error('Incorrect value for property ''numChans'' .');
      end
    case 'accelsampfq'
      posSampFq = varargin{i+1};
      if ~isa(posSampFq,'numeric')
        error('Incorrect value for property ''posSampFq'' .');
      end
    otherwise
      error(['Unknown property ''' num2str(varargin{i}) ''' .']);
  end
end

if ~exist(accelFile)
    error(['couldnt find ' accelFile ' file for accelerometer data...'])
end
    % read data into workspace
    fid = fopen(accelFile);
    dat = fread(fid,[numChans inf],'uint16=>uint16');
    fclose(fid);

%     fid = fopen('time.dat');
%     ts = fread(fid,[1,inf],'int32=>int32');  % we can do this without IO
    ts = 0:1/recSampFq:length(dat)-1/recSampFq;
%     fclose(fid);

    % downsample to 1250 Hz
    behav.accel.x = uint16(resample((double(dat(1,:))),accelSampFq,recSampFq));
    behav.accel.y = uint16(resample((double(dat(2,:))),accelSampFq,recSampFq));
    behav.accel.z = uint16(resample((double(dat(3,:))),accelSampFq,recSampFq));
    behav.accel.ts = ts;
end


text:
function n = p32_n ( )

%*****************************************************************************80
%
%% P32_N returns the number of variables for problem 32.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer N.  If N is positive, N represents the
%    only legal value for N for this problem.  If N is
%    negative, the absolute value of N is the least legal
%    value of N, but other values are allowable.
%
  n = 4;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function cascade = separable_wavelet_factory_2d(N, filter_options, scat_options, M)
	filters{1} = filter_bank(N(1),filter_options{1});
	filters{2} = filter_bank(N(2),filter_options{2});

	for m = 0:M
		filt1_ind = min(numel(filters{1}), m+1);
		filt2_ind = min(numel(filters{2}), m+1);
		cascade{m+1} = @(X)(separable_wavelet_layer_2d(X, ...
			{filters{1}{filt1_ind} filters{2}{filt2_ind}}, scat_options));
	end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [h, array] = display_network(A, opt_normalize, opt_graycolor, cols, opt_colmajor)
% This function visualizes filters in matrix A. Each column of A is a
% filter. We will reshape each column into a square image and visualizes
% on each cell of the visualization panel.
% All other parameters are optional, usually you do not need to worry
% about it.
% opt_normalize: whether we need to normalize the filter so that all of
% them can have similar contrast. Default value is true.
% opt_graycolor: whether we use gray as the heat map. Default is true.
% cols: how many columns are there in the display. Default value is the
% squareroot of the number of columns in A.
% opt_colmajor: you can switch convention to row major for A. In that
% case, each row of A is a filter. Default value is false.
oldwarning = warning;
warning off all

figure();

if ~exist('opt_normalize', 'var') || isempty(opt_normalize)
    opt_normalize= true;
end

if ~exist('opt_graycolor', 'var') || isempty(opt_graycolor)
    opt_graycolor= true;
end

if ~exist('opt_colmajor', 'var') || isempty(opt_colmajor)
    opt_colmajor = false;
end

% rescale
A = A - mean(A(:));

if opt_graycolor, colormap(gray); end

% compute rows, cols
[L M]=size(A);
sz=sqrt(L);
buf=1;
if ~exist('cols', 'var')
    if floor(sqrt(M))^2 ~= M
        n=ceil(sqrt(M));
        while mod(M, n)~=0 && n<1.2*sqrt(M), n=n+1; end
        m=ceil(M/n);
    else
        n=sqrt(M);
        m=n;
    end
else
    n = cols;
    m = ceil(M/n);
end

array=-ones(buf+m*(sz+buf),buf+n*(sz+buf));

if ~opt_graycolor
    array = 0.1.* array;
end

if ~opt_colmajor
    k=1;
    for i=1:m
        for j=1:n
            if k>M,
                continue;
            end
            clim=max(abs(A(:,k)));
            if opt_normalize
                array(buf+(i-1)*(sz+buf)+(1:sz),buf+(j-1)*(sz+buf)+(1:sz))=reshape(A(:,k),sz,sz)/clim;
            else
                array(buf+(i-1)*(sz+buf)+(1:sz),buf+(j-1)*(sz+buf)+(1:sz))=reshape(A(:,k),sz,sz)/max(abs(A(:)));
            end
            k=k+1;
        end
    end
else
    k=1;
    for j=1:n
        for i=1:m
            if k>M,
                continue;
            end
            clim=max(abs(A(:,k)));
            if opt_normalize
                array(buf+(i-1)*(sz+buf)+(1:sz),buf+(j-1)*(sz+buf)+(1:sz))=reshape(A(:,k),sz,sz)/clim;
            else
                array(buf+(i-1)*(sz+buf)+(1:sz),buf+(j-1)*(sz+buf)+(1:sz))=reshape(A(:,k),sz,sz);
            end
            k=k+1;
        end
    end
end

if opt_graycolor
    % h=imagesc(array,'EraseMode','none',[-1 1]);
    h=imagesc(array, [-1 1]);
else
    % h=imagesc(array,'EraseMode','none',[-1 1]);
    h=imagesc(array, [-1 1]);
end
axis image off
drawnow;

% warning on all
warning(oldwarning);

author: zellyn


text:
function coef = plotnsdgtreal(coef,a,varargin)
%PLOTNSDGTREAL Plot NSDGTREAL coefficients
%   Usage:  plotnsdgtreal(c,a,fs,dynrange);
%
%   Input parameters:
%         coef     : Cell array of coefficients.
%         a        : Vector of time positions of windows.
%         fs       : signal sample rate in Hz (optional).
%         dynrange : Colorscale dynamic range in dB (optional).
%
%   `plotnsdgtreal(coef,a)` plots coefficients computed using |nsdgtreal| or
%   |unsdgtreal|. For more details on the format of the variables *coef* and *a*,
%   please read the function help for these functions.
%
%   `plotnsdgtreal(coef,a,fs)` does the same assuming a sampling rate of
%   *fs* Hz of the original signal.
%
%   `plotnsdgtreal(coef,a,fs,dynrange)` additionally limits the dynamic range.
%
%   `C=plotnsdgtreal(...)` returns the processed image data used in the
%   plotting. Inputting this data directly to `imagesc` or similar
%   functions will create the plot. This is useful for custom
%   post-processing of the image data.
%
%   `plotnsdgtreal` supports all the optional parameters of |tfplot|. Please
%   see the help of |tfplot| for an exhaustive list. In addition, the
%   following parameters may be specified:
%
%     'xres',xres  Approximate number of pixels along x-axis /time.
%                  Default value is 800
%
%     'yres',yres  Approximate number of pixels along y-axis / frequency
%                  Default value is 600
%
%   See also: tfplot, nsdgt, nsdgtreal

%   AUTHOR : Florent Jaillet & Peter L. Søndergaard
%   TESTING: OK
%   REFERENCE: NA

if nargin<2
  error('%s: Too few input parameters.',upper(mfilename));
end;

definput.import={'ltfattranslate','tfplot'};

definput.keyvals.xres=800;
definput.keyvals.yres=600;

[flags,kv,fs]=ltfatarghelper({'fs','dynrange'},definput,varargin);

timepos=cumsum(a)-a(1);

N=length(a);
cwork=zeros(kv.yres,N);

%% -------- Interpolate in frequency ---------------------

for ii=1:N
  column=coef{ii};
  M=length(column);
  cwork(:,ii)=interp1(linspace(0,1,M),column,linspace(0,1,kv.yres),'nearest');
end;

%% --------  Interpolate in time -------------------------

% Time step in next equidistant spacing on the x-axis (in samples)
aplot=timepos(end)/kv.xres;

% Time positions where we want our pixels plotted (in samples)
xr=(0:kv.xres-1)*aplot;

coef=zeros(kv.yres,kv.xres);
for ii=1:kv.yres
  data=interp1(timepos,cwork(ii,:).',xr,'nearest').';
  coef(ii,:)=data;
end;

yr=[0,1];

coef=tfplot(coef,aplot,yr,'argimport',flags,kv);

if nargout<1
    clear coef;
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [ Pos ] = getEEFCartesianPosition( t )
%% This function is used to get the endeffector cartizian position of the KUKA iiwa 7 R 800.

%% Syntax:
% [ Pos ] = getEEFCartesianPosition( t )

%% About
% This function is used to get the endeffector Cartesian position of the KUKA iiwa 7 R 800.
% The position is of the media flange of the robot, and it is measured
% relative to the robot base frame.

%% Arreguments:
% t: is the TCP/IP connection
% Pos: is 1x3 cell array. Representing the  X,Y and Z positions of the
% endeffector. The returned values are in (millimeters)

% Copyright, Mohammad SAFEEA, 3rd of May 2017

theCommand='Eef_pos';
fprintf(t, theCommand);
message=fgets(t);
%disp(message)
[P1,N]=getDoubleFromString(message);

for i=1:3
    Pos{i}=P1{i};
end

end

function [jPos,j]=getDoubleFromString(message)
n=max(max(size(message)));
j=0;
numString=[];
for i=1:n
    if message(i)=='_'
        j=j+1;
        jPos{j}=str2num(numString);
        numString=[];
    else
        numString=[numString,message(i)];
    end
end
end

repo: KST-Kuka-Sunrise-Toolbox


text:
%% test modulate for all oasis functions.
col = {[0 114 178],[0 158 115], [213 94 0],[230 159 0],...
    [86 180 233], [204 121 167], [64 224 208], [240 228 66]}; % colors
plot_cvx = false;

%% example 2: foopsi, AR2 model
g = [1.7, -0.712];         % AR coefficient
noise = 1;
T = 3000;
framerate = 30;
firerate = 0.5;
b = 0;              % baseline
N = 20;              % number of trials
seed = 3;          % seed for genrating random variables
[Y, trueC, trueS] = gen_data(g, noise, T, framerate, firerate, b, N, seed);
y = Y(1,:);
true_c = trueC(1,:);  %#ok<*NASGU>
true_s = trueS(1,:);
% case 1: all parameters are known
lambda = 25;
[c_oasis, s_oasis] = deconvolveCa(y, 'ar2', g, 'foopsi', 'lambda', lambda);  %#ok<*ASGLU>

figure('name', 'FOOPSI, AR2, known: g, lambda', 'papersize', [15, 4]);
show_results;

% case 2: know lambda
lambda = 2.5;
[c_oasis, s_oasis, options] = deconvolveCa(y, 'ar2', 'sn', noise, 'foopsi', 'lambda',...
    lambda);
fprintf('true gamma:        %.3f\t %.3f\n', g(1), g(2));
fprintf('estimated gamma:   %.3f\t %.3f\n', options.pars(1),  options.pars(2));

figure('name', 'FOOPSI, AR2, known:lambda, estimated: g', 'papersize', [15, 4]);
show_results;

%%%%%%%%%%%%%%  END %%%%%%%%%%%%%%%%%%

author: flatironinstitute

repo: CaImAn-MATLAB