text:
function iceweb2017(thissubnet, ds, ChannelTagList, ...
    snum, enum, nummins, products, PARAMS, paths)
    debug.printfunctionstack('>');

    % load state
    statefile = sprintf('iceweb_%s_state.mat',thissubnet);
    if exist(statefile, 'file') && ~strcmp(PARAMS.runmode, 'test')
        load(statefile)
    end

    % end time
    if enum==0
        enum = utnow - delaymins/1440;
    end

    % loop one day at a time
    for dnum = floor(snum):ceil(enum)
        disp(datestr(dnum))

%         % figure out which channeltags are active for this day?
%         ctags_thisday = get_channeltags_thisday(ds, ctags, dnum); % subset to channeltags valid
%         if isempty(todaysites)
%             continue;
%         end

%         % subset to sites for which the files pointed to by the
%         % wfdisc table, actually exist
%         m = listMiniseedFiles(ds, chantag, dnum, dnum+1);
%         ctags_thisday = ctags_thisday([m.exists]==2);

        % get timewindow for this day
        timewindows = get_timewindow(min([enum dnum+1]), nummins, max([dnum snum]));

        if ~strcmp(PARAMS.runmode,'test')
            % loop over timewindows
            for count = 1:length(timewindows.start)
                this_timewindow.start = timewindows.start(count);
                this_timewindow.stop = timewindows.stop(count);
                iceweb_helper(paths, PARAMS, newsubnets, this_timewindow, ds, products);
            end
        end
    end
    debug.printfunctionstack('<');
end

function iceweb_helper(paths, PARAMS, subnets, tw, ds, products)
    debug.printfunctionstack('>');

    MILLISECOND_IN_DAYS = (1 / 86400000);

    makeSamFiles = false;
    makeSoundFiles = true;

    snum = tw.start;
    enum = tw.stop - MILLISECOND_IN_DAYS; % try to skip last sample

    % load state
    statefile = sprintf('iceweb_%s_state.mat',subnet);
    if exist(statefile, 'file')
        load(statefile)
        if snum < snum0
            continue
        end
    end

    % save state
    ds0=ds; sites0=sites; snum0=snum; enum0=enum; subnet0 = subnet;
    save(statefile, 'ds0', 'sites0', 'snum0', 'enum0', 'subnet0');
    clear ds0 sites0 snum0 enum0 subnet0

    %% Save raw waveform data to 1 hour MAT file
    dv = datevec(snum);
    yyyy = dv(1);
    jjj = floor(snum - datenum(yyyy,1,1) + 1);
    hh = dv(4);
    wavematdir = 'wavemat';
    wavrawmat = sprintf('%s/%s/%04d/%03d/%s.%04d.%03d.%02d.raw.mat',wavematdir,subnet,yyyy,jjj,subnet,yyyy,jjj,hh);
    if ~exist(wavrawmat,'file')
        %% Get waveform data
        debug.print_debug(0, sprintf('%s %s: Getting waveforms for %s from %s to %s at %s',mfilename, datestr(utnow), subnet , datestr(snum), datestr(enum)));
        w = waveform_wrapper(ds, ChannelTagList, snum, enum);
        mkdir(fileparts(wavmatfile));
        disp(sprintf('Saving waveform data to %s',wavrawmat));
        save(wavrawmat);
    end

    % continue

    % Save the cleaned waveform data to MAT file
    wavcleanmat = sprintf('%s/%s/%04d/%03d/%s.%04d.%03d.%02d.clean.mat',wavematdir,subnet,yyyy,jjj,subnet,yyyy,jjj,hh);
    if ~exist(wavcleanmat,'file')

        % Eliminate empty waveform objects
        w = waveform_remove_empty(w);
        if numel(w)==0
            debug.print_debug(0, 'No waveform data returned - skipping');
            continue
        end

        % Clean the waveforms
        w = fillgaps(w, 'interp');
        w = detrend(w);

        % Apply calibs which should be stored within sites structure to
        % waveform objects to convert from counts to real physical
        % units
        w = apply_calib(w, sites);

        % Pad all waveforms to same start/end
        [wsnum wenum] = gettimerange(w); % assume gaps already filled, signal
        w = pad(w, min([snum wsnum]), max([enum wenum]), 0);

        % Apply filter to all signals
        w = apply_filter(w, PARAMS);

        mkdir(fileparts(wavmatfile));
        disp(sprintf('Saving waveform data to %s',wavcleanmat));
        save(wavcleanmat);
    end

    %% ICEWEB PRODUCTS

    % WAVEFORM PLOT
    if products.waveform_plot.doit
        close all
        plot_panels(w)
        [spdir,spbase,spext] = fileparts(spectrogramFilename);
        fname = fullfile('plots', 'waveforms', subnet, sprintf('%s.png',datestr(snum,31)) );
        orient tall;
        saveImageFile(fname, 72); % this should make directory tree too
    end

    % RSAM
    if products.rsam.doit
        for measure = products.rsam.measures
            rsamobj = waveform2rsam(w);
            rsamobj.save_to_bob_file(fullfile('data', 'rsam', subnet, 'SSSS.CCC.YYYY.MMMM.bob'));
        end
    end

    % CREATE & SAVE HELICORDER PLOT
    if products.helicorders.doit
        close all
        for wi=1:numel(w)
            heliplot = helicorder(w(wi),'mpl',3);
            build(heliplot)
            ct = get(w(wi),'channeltag');
            helicorderFilename = fullfile(spdir, sprintf('heli_%s_%s.%s.%s',spbase,ct.station,ct.channel,spext));
            orient tall;
            saveImageFile(helicorderFilename, 72);
            clear ct helicorderFilename
        end
        clear wi
    end

    % PLOT SPECTROGRAMS
    if products.spectrograms.doit
        debug.print_debug(1, sprintf('Creating %s',spectrogramFilename))
        waveSp = w;

		% restrict spectrogram to Z channels only
		w2=[];
		for wi=1:numel(w)
			ct=get(w(wi),'channeltag');
			if strfind(ct.channel,'Z')
				w2 = [w2 w(wi)];
			end
        end
        waveSp = w2;

        % split into 10-minute timewindows
        for someloopover10mintimewindows
            % make 10-minute spectrogram
            close all
            spectrogramFraction = 0.75;
            [sgresult, Tcell, Fcell, Ycell] = spectrogram_iceweb(PARAMS.spectralobject, waveSp, spectrogramFraction, extended_spectralobject_colormap);
            clear w2 wi ct

            if sgresult > 0 % sgresult = number of waveforms for which a spectrogram was successfully plotted
                % SAVE SPECTROGRAM PLOT TO IMAGE FILE AND CREATE THUMBNAIL
                orient tall;

                if saveImageFile(spectrogramFilename, 72)

                    fileinfo = dir(spectrogramFilename); % getting a weird Index exceeds matrix dimensions error here.
                    debug.print_debug(0, sprintf('%s %s: spectrogram PNG size is %d',mfilename, datestr(utnow), fileinfo.bytes));

                    % make thumbnails
                    makespectrogramthumbnails(spectrogramFilename, spectrogramFraction);

                end

                % save spectral data
                if products.spectral_data.doit
                    frequency_index_divider = 5.0; % Hz
                    for spi = 1:numel(Fcell)
                        thisY = Ycell{spi};
                        thisF = Fcell{spi};
                        thisT = Tcell{spi};
                        sta = get(w(spi),'station');
                        chan = get(w(spi),'channel');

                        % peakF and meanF for each spectrogram window
                        [Ymax,imax] = max(thisY);
                        PEAK_F = thisF(imax);
                        MEAN_F = (thisF' * thisY)./sum(thisY);

                        % frequency index for each spectrogram window
                        fUpperIndices = find(thisF > frequency_index_divider);
                        fLowerIndices = find(thisF < frequency_index_divider);
                        fupper = sum(thisY(fUpperIndices,:));
                        flower = sum(thisY(fLowerIndices,:));
                        F_INDEX = log2(fupper ./ flower);

                        % Peak spectral value in each frequency bin - or in
                        % each minute?

                        % Now downsample to 1 sample per minute
                        dnum = unique(floorminute(thisT));
                        for k=1:length(dnum)
                            p = find(floorminute(thisT) == dnum(k));
                            downsampled_peakf(k) = nanmean(PEAK_F(p));
                            downsampled_meanf(k) = nanmean(MEAN_F(p));
                            downsampled_findex(k) = nanmean(F_INDEX(p));
                            suby = thisY(:,p);
                            if size(suby,2) == 1
                                            max_in_each_freq_band(:,k) = suby;
                            else
                                            max_in_each_freq_band(:,k) = max(suby');
                            end
                        end

        %                 % Plot for verification
        %                 close all
        %                 subplot(2,1,1),plot(dnum,downsampled_peakf,':');datetick('x');
        %                 hold on
        %                 plot(dnum,downsampled_meanf);datetick('x');
        %                 subplot(2,1,2),plot(dnum,downsampled_findex);datetick('x');
        %                 anykey = input('Press any key to continue');

                        % Save data
                        r1 = rsam(dnum, downsampled_peakf, 'sta', sta, ...
                            'chan', chan, 'measure', 'peakf', ...
                            'units', 'Hz', 'snum', min(dnum), 'enum', max(dnum));
                        r1.save_to_bob_file(fullfile('spectrograms', subnet, 'SSSS.CCC.YYYY.peakf'))

                        r2 = rsam(dnum, downsampled_meanf, 'sta', sta, ...
                            'chan', chan, 'measure', 'meanf', ...
                            'units', 'Hz', 'snum', min(dnum), 'enum', max(dnum));
                        r2.save_to_bob_file(fullfile('spectrograms', subnet, 'SSSS.CCC.YYYY.meanf'))

                        r3 = rsam(dnum, downsampled_findex, 'sta', sta, ...
                            'chan', chan, 'measure', 'findex', ...
                            'units', 'none', 'snum', min(dnum), 'enum', max(dnum));
                        r3.save_to_bob_file(fullfile('spectrograms', subnet, 'SSSS.CCC.YYYY.findex'))

                        specdatafilename = fullfile('data', 'spectral', subnet, datestr(min(dnum),'yyyy/mm/dd'), sprintf( '%s_%s_%s.mat', datestr(min(dnum),30), sta, chan) );
                        specdatadir = fileparts(specdatafilename); % make the directory in case it does not exist
                        mkdir(specdatadir); % make the directory in case it does not exist
                        save(specdatafilename, 'dnum', 'max_in_each_freq_band')
                        clear r1 r2 r3 k p   downsampled_peakf downsampled_meanf ...
                            downsampled_findex fUpperIndices fLowerIndices flower ...
                            fupper F_INDEX PEAK_F MEAN_F Ymax imax thisY thisF thisT ...
                            suby max_in_each_freq_band

                    end
                end
            end
        end
    end

    % SOUND FILES
    if products.soundfiles.doit
        try
            % 20120221 Added a "sound file" like 201202211259.sound which simply records order of stachans in waveform object so
            % php script can match spectrogram panel with appropriate wav file
            % 20121101 GTHO COmment: Could replace use of bnameroot below with strrep, since it is just used to change file extensions
            % e.g. strrep(spectrogramFilename, '.png', sprintf('_%s_%s.wav', sta, chan))
            soundfileroot = fullfile('soundfiles', subnet, datestr(get(w(c),'start'),31),
            [dname, bnameroot, bnameext] = fileparts(soundfileroot);
            soundfilelist = sprintf('%s.sound',soundfileroot);
            fsound = fopen(soundfilelist,'a');
            for c=1:length(w)
                soundfilename = fullfile('soundfiles', sprintf('%s_%s_%s.wav',bnameroot, get(w(c),'station'), get(w(c), 'channel')  ) );
                fprintf(fsound,'%s\n', soundfilename);
                debug.print_debug(0, sprintf('Writing to %s',soundfilename));
                isSuccessful = waveform2sound(w(c), 60, soundfilename)
            end
            fclose(fsound);
        end
    end

    debug.printfunctionstack('<');
end

author: geoscience-community-codes

repo: GISMO


text:
function [MSE, sequences, indexes, mx, sds, mn, sd] = motionCompDetectMotionMSE(view,ROI,scans,plotMSE,histogramClip)

%    [MSE, sequences, indexes, mx, sds, mn, sd] = motionCompDetectMotionMSE(view, [ROI], [scans], [plotMSE], [histogramClip])
%
% gb 01/26/05
%
% Remove outliers within a scan and returns an array 'sequences' whose rows
% are the sequences between the outliers. The error computed is a Mean
% Squared Error.
%
% Input arguments :
%       - view : current inplane
%       - ROI : Region of interest where the error has to be computed
%           Default = ones
%       - scans : scans to analyse.
%           Default = all scans
%       - plotMSE : indicates if the user wants a plot of the MSE error.
%       Should be the value '1' if so.
%           Default = 1
%       - histogramClip : sets the clipping value of the histogram
%           Default = 0.9
%
% It returns :
%       - MSE = normalized squared sum of the difference of two
%       consecutive frames
%       - sequences = array representing the resulting sequences between the outliers
%       - indexes = indexes correponding to the outliers
%       - mx = maximum of MSE
%       - sds = array representing the standard deviations of each sequence
%       - mn = mean of MSE after having been clipped
%       - sd = standard deviation of MSE after having been clipped

% Initializes arguments and variables
if ieNotDefined('histogramClip')
    histogramClip = .9;
end

if ieNotDefined('scans')
    scans = selectScans(view);
end

if ieNotDefined('ROI')
    ROI = '';
end

if ieNotDefined('resize')
    resize = 0;
end

% Starts the main calculation
% For each scan, it computes the difference between the first frame
% of this scan and the last frame of the preceding scan before the error
% between consecutive frames.

MSE = [];
lastFrame = 0;
h = mrvWaitbar(0,'Computing Mean Squared Error...');

for scan = scans

    mrvWaitbar((scan - 1)/scans(end),h,['Computing Mean Squared Error for scan ' num2str(scan)]);

    % Loads the tSeries
    clear tSeriesAllSlices;
    tSeriesAllSlices = motionCompLoadImages(view,scan);
    nFrames  = size(tSeriesAllSlices,1);
    nVoxelsX = size(tSeriesAllSlices,2);
    nVoxelsY = size(tSeriesAllSlices,3);
    nSlices  = size(tSeriesAllSlices,4);

    % Computes the Mean Squared Error when changing scans
    firstFrame = tSeriesAllSlices(1,:,:,:);

    if isequal(lastFrame,0)
        MSE = [MSE;0];
    else
        MSE = [MSE;motionCompMSE(firstFrame,ROI,lastFrame)];
    end

    lastFrame = tSeriesAllSlices(nFrames,:,:,:);

    % Computes the difference between two consecutive frames
    tSeriesAllSlices = diff(tSeriesAllSlices,1,1);

	% Computes the Mean Square Error
	MSE = [MSE;motionCompMSE(tSeriesAllSlices,ROI)];

end

mrvWaitbar(1,h);

% Clips the histogram to calculate the mean and the standard deviation
% of MSE without taking into account the outliers
histClip = mrAnatHistogramClip(MSE,0,histogramClip,0);

% Computes the mean and the standard deviation of the clipped MSE
mn = mean(histClip);
sd =  std(histClip);

% Finds the outliers. Are considered outliers all the values above the mean
% + two standard deviations
indexes = (find(MSE > mn + 2*sd));
mx = max(MSE);

% Finds the resulting sequences between the outliers.
sequences = horzcat([1;indexes + 2],[indexes - 2;nFrames - 1]);
sequences = sequences(find(sequences(:,2) - sequences(:,1) > 5),:);

% Computes all the standard deviations
sds = zeros(size(sequences,1),1);
for i = 1:size(sequences,1)
    sds(i) = std(MSE(sequences(i,1):sequences(i,2)));
end

% Plots the curves
if ieNotDefined('plotMSE')
    plotMSE = 0;
end

if plotMSE > 0
    figure
    title(['MSE for ' view.sessionCode ]);

    plot(MSE)
    hold on
    plot(1:length(MSE),mn*ones(1,length(MSE)),'-.')
    plot([1 length(MSE)],[mn + 2*sd, mn + 2*sd],'r')

    % Another measure can be used. Sometimes we can assume that the signal
    % is equally distributed around its mean. The explicit measure is :
    % mn + (mn - min(MSE))
    %
    % plot([1 length(MSE)],[2*mn - min(MSE),2*mn - min(MSE)],'black')

    for scanIndex = 1:(length(scans) - 1)
        plot([nFrames*scanIndex nFrames*scanIndex],[min(MSE) mx],'-.','Color','m');
    end

end

close(h)

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%%******************************************************************
%% geometric_mean: an example with rotated cone variables
%%
%%  max {prod(d+B*x) : d+B*x > 0, x <= 10}
%%
%%  where B = 4xn matrix,
%%        d = 4x1 vector
%%
%% [blk,At,C,b,xx] = geometric_mean(B,d,solve);
%%
%% E.g. p = 6; m = 4; B = [rand(2,p); -rand(2,p)]; d = rand(m,1);
%%
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

   function [blk,At,C,b,xx] = geometric_mean(B,d,solve);

   if (nargin == 2); solve = 0; end

   n = size(B,2);
   zz = zeros(1,n);
   r2 = sqrt(2);
   blk{1,1} = 'r'; blk{1,2} = [3,3,3];
   blk{2,1} = 'l'; blk{2,2} = [2];
   At{1,1} = -[B(1,:),0,0,0; B(2,:),0,0,0; zz,r2,0,0; ...
               B(3,:),0,0,0; B(4,:),0,0,0; zz,0,r2,0; ...
	       zz, 1,0,0; zz, 0,1,0; zz, 0,0,r2];
   At{2,1} = -[B(1,:),0,0,0; B(3,:),0,0,0];
   C{1,1} = [d(1:2); 0; d(3:4); 0; 0;0;0];
   C{2,1} = [d(1); d(3)];
   b = [zeros(n,1); 0;0;1];
   blk{3,1} = 'l'; blk{3,2} = n;
   At{3,1} = [eye(n), zeros(n,3)]; C{3,1} = 10*ones(n,1);
   if (solve)
      [bblk,AAt,CC,bb,T] = convertRcone(blk,At,C,b);
      [obj,X,y,Z] = sqlp(bblk,AAt,CC,bb);
      xx = y(1:n);
   end
%%******************************************************************

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


% resolve conflict between trackers
function trackers = handle_conflicting_trackers(trackers, bboxes_det, opt)

% collect bboxes from trackers
bboxes_track = [];
for i = 1:numel(trackers)
    tracker = trackers{i};
    bboxes = sub_bboxes(tracker.bboxes, numel(tracker.bboxes.fr));

    if tracker.state == opt.STATE_TRACKED || tracker.state == opt.STATE_ACTIVATED
        if isempty(bboxes_track)
            bboxes_track = bboxes;
        else
            bboxes_track = concatenate_bboxes(bboxes_track, bboxes);
        end
    end

    if tracker.state == opt.STATE_STOP
        tracker.eco = [];
    end
end

% compute overlaps
num_det = numel(bboxes_det.fr);
if isempty(bboxes_track)
    num_track = 0;
else
    num_track = numel(bboxes_track.fr);
end

flag = zeros(num_track, 1);
for i = 1:num_track
    if flag(i) == 1
        continue;
    end
    [~, o] = calc_overlap(bboxes_track, i, bboxes_track, 1:num_track);
    o(i) = 0;
    o(flag == 1) = 0;
    [max_ov, j] = max(o);
    if max_ov > opt.overlap_sup
        num1 = trackers{bboxes_track.id(i)}.num_tracked;
        num2 = trackers{bboxes_track.id(j)}.num_tracked;
        o1 = max(calc_overlap(bboxes_track, i, bboxes_det, 1:num_det));
        o2 = max(calc_overlap(bboxes_track, j, bboxes_det, 1:num_det));

        if num1 > num2
            suppressed_idx = j;
            winner_idx = i;
        elseif num1 < num2
            suppressed_idx = i;
            winner_idx = j;
        else
            if o1 > o2
                suppressed_idx = j;
                winner_idx = i;
            else
                suppressed_idx = i;
                winner_idx = j;
            end
        end

        if numel(trackers{bboxes_track.id(suppressed_idx)}.bboxes.fr) == 1
            trackers{bboxes_track.id(suppressed_idx)}.state = opt.STATE_STOP;
            trackers{bboxes_track.id(suppressed_idx)}.bboxes.state(end) = opt.STATE_STOP;
        else
            trackers{bboxes_track.id(suppressed_idx)}.state = opt.STATE_LOST;
            trackers{bboxes_track.id(suppressed_idx)}.bboxes.state(end) = opt.STATE_LOST;
        end
        if opt.is_print
            fprintf('target %d suppressed by %d\n', bboxes_track.id(suppressed_idx), bboxes_track.id(winner_idx));
        end
        flag(suppressed_idx) = 1;
    end
end


text:
function svgGaussianBlur(s, source, deviation, result)
% Adds a feGaussianBlur SVG filter
% PRELIMINARY IMPLEMENTATION (Parameters may change)
%
% svgGaussianBlur(s, source, deviation, result)
% Parameters:
%   s : Array of plot object handles
%   source : Any previous defined filter result string, 'SourceGraphic',
%            or 'SourceAlpha'.
%   deviation : Blur strength
%   result : String that identifies the filter result for following filter
%            stages.
for i = 1:length(s)
    userdata = get(s(i),'UserData');
    if isfield(userdata, 'svg') && isfield(userdata.svg, 'Filter')
        next = length(userdata.svg.Filter) + 1;
    else
        next = 1;
    end
    userdata.svg.Filter(next).Subfilter.Type = 'feGaussianBlur';
    userdata.svg.Filter(next).Subfilter.Deviation = deviation;
    userdata.svg.Filter(next).Subfilter.Source = source;
    userdata.svg.Filter(next).Subfilter.Result = result;
    set(s(i),'UserData', userdata);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function indices = argmin(v)
% ARGMIN Return as a subscript vector the location of the smallest element of a multidimensional array v.
% indices = argmin(v)
%
% Returns the first minimum in the case of ties.
% Example:
% X = [2 8 4; 7 3 9];
% argmin(X) = [1 1], i.e., row 1 column 1

[m i] = min(v(:));
indices = ind2subv(mysize(v), i);
%indices = ind2subv(size(v), i);


text:
function [EEG com] = pop_modprobplot(EEG,models2plot,smoothed,smoothlength,store,events2plot)
if nargin<1
    help pop_modprobplot
    return
end

if ~isfield(EEG.etc,'amica') || isempty(EEG.etc.amica)
    error('No AMICA solution found. You should first load AMICA components');
end
com = '';
if nargin<2
    plotopts = 'Smoothed probabilities|Unsmoothed Probabilities';
    guititle = 'Plot Model Probabilities -- pop_modprobplot()';
    if isfield(EEG.etc.amica,'smooth_length')
        smoothlength = num2str(EEG.etc.amica.smooth_length);
    else
        smoothlength = '2';
    end
    comsmooth = ['if get(gcbo,''value'') == 2, set(findobj(''Tag'',''smoothlengthTag''),''enable'',''off'');'...
        'set(findobj(''Tag'',''saveTag''),''enable'',''off'');else, set(findobj(''Tag'',''smoothlengthTag''),''enable'',''on'');' ...
        'set(findobj(''Tag'',''saveTag''),''enable'',''on'');end;'];

    comtype = [ 'if ~isfield(EEG.event, ''type'')' ...
        '   errordlg2(''No type field'');' ...
        'else' ...
        '   if isnumeric(EEG.event(1).type),' ...
        '        [tmps,tmpstr] = pop_chansel(unique([ EEG.event.type ]));' ...
        '   else,' ...
        '        [tmps,tmpstr] = pop_chansel(unique({ EEG.event.type }));' ...
        '   end;' ...
        '   if ~isempty(tmps)' ...
        '       set(findobj(''parent'', gcbf, ''tag'', ''eventtype''), ''string'', tmpstr);' ...
        '   end;' ...
        'end;' ...
        'clear tmps tmpstr;' ];

    allmodels = ['1:' num2str(EEG.etc.amica.num_models)];

    uilist = {{'style' 'text' 'string' 'Probabilities of which models'} ...
        {'style' 'edit' 'string' allmodels } ...
        {'style' 'text' 'string' 'Plot Option: '} ...
        {'style' 'popupmenu' 'string' plotopts 'value' 1 'callback' comsmooth} ...
        {'style' 'text' 'string' 'Smoothing Hanning Window Length (sec)'} ...
        {'style' 'edit' 'tag' 'smoothlengthTag' 'string' smoothlength 'enable' 'on'} ...
        {'style' 'checkbox' 'string' 'Save new smoothed probabilites as EEG.etc.amica.v_smooth' 'value' 0 'tag' 'saveTag'} ...
        {'style' 'pushbutton' 'string' 'Mark events of type..' 'callback' comtype} ...
        {'style' 'edit' 'string' '' 'tag' 'eventtype'}};
    uigeom = {[1 1] [1 1] [1 1] [1] [1 1]};

    result = inputgui(uigeom, uilist, 'pophelp(''pop_modprobplotl'')', guititle, [], 'normal');
    if isempty(result)
        return
    end

    models2plot = eval(['[' result{1} ']']);
    if min(models2plot)<1 && max(models2plot) > EEG.etc.amica.num_models
        error('Models chosen are not in the proper range');
    end

    smoothed = (result{2} == 1);
    store = result{4};
    events = result{5};

    if smoothed
        smoothlength = str2num(result{3});
    end

    if ~isempty(events)
        events2plot = parsetxt(events);

        if isnumeric(EEG.event(1).type)
            eeg_events = unique([EEG.event.type]);
            for i = 1:length(events2plot)
                if isempty(find(eeg_events == str2num(events2plot{i}), 1))
                    disp(['No event called ''' events2plot{i} '''' 'is found, ignoring that event'])
                    events2plot(i) = [];
                end
            end
        else
            eeg_events = unique({EEG.event.type});
            events2plot(ismember(events2plot,eeg_events) == 0) = [];

        end
    else
        events2plot = {};
    end
end

if ~smoothed
    smoothlength = 0;
end

if ~isfield(EEG.etc.amica,'LLt') || size(EEG.data,2) ~= size(EEG.etc.amica.LLt,2) || size(EEG.data,3) ~= size(EEG.etc.amica.LLt,3)
    recompute_prob = questdlg('Data size and AMICA probabilities do not match. Would you like to recompute the probabilities?', ...
        'Recompute probabilities?','Yes','No','No');
    if strcmpi(recompute_prob,'Yes')
        EEG.etc.amica.LLt = findLLt(EEG,EEG.etc.amica);
        EEG.etc.amica.v = LLt2v(EEG.etc.amica.LLt);
        [~,EEG.etc.amica.LLt_smooth] = smooth_amica_prob(EEG.srate,EEG.etc.amica.LLt,smoothlength);
        EEG = eeg_checkamica(EEG);
    end
end

[v2plot llt2plot] = modprobplot(EEG,models2plot,smoothlength,events2plot);

if store && smoothed
    EEG.etc.amica.v_smooth = v2plot;
    EEG.etc.amica.LLt_smooth = llt2plot;
    EEG.etc.amica.smooth_length = smoothlength;
end

com = sprintf('%s = pop_modprobplot(%s,%s,%d,%d,%d,%s)',inputname(1),inputname(1),vararg2str(models2plot),smoothed,smoothlength,store,vararg2str(events2plot));

author: PatternRecognition

repo: OpenBMI


text:
% PERTURBC - Perturbs the current solution to the solution valid for the
%            given regularization parameters.
%
% Syntax: [a,b,g,ind,X_mer,y_mer,Rs,Q] = perturbc(C)
%
%      a: alpha coefficients
%      b: bias
%      g: partial derivatives of cost function w.r.t. alpha coefficients
%    ind: cell array containing indices of margin, error and reserve vectors
%         ind{1}: indices of margin vectors
%         ind{2}: indices of error vectors
%         ind{3}: indices of reserve vectors
%  X_mer: matrix of margin, error and reserve vectors stored columnwise
%  y_mer: column vector of class labels (-1/+1) for margin, error and reserve vectors
%     Rs: inverse of extended kernel matrix for margin vectors
%      Q: extended kernel matrix for all vectors
%      C: soft-margin regularization parameter(s)
%         dimensionality of C       assumption
%         1-dimensional vector      universal regularization parameter
%         2-dimensional vector      class-conditional regularization parameters (-1/+1)
%         n-dimensional vector      regularization parameter per example
%         (where n = # of examples)
%
% Version 3.22e -- Comments to diehl@alumni.cmu.edu
%

function [a,b,g,ind,X,y,Rs,Q] = perturbc(C_new)

% flags for example state
MARGIN    = 1;
ERROR     = 2;
RESERVE   = 3;
UNLEARNED = 4;

% define global variables
global a;               % alpha coefficients
global b;               % bias
global C;               % regularization parameters
global deps;            % jitter factor in kernel matrix
global g;               % partial derivatives of cost function w.r.t. alpha coefficients
global ind;             % cell array containing indices of margin, error, reserve and unlearned vectors
global perturbations;   % number of perturbations
global Q;               % extended kernel matrix for all vectors
global Rs;              % inverse of extended kernel matrix for margin vectors
global scale;           % kernel scale
global type;            % kernel type
global X;               % matrix of margin, error, reserve and unlearned vectors stored columnwise
global y;               % column vector of class labels (-1/+1) for margin, error, reserve and unlearned vectors

kernel_evals_begin = kevals;

% create a vector containing the regularization parameter
% for each example if necessary
if (length(C_new) == 1)              % same regularization parameter for all examples
   C_new = C_new*ones(size(y));
elseif (length(C_new) == 2)          % class-conditional regularization parameters
   flags = (y == -1);
   C_new = C_new(1)*flags + C_new(2)*(~flags);
end;

% compute the regularization sensitivities
lambda = C_new-C;

% if there are no error vectors initially...
if (length(ind{ERROR}) == 0)

   % find all the examples that have changing regularization parameters
   inde = find(lambda ~= 0);

   % find the subset of the above examples that could become error vectors
   delta_p = (a(inde)-C(inde))./lambda(inde);
   i = find(delta_p > 0);

   % determine the minimum acceptable change in p and adjust the regularization parameters
   p = min([delta_p(i) ; 1]);
   C = C + lambda*p;

   % if one example becomes an error vector, perform the necessary bookkeeping
   if (p < 1)
      i = find(delta_p == p);
      indco = bookkeeping(inde(i),MARGIN,ERROR);
      updateRQ(indco);
   end;

else
   p = 0;
end;

% if there are error vectors to adjust...
if (p < 1)

   % compute sum{k in E} Qik lambda k and sum{k in E} yk lambda k
   SQl = ((y*y(ind{ERROR})').*kernel(X,X(:,ind{ERROR}),type,scale))*lambda(ind{ERROR});
   SQl(ind{ERROR}) = SQl(ind{ERROR}) + deps*lambda(ind{ERROR});
   Syl = y(ind{ERROR})'*lambda(ind{ERROR});

end;

s = sprintf('p = %.2f',p);
disp(s);

% change the regularization parameters incrementally
disp_p_delta = 0.2;
disp_p_count = 1;
num_MVs = length(ind{MARGIN});
perturbations = 0;
while (p < 1)

   perturbations = perturbations + 1;

   % compute beta and gamma
   if (num_MVs > 0)

      v = zeros(num_MVs+1,1);
      if (p < 1-eps)
         v(1) = -Syl - sum(y.*a)/(1-p);
      else
         v(1) = -Syl;
      end;
      v(2:num_MVs+1) = -SQl(ind{MARGIN});
      beta = Rs*v;
      gamma = zeros(size(Q,2),1);
      ind_temp = [ind{ERROR} ind{RESERVE} ind{UNLEARNED}];
      if (length(ind_temp) > 0)
         gamma(ind_temp) = Q(:,ind_temp)'*beta + SQl(ind_temp);
      end;

   else

      beta = 0;
      gamma = SQl;

   end;

   % minimum acceptable parameter change
   [min_delta_p,indss,cstatus,nstatus] = min_delta_p_c(p,gamma,beta,lambda);

   % update a, b, g and p
   if (length(ind{ERROR}) > 0)
      a(ind{ERROR}) = a(ind{ERROR}) + lambda(ind{ERROR})*min_delta_p;
   end;
   if (num_MVs > 0)
      a(ind{MARGIN}) = a(ind{MARGIN}) + beta(2:num_MVs+1)*min_delta_p;
   end;
   b = b + beta(1)*min_delta_p;
   g = g + gamma*min_delta_p;
   p = p + min_delta_p;
   C = C + lambda*min_delta_p;

   % perform bookkeeping
   indco = bookkeeping(indss,cstatus,nstatus);

   % update SQl and Syl when the status of indss changes from MARGIN to ERROR
   if ((cstatus == MARGIN) & (nstatus == ERROR))
      SQl = SQl + Q(indco,:)'*lambda(indss);
      Syl = Syl + y(indss)*lambda(indss);
   end;

   % set g(ind{MARGIN}) to zero
   g(ind{MARGIN}) = 0;

   % update Rs and Q if necessary
   if (nstatus == MARGIN)

      num_MVs = num_MVs + 1;
      if (num_MVs > 1)

         % compute beta and gamma for indss
         beta = -Rs*Q(:,indss);
         gamma = kernel(X(:,indss),X(:,indss),type,scale) + deps + Q(:,indss)'*beta;

      end;

      % expand Rs and Q
      updateRQ(beta,gamma,indss);

   elseif (cstatus == MARGIN)

      % compress Rs and Q
      num_MVs = num_MVs - 1;
      updateRQ(indco);

   end;

   % update SQl and Syl when the status of indss changes from ERROR to MARGIN
   if ((cstatus == ERROR) & (nstatus == MARGIN))
      SQl = SQl - Q(num_MVs+1,:)'*lambda(indss);
      Syl = Syl - y(indss)*lambda(indss);
   end;

   if (p >= disp_p_delta*disp_p_count)
      disp_p_count = disp_p_count + 1;
      s = sprintf('p = %.2f',p);
      disp(s);
   end;

end;
disp('Perturbation complete!');

% summary statistics
s = sprintf('\nMargin vectors:\t\t%d',length(ind{MARGIN}));
disp(s);
s = sprintf('Error vectors:\t\t%d',length(ind{ERROR}));
disp(s);
s = sprintf('Reserve vectors:\t%d',length(ind{RESERVE}));
disp(s);
s = sprintf('Kernel evaluations:\t%d\n',-kernel_evals_begin+kevals);
disp(s);

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%BORDERINTERPOLATE  Computes the source location of an extrapolated pixel
%
%     loc = cv.borderInterpolate(p, len)
%     [...] = cv.borderInterpolate(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __p__ 0-based coordinate of the extrapolated pixel along one of the axes,
%   likely `p < 0` or `p >= len`.
% * __len__ Length of the array along the corresponding axis.
%
% ## Output
% * __loc__ Pixel coordinate along the specified axis.
%
% ## Options
% * __BorderType__ Border type, one of:
%   * __Constant__ `iiiiii|abcdefgh|iiiiiii` with some specified `i`
%   * __Replicate__ `aaaaaa|abcdefgh|hhhhhhh`
%   * __Reflect__ `fedcba|abcdefgh|hgfedcb`
%   * __Reflect101__ `gfedcb|abcdefgh|gfedcba`
%   * __Wrap__ `cdefgh|abcdefgh|abcdefg`
%   * __Default__ same as 'Reflect101' (default)
%
%   When `BorderType=='Constant'`, the function always returns -1, regardless
%   of `p` and `len`.
%
% The function computes and returns the coordinate of a donor pixel
% corresponding to the specified extrapolated pixel when using the specified
% extrapolation border mode.
%
% Normally, the function is not called directly. It is used inside filtering
% functions and also in cv.copyMakeBorder.
%
% See also: cv.copyMakeBorder
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function numVax = freadVAXG(fid, number, method)
% FREADVAXG Converts a number read in IEEE-LE format to its VAXG
% floating point representation.
%
% Usage:
%    fid = fopen('file.vaxg', 'r', 'ieee-le');
%    num = freadVAXG(fid, numElements, 'format')
%
% 'format' options are the same as FREAD.
%
%  The function is intended to be called by the user.
%
%   See also UINT32LE_TO_VAXF, UINT64LE_TO_VAXG, FREADVAXD

%   Copyright 2009 The MathWorks, Inc.

switch method
    case {'float32', 'single'}
      rawUINT32 = fread(fid,number,'uint32=>uint32');
      numVax = uint32le_to_VAXF( rawUINT32 );

    case {'float64', 'double'}
      rawUINT32 = fread(fid,2*number,'uint32=>uint32');%read 2 32bit numbers
      numVax = uint64le_to_VAXG(rawUINT32);

    case {'float'}
      if intmax == 2147483647 %32bit OS float is 32 bits
	     rawUINT32 = fread(fid,  number,'uint32=>uint32');
         numVax = uint32le_to_VAXF(rawUINT32);
      else
         rawUINT32 = fread(fid,2*number,'uint32=>uint32');%read 2 32bit numbers
         numVax = uint64le_to_VAXG(rawUINT32);
      end

    otherwise
      numVax    = fread(fid, number, method);

end

end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function mtrPlotCorr(ccFilenameBase,threshVec,paramNames,midP,whichParams,strLineProp)
%
% paramNames = {'kLength','kSmooth','kMidSD'};
% midP = [0 18 0.175];

for pp = 1:length(whichParams)
    strThreshVec = {};
    parVecs = [];
    corrVecs = [];
    ovrVecs = [];
    for tt = 1:length(threshVec)
        cc = load([ccFilenameBase '_thresh_' num2str(threshVec(tt)) '.mat']);
        ccGrid = mtrCCMatrix2Grid(cc.ccMatrix,cc.paramData,paramNames);
        indGrid = ones(size(ccGrid,1),1);
        for ss = 1:length(paramNames)
            if ss ~= whichParams(pp)
                indGrid = indGrid & ccGrid(:,ss) == midP(ss);
            end
        end
        subGrid = ccGrid(indGrid,:);
        [foo, sortI] = sort(subGrid(:,whichParams(pp)));
        parVecs(:,tt) = subGrid(sortI(:),whichParams(pp));
        corrVecs(:,tt) = subGrid(sortI(:),4);
        strThreshVec{tt} = ['Top ' num2str(threshVec(tt))];
    end
    if( pp>1 ) figure; end
    if strcmp(paramNames{whichParams(pp)},'kSmooth')
        parVecs = asin(sqrt(1./parVecs)).*180./pi;
        plot(parVecs,corrVecs,strLineProp);
    else
        plot(parVecs,corrVecs,strLineProp);
    end
    %legend(strThreshVec);
    xlabel([paramNames{whichParams(pp)} ' parameter']);
    ylabel('Corr. Coef.');
end


text:
function [VAR, VARopt] = VARmodel(ENDO,nlag,const,EXOG,nlag_ex)
% =======================================================================
% Perform vector autogressive (VAR) estimation with OLS
% =======================================================================
% [VAR, VARopt] = VARmodel(ENDO,nlag,const,EXOG,nlag_ex)
% -----------------------------------------------------------------------
% INPUT
%	- ENDO: an (nobs x nvar) matrix of y-vectors
%	- nlag: lag length
% -----------------------------------------------------------------------
% OPTIONAL INPUT
%	- const: 0 no constant; 1 constant; 2 constant and trend; 3 constant,
%       trend, and trend^2 [dflt = 0]
%	- EXOG: optional matrix of variables (nobs x nvar_ex)
%	- nlag_ex: number of lags for exogeonus variables [dflt = 0]
% -----------------------------------------------------------------------
% OUTPUT
%   - VAR: structure including VAR estimation results
%   - VARopt: structure including VAR options (see VARoption)
% =======================================================================
% Ambrogio Cesa Bianchi, March 2015
% ambrogio.cesabianchi@gmail.com

% Note: this code is a modified version of of the vare.m function of James
% P. LeSage

% Representation -->  Y = Y(-1)*F' + u

% Note: compared to Eviews, there is a difference in the estimation of the
% constant when lag is > 2. This is because Eviews initialize the trend
% with the number of lags (i.e., when lag=2, the trend is [2 3 ...T]),
% while VARmakexy.m initialize the trend always with 1.

% I thank Jan Capek for spotting and addressing a compatibility issue with
% Matlab R2014a

%% Check inputs
%===============================================
[nobs, nvar] = size(ENDO);

% Create VARopt and update it
VARopt = VARoption;
VAR.ENDO = ENDO;
VAR.nlag = nlag;

% Check if ther are constant, trend, both, or none
if ~exist('const','var')
    const = 1;
end
VAR.const = const;

% Check if there are exogenous variables
if exist('EXOG','var')
    [nobs2, nvar_ex] = size(EXOG);
    % Check that ENDO and EXOG are conformable
    if (nobs2 ~= nobs)
        error('var: nobs in EXOG-matrix not the same as y-matrix');
    end
    clear nobs2
    % Check if there is lag order of EXOG, otherwise set it to 0
    if ~exist('nlag_ex','var')
        nlag_ex = 0;
    end
    VAR.EXOG = EXOG;
else
    nvar_ex = 0;
    nlag_ex = 0;
    VAR.EXOG = [];
end

%% Save some parameters and create data matrices
%===============================================
    nobse         = nobs - max(nlag,nlag_ex);
    VAR.nobs      = nobse;
    VAR.nvar      = nvar;
    VAR.nvar_ex   = nvar_ex;
    VAR.nlag      = nlag;
    VAR.nlag_ex   = nlag_ex;
    ncoeff        = nvar*nlag;
    VAR.ncoeff    = ncoeff;
    ncoeff_ex     = nvar_ex*(nlag_ex+1);
    ntotcoeff     = ncoeff + ncoeff_ex + const;
    VAR.ntotcoeff = ntotcoeff;
    VAR.const     = const;

% Create independent vector and lagged dependent matrix
[Y, X] = VARmakexy(ENDO,nlag,const);

% Create (lagged) exogeanous matrix
if nvar_ex>0
    X_EX  = VARmakelags(EXOG,nlag_ex);
    if nlag == nlag_ex
        X = [X X_EX];
    elseif nlag > nlag_ex
        diff = nlag - nlag_ex;
        X_EX = X_EX(diff+1:end,:);
        X = [X X_EX];
    elseif nlag < nlag_ex
        diff = nlag_ex - nlag;
        Y = Y(diff+1:end,:);
        X = [X(diff+1:end,:) X_EX];
    end
end

%% OLS estimation equation by equation
%===============================================

for j=1:nvar;
    Yvec = Y(:,j);
    OLSout = OLSmodel(Yvec,X,0);
    aux = ['eq' num2str(j)];
    eval( ['VAR.' aux '.beta  = OLSout.beta;'] );  % bhats
    eval( ['VAR.' aux '.tstat = OLSout.tstat;'] ); % t-stats
    % compute t-probs
    tstat = zeros(ncoeff,1);
    tstat = OLSout.tstat;
    tout = tdis_prb(tstat,nobse-ncoeff);
    eval( ['VAR.' aux '.tprob = tout;'] );        % t-probs
    eval( ['VAR.' aux '.resid = OLSout.resid;'] );% resids
    eval( ['VAR.' aux '.yhat  = OLSout.yhat;'] ); % yhats
    eval( ['VAR.' aux '.y     = Yvec;'] );        % actual y
    eval( ['VAR.' aux '.rsqr  = OLSout.rsqr;'] ); % r-squared
    eval( ['VAR.' aux '.rbar  = OLSout.rbar;'] ); % r-adjusted
    eval( ['VAR.' aux '.sige  = OLSout.sige;'] ); % standard error
end

%% Compute the matrix of coefficients & VCV
%===============================================
Ft = (X'*X)\(X'*Y);
VAR.Ft = Ft;
SIGMA = (1/(nobse-ntotcoeff))*(Y-X*Ft)'*(Y-X*Ft); % adjusted for # of estimated coeff per equation
VAR.sigma = SIGMA;
VAR.residuals = Y - X*Ft;
VAR.X = X;
VAR.Y = Y;
if nvar_ex > 0
    VAR.X_EX = X_EX;
end

%% Companion matrix of Ft' and max eigenvalue
%===============================================
F = Ft';
Fcomp = [F(:,1+const:nvar*nlag+const); eye(nvar*(nlag-1)) zeros(nvar*(nlag-1),nvar)];
VAR.Fcomp = Fcomp;
VAR.maxEig = max(abs(eig(Fcomp)));

%% Initialize other results
%===============================================
VAR.invA = [];  % inverse of teh A matrix (need identification: see VARir/VARfevd)
VAR.S    = [];  % Orthonormal matrix (need identification: see SR)

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
% PCA with automatic selection of the method: covariance or gram matrix
% Usage: [X, eigvec, eigval, Xm] = pca (X, dout, center, verbose)
%   X       input vector set (1 vector per column)
%   dout    number of principal components to be computed
%   center  need to center data?
%
% Note: the eigenvalues are given in decreasing order of magnitude
%
% Author: Herve Jegou, 2011.
% Last revision: 08/10/2013
function [X, eigvec, eigval, Xm] = yael_pca (X, dout, center, verbose)

if nargin < 3,         center = true; end
if ~exist ('verbose', 'var'), verbose = false; end

X = double (X);
d = size (X, 1);
n = size (X, 2);

if nargin < 2
  dout = d;
end

if center
  Xm = mean (X, 2);
  X = bsxfun (@minus, X, Xm);
else
  Xm = zeros (d, 1);
end

opts.issym = true;
opts.isreal = true;
opts.tol = eps;
opts.disp = 0;

% PCA with covariance matrix
if n > d
  if verbose, fprintf ('PCA with covariance matrix: %d -> %d\n', d, dout); end
  Xcov = X * X';
  Xcov = (Xcov + Xcov') / (2 * n);

  if dout < d
    [eigvec, eigval] = eigs (Xcov, dout, 'LM', opts);
  else
    [eigvec, eigval] = eig (Xcov);
  end
else
  % PCA with gram matrix
  if verbose, fprintf ('PCA with gram matrix: %d -> %d\n', d, dout); end
  Xgram = X' * X;
  Xgram = (Xgram + Xgram') / 2;
  if dout < d
    [eigvec, eigval] = eigs (Xgram, dout, 'LM', opts);
  else
    [eigvec, eigval] = eig (Xgram);
  end
  eigvec = single (X * eigvec);
  eigvec = utls.yael_vecs_normalize (eigvec);
end

X = eigvec' * X;
X = single (X);
eigval = diag(eigval);

% We prefer a consistent order
[~, eigord] = sort (eigval, 'descend');
eigval = eigval (eigord);
eigvec = eigvec (:, eigord);
X = X(eigord, :);

author: hpatches

repo: hpatches-benchmark


text:
%% ggremesh
% Below is a demonstration of the features of the |ggremesh| function

%% Syntax
% |[Fn,Vn]=ggremesh(F,V,optionStruct)|

%% Description
% This function uses the external library Geogram to remesh the input
% triangulation defined by the faces F and the vertices V. In particular
% the code "vorpalite" is used. An additional option structure may be
% provided where users can set particular parameters for Geogram.
%
% Below the options and defaults are provided:
% optionStruct.nb_pts=size(V,1); %number of points
% optionStruct.anisotropy=0; %Use anisotropy (~=0) to capture geometry or favour isotropic triangles (=0)
% optionStruct.pre.max_hole_area=100; %Max hole area for pre-processing step
% optionStruct.pre.max_hole_edges=0; %Max number of hole edges for pre-processing step
% optionStruct.post.max_hole_area=100; %Max hole area for post-processing step
% optionStruct.post.max_hole_edges=0; %Max number of hole edges for post-processing step
% optionStruct.disp_on=1; %Turn on/off displaying of Geogram text
%
% Instead of nb_pts users can also specify a pointSpacing to be used
% instead of nb_pts. This is not a Geogram feature but a GIBBON option
% which is translated to the number of points for Geogram remeshing. This
% is and example for a desired point spacing of 4:
% optionStruct.pointSpacing=4
%
% Geogram website:
% http://alice.loria.fr/index.php/software/4-library/75-geogram.html
%
% Geogram license:
% http://alice.loria.fr/software/geogram/doc/html/geogram_license.html
%
% Lévy B., Bonneel N. (2013) Variational Anisotropic Surface Meshing with
% Voronoi Parallel Linear Enumeration. In: Jiao X., Weill JC. (eds)
% Proceedings of the 21st International Meshing Roundtable. Springer,
% Berlin, Heidelberg. https://doi.org/10.1007/978-3-642-33573-0_21
%
% See also:
% http://alice.loria.fr/publications/papers/2012/Vorpaline_IMR/vorpaline.pdf
% https://www.ljll.math.upmc.fr/hecht/ftp/ff++days/2013/BrunoLevy.pdf

%% Examples

%%
clear; close all; clc;

%%
% Plot settings
fontSize=15;
faceColor='b';
faceAlpha=1;
edgeColor='k';
edgeWidth=0.5;

%% Example 1: Remeshing a triangulated surface isotropically

%%
% Get example geometry
[F,V]=graphicsModels(5); % Get surface

%%
% Remesh using ggremesh

[Fn,Vn]=ggremesh(F,V);

%%
% Visualiza patch data

cFigure;
subplot(1,2,1); hold on;
title('Input mesh');
gpatch(F,V,'w','k');
axisGeom;
view(-75,-36);
camlight headlight; axis off;

subplot(1,2,2); hold on;
title('Geogram remeshed');
gpatch(Fn,Vn,'gw','k',1,1);
axisGeom;
view(-75,-36);
camlight headlight; axis off;

gdrawnow;

%% Example 2: Remeshing a triangulated surface with desired number of points

%%
% Get example geometry
[F,V]=graphicsModels(1); % Get surface

%%
% Remesh using ggremesh

optionStruct1.nb_pts=500; %Set desired number of points
[Fn,Vn]=ggremesh(F,V,optionStruct1);

%%
% Visualiza patch data

cFigure;
subplot(1,2,1); hold on;
title('Input mesh');
gpatch(F,V,'w','k');
axisGeom;
camlight headlight; axis off;

subplot(1,2,2); hold on;
title('Geogram remeshed');
gpatch(Fn,Vn,'gw','k',1,1);
axisGeom;
camlight headlight; axis off;

gdrawnow;

%% Example 3: Remeshing a triangulated surface with desired point spacing

%%
% Get example geometry
[F,V]=graphicsModels(11); % Get surface

%%
% Remesh using ggremesh

optionStruct2.pointSpacing=4; %Set desired point spacing
optionStruct2.disp_on=1; % Turn off command window text display
[Fn,Vn]=ggremesh(F,V,optionStruct2);

%%
% Visualiza patch data

cFigure;
subplot(1,2,1); hold on;
title('Input mesh');
gpatch(F,V,'w','k');
axisGeom;
camlight headlight;

subplot(1,2,2); hold on;
title('Geogram remeshed');
gpatch(Fn,Vn,'gw','k',1,1);
axisGeom;
camlight headlight;

gdrawnow;

%% Example 4: Setting pre- and prost-processing settings e.g. to close holes

%%
% Get example geometry

inputStruct.cylRadius=1;
inputStruct.numRadial=15;
inputStruct.cylHeight=3;
inputStruct.numHeight=11;
inputStruct.meshType='tri';

% Derive patch data for a cylinder
[F,V]=patchcylinder(inputStruct);

%%
% Remesh using ggremesh
optionStruct3.nb_pts=size(V,1); %Set desired number of points
optionStruct3.disp_on=0; % Turn off command window text display
optionStruct3.pre.max_hole_area=10; %Max hole area for pre-processing step
optionStruct3.pre.max_hole_edges=20; %Max number of hole edges for pre-processing step
% optionStruct3.post.max_hole_area=10; %Max hole area for post-processing step
% optionStruct3.post.max_hole_edges=20; %Max number of hole edges for post-processing step

[Fn,Vn]=ggremesh(F,V,optionStruct3);

% Visualiza patch data
Eb=patchBoundary(F);
cFigure;
subplot(1,2,1); hold on;
title('Input mesh with holes');
gpatch(F,V,'w','k');
gpatch(Eb,V,'none','b',1,2);
axisGeom;
camlight headlight;

subplot(1,2,2); hold on;
title('Geogram remeshed and closed');
gpatch(Fn,Vn,'gw','k',1,1);
axisGeom;
camlight headlight;

gdrawnow;

%% Example 5: Setting pre- and prost-processing settings e.g. to avoid closure of holes

%%
% Get example geometry

inputStruct.cylRadius=1;
inputStruct.numRadial=15;
inputStruct.cylHeight=3;
inputStruct.numHeight=11;
inputStruct.meshType='tri';

% Derive patch data for a cylinder
[F,V]=patchcylinder(inputStruct);

%%
% Remesh using ggremesh
optionStruct3.nb_pts=size(V,1); %Set desired number of points
optionStruct3.disp_on=1; % Turn off command window text display
optionStruct3.pre.max_hole_area=100; %Max hole area for pre-processing step
optionStruct3.pre.max_hole_edges=0; %Max number of hole edges for pre-processing step

[Fn,Vn]=ggremesh(F,V,optionStruct3);

% Visualiza patch data
Eb=patchBoundary(F);
cFigure;
subplot(1,2,1); hold on;
title('Input mesh with holes');
gpatch(F,V,'w','k');
gpatch(Eb,V,'none','b',1,2);
axisGeom;
camlight headlight;

subplot(1,2,2); hold on;
title('Geogram remeshed with holes');
gpatch(Fn,Vn,'gw','k',1,1);
axisGeom;
camlight headlight;

gdrawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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 3 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, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [h] = fscatter3(X,Y,Z,C,cmap);
% [h] = fscatter3(X,Y,Z,C,cmap);
% Plots point cloud data in cmap color classes and 3 Dimensions,
% much faster and very little memory usage compared to scatter3 !
% X,Y,Z,C are vectors of the same length
% X,Y,Z,C might be put in as structure points.x,points.y,points.z,points.int
% with C being used as index into colormap (can be any values though)
% cmap is optional colourmap to be used
% h are handles to the line objects

% Felix Morsdorf, Jan 2003 (last update Oct. 2010), Remote Sensing Laboratory Zuerich

  if nargin == 1
    if isfield(X,'int')
      C = X.int;
      Z = X.z;
    else
      C = X.z;
      Z = C;
    end
    Y = X.y;
    NX = X.x;
    clear X;
    X = NX;clear NX;
    numclass = 256; % Number of color classes
    cmap = myspecmap(256);
    siz = 5;
  elseif nargin == 4
    numclass = 256; % Number of color classes
    cmap = hsv(256);
    siz = 5;
  elseif nargin == 5
    numclass = max(size(cmap));
    siz = 5;
    if numclass == 1
      siz = cmap;
      cmap = hsv(256);
      numclass = 256;
    end
  elseif nargin == 6
    numclass = max(size(cmap));
    if numclass == 1
      siz = cmap;
      cmap = hsv(256);
      numclass = 256;
    end
  end

% avoid too many calculations

mins = min(C);
maxs = max(C);
minz = min(Z);
maxz = max(Z);
minx = min(X);
maxx = max(X);
miny = min(Y);
maxy = max(Y);

% construct colormap :

col = cmap;

% determine index into colormap

ii = floor( (C - mins ) * (numclass-1) / (maxs - mins) );
ii = ii + 1;

colormap(cmap);

hold on
k = 0;o = k;
for j = 1:numclass
  jj = (ii(:)== j);
  if ~isempty(jj)
    k = k + 1;
    h = plot3(X(jj),Y(jj),Z(jj),'.','color',col(j,:),'markersize',siz);
    if ~isempty(h)
      o = o+1;
        hp(o) = h;
    end
  end
end
caxis([min(C) max(C)])
axis equal;rotate3d on;view(3);
box on
hcb = colorbar('location','east');

author: Sable

repo: mcbench-benchmarks


text:
function output = callpennlp(interfacedata)

% Retrieve needed data
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
x0      = interfacedata.x0;
Q       = interfacedata.Q;
lb      = interfacedata.lb;
ub      = interfacedata.ub;
monomtable = interfacedata.monomtable;

switch options.verbose
    case 0
        options.fmincon.Display = 'off';
    case 1
        options.fmincon.Display = 'final';
    otherwise
        options.fmincon.Display = 'iter';
end

% Do some pre-calc to be used in calls from pennlp
temp = sum(interfacedata.monomtable,2)>1;
temp = (sum(interfacedata.monomtable,2)==1 & sum(interfacedata.monomtable~=0,2)==1);
nonlinearindicies = find(full(~temp));
linearindicies = find(full(temp));
interfacedata.nonlinearindicies = nonlinearindicies;
interfacedata.linearindicies = linearindicies;

Infinity = 1.0E38;
if isempty(lb)
    lb = -Infinity*ones(length(linearindicies),1);
else
    lb(isinf(lb)) =  -Infinity;
end
if isempty(ub)
    ub = Infinity*ones(length(linearindicies),1);
else
    ub(isinf(ub)) =  Infinity;
end

% Numer of constraints
m = K.l+K.f;
n = length(linearindicies);

lb = lb(linearindicies);
ub = ub(linearindicies);

lbc =  0*ones(m,1);
ubc =  Infinity*ones(m,1);
ubc(1:K.f) = 0;

if isempty(x0)
    x0 = zeros(length(linearindicies),1);
else
    x0 = x0(linearindicies);
end

linearConstraints = ~any(F_struc(1:K.f+K.l,1+nonlinearindicies),2);
if ~isempty(linearConstraints)
    F = F_struc(1:K.f+K.l,:);
    F_struc = [F(~linearConstraints,:);F(linearConstraints,:);F_struc(K.f+K.l+1:end,:)];
    lbc = [lbc(~linearConstraints);lbc(linearConstraints)];
    ubc = [ubc(~linearConstraints);ubc(linearConstraints)];
end

functiondata.nonlinearindicies = nonlinearindicies;
functiondata.linearindicies = linearindicies;
functiondata.monomtable = monomtable;
functionGradientdata.nonlinearindicies = nonlinearindicies;
functionGradientdata.linearindicies = linearindicies;
functionGradientdata.monomtable = monomtable;
functionHessiandata.nonlinearindicies = nonlinearindicies;
functionHessiandata.linearindicies = linearindicies;
functionHessiandata.monomtable = monomtable;

constraintdata.nonlinearindicies = nonlinearindicies;
constraintdata.linearindicies = linearindicies;
constraintdata.monomtable = monomtable;
constraintGradientdata.nonlinearindicies = nonlinearindicies;
constraintGradientdata.linearindicies = linearindicies;
constraintGradientdata.monomtable = monomtable;
constraintHessiandata.nonlinearindicies = nonlinearindicies;
constraintHessiandata.linearindicies = linearindicies;
constraintHessiandata.monomtable = monomtable;

functiondata.F_struc = [interfacedata.f c(:)'];
constraintdata.F_struc = F_struc(1:K.l+K.f,:);start = K.f + K.l + 1;
functionGradientdata.F_struc   = F_struc(start:start+n-1,:);start = start+n;
functionHessiandata.F_struc    = F_struc(start:start+n^2-1,:);start = start+n^2;
constraintGradientdata.F_struc = F_struc(start:start+n*m-1,:);start = start+n*m;
constraintHessiandata.F_struc  = F_struc(start:end,:);start = start+n*m;

datasaver(functiondata,functionGradientdata,functionHessiandata,constraintdata,constraintGradientdata,constraintHessiandata)

Infinity = 1.0E38;
pen.nvars = n;
pen.nlin = length(find(linearConstraints));
pen.nconstr = m;
pen.nnz_gradient = n;
pen.nnz_hessian = (n+1)*n/2;
pen.lbv = lb(:)';
pen.ubv = ub(:)';
pen.lbc = lbc(:)';
pen.ubc = ubc(:)';
pen.xinit= x0(:)';
pen.my_f = 'pennlp_fun';
pen.my_f_gradient = 'pennlp_fungrad';
pen.my_f_hessian = 'pennlp_funhess';
pen.my_g = 'pennlp_con';
pen.my_g_gradient = 'pennlp_congrad';
pen.my_g_hessian =  'pennlp_conhess';

ops = struct2cell(options.pennlp);ops = [ops{1:end}];
pen.ioptions = [ops(1:2) options.verbose ops(3:15)];
if options.verbose>0
    pen.ioptions(3) = pen.ioptions(3)+1;
end
pen.doptions = ops(16:end);

%pen.ioptions = ops(1:15);
%pen.foptions = ops(16:end);
%pen.ioptions(4) = max(0,min(3,options.verbose+1));

%pen.ioptions = [100 100 2 0 0 0 1 0 0 2 0 0 0 -1 0 1];
%pen.doptions = [1.0E-7 1.0E0 1.0E-0 1.0E-2 5.0E-1 1.0E-1 1.0E-6 1.0E-12 1.0e-1 0.05 1.0 1.0 1.0];

if options.savedebug
    save pennlpdebug pen
end

showprogress('Calling PENNLP',options.showprogress);

solvertime = tic;
[obj,xout,duals,flag] = pennlpm(pen);
solvertime = toc(solvertime);

if isempty(nonlinearindicies)
    x = xout(:);
else
    x = zeros(length(interfacedata.c),1);
    for i = 1:length(linearindicies)
        x(linearindicies(i)) = xout(i);
    end
end

switch flag
    case 0
        problem = 0;
    case {1,6,7}
        problem = 4;
    otherwise
        problem = 9;
end

% Internal format for duals
D_struc = [];

% Save all data sent to solver?
if options.savesolverinput
    solverinput.pen = pen;
else
    solverinput = [];
end

% Save all data from the solver?
if options.savesolveroutput
    solveroutput.obj = obj;
    solveroutput.xout = xout;
    solveroutput.duals = duals;
    solveroutput.flag = flag;
else
    solveroutput = [];
end

% Standard interface
output = createOutputStructure(x,D_struc,[],problem,interfacedata.solver.tag,solverinput,solveroutput,solvertime);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
classdef Symmetrizer < handle

    properties (Access = public)
        symmetricMesh
    end

    properties (Access = private)
        symmetricPoint
        meshMerger
    end

    properties (Access = private)
        mesh
        symmetrizedMesh
        symmetricLine
    end

    methods (Access = public)

        function obj = Symmetrizer(cParams)
            obj.init(cParams);
            obj.createSymmetricPointComputer();
            obj.createSymmetrizedMesh();
            obj.createMeshMerger();
        end

        function m = computeSymmetricMesh(obj)
            m = obj.meshMerger.compute();
        end

        function sF = symmetrizeScalarField(obj,f)
            nodes = obj.meshMerger.computeRemainingNodes();
            fS = f;
            fD = [f;fS];
            sF   = fD(nodes,:);
        end

        function sF = symmetrizeVectorField(obj,f)
            nodes = obj.meshMerger.computeRemainingNodes();
            fS = obj.computeFsymmetric(f);
            fD   = [f;fS];
            sF   = fD(nodes,:);
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.mesh          = cParams.mesh;
            obj.symmetricLine = cParams.symmetricLine;
        end

        function fSym = computeFsymmetric(obj,f)
            s.vector = f;
            s.line  = obj.symmetricLine;
            s.line.point = [0,0];
            sP = SymmetricPointComputer(s);
            fSym = sP.computeSymmetricVector();
        end

        function createSymmetricPointComputer(obj)
            s.vector = obj.mesh.coord;
            s.line   = obj.symmetricLine;
            sP = SymmetricPointComputer(s);
            obj.symmetricPoint = sP;
        end

        function createSymmetrizedMesh(obj)
            sConnec(:,[1 3 2]) = obj.mesh.connec(:,[1 2 3]);
            sCoord = obj.symmetricPoint.computeSymmetricVector();
            s.connec = sConnec;
            s.coord  = sCoord;
            m = Mesh(s);
            obj.symmetrizedMesh = m;
        end

        function createMeshMerger(obj)
            isInSymLine = obj.symmetricPoint.isNodeInLine;
            s.meshA  = obj.mesh;
            s.meshB  = obj.symmetrizedMesh;
            s.isMergedNodeA = isInSymLine;
            s.isMergedNodeB = isInSymLine;
            obj.meshMerger = MeshMerger(s);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function tapas_uniqc_new_method(varargin)
% tapas_uniqc_new_method creates a new method assuming that you are in the corresponding class
% folder; creates header using a template.
% tapas_uniqc_new_method(funname) opens the editor and pastes the content
% of a user-defined template into the file funname.m.
%
%   Example
%       tapas_uniqc_new_method myfunfun
%           OR
%       tapas_uniqc_new_method myfunfun author
%
%   opens the editor and pastes the following
%
% 	function output = myfunfun(input)
% 	%MYFUNFUN  One-line description here, please.
% 	%   output = myfunfun(input)
% 	%
% 	%   Example
% 	%   myfunfun
% 	%
% 	%   See also
%
% 	% Author: Your name
% 	% Created: 2005-09-22
% 	% Copyright 2005 Your company.
%
%   See also edit, mfiletemplate

% Author: Peter (PB) Bodin
% Created: 2005-09-22
% Modified: 2014-04-15 (Saskia Klein and Lars Kasper, IBT Zurich)

% See the variables repstr, repwithstr and tmpl to figure out how
% to design your own template.
% Edit tmpl to your liking, if you add more tokens in tmpl, make
% sure to add them in repstr and repwithstr as well.

% I made this function just for fun to check out some java handles to
% the editor. It would probably be better to fprintf the template
% to a new file and then call edit, since the java objects might change
% names between versions.

switch nargin
    case 0
        edit
        warning('new_method without argument is the same as edit')
        return;
    case 1
        fname=varargin{1};
        edit(fullfile(pwd,fname));
        authors = 'Saskia Bollmann & Lars Kasper'; %defuaults authors, set down in function authors
    case 2
        fname = varargin{1};
        authors = varargin{2};
    otherwise
        error('tapas:uniqc:TooManyInputArguments', ...
            'too many input arguments')
end

% determine class name from directory name
[~, classname] = fileparts(pwd);
classname = regexprep(classname, '@', '');

try lasterror
    edhandle=com.mathworks.mlservices.MLEditorServices;

    % R2009a => 2009.0, R2009b = 2009.5
    vs = version('-release');
    v = str2double(vs(1:4));
    if vs(5)=='b'
        v = v + 0.5;
    end

    if v < 2009.0
        edhandle.builtinAppendDocumentText(strcat(fname,'.m'),...
            parse(fname,authors, classname));
    else
        edhandle.getEditorApplication.getActiveEditor.appendText(...
            parse(fname, authors, classname));
    end
catch
    rethrow(lasterror)
end

    function out = parse(func, authors, classname)

        tmpl={ ...
            'function this = $filename(this)'
            '%ONE_LINE_DESCRIPTION'
            '%'
            '%   Y = $classname()'
            '%   Y.$filename(inputs)'
            '%'
            '% This is a method of class $classname.'
            '%'
            '% IN'
            '%'
            '% OUT'
            '%'
            '% EXAMPLE'
            '%   $filename'
            '%'
            '%   See also $classname'
            ' '
            '% Author:   $author'
            '% Created:  $date'
            '% Copyright (C) $year $institute'
            '%                    $company'
            '%'
            '% This file is part of the TAPAS UniQC Toolbox, which is released'
            '% under the terms of the GNU General Public License (GPL), version 3. '
            '% You can redistribute it and/or modify it under the terms of the GPL'
            '% (either version 3 or, at your option, any later version).'
            '% For further details, see the file COPYING or'
            '%  <http://www.gnu.org/licenses/>.'
            };

        repstr={...
            '$classname'
            '$filename'
            '$FILENAME'
            '$date'
            '$year'
            '$author'
            '$institute'
            '$company'};

        repwithstr={...
            classname
            func
            upper(func)
            datestr(now,29)
            datestr(now,10)
            authors
            'Institute for Biomedical Engineering'
            'University of Zurich and ETH Zurich'};

        for k = 1:numel(repstr)
            tmpl = strrep(tmpl,repstr{k},repwithstr{k});
        end
        out = sprintf('%s\n',tmpl{:});
    end
end

author: translationalneuromodeling

repo: tapas


text:
function [mappedX, mapping] = kernel_pca(X, no_dims, varargin)
%KERNEL_PCA Perform the kernel PCA algorithm
%
%   [mappedX, mapping] = kernel_pca(X, no_dims)
%   [mappedX, mapping] = kernel_pca(X, no_dims, kernel)
%   [mappedX, mapping] = kernel_pca(X, no_dims, kernel, param1)
%   [mappedX, mapping] = kernel_pca(X, no_dims, kernel, param1, param2)
%
% The function runs kernel PCA on a set of datapoints X. The variable
% no_dims sets the number of dimensions of the feature points in the
% embedded feature space (no_dims >= 1, default = 2).
% For no_dims, you can also specify a number between 0 and 1, determining
% the amount of variance you want to retain in the PCA step.
% The value of kernel determines the used kernel. Possible values are 'linear',
% 'gauss', 'poly', 'subsets', or 'princ_angles' (default = 'gauss'). For
% more info on setting the parameters of the kernel function, type HELP
% GRAM.
% The function returns the locations of the embedded trainingdata in
% mappedX. Furthermore, it returns information on the mapping in mapping.
%
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

    if ~exist('no_dims', 'var')
        no_dims = 2;
    end
    kernel = 'gauss';
    param1 = 1;
	param2 = 3;
    if nargin > 2
		kernel = varargin{1};
		if length(varargin) > 1 & strcmp(class(varargin{2}), 'double'), param1 = varargin{2}; end
		if length(varargin) > 2 & strcmp(class(varargin{3}), 'double'), param2 = varargin{3}; end
    end

    % Store the number of training and test points
    ell = size(X, 1);

    if size(X, 1) < 2000

        % Compute Gram matrix for training points
        disp('Computing kernel matrix...');
        K = gram(X, X, kernel, param1, param2);

        % Normalize kernel matrix K
        mapping.column_sums = sum(K) / ell;                       % column sums
        mapping.total_sum   = sum(mapping.column_sums) / ell;     % total sum
        J = ones(ell, 1) * mapping.column_sums;                   % column sums (in matrix)
        K = K - J - J';
        K = K + mapping.total_sum;

        % Compute first no_dims eigenvectors and store these in V, store corresponding eigenvalues in L
        disp('Eigenanalysis of kernel matrix...');
        K(isnan(K)) = 0;
        K(isinf(K)) = 0;
        [V, L] = eig(K);
    else
        % Compute column sums (for out-of-sample extension)
        mapping.column_sums = kernel_function([], X', 1, kernel, param1, param2, 'ColumnSums') / ell;
        mapping.total_sum   = sum(mapping.column_sums) / ell;

        % Perform eigenanalysis of kernel matrix without explicitly
        % computing it
        disp('Eigenanalysis of kernel matrix (using slower but memory-conservative implementation)...');
        options.disp = 0;
        options.isreal = 1;
        options.issym = 1;
        [V, L] = eigs(@(v)kernel_function(v, X', 1, kernel, param1, param2, 'Normal'), size(X, 1), no_dims, 'LM', options);
        disp(' ');
    end

    % Sort eigenvalues and eigenvectors in descending order
    [L, ind] = sort(diag(L), 'descend');
    L = L(1:no_dims);
	V = V(:,ind(1:no_dims));

    % Compute inverse of eigenvalues matrix L
	disp('Computing final embedding...');
    invL = diag(1 ./ L);

    % Compute square root of eigenvalues matrix L
    sqrtL = diag(sqrt(L));

    % Compute inverse of square root of eigenvalues matrix L
    invsqrtL = diag(1 ./ diag(sqrtL));

    % Compute the new embedded points for both K and Ktest-data
    mappedX = sqrtL * V';                     % = invsqrtL * V'* K

    % Set feature vectors in original format
    mappedX = mappedX';

    % Store information for out-of-sample extension
    mapping.X = X;
    mapping.V = V;
    mapping.invsqrtL = invsqrtL;
    mapping.kernel = kernel;
    mapping.param1 = param1;
    mapping.param2 = param2;

repo: cluster


text:
% Test: Linear subdivision for triangle meshes
%
% Author: Jesus Mena

% Example: Box
vertices = [10 10 10; -10 10 10; 10 -10 10; -10 -10 10; 10 10 -10; -10 10 -10; 10 -10 -10; -10 -10 -10]';
faces = [1 2 3; 4 3 2; 1 3 5; 7 5 3; 1 5 2; 6 2 5; 8 6 7; 5 7 6; 8 7 4; 3 4 7; 8 4 6; 2 6 4]';

figure(1);
subplot(1,4,1);
plotMesh(vertices, faces);
for i=2:4
 subplot(1,4,i);
 [vertices, faces] = linearSubdivision(vertices, faces);
 plotMesh(vertices, faces);
end

% Example: Tetrahedron
vertices = [10 10 10; -100 10 -10; -100 -10 10; 10 -10 -10]';
faces = [1 2 3; 1 3 4; 1 4 2; 4 3 2]';

figure(2);
subplot(1,4,1);
plotMesh(vertices, faces);
for i=2:4
 subplot(1,4,i);
 [vertices, faces] = linearSubdivision(vertices, faces);
 plotMesh(vertices, faces);
end

% Example: Cylinder
vertices = [0 -5 0; 0 5 0; 10 -5 0; 9.65 -5 2.58; 8.66 -5 5; 7.07 -5 7.07; 5 -5 8.66; 2.58 -5 9.65; 0 -5 10; -2.58 -5 9.65; -5 -5 8.66; -7.07 -5 7.07; -8.66 -5 5; -9.65 -5 2.58; -10 -5 0; -9.65 -5 -2.58; -8.66 -5 -5; -7.07 -5 -7.07; -5 -5 -8.66; -2.58 -5 -9.65; -0 -5 -10; 2.58 -5 -9.65; 5 -5 -8.66; 7.07 -5 -7.07; 8.66 -5 -5; 9.65 -5 -2.58; 10 5 0 ; 9.65 5 2.58; 8.66 5 5; 7.07 5 7.07; 5 5 8.66; 2.58 5 9.65; 0 5 10; -2.58 5 9.65; -5 5 8.66; -7.07 5 7.07; -8.66 5 5; -9.65 5 2.58; -10 5 0; -9.65 5 -2.58; -8.66 5 -5; -7.07 5 -7.07; -5 5 -8.66; -2.58 5 -9.65; -0 5 -10; 2.58 5 -9.65; 5 5 -8.66; 7.07 5 -7.07; 8.66 5 -5; 9.65 5 -2.58]';
faces = [1 3 4; 1 4 5; 1 5 6; 1 6 7; 1 7 8; 1 8 9; 1 9 10; 1 10 11; 1 11 12; 1 12 13; 1 13 14; 1 14 15; 1 15 16; 1 16 17; 1 17 18; 1 18 19; 1 19 20; 1 20 21; 1 21 22; 1 22 23; 1 23 24; 1 24 25; 1 25 26; 1 26 3; 2 28 27; 2 29 28; 2 30 29; 2 31 30; 2 32 31; 2 33 32; 2 34 33; 2 35 34; 2 36 35; 2 37 36; 2 38 37; 2 39 38; 2 40 39; 2 41 40; 2 42 41; 2 43 42; 2 44 43; 2 45 44; 2 46 45; 2 47 46; 2 48 47; 2 49 48; 2 50 49; 2 27 50; 3 27 28; 3 28 4; 4 28 29; 4 29 5; 5 29 30; 5 30 6; 6 30 31; 6 31 7; 7 31 32; 7 32 8; 8 32 33; 8 33 9; 9 33 34; 9 34 10; 10 34 35; 10 35 11; 11 35 36; 11 36 12; 12 36 37; 12 37 13; 13 37 38; 13 38 14; 14 38 39; 14 39 15; 15 39 40; 15 40 16; 16 40 41; 16 41 17; 17 41 42; 17 42 18; 18 42 43; 18 43 19; 19 43 44; 19 44 20; 20 44 45; 20 45 21; 21 45 46; 21 46 22; 22 46 47; 22 47 23; 23 47 48; 23 48 24; 24 48 49; 24 49 25; 25 49 50; 25 50 26; 26 50 27; 26 27 3]';

figure(3);
subplot(1,4,1);
plotMesh(vertices, faces);
for i=2:4
 subplot(1,4,i);
 [vertices, faces] = linearSubdivision(vertices, faces);
 plotMesh(vertices, faces);
end

% Example: Grid
vertices = [-4 -4 0; -2 -4 0; 0 -4 0; 2 -4 0; 4 -4 0; -4 -2 0; -2 -2 0; 0 -2 0; 2 -2 0; 4 -2 0; -4 0 0; -2 0 0; 0 0 0; 2 0 0; 4 0 0; -4 2 0; -2 2 0; 0 2 0; 2 2 0; 4 2 0; -4 4 0; -2 4 0; 0 4 0; 2 4 0; 4 4 0]';
faces = [7 2 1; 1 6 7; 8 3 2; 2 7 8; 9 4 3; 3 8 9; 10 5 4; 4 9 10; 12 7 6; 6 11 12; 13 8 7; 7 12 13; 14 9 8; 8 13 14; 15 10 9; 9 14 15; 17 12 11; 11 16 17; 18 13 12; 12 17 18; 19 14 13; 13 18 19; 20 15 14; 14 19 20; 22 17 16; 16 21 22; 23 18 17; 17 22 23; 24 19 18; 18 23 24; 25 20 19; 19 24 25]';

figure(4);
subplot(1,4,1);
plotMesh(vertices, faces);
for i=2:4
 subplot(1,4,i);
 [vertices, faces] = linearSubdivision(vertices, faces);
 plotMesh(vertices, faces);
end


text:
function [adfstat,pval,critval,resid,lags,ICs]=augdfautolag(y,p,maxlags,IC)
% Dickey-Fuller and Augmented Dickey Fuller with automatic lag selection
%
% USAGE:
%  [ADFSTAT,PVAL,CRITVAL] = augdfautolag(Y,P,LAGS,IC)
%  [ADFSTAT,PVAL,CRITVAL,RESID,LAGS] = augdfautolag(Y,P,LAGS,IC)
%
% INPUTS:
%  Y         - A T by 1 vector of data
%  P         - Order of the polynomial of include in the ADF regression:
%                0 : No deterministic terms
%                1 : Constant
%                2 : Time Trend
%                3 : Constant, DGP assumed to have a time trend
%  MAXLAGS   - The maximum number of lags to include in the ADF test
%  IC        - [OPTIONAL] String, either 'AIC' (default) or 'BIC' to choose the criteria to select
%                the model
%
% OUTPUTS:
%  ADFSTAT   - Dickey-Fuller statistic
%  PVAL      - Probability the series is a unit root
%  CRITVALS  - A 6 by 1 vector with the [.01 .05 .1 .9 .95 .99] values from the DF distribution
%  LAGS      - The selected number of lags
%  IC        - The value at all lags of the selected IC
%
% COMMENTS:
%
% See also AUGDF

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3.0.1    Date: 1/1/2007

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<3 || nargin>4
    error('3 or 4 inputs required.')
end
T=length(y);
if T<=(maxlags+1)
    error('Length of data must be larger than LAGS')
end
if size(y,1)~=T,
    y=y';
end
if size(y,2)~=1
    error('Y must be a column vector')
end
if ~isscalar(maxlags) && maxlags>0 && floor(maxlags)==maxlags
    error('LAGS must be a positive integer')
end
if ~isscalar(p)
    error('P must be a scalar integer in {0, 1, 2, 3}')
end
if ~ismember(p,[0 1 2 3])
    error('P must be a scalar integer in {0, 1, 2, 3}')
end
if nargin==3
    IC='AIC';
else
    if ~ischar(IC)
        error('IC must be a string, either ''AIC'' or ''BIC''');
    elseif ~ismember(IC,{'AIC','BIC'})
        error('IC must be a string, either ''AIC'' or ''BIC''');
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Setup common to all problems
%y=y-y(1);
ydiff=diff(y);
[ydiffcurr, ydifflags]=newlagmatrix(ydiff,maxlags); %#ok<ASGLU>
T=length(y);
Y=y(maxlags+2:T);
tau=length(Y);
%
switch p
    case 0
        % Case 1
        i=0;
        X=y(maxlags+1:T-1);
        rho = X\Y;
        % Compute the errors
        e= Y-X*rho;
        s2(i+1)=e'*e/tau;
        K(i+1)=size(X,2);

        % Loop
        for i=1:maxlags
            X=[y(maxlags+1:T-1) ydifflags(:,1:i)];
            rho = X\Y;
            % Compute the errors
            e= Y-X*rho;
            s2(i+1)=e'*e/tau;
            K(i+1)=size(X,2);
        end

    case {1,3}
        %Case 2
        i=0;
        X=[ones(size(Y)) y(maxlags+1:T-1)];
        rho = X\Y;
        % Compute the errors
        e= Y-X*rho;
        s2(i+1)=e'*e/tau;
        K(i+1)=size(X,2);

        % Loop
        for i=1:maxlags
            X=[ones(size(Y)) y(maxlags+1:T-1) ydifflags(:,1:i)];
            rho = X\Y;
            % Compute the errors
            e= Y-X*rho;
            s2(i+1)=e'*e/tau;
            K(i+1)=size(X,2);
        end

    case 2
        %Case 4
        i=0;
        X=[ones(size(Y)) y(maxlags+1:T-1) (1:tau)'];
        rho = X\Y;
        % Compute the errors
        e= Y-X*rho;
        s2(i+1)=e'*e/tau;
        K(i+1)=size(X,2);

        % Loop
        for i=1:maxlags
            X=[ones(size(Y)) y(maxlags+1:T-1) (1:tau)' ydifflags(:,1:i)];
            rho = X\Y;
            % Compute the errors
            e= Y-X*rho;
            s2(i+1)=e'*e/tau;
            K(i+1)=size(X,2);
        end
end

if strcmp(IC,'AIC')
    ICs=log(s2) + 2*K/tau;
else
    ICs=log(s2) + K*log(tau)/tau;
end
[~,lags]=min(ICs);
lags=lags-1;
[adfstat,pval,critval,resid]=augdf(y,p,lags);

author: bashtage

repo: mfe-toolbox


text:
function hdr=cosmo_map2fmri(dataset, fn, varargin)
% maps a dataset structure to a NIFTI, AFNI, or BV structure or file
%
% Usage 1: hdr=cosmo_map2fmri(dataset, '-{FMT}) returns a header structure
% Usage 2: cosmo_map2fmri(dataset, fn) saves dataset to a volumetric file.
%
% In Usage 1, {FMT} can be one of 'nii','bv_vmp',bv_vmr','bv_msk','afni'.
% In Usage 2, fn should end with '.nii.gz', '.nii', '.hdr', '.img', '.vmp',
%             '.vmr', '.msk', '+orig','+orig.HEAD','+orig.BRIK',
%             '+orig.BRIK.gz','+tlrc','+tlrc.HEAD','+tlrc.BRIK', or
%             '+tlrc.BRIK.gz'.
%
% Use cosmo_map2fmri(...,'deoblique',true) to de-oblique the dataset. For
% information on this option, see cosmo_fmri_deoblique
%
% - for NIFTI files, it requires the following toolbox:
%   http://www.mathworks.com/matlabcentral/fileexchange/8797-tools-for-nifti-and-analyze-image
%   (note that his toolbox is included in CoSMoMVPA in /externals)
% - for Brainvoyager files (.vmp and .vtc), it requires the NeuroElf
%   toolbox, available from: http://neuroelf.net
% - for AFNI files (+{orig,tlrc}.{HEAD,BRIK[.gz]}) it requires the AFNI
%   Matlab toolbox, available from: http://afni.nimh.nih.gov/afni/matlab/
%
% See also: cosmo_fmri_deoblique
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    defaults=struct();
    defaults.deoblique=false;
    defaults.bv_force_fit=[];
    opt=cosmo_structjoin(defaults,varargin{:});

    cosmo_check_dataset(dataset, 'fmri');

    img_formats=get_img_formats();
    sp=cosmo_strsplit(fn,'-');
    save_to_file=~isempty(sp{1});

    if save_to_file
        fmt=get_format(img_formats, fn);
    else
        if numel(sp)~=2
            error('expected -{FORMAT}');
        end
        fmt=sp{2};
    end

    if ~isfield(img_formats,fmt)
        error('Unsupported format %s', fmt);
    end

    methods=img_formats.(fmt);
    externals=methods.externals;
    cosmo_check_external(externals);

    creator=methods.creator;

    % preprocess the data (de-oblique if necessary)
    dataset=preprocess(dataset,opt);

    % build header structure
    hdr=creator(dataset);

    if nargout==0
        cleaner=onCleanup(get_cleaner(methods,hdr));
    end

    if save_to_file
        writer=methods.writer;
        writer(fn, hdr);
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% general helper functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f=get_cleaner(methods,hdr)
    if isfield(methods,'cleaner')
        f=@()methods.cleaner(hdr);
    else
        f=@do_nothing;
    end

function do_nothing
    % do nothing

function ds=preprocess(ds,opt)
    if opt.deoblique
        ds=cosmo_fmri_deoblique(ds);
    end

    if opt.bv_force_fit
        ds=bv_force_fit(ds);
    end

function ds=bv_force_fit(ds)

    mat=ds.a.vol.mat;

    dim=[ds.a.vol.dim(:);1];
    rot=mat(1:3,1:3);
    vox_size=sqrt(sum(rot.^2,1)); %.*sign(sum(rot,1));
    resolution=prod(vox_size).^(1/3);

    % put the center defined by the original matrix close to the center
    % defined by the new matrix
    old_center=mat*(dim+1)/2;
    mat(1:3,1:3)=bsxfun(@times,mat(1:3,1:3),resolution./vox_size);

    new_center=mat*(dim+1)/2;
    mat(:,4)=round(mat(:,4)+old_center-new_center)+.5;

    ds.a.vol.mat=mat;

function check_plump_orientation(ds)
    rot=ds.a.vol.mat(1:3,1:3);
    nonzero=rot~=0;
    is_oblique=any(sum(nonzero,1)~=1) || any(sum(nonzero,2)~=1);
    if is_oblique
        error(['Dataset has oblique orientation, and the specified '...
                'output format does not support this. Options are:\n\n'...
                '(1) use %s(...,''deoblique'',true) to deoblique the '...
                'dataset automatically when writing it\n'...
                '(2) use cosmo_fmri_deoblique to deoblique the dataset '...
                'manually\n'...
                '(3) store the dataset in another file format (such as '...
                'NIFTI)\n\n'...
                'For options (1) and (2), the spatial location of the '...
                'voxels will change. It is recommended to inspect the '...
                'results visually when using one of these options.'],...
                mfilename());
    end

function check_isotropic_voxels(ds)
    rot=ds.a.vol.mat(1:3,1:3);
    vox_size=sqrt(sum(rot.^2,1));

    max_delta=1e-5;
    delta=max(vox_size)-min(vox_size);

    if delta>max_delta
        error(['Dataset has non-isotropic voxels, and the specified '...
                'output format does not support this. Options are:\n\n'...
                '(1) use %s(...,''bv_force_fit'',true) to set the '...
                'voxel size to be isotropic\n'...
                '(2) store the dataset in another file format (such as '...
                'NIFTI)\n\n'...
                'For option (1), the spatial location of the '...
                'voxels will change. It is recommended to inspect the '...
                'results visually when using one of these options.'],...
                mfilename());
    end

function unfl_ds_timelast=unflatten(ds)
    % puts the time dimension last, instead of first
    unfl_ds=cosmo_unflatten(ds);

    % make time dimension last instead of first
    unfl_ds_timelast=shiftdim(unfl_ds,1);

    % deal with singleton dimension
    sz=size(unfl_ds);
    dim_size=ones(1,4);
    dim_size(4)=sz(1);
    dim_size(1:(numel(sz)-1))=sz(2:end);

    nsamples=size(ds.samples,1);
    assert(isequal([ds.a.vol.dim nsamples],dim_size));

    if ~isequal(size(unfl_ds_timelast),dim_size)
        unfl_ds_timelast=reshape(unfl_ds_timelast,dim_size);
    end

function b=ends_with(end_str, str)
    if iscell(end_str)
        b=any(cellfun(@(x) ends_with(x,str),end_str));
    else
        b=isempty(cosmo_strsplit(str, end_str,-1));
    end

function fmt=get_format(img_formats, file_name)
    fns=fieldnames(img_formats);
    fmt=[];
    for k=1:numel(fns)
        fn=fns{k};
        exts=img_formats.(fn).exts;
        if ends_with(exts, file_name)
            fmt=fn;
        end
    end
    if isempty(fmt)
        error('Not found: format for %s', file_name);
    end

function img_formats=get_img_formats()
    img_formats=struct();

    img_formats.nii.creator=@new_nii;
    img_formats.nii.writer=@write_nii;
    img_formats.nii.externals={'nifti'};
    img_formats.nii.exts={'.nii','.nii.gz','.hdr','.img'};

    img_formats.bv_vmp.creator=@new_bv_vmp;
    img_formats.bv_vmp.writer=@write_bv;
    img_formats.bv_vmp.cleaner=@clean_bv;
    img_formats.bv_vmp.externals={'neuroelf'};
    img_formats.bv_vmp.exts={'.vmp'};

    img_formats.bv_vmr.creator=@new_bv_vmr;
    img_formats.bv_vmr.writer=@write_bv;
    img_formats.bv_vmr.cleaner=@clean_bv;
    img_formats.bv_vmr.externals={'neuroelf'};
    img_formats.bv_vmr.exts={'.vmr'};

    img_formats.bv_msk.creator=@new_bv_msk;
    img_formats.bv_msk.writer=@write_bv;
    img_formats.bv_msk.cleaner=@clean_bv;
    img_formats.bv_msk.externals={'neuroelf'};
    img_formats.bv_msk.exts={'.msk'};

    img_formats.afni.creator=@new_afni;
    img_formats.afni.writer=@write_afni;
    img_formats.afni.externals={'afni'};
    img_formats.afni.exts={'+orig','+orig.HEAD','+orig.BRIK',...
                           '+orig.BRIK.gz','+tlrc','+tlrc.HEAD',...
                           '+tlrc.BRIK','+tlrc.BRIK.gz'};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% format-specific helper functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Nifti

function ni=new_nii(ds)
    a=ds.a;
    vol=a.vol;
    vol_data=unflatten(ds);
    dim_all=[4 ones(1,7)];
    vol_dim=size(vol_data);
    nvol_dim=numel(vol_dim);
    dim_all(1+(1:nvol_dim))=vol_dim;

    mat=vol.mat;
    mat(1:3,4)=mat(1:3,4)+mat(1:3,1:3)*[1 1 1]';
    pix_dim=sqrt(sum(mat(1:3,1:3).^2,1));
    hdr=struct();

    dime=struct();
    dime.datatype=16; %single
    dime.dim=dim_all;
    dime.pixdim=[0 abs(pix_dim(:))' 0 0 0 0]; % ensure positive values

    fns={'intent_p1','intent_p2','intent_p3','intent_code',...
        'slice_start','slice_duration','slice_end',...
        'scl_slope','scl_inter','slice_code','cal_max',...
        'cal_min','toffset'};

    dime=set_all(dime,fns);
    dime=cosmo_structjoin(dime, cosmo_statcode(ds,'nifti'));
    dime.xyzt_units=10;
    hdr.dime=dime;

    hk=struct();
    hk.sizeof_hdr=348;
    hk.data_type='';
    hk.db_name='';
    hk.extents=0;
    hk.session_error=0;
    hk.regular='r';
    hk.dim_info=0;
    hdr.hk=hk;

    hist=struct();

    if isfield(ds.a.vol,'xform')
        xform=ds.a.vol.xform;
    else
        xform=''; % unknown
    end
    xform_code=cosmo_fmri_convert_xform('nii',xform);

    hist.sform_code=xform_code;
    % hist.originator=[1 1 1 1 0];
    hist=set_all(hist,{'descrip','aux_file','intent_name'},'');
    hist=set_all(hist,{'qform_code','quatern_b',...
                        'quatern_d',...
                        'qoffset_x','qoffset_y','qoffset_z'});
    hist=set_all(hist,{'intent_name'},'');
    hist.srow_x=mat(1,:);
    hist.srow_y=mat(2,:);
    hist.srow_z=mat(3,:);
    hist.quatern_c=0;

    % required to avoid complaints with hdr/img
    hist.originator=round(dim_all(2:5)/2);
    hdr.hist=hist;

    ni.img=single(vol_data);
    ni.hdr=hdr;

function write_nii(fn, hdr)
    save_nii(hdr, fn);

    %% BrainVoyager helper
function mat=neuroelf_bvcoordconv_wrapper(varargin)
    % converts BV bounding box to affine transformation matrix
    % helper functions that deals with both new neuroelf (version 1.0)
    % and older versions.
    % the old version provides a 'bvcoordconv' .m file
    % the new version privides this function in the neuroelf class
    has_bvcoordconv=~isempty(which('bvcoordconv'));

    % set function handle
    if has_bvcoordconv
        f=@bvcoordconv;
    else
        n=neuroelf();
        f=@n.bvcoordconv;
    end

    mat=f(varargin{:});

function result=neuroelf_bless_wrapper(arg)
    % deals with recent neuroelf (>v1.1), where bless is deprecated
    s=warning('off','neuroelf:xff:deprecated');
    resetter=onCleanup(@()warning(s));

    result=bless(arg);

    %% Brainvoyager VMP
function [hdr,ds]=new_bv_mat_hdr(ds,bv_type)
    % ensure dataset is plump
    check_plump_orientation(ds);

    if ~strcmp(bv_type,'vmr')
        % require voxels to be isotropic
        check_isotropic_voxels(ds);
    end

    % automatically set orientation to ASR
    if ~strcmp(cosmo_fmri_orientation(ds),'ASR')
        ds=cosmo_fmri_reorient(ds,'ASR');
    end

    % helper to set matrix in header
    mat=ds.a.vol.mat;

    % ensure proper VMP/VMR orientation
    rot_asr=mat(1:3,1:3)*[0 0 -1; -1 0 0; 0 -1 0]';

    if ~isequal(rot_asr,diag(diag(rot_asr)))
        % should never get here because of the ASR orientation
        % forced earlier
        error('Unsupported orientation: need ASR');
    end

    hdr=xff(['new:' bv_type]);
    hdr=neuroelf_bless_wrapper(hdr);

    % Set {X,Y,Z}{Start,End} values based on the transformation matrix
    tal_coords=mat*[1 1 1 1; ds.a.vol.dim+1, 1]';
    bv_coords=neuroelf_bvcoordconv_wrapper(tal_coords(1:3,:), ...
                                        'tal2bvs',hdr.BoundingBox);

    switch bv_type
        case {'vmp','msk'}
            dg=diag(rot_asr);
            resolution=prod(dg)^(1/3);

            labels={'ZStart','ZEnd','XStart','XEnd','YStart','YEnd'};
            for k=1:numel(labels)
                label=labels{k};
                hdr.(label)=bv_coords(k);
            end

            hdr.Resolution=resolution;

        case 'vmr'
            % this *should* be a 256^3 volume
            % XXX check for this?
            labels={'X','Y','Z'};
            dg=diag(rot_asr);
            for k=1:3
                label=labels{k};
                hdr.(['VoXRes' label])=dg(k);
                hdr.(['Dim' label])=round(bv_coords(2,k));
            end

        otherwise
            error('Unsupported type %s', bv_type);
    end

function hdr=new_bv_vmp(ds)
    [hdr,ds]=new_bv_mat_hdr(ds,'vmp');

    % Get data
    nsamples=size(ds.samples,1);
    maps=cell(1,nsamples);

    stats=cosmo_statcode(ds,'bv');

    for k=1:nsamples
        empty_hdr=xff('new:vmp');
        map=empty_hdr.Map;
        ds_k=cosmo_slice(ds,k);
        map.VMPData=unflatten(ds_k);
        if cosmo_isfield(ds_k,'sa.labels')
            map.Name=[ds_k.sa.labels{:}];
        end
        if ~isempty(stats)
            stats_k=stats{k};

            stats_k_fieldnames=fieldnames(stats_k);
            for j=1:numel(stats_k_fieldnames)
                key=stats_k_fieldnames{j};
                map.(key)=stats_k.(key);
            end
        end
        maps{k}=map;
        empty_hdr.ClearObject();
    end

    hdr.Map=cat(2,maps{:});
    hdr.NrOfMaps=nsamples;

    neuroelf_bless_wrapper(hdr);

function write_bv(fn, hdr)
    % general storage function
    hdr.SaveAs(fn);

function clean_bv(hdr)
    hdr.ClearObject();

    %% Brainvoyager VMR
function hdr=new_bv_vmr(ds)
    nsamples=size(ds.samples,1);
    if nsamples~=1,
        error('Unsupported: more than 1 sample');
    end

    [hdr,ds]=new_bv_mat_hdr(ds,'vmr');

    % scale to 0..255
    vol_data=unflatten(ds);
    hdr.VMRData=scale_uint8(vol_data(:,:,:,1));

    %% Brainvoyager mask
function hdr=new_bv_msk(ds)
    nsamples=size(ds.samples,1);
    if nsamples~=1,
        error('Unsupported: more than 1 sample');
    end

    [hdr,ds]=new_bv_mat_hdr(ds,'msk');

    vol_data=unflatten(ds);
    hdr.Mask=scale_uint8(vol_data);

function scaled_data=scale_uint8(data)
    mn=min(data(:));
    mx=max(data(:));
    if mn<0 || mx>255 || ~isequal(round(data),data)
        cosmo_warning(['.samples does not fit integer range 0..255 ',...
                            'data will be rescaled']);
        data=255*(data-mn)/(mx-mn);
    end

    scaled_data=uint8(data);

%% AFNI
function afni_info=new_afni(ds)
    check_plump_orientation(ds);

    a=ds.a;
    mat=a.vol.mat;

    % deal with orientation
    idxs=zeros(1,3);
    m=zeros(1,3);

    % this part should be necessary
    for k=1:3
        % for each spatial dimension, find which row transforms it
        idx=find(mat(1:3,k));
        assert(numel(idx)==1); % ds is plump

        % store index and transformation matrix
        idxs(k)=idx;
        m(k)=mat(idx,k);
    end

    % set voxel size and origin
    % as the header was stored in LPI but AFNI likes RAI,
    % convert coordinates to RAI
    lpi2rai=[-1 -1 1];

    delta=m.*lpi2rai(idxs);
    origin=mat(idxs,:)*[1 1 1 1]'.*lpi2rai(idxs)';

    % set orientation code.
    % thse are neither RAI or LPI, but RPI (when ordered logically)
    % No idea why Bob Cox made that decision.
    lpi2orient=[-1 1 1];
    offset=(1-sign(m).*lpi2orient(idxs))/2;
    orient=(idxs-1)*2+offset;

    [nsamples,nfeatures]=size(ds.samples);

    vol_data=unflatten(ds);
    dim=ds.a.vol.dim(:)';
    vol_dim=[size(vol_data) 1]; % in case of singleton in last dimension
    assert(isequal(vol_dim(1:3),dim));
    assert(prod(dim)>=nfeatures);

    brik_type=1; %functional head
    brik_typestring='3DIM_HEAD_FUNC';
    brik_func=11; % FUNC_BUCK_TYPE

    if isfield(ds.a.vol,'xform')
        xform=ds.a.vol.xform;
    else
        xform=''; % unknown
    end
    xform_code=cosmo_fmri_convert_xform('afni',xform);

    brik_view=xform_code; % default to +orig, but overriden by writer

    afni_info=struct();
    afni_info.SCENE_DATA=[brik_view, brik_func, brik_type];
    afni_info.TYPESTRING=brik_typestring;
    afni_info.BRICK_TYPES=3*ones(1,nsamples); % store in float format
    afni_info.BRICK_STATS=[];                 % ... and thus need no stats
    afni_info.BRICK_FLOAT_FACS=[];            % ... or multipliers
    afni_info.DATASET_RANK=[3 nsamples];      % number of volumes
    afni_info.DATASET_DIMENSIONS=dim(1:3);
    afni_info.ORIENT_SPECIFIC=orient;
    afni_info.DELTA=delta(:)';
    afni_info.ORIGIN=origin(:)';
    afni_info.SCALE=0;
    afni_info.NOTES_COUNT=0;
    afni_info.WARP_TYPE=[0 0];
    afni_info.FileFormat='BRIK';

    [unused, unused, endian_ness] = computer();
    afni_info.BYTEORDER_STRING = sprintf('%sSB_FIRST', endian_ness);

    set_empty={'BRICK_LABS','BRICK_KEYWORDS',...
                'BRICK_STATS','BRICK_FLOAT_FACS',...
                'BRICK_STATAUX','STAT_AUX'};
    for k=1:numel(set_empty)
        fn=set_empty{k};
        afni_info.(fn)=[];
    end

    % if labels for the samples, store them in the header
    if isfield(ds.sa,'labels') && ~isempty(ds.sa.labels)
        afni_info.BRICK_LABS=cosmo_strjoin(ds.sa.labels,'~');
    end

    if isfield(ds.sa,'stats') && ~isempty(ds.sa.stats)
        afni_info=cosmo_structjoin(afni_info,cosmo_statcode(ds,'afni'));
    end

    % store data in non-afni field 'img'
    afni_info.img=vol_data;

function write_afni(fn, hdr)

    hdr.RootName=fn;
    data=hdr.img; % get the data
    hdr=rmfield(hdr,'img'); % remove the non-afni field 'img'

    afniopt=struct();
    afniopt.Prefix=fn; %the second input argument
    afniopt.OverWrite='y';
    afniopt.NoCheck=0;
    afniopt.AppendHistory=false;
    afniopt.verbose=0;
    afniopt.Scale=0;
    afniopt.AdjustHeader='no';

    if ends_with({'+orig','+orig.HEAD','+orig.BRIK',...
                           '+orig.BRIK.gz'},fn)
        hdr.SCENE_DATA(1)=0;
        afniopt.View='+orig';
    elseif ends_with({'+tlrc','+tlrc.HEAD',...
                           '+tlrc.BRIK','+tlrc.BRIK.gz'},fn)
        hdr.SCENE_DATA(1)=2;
        afniopt.View='+tlrc';
    else
        error('Unsupported scene data for %s', fn);
    end

    [err, ErrMessage]=WriteBrik(data, hdr, afniopt);
    if err
        error(ErrMessage);
    end

 function s=set_all(s, fns, v)
    % sets all fields in fns in struct s to v
    % if v is omitted it is set to 0.
    if nargin<3, v=0; end
    n=numel(fns);
    for k=1:n
        fn=fns{k};
        s.(fn)=v;
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function anonTextureS = defineAnonTexture()
%
% function anonTextureS = defineAnonCerroptions()
%
% Defines the anonymous "texture" data structure. The fields containing PHI have
% been left out by default. In order to leave out additional fields, remove
% them from the list below.
%
% The following three operations are allowed per field of the input data structure:
% 1. Keep the value as is: listed as 'keep' in the list below.
% 2. Allow only the specific values: permitted values listed within a cell array.
%    If no match is found, a default anonymous string will be inserted.
% 3. Insert dummy date: listed as 'date' in the list below. Date will be
%    replaced by a dummy value of 11/11/1111.
%
% APA, 1/11/2018

anonTextureS = struct( ...
    'category',  'keep', ...
    'description',  'keep', ...
    'patchSize',  'keep', ...
    'patchUnit',  'keep', ...
    'paramS',  'keep', ...
    'assocScanUID',  'keep', ...
    'assocStructUID',  'keep', ...
    'textureUID',  'keep' ...
           );

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [currSelectedId,selTypeC] = dispSelCriteriaROE(hObj,hEvt,hFig,varargin)
%Display selected constraints
% AI 05/12/21

%cMode = datacursormode(hFig);
%cMode.removeAllDataCursors;

ud = guidata(hFig);
legH = ud.handle.legend;

if isempty(hEvt)  %Prog call

    %Get handles to constraints
    protS = ud.Protocols;
    type = varargin{1};
    idxV = varargin{2};
    pNum = varargin{3};

    %Turn off currently displayed limits
    hCrit = protS(pNum).criteria;
    hGuide = protS(pNum).guidelines;
    for k = 1:numel(hCrit)
        set(hCrit(k),'Visible','Off')
    end
    for k = 1:numel(hGuide)
        set(hGuide(k),'Visible','Off')
    end
    protS(pNum).criteria = hCrit;
    protS(pNum).guidelines = hGuide;

    %Turn on selected limit
    if strcmp(type,'criteria')
        hCrit = [protS(pNum).criteria];
        set(hCrit(idxV),'Visible','On');
        %numElements = [0,cumsum(arrayfun(@(x)numel(x.criteria),protS))];
        %for pNum = 1:numel(protS)
        %protS(pNum).criteria = hCrit(numElements(pNum)+1:numElements(pNum+1));
        %end
        protS(pNum).criteria = hCrit;

        drawnow;
        legH.EntryContainer.NodeChildren(2).Label.Color = [0,0,0];
        legH.EntryContainer.NodeChildren(1).Label.Color = [0.65,0.65,0.65];
    else
        hGuide = [protS(pNum).guidelines];
        set(hGuide(idxV),'Visible','On');
        %numElements = [0,cumsum(arrayfun(@(x)numel(x.guidelines),protS))];
        %for pNum = 1:numel(protS)
        %protS(pNum).guidelines = hGuide(numElements(pNum)+1:numElements(pNum+1));
        %end
        protS(pNum).guidelines = hGuide;

        drawnow;
        legH.EntryContainer.NodeChildren(1).Label.Color = [0,0,0];
        legH.EntryContainer.NodeChildren(2).Label.Color = [0.65,0.65,0.65];
    end

    ud.Protocols = protS;
    set(hFig,'userdata',ud);

else %Checkbox selection

    %Get handles to constraints
    protS = ud.Protocols;

    %Get slelected constraint
    currSelectedId = hEvt.Indices(:,1);
    stateV = cell2mat(hObj.Data(currSelectedId,1));
    stateC = {'Off','On'};

    if currSelectedId==1  %'All'
        for pNum = 1:numel(protS)
            hCrit = protS(pNum).criteria;
            hGuide = protS(pNum).guidelines;
            %Criteria
            for k = 1:numel(hCrit)
                set(hCrit(k),'Visible',stateC{stateV+1});
            end
            %Guidelines
            for k = 1:numel(hGuide)
                set(hGuide(k),'Visible',stateC{stateV+1});
            end
        end

        drawnow;
        legH.EntryContainer.NodeChildren(1).Label.Color = [0,0,0];
        legH.EntryContainer.NodeChildren(2).Label.Color = [0,0,0];

    elseif currSelectedId==2 %'None'
        if stateV == 1
            %Criteria
            for pNum = 1:numel(protS)
                hCrit = protS(pNum).criteria;
                hGuide = protS(pNum).guidelines;
                for k = 1:numel(hCrit)
                    set(hCrit(k),'Visible','Off');
                end
                %Guidelines
                for k = 1:numel(hGuide)
                    set(hGuide(k),'Visible','Off');
                end
            end
        end
        hObj.Data(:,1) = {false};

        drawnow;
        legH.EntryContainer.NodeChildren(1).Label.Color = [0.65,0.65,0.65];
        legH.EntryContainer.NodeChildren(2).Label.Color = [0.65,0.65,0.65];

    else

        protS = ud.Protocols;
        currProtocol = ud.foreground;
        gNum = numel(protS(currProtocol).guidelines);

        currSelectedId = currSelectedId-2;
        type =  get(hObj,'userdata');
        selTypeC = type(currSelectedId);

        for pNum = 1:numel(protS)
            for k = 1:numel(currSelectedId)
                if strcmp(selTypeC(k),'guidelines') %guidelines
                    selNum = currSelectedId(k);
                    if strcmp(stateC{stateV+1},'On')
                        drawnow;
                        legH.EntryContainer.NodeChildren(1).Label.Color = [0,0,0];
                    else
                        drawnow;
                        legH.EntryContainer.NodeChildren(1).Label.Color = [0.65,0.65,0.65];
                    end
                else                               %criteria
                    selNum = currSelectedId(k)- gNum;
                    if strcmp(stateC{stateV+1},'On')
                        drawnow;
                        legH.EntryContainer.NodeChildren(2).Label.Color = [0,0,0];
                    else
                        drawnow;
                        legH.EntryContainer.NodeChildren(2).Label.Color = [0.65,0.65,0.65];
                    end
                end
                %Toggle display on/off
                set(protS(pNum).(selTypeC{k})(selNum),'Visible',stateC{stateV(k)+1});
                %Expand tooltip if on
                if strcmp(stateC{stateV+1},'On')
                    %Check if datatip exists
                    if ~isfield(ud.handle,'datatips')
                        ud.handle.datatips = [];
                    end

                    target = protS(pNum).(selTypeC{k})(selNum);
                    guidata(hFig,ud);
                    createDataTipROE(hFig,target);
                end

            end
        end

    end

end

ud.Protocols = protS;
ud.handle.legend = legH;

end