text:
function canlab_glm_group_levels(varargin)
% Performs group level robust GLM analysis with robfit
%   1. sets up analysis
%   2. runs robfit
%   3. (optionally) make inverse p maps (for FSL viewing)
%   4. (optionally) estimate significant cluster sizes
%   5. publishes analysis with robfit_results_batch
%
% :Usage:
% ::
%
%     canlab_glm_group_levels([options])
%
% :Optional Inputs:
%
%   **'s', subjects:**
%        cell array of filenames of subject-level analysis directories IN
%        modeldir
%        (note: modeldir won't be prepended to absolute paths)
%
%   **'m', modeldir:**
%        filename of directory containing subject level analyses
%
%   **'o', grpmodeldir:**
%        output directory name
%
%   **'c', cov:**
%        a matrix describing group level model
%        (do not include intercept, it is automatically included as first regressor)
%
%        see help robfit
%
%        note: requires specifying an output directory name
%
%        note: ordering of inputs (rows) must match subjects ordering
%
%   **'n', covname:**
%        a cell array of names for the covariates in cov
%
%   **'f', covfile:**
%        a csv file will specify the group level model:
%        first column, header 'subject', contains names of subject
%        directories (to be found in modeldir)
%        subsequent columns have covariates, headers are names of covariates
%        name of covfile will be name of group analysis directory (placed in
%        grpmodeldir)
%
%   **'mask', maskimage:**
%        filename of mask image
%
%   **DSGN:**
%        will use the following fields of the DSGN structure:
%        modeldir    = DSGN.modeldir
%
%        subjects    = DSGN.subjects
%
%        maskimage   = DSGN.mask
%
% :Note:
%   A covfile will cause other specifications of subject, cov, and
%   covnames to be ignored.
%
%   If parameters are defined more than once (e.g., modeldir or subjects),
%   only the last entered option will count.
%
% :Defaults:
%
%   **subjects:**
%        all SPM.mat-containing directories in modeldir
%
%   **modeldir:**
%        pwd
%
%   **grpmodeldir:**
%        modeldir/one_sample_t_test
%
%   **cov:**
%        {} (run 1 sample t-test, see help robfit)
%
%   **covname:**
%        'groupmean'
%
%   **mask:**
%        'brainmask.nii'
%
% :Options:
%
%   **'README':**
%        prints canlab_glm_README, an overview of canlab_glm_{subject,group}_levels
%
%   **'overwrite':**
%        overwrite existing output directories
%
%   **'noresults':**
%        don't run/publish robfit_results_batch
%
%   **'onlyresults':**
%        just run/publish robfit_results_batch, don't run robfit (assumes existing analyses)
%
%   **'whichcons', [which cons]
%       vector of contrasts to analyze (DEFAULT: aall subject level contrasts)
%       see [which cons] in help robfit
%
%   **'invp' [, target_space_image]:**
%       generate inverse p maps and resample to the voxel and image dimensions
%       of target_space_image
%       (viewable on dream with ~ruzicl/scripts/invpview)
%
%   **'nolinks':**
%        do not make directory of named links to robust directories (using contrast names)
%
%   **'dream':**
%        if you're running on the dream cluster, this option will cause
%       all analyses (e.g., lower level contrasts) to be run in parallel
%       (submitted with matlab DCS and the Sun Grid Engine)
%       Note: currently only works with MATLAB R2009a
%
%   **'email', address:**
%       send notification email to address when done running
%
% ..
%     ----------------------------------------------------------------------
%     Copyright (C) 2013  Luka Ruzic
%
%     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/>.
%     ----------------------------------------------------------------------
%
%     Programmers' notes
%     NOT READY YET:
%     'grf'
%       estimate significant cluster sizes for each contrast using GRF theory
%       will run at a=.05, voxelwise thresholds of .05, .01, .005, .001
%       see help estimate_cluster_extent
% ..

STARTTIME = datestr(now,31); % SET UP

STARTINGDIR = pwd;

% path
addpath(genpath('/usr/local/spm/spm8/matlabbatch'));
addpath(genpath('/usr/local/spm/spm8/apriori'));
addpath(genpath('/usr/local/spm/spm8/canonical'));

% set defaults
OPTS.parallel_dream = '';
OPTS.overwrite = false;
OPTS.run_results = true;
OPTS.run_robfit = true;
OPTS.run_printconnames = true;
OPTS.run_invp = false;
OPTS.run_grf = false;
OPTS.run_links = true;
OPTS.nocatch = false;

pthresh = [.05 .01 .005 .001]; % for significant cluster extent estimation

EXPT.cov = [];
EXPT.mask = which('brainmask.nii');

% set up statuses
robfitstatus = 0;
invpstatus = 0;
grfstatus = 0;
linkdirstatus = 0;
reviewstatus = 0;

diarydirname = 'canlab_glm_logs';
diaryfilename = ['group_levels_' regexprep(regexprep(STARTTIME,' ','_'), ':','') '.log'];

%% parse arguments
i=1;
while i<=numel(varargin)
    if isstruct(varargin{i})
        DSGN = varargin{i};
        if isfield(DSGN,'modeldir')
            modeldir = DSGN.modeldir;
        end
        if isfield(DSGN,'subjects')
            for s = 1:numel(DSGN.subjects)
                [dummy subname{s}] = fileparts(DSGN.subjects{s}); %#ok
                sublevs{s} = subname{s}; %#ok
            end
        end
        if isfield(DSGN,'mask')
            EXPT.mask = DSGN.mask;
        end
    elseif ischar(varargin{i})
        switch varargin{i}
            case {'README'}
                if ispc()
                    system(sprintf('type %s',which('canlab_glm_README.txt')));
                else
                    system(sprintf('cat %s',which('canlab_glm_README.txt')));
                end
                return;
            case {'dream'}
                version = ver('matlab');
                switch version.Release
                    case '(R2009a)'
                        OPTS.parallel_dream = '2009a';
                    case '(R2011a)'
                        OPTS.parallel_dream = '2011a';
                    otherwise
                        error('Current version of matlab (%s)',version.Release)
                end
            case {'f'}
                i=i+1;
                covfile = varargin{i};
            case {'c'}
                i=i+1;
                EXPT.cov = varargin{i};
            case {'n'}
                i=i+1;
                EXPT.covnames = ['intercept' varargin{i}];
            case {'s'}
                i=i+1;
                sublevs = varargin{i};
            case {'d' 'm'}
                i=i+1;
                modeldir = varargin{i};
            case {'mask'}
                i=i+1;
                EXPT.mask = varargin{i};
            case {'overwrite'}
                OPTS.overwrite = true;
            case {'noresults'}
                OPTS.run_results = false;
            case {'onlyresults'}
                OPTS.run_robfit = false;
            case {'whichcons'}
                i=i+1;
                includedcons = varargin{i};
            case {'invp'}
                OPTS.run_invp = true;
            case {'grf'}
                OPTS.run_grf = true;
            case {'nolinks'}
                OPTS.run_links = false;
            case {'nocatch'}
                OPTS.nocatch = true;
            case {'o'}
                i=i+1;
                grpmodeldir = varargin{i};
            case {'email'}
                i=i+1;
                address = varargin{i};
            otherwise
                error(['Unrecognized argument: ' varargin{i}])
        end
    else
        disp(varargin{i})
        error('Above argument is unrecognized.')
    end
    i=i+1;
end

% MODELDIR
if ~exist('modeldir','var')
    modeldir = pwd;
end
if ~exist(modeldir,'dir')
    error(['No such directory: ' modeldir])
else
    % make absolute (unsure if necessary)

    % This is a bug on discovery
%     modeldir = filenames(modeldir,'absolute','char');
end

% SUBLEVS
if exist('sublevs','var')
    for s = 1:numel(sublevs)
        if ~strcmp(sublevs{s}(1),'/')
            sublevs{s} = fullfile(modeldir,sublevs{s}); %#ok
        end
    end
end

% GRPMODELDIR
if ~exist('grpmodeldir','var')
    if exist('modeldir','var')
        grpmodeldir = fullfile(modeldir,'one_sample_t_test');
    else
        grpmodeldir = fullfile(pwd,'one_sample_t_test');
    end
end
% make grpmodeldir
if ~exist(grpmodeldir,'dir')
    mkdir(grpmodeldir)
end
% make absolute
grpmodeldir = filenames(grpmodeldir,'absolute','char');

% COV
if exist('covfile','var')
    if ~exist(covfile,'file')
        error(['No such file: ' covfile])
    else
        [d f e] = fileparts(covfile); %#ok

        % grpmodeldir
        grpmodeldir = fullfile(fileparts(grpmodeldir),f);

        % sublevs
        tempsubsfile = fullfile(d,[f '_tempsubs.txt']);
        system(['cut -d, -f1 < ' covfile ' | tail -n +2 > ' tempsubsfile]);
        sublevs = importdata(tempsubsfile); delete(tempsubsfile);
        if isnumeric(sublevs), sublevs = textscan(num2str(sublevs'), '%s'); sublevs=sublevs{1}; end
        for i=1:numel(sublevs)
            sublevs{i} = fullfile(modeldir,sublevs{i});
        end

        % cov, covnames
        tempcovsfile = fullfile(d,[f '_tempcovs.txt']);
        system(['cut -s -d, -f2- < ' covfile ' > ' tempcovsfile]);
        s = importdata(tempcovsfile); delete(tempcovsfile);
        if isfield(s,'data')
            EXPT.cov = s.data;
            EXPT.covnames = ['intercept' s.colheaders];
        else
            EXPT.covnames = {'groupmean'};
        end
    end
else
    if isempty(EXPT.cov)
        EXPT.covnames = {'groupmean'};
    else
        if ~exist('grpmodeldir','var')
            error('cov has been specified: must specify grpmodeldir.')
        end
        if isfield(EXPT,'covnames')
            if numel(EXPT.cov) ~= numel(EXPT.covnames)
                error('Number of covariate names (%d) does not match number of covariates (%d).',numel(EXPT.covnames),numel(EXPT.cov))
            end
        else
            EXPT.covnames{1} = 'intercept';
            for i=2:numel(EXPT.cov)
                EXPT.covnames{i} = sprintf('grpcov%d',i-1);
            end
        end
    end
end

% set up diary - to be on except when calling other functions that do their own logging
diarydir = fullfile(grpmodeldir,diarydirname);
if ~exist(diarydir,'dir'), mkdir(diarydir); end
% diarydir = filenames(diarydir,'absolute','char');
diarydir = filenames(diarydir, 'char');
diaryname = fullfile(diarydir,diaryfilename);
fulldiaryname = regexprep(diaryname,'\.log$','_full.log');
fprintf('Writing to logfile: %s\n',diaryname)
diary(diaryname), fprintf('> STARTED: %s\n',STARTTIME), diary off

% make working directory
wd = regexprep(diaryname,'\.log$','');
mkdir(wd)

%% SET UP AND RUN ROBFIT
diary(diaryname), announce_string('ROBUST REGRESSION'), diary off
if ~OPTS.run_robfit
    diary(diaryname), fprintf('> SKIPPED: turned off in options\n'); diary off
else
    % subject level analyses + 1 SPM.mat (assumed representative)
    if ~exist('sublevs','var')
        diary(diaryname), fprintf('> No subject levels specified: looking for all subject level analyses in:\n\t%s\n',modeldir); diary off
%         spmfiles = filenames(fullfile(modeldir,'*/SPM.mat'),'absolute');
        spmfiles = filenames(fullfile(modeldir,'*/SPM.mat'));
        if isempty(spmfiles)
            diary(diaryname), fprintf('> ERROR: No subject level analyses found.\n'), diary off
            cd(STARTINGDIR)
            return
        end
        for i=1:numel(spmfiles), [sublevs{i}] = fileparts(spmfiles{i}); end %#ok
        sublevs = unique(sublevs);
        spmfile = spmfiles{1};
    else
        spmfile = fullfile(sublevs{1},'SPM.mat');
    end

    diary(diaryname)
    fprintf('> Subject level analyses:\n')
    for i=1:numel(sublevs)
        fprintf('>  %-8d%s\n',i,sublevs{i})
    end
    diary off

    if exist(spmfile,'file')
        load(spmfile)
    else
        diary(diaryname), fprintf('> ERROR: No such file: %s\n',spmfile); diary off
        cd(STARTINGDIR)
        return
    end

    % connames, connums, P
    diary(diaryname), fprintf('> Loading contrast names from:\n> \t%s\n',spmfile); diary off
    tcons = strmatch('T', char(SPM.xCon.STAT));
    EXPT.SNPM.connames = strvcat(SPM.xCon(tcons).name); %#ok

    ncons = size(EXPT.SNPM.connames,1);
    EXPT.SNPM.connums = 1:ncons;
    EXPT.SNPM.P = {};
    for c = 1:ncons
        confiles = {};
        for s = 1:numel(sublevs)
            confilename = fullfile(sublevs{s}, sprintf('con_%04d.nii', tcons(c)));
%             confile = filenames(confilename, 'char', 'absolute');
            confile=confilename;
            if exist(confile,'file')
                confiles{s} = confile; %#ok
            else
                diary(diaryname), fprintf('> ERROR: No such contrast file: %s\n',confilename); diary off
                cd(STARTINGDIR)
                return
            end
        end
        EXPT.SNPM.P{c} = strvcat(confiles); %#ok
    end

    %% RUN ROBFIT
    diary(diaryname), fprintf('> Group level analysis directory:\n> \t%s\n',grpmodeldir); diary off
    cd(grpmodeldir);

    save EXPT EXPT

    % which contrasts should be run
    if ~exist('includedcons','var')
        includedcons = [];
        for c = 1:numel(EXPT.SNPM.connums)
            if OPTS.overwrite || numel(filenames(fullfile(pwd,sprintf('robust%04d', c),'rob_p_*img'))) < size(EXPT.cov,2)+1
                includedcons = [includedcons, c]; %#ok
            end
        end
    end

    if isempty(includedcons)
        diary(diaryname), fprintf('> No robust regressions to run.\n'); diary off
    else
        diary(diaryname)
        fprintf('> Subject level contrasts:\n');
        for i=1:size(EXPT.SNPM.connames,1)
            fprintf('> ');
            if ~sum(find(includedcons == i)), fprintf('*EXCLUDED*  '); end
            fprintf(' %-4d%s\n',i,EXPT.SNPM.connames(i,:));
        end
        fprintf('> Group level predictors:\n');
        for i=1:numel(EXPT.covnames)
            fprintf('>  %-4d%s\n',i,EXPT.covnames{i});
        end
        diary off

        njobs = numel(includedcons) * (1+numel(EXPT.covnames));
        if ~isempty(OPTS.parallel_dream) && njobs>1
            % define scheduler
            sched = findResource('scheduler', 'type', 'generic');
            set(sched,'ClusterSize', 1);
            set(sched,'DataLocation',wd);
            switch OPTS.parallel_dream
                case '2009a'
                    set(sched,'ClusterMatlabRoot','/usr/local/matlab/R2009a');
                    set(sched,'ParallelSubmitFcn',@sgeParallelSubmitFcn);
                    set(sched,'SubmitFcn',@sgeSubmitFcn);
                    set(sched,'DestroyJobFcn', @sgeDestroyJobFcn);
                    set(sched,'DestroyTaskFcn',@sgeDestroyTaskFcn);
                case '2011a'
                    set(sched,'ClusterMatlabRoot','/usr/local/matlab/R2011a');
                    set(sched,'ParallelSubmitFcn',@parallelSubmitFcn);
                    set(sched,'SubmitFcn',@distributedSubmitFcn);
                    set(sched,'DestroyJobFcn', @destroyJobFcn);
                    set(sched,'DestroyTaskFcn',@destroyTaskFcn);
            end
            set(sched,'ClusterOsType','unix');
            set(sched,'HasSharedFileSystem', true);
            nargsout = 1;

            % submit jobs
            diary(diaryname), fprintf('> %s\tSubmitting jobs to cluster\n',datestr(now,31)); date; diary off
            for c = 1:numel(includedcons)
                %                 save(fullfile(wd,sprintf('env_%04d',c)),'EXPT','includedcons','OPTS','grpmodeldir','pthresh');
                save(fullfile(wd,sprintf('env_%04d',c)),'EXPT','includedcons','OPTS','grpmodeldir','STARTINGDIR');
                j = sched.createJob();
                set(j,'PathDependencies',cellstr(path));
                createTask(j, str2func('canlab_glm_group_levels_run1input'), nargsout, {wd c});
                alltasks = get(j, 'Tasks');
                set(alltasks, 'CaptureCommandWindowOutput', true);
                diary(diaryname), fprintf('> Submitting analysis for contrast number %d\n',includedcons(c)); diary off
                submit(j);
            end

            % wait
            diary(diaryname), fprintf('> WAITING for jobs to stop running (they may run for a while)\n'); diary off
            t1 = clock;
            notdone = true;
            fprintf('Elapsed minutes:         ')
            while notdone
                notdone = false;
                for i=1:8, fprintf('\b'); end; fprintf('%8.1f',etime(clock,t1)/60);
                pause(10)
%                 jobstatefiles = filenames(sprintf('%s/Job*.state.mat',wd),'absolute');
                jobstatefiles = filenames(sprintf('%s/Job*.state.mat',wd));
                for s = 1:numel(jobstatefiles)
                    jobstate = textread(jobstatefiles{s},'%s');
                    if strcmp(jobstate,'running') || strcmp(jobstate,'queued'), notdone = true; break; end
                end
            end
            diary(diaryname), fprintf('> \n> %s\tJobs done running\n',datestr(now,31)); diary off

            % handle output arguments/streams
            diary(diaryname), fprintf('> GATHERING job outputs\n'); diary off
%             jobdirs = filenames(sprintf('%s/Job*[0-9]',wd),'absolute');
            jobdirs = filenames(sprintf('%s/Job*[0-9]',wd));
            for i = 1:numel(jobdirs)
                % load state and output
                jobin = load(sprintf('%s/Task1.in.mat',jobdirs{i}));
                jobout = load(sprintf('%s/Task1.out.mat',jobdirs{i}));

                % parse output arguments
                robfitstatuses(jobin.argsin{2}) = jobout.argsout{1}; %#ok
%                 grfstatuses(jobin.argsin{2}) = jobout.argsout{2}; %#ok

                % output stream
                diary(fullfile(wd,sprintf('cmdwnd_%04d.txt',jobin.argsin{2}))), regexprep(jobout.commandwindowoutput,'\b[^\n]*\n',' LINES WITH BACKSPACES OMITTED\n'), diary off
            end
            % merge statuses
            if any(robfitstatuses==-1), robfitstatus = -1; else robfitstatus = 1; end
%             if any(grfstatuses==-1), grfstatus = -1; else grfstatus = 1; end

            % output stream
            if ispc()
                [ignore ignore] = system(sprintf('type %s/cmdwnd_*txt > %s',wd,fulldiaryname));
            else
                [ignore ignore] = system(sprintf('cat %s/cmdwnd_*txt > %s',wd,fulldiaryname));
            end

%             system(sprintf('cat %s/cmdwnd_*txt > %s',wd,fulldiaryname)); %#ok

            % merge diaries
            [ignore ignore] = system(sprintf('grep ''^> '' %s >> %s',fulldiaryname,diaryname)); %#ok
        else
            cmd = 'EXPT = robfit(EXPT, includedcons, 0, EXPT.mask);';
            diary(diaryname), disp(cmd); diary off
            try
                eval(cmd)
                robfitstatus = 1;
            catch exc
                if OPTS.nocatch, cd(STARTINGDIR); rethrow(exc);
                else diary(diaryname), fprintf('> %s\n',getReport(exc,'extended')), diary off; end
                robfitstatus = -1;
            end
        end
    end

    %% post processes
    % write contrast names to textfiles
    if OPTS.run_printconnames
        subconnamefile = fullfile(grpmodeldir,'subconnames.txt');
        grpconnamefile = fullfile(grpmodeldir,'grpconnames.txt');
        diary(diaryname), announce_string('WRITING CONTRAST NAMES FILES'), diary off
        diary(diaryname), fprintf('> %s\n> %s\n',subconnamefile,grpconnamefile), diary off

        fid = fopen(subconnamefile,'w');
        for i = 1:size(EXPT.SNPM.connames,1)
            fprintf(fid,'> %d %s\n',i,deblank(EXPT.SNPM.connames(i,:)));
        end
        fclose(fid);

        fid = fopen(grpconnamefile,'w');
        for i = 1:numel(EXPT.covnames)
            fprintf(fid,'> %d %s\n',i,EXPT.covnames{i});
        end
        fclose(fid);
    end

    % estimate significant cluster sizes
    %     if ~OPTS.parallel_dream && OPTS.run_grf
    if OPTS.run_grf
        diary(diaryname), announce_string('SIGNIFICANT CLUSTER SIZE ESTIMATION using GRF'), diary off
        graymattermask = which('scalped_avg152T1_graymatter_smoothed.img'); %#ok used in eval
        try
%             robdirs = filenames(fullfile(grpmodeldir,'robust[0-9][0-9][0-9][0-9]'),'absolute');
            robdirs = filenames(fullfile(grpmodeldir,'robust[0-9][0-9][0-9][0-9]'));

            for i = 1:numel(robdirs)
                load(fullfile(robdirs{i},'SETUP.mat'));

                cmd = 'sigclext = estimate_cluster_extent(.05, pthresh, SETUP.files, ''mask'', graymattermask);';
                diary(diaryname), fprintf('> %s\n',cmd), diary off
                eval(cmd);
                close all

                fout = fullfile(robdirs{i},'significant_cluster_extents_grf.txt');
                dlmwrite(fout,[pthresh' sigclext(:,1)],'precision','%g','delimiter',' '); %#ok
            end

            grfstatus = 1;
        catch exc
            if OPTS.nocatch, cd(STARTINGDIR); rethrow(exc)
            else diary(diaryname), fprintf('> %s\n',getReport(exc,'extended')), diary off; end
            grfstatus = -1;
        end
    end

    % make invp maps
    if OPTS.run_invp
        cmd = sprintf('robust_results_make_invp_maps(''d'',''%s'');',grpmodeldir);
        diary(diaryname), announce_string('MAKING INVERSE P MAPS'), diary off
        try
            eval(cmd);
            invpstatus = 1;
        catch exc
            if OPTS.nocatch, cd(STARTINGDIR); rethrow(exc);
            else diary(diaryname), fprintf('> %s\n',getReport(exc,'extended')), diary off; end
            invpstatus = -1;
        end
    end

    % make named links to robust directories
    if OPTS.run_links
        diary(diaryname), announce_string('MAKING NAMED LINKS TO ROBUST DIRECTORIES'), diary off
        try
            % get grpmodeldir without /../
            if ismac
                [dummy absgrpmodeldir] = system(sprintf('cd %s; pwd -P',grpmodeldir)); %#ok ~ doesn't work in older matlabs
            else
                [dummy absgrpmodeldir] = system(sprintf('readlink -f %s',grpmodeldir)); %#ok ~ doesn't work in older matlabs
            end
            absgrpmodeldir = deblank(absgrpmodeldir);
            % make output directory
            linkdir = fullfile(absgrpmodeldir,'named_robust_directories');
            mkdir(linkdir);
            % make links
            for i = 1:size(EXPT.SNPM.connames,1)
                % get absolute path of robust directory (with no /../)
                robdir = fullfile(absgrpmodeldir, sprintf('robust%04d',i));

                % clean up contrast name for filesystem
                linkname = deblank(EXPT.SNPM.connames(i,:));
                linkname = regexprep(linkname,'[()]','');
                linkname = regexprep(linkname,'[^0-9A-Za-z-_]','_');
                linkname = fullfile(linkdir,linkname);

                % make link
                diary(diaryname)
                fprintf('> %20s: %s\n','robust directory',robdir)
                fprintf('> %20s: %s\n','named link',linkname)
                eval(sprintf('!ln -s %s %s',robdir,linkname))
                diary off
            end

            linkdirstatus = 1;
        catch exc
            if OPTS.nocatch, cd(STARTINGDIR); rethrow(exc);
            else diary(diaryname), fprintf('> %s\n',getReport(exc,'extended')), diary off; end
            linkdirstatus = -1;
        end
    end

end

%% PUBLISH RESULTS
diary(diaryname), announce_string('PUBLISH RESULTS'), diary off
if ~OPTS.run_results
    diary(diaryname), fprintf('> SKIPPED: turned off in options.\n'); diary off
elseif robfitstatus == -1
    diary(diaryname), fprintf('> SKIPPED: robfit failed.\n'); diary off
else
    try
        canlab_glm_publish('g', grpmodeldir);
        reviewstatus = 1;
    catch exc
        if OPTS.nocatch, cd(STARTINGDIR); rethrow(exc);
        else diary(diaryname), fprintf('> %s\n',getReport(exc,'extended')), diary off; end
        reviewstatus = -1;
    end
end

%% failure report
diary(diaryname)
announce_string('FAILURE REPORT')
fail = 0;
if robfitstatus == -1
    fprintf('> \n> WARNING: robfit failed.\n')
    fail=fail+1;
end
if invpstatus == -1
    fprintf('> \n> WARNING: inverse p map making failed.\n')
    fail=fail+1;
end
if grfstatus == -1
    fprintf('> \n> WARNING: cluster size estimation failed.\n')
    fail=fail+1;
end
if linkdirstatus == -1
    fprintf('> \n> WARNING: named linking to robust directories failed.\n')
    fail=fail+1;
end
if reviewstatus == -1
    fprintf('> \n> WARNING: review publishing failed.\n')
    fail=fail+1;
end

if fail==0, fprintf('> \n> RAN WITH NO PROBLEMS (or at least so it seems).\n'); end
diary off

%% email notification
if exist('address','var')
    diary(diaryname)
    try
        [ignore output] = system(sprintf('printf "canlab_glm_group_levels has finished running.\nDirectory: %s\n\nLog file: %s\nFailed parts: %d\n" | mail -v -s "canlab_glm_group_levels done" %s',grpmodeldir,diaryname,fail,address)); %#ok
    catch %#ok
        fprintf('The email notification job failed.\n');
    end
    diary off
end

%% clean up
diary(diaryname), fprintf('> \n> \n> FINISHED: %s\n',datestr(now,31)), diary off

cd(STARTINGDIR)

end

% -------------------------------------------------------------------------
%  SUBFUNCTIONS  ----------------------------------------------------------
% -------------------------------------------------------------------------

function announce_string(string)

s = sprintf('--  %s  --',string);
l = regexprep(s,'.','-');
fprintf('> \n> \n> \n> %s\n> %s\n> %s\n> \n',l,s,l);

end

%% From Kathy Pearson:
% replace each first instance of SPM-like output backspace with newline;
% ignore additional backspaces found in sequence
%
% function [wrapstr] = nobackspace(str)
%
% wrapstr = str;
% i = strfind(wrapstr, 8);
% if ~isempty(i)
%     k = 0;
%     n = length(str);
%     first8 = 1;
%     for j = 1:n
%         if str(j) == 8
%             if first8
%                 k = k + 1;
%                 wrapstr(k) = 10;
%                 first8 = 0;
%             end
%         else
%             k = k + 1;
%             wrapstr(k) = str(j);
%             first8 = 1;
%         end
%     end
%     wrapstr = wrapstr(1:k);
% end
%
% end


text:
%GETFEATDOM Get feature domain of a dataset
%
%	FEATDOM = GETFEATDOM(A,K)
%
% Returns the cell array with the feature domains of the dataset A.

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function varargout = integral2(varargin)
%INTEGRAL2  Double integral of a DISKFUN over its domain.
%   I = INTEGRAL2(F) returns a value representing the double integral of a
%   DISKFUN.
%
%
%   I = INTEGRAL2(F, [a b c d]) integrates F in polar coordinates
%   over the region [a b] x [c d], where a and b are angular values
%   in the interval [-pi  pi], and  c and d are radial values in the
%   interval [0 1].
%
% See also DISKFUN/INTEGRAL, DISKFUN/SUM2, DISKFUN/QUAD2D.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

[varargout{1:nargout}] = quad2d(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function Y = variantMultiply( psfMatData, X, padsize );
%
%           Y = variantMultiply( psfMatData, X, padsize );
%
%  This function computes the multiplication of a spatially variant
%  point spread function (PSF) times an image:
%                y = A*x
%
%  Here we assume A is made up of a several space invariant PSFs, and
%  use piece-wise constant interpolation of them to define the spatially
%  variant PSF.  That is, A has the form:
%
%       A = D1*A1 + D2*A2 + ... + Dp*Ap
%
%  Input:
%   psfMatData  -  cell array containing the matrix data of each of the
%                  individual PSFs.  This matrix data is usuall computed
%                  from onePsfMatrix.m
%            X  -  array containing the image to which the psfMatrix
%                  is to be multiplied.
%
%  Output:
%            Y  -  contains the result after PSF multiplication.
%

%  J. Nagy 1/7/02

imsize = size( X ) - 2*padsize;

%
%  We partition the image domain into regions of equal sizes,
%  according to the number of PSFs we have ...
%
nregions = size(psfMatData);
rsize = ceil(imsize ./ nregions);

%
%  In order for this to be consistent for 2-D and 3-D images, we need to make
%  sure there is a third dimension ...
%
if length(imsize) == 1
  imsize = [imsize, 1, 1];
  rsize = [rsize, 1, 1];
  nregions = [nregions, 1, 1];
  padsize = [padsize, 0, 0];
elseif length(imsize) == 2
  imsize = [imsize, 1];
  rsize = [rsize, 1];
  nregions = [nregions, 1];
  padsize = [padsize, 0];
end

%
%  Coding the rest of this will be easier if all of the image subregions
%  have the same dimensions.  If it's not, we pad with a few zeros to make
%  it so ...
%
padsize1 = rsize .* nregions - imsize;
if any( padsize1 < 0 )
  error('Something is wrong here ...')
end
X = padarray(X, padsize1, 'post');

%
%  Now we get information about beginning and ending indices of subregions
%  so we can "put" and "get" subregions correctly ...
%
[RIidx, RJidx, RKidx] = region_indices( nregions, rsize );
[EIidx, EJidx, EKidx] = eregion_indices( RIidx, RJidx, RKidx, 2*padsize );

%
%  Now loop over all the subregions ...
%
Y = zeros(imsize);
for k = 1:nregions(3)
  for j = 1:nregions(2)
    for i = 1:nregions(1)
      Xt = X(EIidx(i,1):EIidx(i,2), EJidx(j,1):EJidx(j,2), EKidx(k,1):EKidx(k,2));
      Yt = invariantMultiply( psfMatData{i,j,k}, Xt, padsize(1:length(size(X))) );
      Y(RIidx(i,1):RIidx(i,2), RJidx(j,1):RJidx(j,2), RKidx(k,1):RKidx(k,2)) = Yt;
    end
  end
end

Y = Y(1:imsize(1), 1:imsize(2), 1:imsize(3));

repo: IRtools


text:
%GOBJ_NACA Create NACA 4-series wing profile geometry object.
%
%   [ GOBJ ] = GOBJ_NACA( SERIES, P0, L, TH, N, TAG ) Creates a NACA 4-series
%   wing profile geometry object. Accepts the following input parameters.
%
%       Parameter   Value/{Default}           Description
%       -----------------------------------------------------------------------------------
%       series      string  {0012}             NACA 4-series identifier
%       p0          array   {[0,0]}            Coordinates of leading edge/point
%       l           length  {1}                Length of wing chord
%       th          scalar  {0}                Angle of attack (degrees)
%       n           integer {20}               Number of points to define shape
%       tag         string  {N1}               Geometry object tag/name

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function [W V] = decomp(U, wtype)
%DECOMP Decomposing a matrix with orthonormal columns
%	[W V] = DECOMP(U, wtype)
%
%	U     - matrix with orthogonal column vectors ( U' * U = I )
%	wtype - decomposition type (weight type) (default: 'snnn')
%
%	W * V = U
%	W     - "weight" matrix
%	V     - "vertex" matrix
%
%	wtype parameter controls the type of the column vectors in the
%	resulting W matrix.
%	If W is a "hull" type matrix it means that its row sums are 1 and each
%	element is positive ie.: W*ones == ones && all(all(W >= 0)), in this
%	case the polytope given by the rows of V contains the rows of U (the
%	convex hull of rows(V) contains the convex hull of rows(U)).
%
%	possible values of wtype:
%		'eye': hull type, W=eye(size(U,1)), V=U
%		'ortho': W orthogonal, W=U, V=eye(size(U,2))
%		'snnn' or 'hull': hull type (uses the least possible vertices)
%		'cno': 'snnn' and each column of W comes close to 1 (tight hull)
%		'irno': 'snnn' and each column of W comes close to 0
%		'box': perform box_decomp
%
%	See also HOSVD, SNNN_DECOMP, BOX_DECOMP.

% TODO: fix cno,..
% TODO: different types..

if nargin <= 1
	wtype = 'close';
end

if ~strcmp(wtype,'eye') && ~strcmp(wtype,'ortho')
	% using an affine subspace for convex decompositions
	[n1 n2] = size(U);
	Ushift = mean(U);
	[Uu Su Vu] = svd(U - ones(n1,1)*Ushift, 'econ');
	sv = diag(Su);
	ns = sum(sv > 1e-5); % TODO
	Su = Su(1:ns,1:ns);
	Uu = Uu(:,1:ns);
	Vu = Vu(:,1:ns);

	if ns < n2
		U = Uu;
	end
end

% TODO: simplify
if nargout < 2
	switch wtype
		case 'close'
			W = close_decomp(U);
		case {'snnn', 'hull'}
			W = snnn_decomp(U);
		case 'cnoy'
			W = snnn_decomp(U);
			W = no(W);
		case 'inov'
			W = snnn_decomp(U);
			W = ino(W);
		case 'irno'
			W = snnn_decomp(U);
			W = rnoino(W);
		case 'cno'
			W = snnn_decomp(U);
			% W = cno(W, 0.5, 4, 20, 3, 3);
%			W = cno(W, 0.5, 5, 50, 10, 10);
			W = cno(W, 1, 5, 50, 15, 10);
		case 'eye'
			W = eye(size(U,1));
		case 'ortho'
			W = U;
		case 'box'
			W = box_decomp(U);
		otherwise
			error('unknown wtype');
	end
else
	switch wtype
		case 'close'
			[W V] = close_decomp(U);
		case {'snnn', 'hull'}
			[W V] = snnn_decomp(U);
		case 'cnoy'
			[W V1] = snnn_decomp(U);
			[W V2] = no(W);
			V = V2*V1;
		case 'inov'
			[W V1] = snnn_decomp(U);
			[W V2] = ino(W);
			V = V2*V1;
		case 'irno'
			[W V1] = snnn_decomp(U);
			[W V2] = rnoino(W);
			V = V2*V1;
		case 'cno'
			[W V1] = snnn_decomp(U);
			% [W V2] = cno(W, 0.5, 4, 20, 3, 3);
			[W V2] = cno(W, 1, 5, 50, 15, 10);
			V = V2*V1;
		case 'eye'
			W = eye(size(U,1));
			V = U;
		case 'ortho'
			W = U;
			V = eye(size(U,2));
		case 'box'
			[W V] = box_decomp(U);
		otherwise
			error('unknown wtype');
	end

	% shift back (for convex decompositions)
	if ~strcmp(wtype,'eye') && ~strcmp(wtype,'ortho') && ns < n2
		V = V*Su*Vu' + ones(size(V,1),1)*Ushift;
	end
end

author: Sable

repo: mcbench-benchmarks


text:
function dfdx = shadowjacobian(f,x)
% See SDPVAR/jacobian

if isa(f,'double')
    dfdx = zeros(size(f,1),length(x));
    return
end

if ~isempty(intersect(deepdepends(f),depends(x)))
    % Under development
end

if nargin==1
    if isa(f,'sdpvar')
        x = recover(depends(f));
    else
        x = 0;
    end
else
    if length(getvariables(x))<length(x)
      error('x should be a vector of scalar independant variables');
    end
end

[n,m]=size(f);

if m>1
   error('Jacobian only defined for column vectors.')
end

if n*m==1
    dfdx = scalar_jacobian(f,x);
    % Argh, fix this (sorts inside scalar_jacobian
    for i = 1:length(x)
        var(i,1)=getvariables(x(i));
    end
    [i,j]=sort(var);
    dfdx = dfdx(1,j);
    return

else
    dfdx = [];
    AllVars = recover(unique([depends(f) getvariables(x)]));

    for i = 1:length(f)
        dfdx = [dfdx;scalar_jacobian(f(i),x,AllVars)];
    end
    % Argh, fix this (sorts inside scalar_jacobian
    for i = 1:length(x)
        var(i,1)=getvariables(x(i));
    end
    [i,j]=sort(var);
    dfdx = dfdx(:,j);
end

function [dfdx,dummy] = scalar_jacobian(f,x,AllVars)

if isa(f,'double')
    dfdx = zeros(1,length(x));
    return
end

if nargin==2
    AllVars = recover(uniquestripped([depends(f) getvariables(x)]));
    %AllVars = recover(uniquestripped([deepdepends(f) depends(f) getvariables(x)]));
end

exponent_p = exponents(f,AllVars);

coefficients = getbase(f);
coefficients = coefficients(2:end);
coefficients = coefficients(:);
if nnz(exponent_p(1,:))==0
    exponent_p=exponent_p(2:end,:);
end

x_variables = getvariables(x);
AllVars_variables = getvariables(AllVars);
%AllDeriv = [];
AllDeriv2 = [];
for k = 1:length(x)
    wrt = find(ismembc(AllVars_variables,x_variables(k)));
    deriv = exponent_p;
    deriv(:,wrt) = deriv(:,wrt)-1;
    keep{k} = find(deriv(:,wrt)~=-1);
    AllDeriv2 = [AllDeriv2;deriv(keep{k},:)];
end

if size(AllDeriv2,1)==0
    dummy = 1;
else
    dummy = recovermonoms(AllDeriv2,AllVars);
end

top = 1;
dfdx=[];
for k = 1:length(x)
    wrt = find(ismembc(AllVars_variables,x_variables(k)));
    m = length(keep{k});
    if m>0
        poly = sum((coefficients(keep{k}(:)).*exponent_p(keep{k}(:),wrt)).*dummy(top:top+m-1),1);
        top = top + m;
    else
        poly = 0;
    end
    dfdx = [dfdx ; poly];
end
dfdx = dfdx';

author: shamilmamedov


text:
function F = stoploop(str)
% STOPLOOP - creates stop button to have a user interrupt a loop
%
%  FS = STOPLOOP creates a message box window and returns a structure FS that
%  holds two functions, called FS.Stop and FS.Clear. The function FS.Stop()
%  will return true, if the OK button has been clicked (or the message box
%  has been removed), so that a loop can be interrupted.
%  The function FS.Clear() can be used to remove the message box, if a loop
%  has ended without user interruption.
%
%  FS = STOPLOOP(STR) uses the string STR to display instead of the default
%  'Stop the Loop'.
%
%  Example:
%    tic ;         % We will measure elapsed time in a loop
%                  % Set up the stop box:
%    FS = stoploop({'Stop me before', '5 seconds have elapsed'}) ;
%                  % Display elapsed time
%    fprintf('\nSTOPLOOP: elapsed time (s): %5.2f\n',toc)
%                  % start the loop
%    while(~FS.Stop() && toc < 5),       % Check if the loop has to be stopped
%        fprintf('%c',repmat(8,6,1)) ;   % clear up previous time
%        fprintf('%5.2f\n',toc) ;        % display elapsed time
%    end
%    FS.Clear() ;  % Clear up the box
%    clear FS ;    % this structure has no use anymore
%
%  Notes:
%  - The function call F.Stop() issues a drawnow command.
%  - The function call F.Clear() does nothing when the message box already
%    has been cleared.
%
%  See also MSGBOX, KEYBOARD, DBSTOP, WHILE, FOR
%           GETKEYNOW, GETKEYWAIT, GETKEY (on the File Exchange)

% for Matlab R14+
% version 1.0 (jun 2008)
% (c) Jos van der Geest
% email: jos@jasen.nl

% History:
% 1.0 (jun 2008) Inspired by several requests over the years on CSSM

error(nargoutchk(1,1,nargout)) ;

if nargin,
    if ~ischar(str) && ~iscellstr(str),
        error([mfilename ':InputString'],...
            'Input should be a string, or a cell array of strings.') ;
    end
else
    % default message string
    str = 'Stop the Loop' ;
end

% create a msgbox displaying the string
H = msgbox(str,'STOPLOOP') ;

% create the two anonymous functions
F.Stop = @() stopfun(H) ; % false if message box still exists
F.Clear = @() clearfun(H) ; % delete message box

function r = stopfun(H)
drawnow ;          % ensure that button presses are recorded
r = ~ishandle(H) ; % false if message box still exists

function clearfun(H)
% clear the message box if it still exists
if ishandle(H),
    delete(H) ;
end


text:
function variance = logistic_variance ( a, b )

%*****************************************************************************80
%
%% LOGISTIC_VARIANCE returns the variance of the Logistic PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < B.
%
%    Output, real VARIANCE, the variance of the PDF.
%
  variance = ( pi * b )^2 / 3.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function B=nancumsum(A,dim,nmode)
% NANCUMSUM: Cumulative sum of a matrix, with user-specified treatment of NaNs.
%     Computes the cumulative sum of matrix A along dimension DIM, allowing
%     the user to replace NaNs with zeros, to skip over them, or to reset
%     on NaNs, maintaining NaNs as placeholders.
%
% USAGE: B = nancumsum(A, DIM, NMODE)
%
% ARGUMENTS:
%
% A:    Input matrix.
%
% B:    Output cumulative sum matrix, treating NaNs as determined by nmode.
%
% DIM:  B = nancumsum(A, DIM) returns the nan-cumulative sum of the elements
%       along the dimension of A specified by scalar DIM. For example,
%       nancumsum(A,1) works down the columns, nancumsum(A,2) works
%       across the rows. If DIM is not specified, it defaults to the first
%       non-singleton dimension of A.
%
% NMODE: specifies how NaNs should be treated. Acceptable values are:
%       1: REPLACE NaNs with zeros (default).
%       2: MAINTAIN NaNs as position holders in B. (Skip NaNs without reset.)
%       3: RESET sum on NaNs, replacing NaNs with zeros.
%       4: RESET sum on NaNs, maintaining NaNs as position holders.
%
% EXAMPLES:
%
% 1) a = [NaN,2:5];
%
% nancumsum(a)
% ans =
%     0     2     5     9    14
%
% nancumsum(a,[],2)
% ans =
%   NaN     2     5     9    14
%
% nancumsum(a,[],3)
% ans =
%     2     5     9    14
%
% 2) a = magic(3); a(5)=NaN;
%
% b = nancumsum(a,2) % (Default NMode = 1)
% b =
%     8     9    15
%     3     3    10
%     4    13    15
%
% b = nancumsum(a,2,2)
% b =
%     8     9    15
%     3   NaN    10
%     4    13    15
%
% b = nancumsum(a,2,3)
% b =
%     8     9    15
%     3     0     7
%     4    13    15
%
% b = nancumsum(a,2,4)
% b =
%     8     9    15
%     3   NaN     7
%     4    13    15

% See also: cumsum, nansum, nancumprod, nanmean, nanmedian, ...
% (nancumprod is available from the FEX. Other nan* may require Toolboxes)

% Brett Shoelson
% brett.shoelson@mathworks.com
% 05/04/07
%
% Revision: 08/28/11
% Fixed bug in option 2 (faulty reset). Thanks to Andrew Stevens and Rick
% Patterson for reporting it. Also, eliminated old option 3 (deleting NaNs
% in a vector) as a trivial case and added two new options.
% Revision: 09/15/11
% Fixed a bug with multiple NaNs. Thanks to Tim Yates.
%
% Copyright The MathWorks, Inc. 2011

% Set defaults, check and validate inputs
if nargin < 3
    nmode = 1;
end

if ~ismember(nmode,1:4)
    error('NANCUMSUM: unacceptable value for nmode parameter.');
end

if nargin < 2 || isempty(dim)
    if ~isscalar(A)
        dim = find(size(A)>1);
        dim = dim(1);
    else
        % For scalar inputs (no nonsingleton dimension)
        dim = 1;
    end
end

% Calculate cumulative sum, depending on selection of nmode
switch nmode
    case 1
        % TREAT NaNs as 0's
        B = A;
        B(B~=B) = 0;
        B = cumsum(B, dim);
    case 2
        % DO NOT INCREMENT, BUT USE NaNs AS PLACEHOLDERS.
        B = nancumsum(A, dim, 1);
        B(A~=A) = NaN;
     case 3
        % RESET sum on NaNs, replacing NaNs with zeros.
        naninds = find(A~=A);
        for ii = 1:numel(naninds)
            B = nancumsum(A, dim, 1);
            A(naninds(ii)) = -B(naninds(ii));
        end
        B = cumsum(A,dim);
    otherwise %case 4
        % RESET sum on NaNs, maintaining NaNs as position holders.
        naninds = find(A~=A);
        B = nancumsum(A,dim,3);
        B(naninds)= NaN;
end


text:
function i4_rise_values_test ( )

%*****************************************************************************80
%
%% I4_RISE_VALUES_TEST demonstrates the use of I4_RISE_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 December 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4_RISE_VALUES_TEST:\n' );
  fprintf ( 1, '  I4_RISE_VALUES returns some exact values\n' );
  fprintf ( 1, '  of the integer rising factorial function:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     M     N      I4_RISE(M,N)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, m, n, fmn ] = i4_rise_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %4d  %4d  %10d\n', m, n, fmn );

  end

  return
end


text:
function X = expanded(F,state)

F = flatten(F);
if nargin == 1
    if length(F.LMIid) == 0
        X = [];
    else
        for i = 1:length(F.clauses)
            X(i,1) = F.clauses{i}.expanded;
        end
    end
else
    X = F;
    for i = 1:length(F.clauses)
        X.clauses{i}.expanded = state;
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function display(X)
%DISPLAY (overloaded)

switch(X.typeflag)
    case {0,9}
        n = X.dim(1);
        m = X.dim(2);
        if (n*m==1)

            linearbilinearquadraticsigmonial = is(X,'LBQS');
            if linearbilinearquadraticsigmonial(1)
                classification = 'Noncommuting linear scalar ';
            elseif linearbilinearquadraticsigmonial(4)
                classification = 'Noncommuting sigmonial scalar ';
            elseif linearbilinearquadraticsigmonial(2)
                classification = 'Noncommuting bilinear scalar ';
            elseif linearbilinearquadraticsigmonial(3)
                classification = 'Noncommuting quadratic scalar ';
            else
                classification = 'Noncommuting polynomial scalar ';
            end

            if ~isreal(X.basis)
                classification = [classification '(complex'];
            else
                classification = [classification '(real'];
            end

            if is(X,'compound')
                classification = [classification ', derived'];
            end

            if ~islinear(X)
                variables = getvariables(X);
                monomtable = yalmip('monomtable');
                if ((nnz(getbasematrix(X,0))==0) ) & (sum(diff(sum(monomtable(variables,:),2)))==0)
                    classification = [classification ', homogeneous'];
                end
            end
            if any(ismember(depends(X),yalmip('intvariables')))
                classification = [classification ', integer'];
            else
                if any(ismember(depends(X),yalmip('binvariables')))
                    classification = [classification ', binary'];
                else
                    if any(ismember(depends(X),yalmip('uncvariables')))
                        classification = [classification ', uncertain'];
                    end
                end
            end

            nvars = length(depends(X));
            if nvars == 1
                classification = [classification ', ' num2str(nvars) ' variable'];
            else
                classification = [classification ', ' num2str(nvars) ' variables'];
            end

            if any(ismember(depends(X),yalmip('parvariables')))
                classification = [classification ', parametric'];
            end
            if ~isnan(double(X))
                classification = [classification ', current value : ' num2str(double(X))];
            end
            classification = [classification ')'];

            disp([classification]);
        else

            if islinear(X)
                classification = 'Noncommuting linear matrix variable ';
            elseif is(X,'sigmonial')
                classification = 'Noncommuting sigmonial matrix variable ';
            elseif is(X,'bilinear')
                classification = 'Noncommuting bilinear matrix variable ';
            elseif is(X,'quadratic')
                classification = 'Noncommuting quadratic matrix variable ';
            else
                classification = 'Noncommuting polynomial matrix variable ';
            end

            if isreal(X.basis)
                if issymmetric(X)
                    info = ' (symmetric, ';
                else
                    info = ' (full, ';
                end;
                info = [info 'real'];
            else
                if issymmetric(X)
                    info = ' (symmetric, ';
                elseif ishermitian(X)
                    info = ' (hermitian, ';
                else
                    info = ' (full, ';
                end;
                info = [info 'complex'];
            end;

            if is(X,'compound')
                info = [info ', derived'];
            end

            if X.typeflag==9
                info = [info ', KYP'];
            end
            if any(ismember(depends(X),yalmip('intvariables')))
                info = [info ', integer'];
            else
                if any(ismember(depends(X),yalmip('binvariables')))
                    info = [info ', binary'];
                else
                    if any(ismember(depends(X),yalmip('uncvariables')))
                        info = [info ', uncertain'];
                    end
                end
            end
            if any(ismember(depends(X),yalmip('parvariables')))
                info = [info ', parametric'];
            end

            xvars = depends(X);
            nvars = length(xvars);
            if nvars == 1
                info = [info ', ' num2str(nvars) ' variable'];
            else
                info = [info ', ' num2str(nvars) ' variables'];
            end

            n = X.dim(1);
            m = X.dim(2);

            if (n<100) & (n==m)
                x = recover(xvars);
                if ~any(any(isnan(double(x))))
                    doubleX = double(X);
                    try
                    eigX = eig(doubleX);
                    info = [info ', eigenvalues between [' num2str(min(eigX)) ',' num2str(max(eigX)) ']'];
                    catch
                    end
                end
            elseif n~=m
                x = recover(xvars);
                if ~any(any(isnan(double(x))))
                    doubleX = double(X);
                    try
                        info = [info ', values in range [' num2str(min(min(doubleX))) ',' num2str(max(max(doubleX))) ']'];
                    catch
                    end
                end
            end

            info = [info ')'];
            disp([classification num2str(n) 'x' num2str(m) info]);
        end;
    case 1
        disp('Relational object');
    case 3
        disp('Relational object');
    case 4
        disp('Relational object');
    case 5
        disp('Cone object');
    case 6
        disp('logdet object');
    case 7
        disp('Integrality constraint');
    case 10
        disp('Eigenvalue object');
    case 11
        disp('SOS object');
    otherwise
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% Eric, please check this out.
%
%

function [signal, time, chanlabels, pnts, nchan, serror] = loadtxterp(fullname, transpose, timex, elabel)

signal=[]; time=[]; chanlabels=''; pnts=[]; nchan=[];
serror = 0; % no errors
delimiter  = '\t';

%if transpose == 0; % data matrix is point x elec  (normal, no transpose)
if elabel==1
        nheaderlines = 1;
else
        nheaderlines = 0;
end
%else   % data matrix is elec x points  (transposed. should be corrected)
%    nheaderlines = 0;  %changed from 0 to 3
%end

% try
%  if transpose==1
%         fid = fopen(fullname,'rt');
%         nLines = 0;
%         while (fgets(fid) ~= -1)
%             nLines = nLines+1;
%         end
%         fclose(fid);
%         fid = fopen(fullname,'rt');
%         [A, count]=fscanf(fid, '%s');
%         fclose(fid);
%
%         if timex==1
%         columns_of_header=count/nLines; %counts the number of channels +1 (time header)
%         num_of_col=columns_of_header-1; %stores number of columns (time points)
%         signal_trans=signal(:,2:end);
%         signal=signal_trans;
%         time= dlmread(fullname, '', [0 1 0 num_of_col]);
%         else
%             signal=dlmread(fullname, '', [0 1 nLines num_of_col]);
%             time=[];
%         end
%  end

if nheaderlines==1
        values     = importdata(fullname);
        signal     = values.data;

        if transpose==1
                signal=signal';
                if timex==1
                        time=signal(:,1);
                        signal=signal(:,2:end);
                else
                        time=[];
                end
        else
                if timex==1
                        time = signal(:,1);
                        signal = signal(:,2:end);
                else
                        time = [];
                end
        end
else
        values=importdata(fullname);
        signal=values;
        if transpose==1
                signal=signal';
                if timex==1
                        time=signal(:,1);
                        signal=signal(:,2:end);
                else
                        time=[];
                end
        else
                signal = dlmread(fullname);

                if timex==1
                        time=signal(:,1);
                        signal=signal(:,2:end);
                else
                        time=[];
                end
                if transpose==1
                        signal=signal';
                end
        end
end

%  if nheaderlines==1
%      values     = importdata(fullname);
%      signal     = values.data;
%
%
%      if timex==1
%          time = signal(:,1);
%          signal = signal(:,2:end);
%      else
%          time = [];
%      end
%
%
%      if transpose==1
%          signal = signal'; % corrected
%      end
%  else
%      signal = dlmread(fullname);
%
%      if timex==1
%          time=signal(:,1);
%          signal=signal(:,2:end);
%      else
%          time=[];
%      end
%
%      if transpose==1
%          signal=signal';
%      end
%
%  end

pnts  = size(signal,1);
nchan = size(signal,2);

%       try
if elabel==1
        chanlabels = values.textdata;
        if size(chanlabels,1)==1 && size(chanlabels,2)==1
                %                         try
                chanlabels=regexp(chanlabels,'.*?\s+', 'match');
                chanlabels = [chanlabels{:}];
                chanlabels=strtrim(chanlabels);
                %                         catch
                %                               fprintf('Oops...Please check ERP.chanlocs.labels. Is it fine?');
                %                         end
        end
        if timex==1
                chanlabels = chanlabels(~ismember_bc2(chanlabels,{'time'}));
        end
else
        chanlabels = '';
end

%       catch
%             chanlabels = '';
%       end

if isempty(char(chanlabels))
        chanlabels = cell(1);
        for e=1:nchan
                chanlabels{e} = ['Ch' num2str(e)];
        end
end

% catch
%       serror = 1; % error found
% end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function varargout = detectChange(alg_cd, config_alg_cd, alg_thre, config_alg_thre, varargin)
%detectChange Perform change detection on a given image pair or data folder.
%
% [CM, DI] = detectChange(alg_cd, config_alg_cd, alg_thre, config_alg_thre, path_t1, path_t2)
%
% returns the change map CM and the difference image DI. alg_cd specifies the name of the CD
% algorithm to use, and config_alg_cd is a cell array that contains the arguments to construct
% the CDAlg object. alg_thre specifies the name of the thresholding algorithm, with
% config_alg_thre being its corresponding configs. path_t1 and path_t2 are the absolute or
% relative paths to the input bi-temporal images.
%
% [CMs, DIs] = detectChange(alg_cd, config_alg_cd, alg_thre, config_alg_thre, data_dir)
%
% handles an image folder, which is specified by data_dir. The image folder
% should be structured as follows:
% |- DATA_DIR
% |     |- t1
% |     |- t2
% |     |_ gt
% where the gt subfolder is optional. t1 contains images acquired at the first time instance
% while t2 at the second time instance. see Datasets.Folder for more details. Note that the CMs
% and the DIs will be returned in cell arrays.
%
% [CM, DI, results] = detectChange(alg_cd, config_alg_cd, alg_thre, config_alg_thre, path_t1,
% path_t2, path_ref, metrics, configs_metrics)
%
% returns evaluation results, additionally. To enable the evaluation metrics, you will have
% to specify the path to the reference image, path_ref, and the names of the metrics to use.
% All names ought to be contained in the cell array metrics, and there should also be
% configuration cells for each metric, which are nested in the cell configs_metrics.
%
% [CMs, DIs, results] = detectChange(alg_cd, config_alg_cd, alg_thre, config_alg_thre, data_dir,
% metrics, configs_metrics)
%
% evaluates each image pair in a folder.In this case, the gt subfolder is required.
%

% Create instances
algCD = Algorithms.(alg_cd)(config_alg_cd{:});
algThre = ThreAlgs.(alg_thre)(config_alg_thre{:});
% Choose the right entrance
if nargin == 9
    % varargin being path_t1, path_t2, path_ref, metrics, and
    % config_metrics
    varargout = cell(1,3);
    [varargout{:}] = cd4Pair_(algCD, algThre, varargin{:});
elseif nargin == 7
    % varargin being dir_, metrics, and config_metrics
    varargout = cell(1,3);
    [varargout{:}] = cd4Dir_(algCD, algThre, varargin{:});
elseif nargin == 6
    % varargin being path_t1 and path_t2
    varargout = cell(1,2);
    [varargout{:}] = cd4Pair_(algCD, algThre, varargin{:}, '', {}, {});
elseif nargin == 5
    % varargin being dir_
    [varargout{1}, varargout{2}] = cd4Dir_(algCD, algThre, varargin{:}, {}, {});
else
    error('Wrong number of input arguments');
end
end

function [CMs, DIs, results] = cd4Dir_(alg_cd, alg_thre, dir_, metrics, configs)
dataset = Datasets.Folder(dir_);
iterator = Datasets.CDDIterator(dataset);

% Create metric instances
nMetrics = length(metrics);
meters = cell(1, nMetrics);
for ii = 1:nMetrics
    meters{ii} = Metrics.(metrics{ii})(configs{ii}{:});
end

DIs = {}; CMs = {};
% Do loop
while(iterator.hasNext())
    % Fetch data
    [t1, t2, ref] = iterator.nextChunk();
    % Check rank
    getRank = @(x) rank(reshape(double(x), size(x, 1)*size(x, 2), size(x, 3)));
    if getRank(t1) <= 1 || getRank(t2) <= 1
        warning('The rank of the input image is detected to be lower than 2.');
    end
    DI = alg_cd.detectChange(t1, t2);
    CM = alg_thre.segment(DI);
    if nMetrics > 0 %&& ~isempty(ref)
        % Update value
        cellfun(@(obj) obj.update(CM, ref, DI), meters);
        % Show logs
        for ii = 1:nMetrics
            m = meters{ii};
            fprintf('type: %s\n', metrics{ii});
            fprintf('\tnewest: %f\n', m.val(end));
            fprintf('\taverage: %f\n', m.avg);
        end
        fprintf('\n');
    end
    DIs = [DIs, DI];
    CMs = [CMs, CM];
end

results = struct();
for ii = 1:nMetrics
    results.(metrics{ii}) = meters{ii}.avg;
end
end

function [CM, DI, results] = cd4Pair_(alg_cd, alg_thre, path_t1, path_t2, path_ref, metrics, configs)
% Load data
imT1 = readFile__(path_t1);
imT2 = readFile__(path_t2);

% Check data redundancy
getRank = @(x) rank(reshape(double(x), size(x, 1)*size(x, 2), size(x, 3)));
if getRank(imT1) <= 1 || getRank(imT2) <= 1
    warning('The rank of the input image is detected to be lower than 2.');
end

DI = alg_cd.detectChange(imT1, imT2);
CM = alg_thre.segment(DI);

if  ~isempty(metrics)% && ~isempty(path_ref)
    ref = readFile__(path_ref);
    % Create metric instances
    nMetrics = length(metrics);
    meters = cell(1, nMetrics);
    for ii = 1:nMetrics
        meters{ii} = Metrics.(metrics{ii})(configs{ii}{:});
        meters{ii}.update(CM, ref, DI);
        results.(metrics{ii}) = meters{ii}.val;
    end
end

    function im = readFile__(path_)
        try
            % Try to read via default loader
            im = Datasets.Loaders.defaultLoader(path_);
        catch ME
            % Otherwise, read as envi file
            if strcmp(ME.identifier, 'MATLAB:imagesci:imread:fileFormat')
                im = Datasets.Loaders.enviLoader(path_);
            else
                rethrow(ME);
            end
        end
    end
end

author: Bobholamovic

repo: ChangeDetectionToolbox


text:
function [Fconv,no_changed,infeasible,Forigquad] = convertquadratics(F)
%CONVERTQUADRATICS Internal function to extract quadratic constraints

% ******************************
% LINEAR?
% ******************************

infeasible = 0;
Fconv = F;
no_changed = 0;
Forigquad = [];

if islinear(F)
    return
end

[monomtable,variabletype] = yalmip('monomtable');

if any(variabletype == 4)
    if issigmonial(F)
        return
    end
end

Fconv = lmi;
no_changed = 0;
i_changed = [];
% Make sure bounds a pre-processed, in case we look for rotated cone which
% requires us to know lower bounds
nv = yalmip('nvars');
LU = yalmip('getbounds',1:nv);
LUhere = getbounds(F,[],LU);
for i = 1:1:length(F)
    if max(variabletype(getvariables(F(i)))) <= 1
        % Definitely no quadratic to model as all variables are bilinear at
        % most
         Fconv = Fconv + F(i);
    elseif is(F(i),'element-wise') & ~is(F(i),'linear') & ~is(F(i),'sigmonial')
        % f-c'*x-x'*Q*x>0
        fi = sdpvar(F(i));fi = fi(:);
        %[Qs,cs,fs,x,info] = vecquaddecomp(fi);
        for j = 1:length(fi)
            fij = fi(j);
            if isa(fij,'double')
                if fij < 0
                    infeasible = 1;
                    warning('The problem is trivially infeasible. You have a term POSITIVE NUMBER < 0.');
                    return
                end
            end
           % Q = Qs{j};c = cs{j};f = fs{j};
           if max(variabletype(getvariables(fij))) <= 1
               % There are at most bilinear terms, so it cannot be convex
               info = 1;
           else
               [Q,c,f,x,info] = quaddecomp(fij);
           end
            if info==0
                if nnz(Q)==0
                    % Oh, linear,...
                    Fconv = Fconv + (fi(j)>=0);
                else
                    % Yes, quadratic, but convex?
                    % Change sign definitions
                    Q = -Q;
                    c = -c;
                    f = -f;
                    % Semi-definite case when only part of x in Q
                    % Occurs, e.g, in constraints like y'*Q*y < t
                    used = find(any(Q));Qred=Q(:,used);Qred = Qred(used,:);xred = x(used);
                    [R,p]=chol(Qred);
                    done = 0;
                    if p
                        % Safety check to account for low rank problems
                        if all(eig(full(Qred))>=-1e-12)
                            [u,s,v]=svd(full(Qred));
                            r=find(diag(s)>1e-12);
                            R=(u(:,r)*sqrt(s(r,r)))';
                            p=0;
                        else
                            % Try to detect rotated SOCP (x-d)'*A*(x-d)+k<= C*y*z
                            yzCandidates = find(~diag(Qred));
                            if length(yzCandidates) == 2 && nnz(c(yzCandidates))==0
                                %LU = yalmip('getbounds',getvariables(xred(yzCandidates)));
                                LU = LUhere(getvariables(xred(yzCandidates)),:);
                                if all(LU(:,1)>=0)
                                    yzSubQ = Qred(yzCandidates,yzCandidates);
                                    if yzSubQ(1,2) < 0
                                        C = -2*yzSubQ(1,2);
                                        xCandidates = setdiff(1:length(used),yzCandidates);
                                        A = Qred(xCandidates,xCandidates);
                                        [B,p]=chol(A);
                                        if ~p
                                            y = xred(yzCandidates(1));
                                            z = xred(yzCandidates(2));
                                            d = -A\(c(xCandidates)/2);
                                            k = f-d'*A*d;
                                            if abs(k)<= 1e-12
                                                Fconv=Fconv + lmi(rcone(B*(xred(xCandidates)-d),.5*C*y,z));
                                                no_changed = no_changed + 1;
                                                i_changed = [i_changed i];
                                                done = 1;
                                            elseif k >= -1e-12
                                                Fconv=Fconv + lmi(rcone([B*(xred(xCandidates)-d);sqrt(abs(k))],.5*C*y,z));
                                                no_changed = no_changed + 1;
                                                i_changed = [i_changed i];
                                                done = 1;
                                            else
                                                p = 1;
                                            end
                                        end
                                    end
                                end
                            end
                        end
                    end
                    if p==0 && ~done
                        % Write as second order cone
                        d = -c'*x-f;
                        if isa(d,'double')
                            if d<0
                                infeasible = 1;
                                return
                            else
                                Fconv=Fconv + lmi(cone([R*xred],sqrt(d)));
                            end
                        else
                            if length(c) == length(xred)
                                ctilde = -(R')\(c/2);
                                Fconv=Fconv + lmi(cone([R*xred;.5*(1-d)],.5*(1+d)));
                            else
                                Fconv=Fconv + lmi(cone([R*xred;.5*(1-d)],.5*(1+d)));
                            end
                        end
                        no_changed = no_changed + 1;
                        i_changed = [i_changed i];
                    elseif ~done
                        Fconv = Fconv + lmi(fi(j));
                    end
                end
            else
                Fconv = Fconv + lmi(fi(j));
            end
        end
    else
        Fconv = Fconv + F(i);
    end
end
if ~isempty(i_changed)
    Forigquad = F(i_changed);
else
    Fconv = F;
end


text:
clear; clf;
		% Incarcarea vetorului x
x=-2.5:0.01:2.5;
		% Evaluarea functiei date
y=funct2(x);
		% Calcularea aproximativa a derivatei
df=diff(y)./diff(x);
		% Reprezentarea grafica a functiei de derivat
plot(x,y,'g:')
hold on
		% Incarcarea vectorului xd cu n-1 elemente
		% din vectorul x
xd=x(2:length(x));
		% Reprezentarea grafica a derivatei
plot(xd,df,'m-')
		% Calcularea elementelor vectorului produs
produs=df(1:length(df)-1).*df(2:length(df));
		% Incarcarea punctelor critice in vectorul minmax
minmax=xd(find(produs<0));
		% Evaluarea functiei in punctele critice
fminmax=y(find(produs<0)+1);
		% Reprezentarea grafica a maximelor si minimelor
        % locale ale functiei initiale
plot(minmax,fminmax,'kx','MarkerSize',14)
grid on
        % Trasarea liniilor de grid prin punctele gasite
        % - obtinerea valorilor la care se plaseaza
        %   liniutele de divizare implicite
xtick=get(gca,'XTick');
        % - gasirea valorilor minime si maxime ale acestora
xtickmin=min(xtick); xtickmax=max(xtick);
        % - generarea vectorului noilor linii de divizare pe axa x
xtick=[xtickmin,minmax,xtickmax];
        % - setarea noilor liniute de divizare
set(gca,'XTick',xtick);
set(gca,'YTick',0);
    	% Plasarea unei legende
legend('Functia','Derivata');


text:
function [interp] = ft_sourceinterpolate(cfg, functional, anatomical)

% FT_SOURCEINTERPOLATE interpolates source activity or statistical maps onto the
% voxels or vertices of an anatomical description of the brain.  Both the functional
% and the anatomical data can either describe a volumetric 3D regular grid, a
% triangulated description of the cortical sheet or an arbitrary cloud of points.
%
% The functional data in the output data will be interpolated at the locations at
% which the anatomical data are defined. For example, if the anatomical data was
% volumetric, the output data is a volume-structure, containing the resliced source
% and the anatomical volume that can be visualized using FT_SOURCEPLOT or written to
% file using FT_SOURCEWRITE.
%
% The following scenarios can be considered:
%
% - Both functional data and anatomical data are defined on 3D regular grids, for
%   example with a low-res grid for the functional data and a high-res grid for the
%   anatomy.
%
% - The functional data is defined on a 3D regular grid and the anatomical data is
%   defined on an irregular point cloud, which can be a 2D triangulated surface mesh.
%
% - The functional data is defined on an irregular point cloud, which can be a 2D
%   triangulated surface mesh, and the anatomical data is defined on a 3D regular grid.
%
% - Both the functional and the anatomical data are defined on an irregular point
%   cloud, which can be a 2D triangulated mesh.
%
% - The functional data is defined on a low-resolution 2D triangulated surface mesh and the
%   anatomical data is defined on a high-resolution 2D triangulated surface mesh, where the
%   low-res vertices form a subset of the high-res vertices. This allows for mesh-based
%   interpolation. The algorithm currently implemented is so-called 'smudging' as it is
%   also applied by the MNE-suite software.
%
% Use as
%   [interp] = ft_sourceinterpolate(cfg, source, anatomy)
%   [interp] = ft_sourceinterpolate(cfg, stat,   anatomy)
% where
%   source  is the output of FT_SOURCEANALYSIS
%   stat    is the output of FT_SOURCESTATISTICS
%   anatomy is the output of FT_READ_MRI, or one of the FT_VOLUMExxx functions,
%           or a cortical sheet that was read with FT_READ_HEADSHAPE,
%           or a regular 3D grid created with FT_PREPARE_SOURCEMODEL.
%
% The configuration should contain:
%   cfg.parameter     = string or cell-array with the functional parameter(s) to be interpolated
%   cfg.downsample    = integer number (default = 1, i.e. no downsampling)
%   cfg.interpmethod  = string, can be 'nearest', 'linear', 'cubic',  'spline', 'sphere_avg', 'sphere_weighteddistance', or 'smudge' (default = 'linear for interpolating two 3D volumes, 'nearest' for all other cases)
%
% For interpolating two 3D regular grids or volumes onto each other the supported
% interpolation methods are 'nearest', 'linear', 'cubic' or 'spline'. For all other
% cases the supported interpolation methods are 'nearest', 'sphere_avg',
% 'sphere_weighteddistance' or 'smudge'.
%
% The functional and anatomical data should be expressed in the same
% coordinate sytem, i.e. either both in MEG headcoordinates (NAS/LPA/RPA)
% or both in SPM coordinates (AC/PC).
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_READ_MRI, FT_READ_HEADSHAPE, FT_SOURCEPLOT, FT_SOURCEANALYSIS,
% FT_SOURCEWRITE

% Copyright (C) 2003-2007, Robert Oostenveld
% Copyright (C) 2011-2014, Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar functional anatomical
ft_preamble provenance functional anatomical

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% this is not supported any more as of 26/10/2011
if ischar(anatomical)
  ft_error('please use cfg.inputfile instead of specifying the input variable as a sting');
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'unused',     {'keepinside' 'voxelcoord'});
cfg = ft_checkconfig(cfg, 'deprecated', {'sourceunits', 'mriunits'});
cfg = ft_checkconfig(cfg, 'required',   'parameter');
cfg = ft_checkconfig(cfg, 'renamedval', {'parameter', 'avg.pow', 'pow'});
cfg = ft_checkconfig(cfg, 'renamedval', {'parameter', 'avg.coh', 'coh'});
cfg = ft_checkconfig(cfg, 'renamedval', {'parameter', 'avg.mom', 'mom'});

% set the defaults
cfg.downsample   = ft_getopt(cfg, 'downsample', 1);
cfg.feedback     = ft_getopt(cfg, 'feedback', 'text');
cfg.interpmethod = ft_getopt(cfg, 'interpmethod', []);   % cfg.interpmethod depends on how the interpolation should be done and actual defaults will be specified below

% replace pnt by pos
anatomical = fixpos(anatomical);
functional = fixpos(functional);

% ensure the functional data to be in double precision, the maxdepth parameter ensure double precision up to the content of functional.avg.mom{:}, avoiding too much recursion
functional = ft_struct2double(functional, 3);

if (strcmp(cfg.interpmethod, 'nearest') || strcmp(cfg.interpmethod, 'mode')) && (ft_datatype(functional, 'volume+label') || ft_datatype(functional, 'source+label') || ft_datatype(functional, 'mesh+label'))
  % the first input argument describes a parcellation or segmentation with tissue labels
  isAtlasFun = true;
else
  isAtlasFun = false;
end

if isfield(anatomical, 'transform') && isfield(anatomical, 'dim')
  % anatomical volume
  isUnstructuredAna  = false;
elseif isfield(anatomical, 'pos') && isfield(anatomical, 'dim')
  % positions that can be mapped onto a 3D regular grid
  isUnstructuredAna  = false;
elseif isfield(anatomical, 'pos')
  % anatomical data that consists of a mesh, but no smudging possible
  isUnstructuredAna  = true;
end

if isfield(functional, 'transform') && isfield(functional, 'dim')
  % functional volume
  isUnstructuredFun  = false;
elseif isfield(functional, 'pos') && isfield(functional, 'dim')
  % positions that can be mapped onto a 3D regular grid
  isUnstructuredFun  = false;
else
  isUnstructuredFun  = true;
end

if isUnstructuredAna
  anatomical = ft_checkdata(anatomical, 'datatype', {'source', 'source+label', 'mesh'}, 'insidestyle', 'logical', 'feedback', 'yes', 'hasunit', 'yes');
else
  anatomical = ft_checkdata(anatomical, 'datatype', {'volume', 'volume+label'}, 'insidestyle', 'logical', 'feedback', 'yes', 'hasunit', 'yes');
end

if isUnstructuredFun
  functional = ft_checkdata(functional, 'datatype', 'source', 'insidestyle', 'logical', 'feedback', 'yes', 'hasunit', 'yes');
else
  functional = ft_checkdata(functional, 'datatype', 'volume', 'insidestyle', 'logical', 'feedback', 'yes', 'hasunit', 'yes');
end

% select the parameters from the data, this needs to be done here, because after running checkdata, the parameterselection fails if the numeric data has nfreq/ntime/etc>1
cfg.parameter = parameterselection(cfg.parameter, functional);

% ensure that the functional data has the same unit as the anatomical data
functional = ft_convert_units(functional, anatomical.unit);

if isfield(functional, 'coordsys') && isfield(anatomical, 'coordsys') && ~isequal(functional.coordsys, anatomical.coordsys)
  % FIXME is this different when smudged or not?
  % ft_warning('the coordinate systems are not aligned');
  % ft_error('the coordinate systems are not aligned');
end

if ~isUnstructuredAna && cfg.downsample~=1
  % downsample the anatomical volume
  tmpcfg = keepfields(cfg, {'downsample', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
  tmpcfg.parameter = 'anatomy';
  anatomical = ft_volumedownsample(tmpcfg, anatomical);
  % restore the provenance information and put back cfg.parameter
  tmpparameter = cfg.parameter;
  [cfg, anatomical] = rollback_provenance(cfg, anatomical);
  cfg.parameter = tmpparameter;
end

% collect the functional volumes that should be converted
dat_name = {};
dat_array = {};
for i=1:length(cfg.parameter)
  if ~iscell(getsubfield(functional, cfg.parameter{i}))
    dat_name{end+1} = cfg.parameter{i};
    dat_array{end+1} = getsubfield(functional, cfg.parameter{i});
  else
    fprintf('not interpolating %s, since it is represented in a cell-array\n', cfg.parameter{i});
  end
end

% hmmmm, if the input data contains a time and/or freq dimension, then the output
% may be terribly blown up; most convenient would be to output only the
% smudging matrix, and project the data when plotting

if isUnstructuredFun && isUnstructuredAna && isfield(anatomical, 'orig') && isfield(anatomical.orig, 'pos') && isfield(anatomical.orig, 'tri')
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % functional data defined on subset of vertices in an anatomical mesh
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % FIXME this should not be decided on the basis of the data structures but on the basis of the cfg.interpmethod option
  % FIXME the distribution of 3 geometries over the 2 structures is weird
  % FIXME a (perhaps extreme) application of this would be to interpolate data from parcels on the sheet, i.e. an inverse parcellation

  % anatomical data consists of a decimated triangulated mesh, containing
  % the original description, allowing for smudging.

  % smudge the low resolution functional data according to the strategy in
  % MNE-suite (chapter 8.3 of the manual)

  interpmat = interp_ungridded(anatomical.pos, anatomical.orig.pos, 'projmethod', 'smudge', 'triout', anatomical.orig.tri);
  interpmat(~anatomical.inside(:), :) = 0;

  % start with an empty structure, keep only some fields
  interp = keepfields(functional, {'time', 'freq'});
  interp = copyfields(anatomical, interp, {'coordsys', 'unit'});
  interp = copyfields(anatomical.orig, interp, {'pos', 'tri'});

  % identify the inside voxels after interpolation
  nzeros     = sum(interpmat~=0,2);
  newinside  = (nzeros~=0);
  newoutside = (nzeros==0);

  interp.inside = false(size(anatomical.pos,1),1);
  interp.inside(newinside) = true;

  % interpolate all functional data
  for i=1:length(dat_name)
    fprintf('interpolating %s\n', dat_name{i});

    dimord = getdimord(functional, dat_name{i});
    dimtok = tokenize(dimord, '_');
    dimf   = getdimsiz(functional, dat_name{i});
    dimf(end+1:length(dimtok)) = 1; % there can be additional trailing singleton dimensions

    % should be 3-D array, can have trailing singleton dimensions
    if numel(dimf)<2
      dimf(2) = 1;
    end
    if numel(dimf)<3
      dimf(3) = 1;
    end

    allav = zeros([size(anatomical.orig.pos,1), dimf(2:end)]);
    for k=1:dimf(2)
      for m=1:dimf(3)
        fv     = double_ifnot(dat_array{i}(:,k,m));
        av     = interpmat*fv;
        av(newoutside) = nan;
        allav(:,k,m)   = av;
      end
    end
    interp = setsubfield(interp, dat_name{i}, allav);
  end

elseif isUnstructuredFun && isUnstructuredAna
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % functional data defined on a point cloud/mesh, anatomy on a point cloud/mesh
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % set default interpmethod for this situation
  cfg.interpmethod = ft_getopt(cfg, 'interpmethod', 'nearest');
  cfg.sphereradius = ft_getopt(cfg, 'sphereradius', 0.5);
  cfg.power        = ft_getopt(cfg, 'power',        1);

  interpmat = interp_ungridded(functional.pos, anatomical.pos, 'projmethod', cfg.interpmethod, 'sphereradius', cfg.sphereradius, 'power', cfg.power); % FIXME include other key-value pairs as well
  interpmat(~anatomical.inside(:), :) = 0;

  % start with an empty structure, keep only some fields
  interp = keepfields(functional, {'time', 'freq'});
  interp = copyfields(anatomical, interp, {'pos', 'tri', 'dim', 'transform', 'coordsys', 'unit'});

  % identify the inside voxels after interpolation
  nzeros     = sum(interpmat~=0,2);
  newinside  = (nzeros~=0);
  newoutside = (nzeros==0);

  interp.inside = false(size(anatomical.pos,1),1);
  interp.inside(newinside) = true;

  % interpolate all functional data
  for i=1:length(dat_name)
    fprintf('interpolating %s\n', dat_name{i});

    dimord = getdimord(functional, dat_name{i});
    dimtok = tokenize(dimord, '_');
    dimf   = getdimsiz(functional, dat_name{i});
    dimf(end+1:length(dimtok)) = 1; % there can be additional trailing singleton dimensions

    % should be 3-D array, can have trailing singleton dimensions
    if numel(dimf)<2
      dimf(2) = 1;
    end
    if numel(dimf)<3
      dimf(3) = 1;
    end

    allav = zeros([size(anatomical.pos,1), dimf(2:end)]);
    for k=1:dimf(2)
      for m=1:dimf(3)
        fv     = double_ifnot(dat_array{i}(:,k,m));
        av     = interpmat*fv;
        av(newoutside) = nan;
        allav(:,k,m)   = av;
      end
    end
    interp = setsubfield(interp, dat_name{i}, allav);
  end

elseif isUnstructuredFun && ~isUnstructuredAna
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % functional data defined on a point cloud/mesh, anatomy on a volume
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % set default interpmethod for this situation
  cfg.interpmethod = ft_getopt(cfg, 'interpmethod', 'nearest');
  cfg.sphereradius = ft_getopt(cfg, 'sphereradius', 0.5);
  cfg.power        = ft_getopt(cfg, 'power',        1);

  [ax, ay, az] = voxelcoords(anatomical.dim, anatomical.transform);
  anatomical.pos = [ax(:) ay(:) az(:)];
  clear ax ay az

  tmp = interp_ungridded(functional.pos, anatomical.pos(anatomical.inside,:), 'projmethod', cfg.interpmethod, 'sphereradius', cfg.sphereradius, 'power', cfg.power); % FIXME include other key-value pairs as well
  interpmat( anatomical.inside(:), :) = tmp;
  interpmat(~anatomical.inside(:), :) = 0;

  % start with an empty structure, keep only some fields
  interp = keepfields(functional, {'time', 'freq'});
  interp = copyfields(anatomical, interp, {'pos', 'tri', 'dim', 'transform', 'coordsys', 'unit', 'anatomy'});

  % identify the inside voxels after interpolation
  nzeros     = sum(interpmat~=0,2);
  newinside  = (nzeros~=0);
  newoutside = (nzeros==0);

  interp.inside = false(anatomical.dim);
  interp.inside(newinside) = true;

  % interpolate all functional data
  for i=1:length(dat_name)
    fprintf('interpolating %s\n', dat_name{i});

    dimord = getdimord(functional, dat_name{i});
    dimtok = tokenize(dimord, '_');
    dimf   = getdimsiz(functional, dat_name{i});
    dimf(end+1:length(dimtok)) = 1; % there can be additional trailing singleton dimensions

    % should be 3-D array, can have trailing singleton dimensions
    if numel(dimf)<2
      dimf(2) = 1;
    end
    if numel(dimf)<3
      dimf(3) = 1;
    end

    av    = zeros([anatomical.dim            ]);
    allav = zeros([anatomical.dim dimf(2:end)]);

    for k=1:dimf(2)
      for m=1:dimf(3)
        fv     = double_ifnot(dat_array{i}(:,k,m));
        av(:)  = interpmat*fv;
        av(newoutside)   = nan;
        allav(:,:,:,k,m) = av;
      end
    end
    if isfield(interp, 'freq') || isfield(interp, 'time')
      % the output should be a source representation, not a volume
      allav = reshape(allav, prod(anatomical.dim), dimf(2), dimf(3));
    end
    interp = setsubfield(interp, dat_name{i}, allav);
  end

  if ~isfield(interp, 'freq') && ~isfield(interp, 'time')
      % the output should be a volume representation, not a source
    interp = rmfield(interp, 'pos');
  end

elseif ~isUnstructuredFun && isUnstructuredAna
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % functional data defined on a volume, anatomy on a point cloud/mesh
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % set default interpmethod for this situation
  cfg.interpmethod = ft_getopt(cfg, 'interpmethod', 'nearest');
  cfg.sphereradius = ft_getopt(cfg, 'sphereradius', []);
  cfg.power        = ft_getopt(cfg, 'power',        1);

  % interpolate the 3D volume onto the anatomy
  if ~strcmp(cfg.interpmethod, 'project')
    % use interp_gridded
    [interpmat, dummy] = interp_gridded(functional.transform, zeros(functional.dim), anatomical.pos, 'projmethod', cfg.interpmethod, 'sphereradius', cfg.sphereradius, 'inside', functional.inside, 'power', cfg.power);

    % use interp_ungridded
    % interpmat = interp_ungridded(functional.pos, anatomical.pos, 'projmethod', cfg.interpmethod, 'sphereradius', cfg.sphereradius, 'inside', functional.inside, 'power', cfg.power);
  else
    % do the interpolation below, the current implementation of the
    % 'project' method does not output an interpmat (and is therefore quite
    % inefficient

    % set the defaults
    cfg.projvec        = ft_getopt(cfg, 'projvec',       1);
    cfg.projweight     = ft_getopt(cfg, 'projweight',    ones(size(cfg.projvec)));
    cfg.projcomb       = ft_getopt(cfg, 'projcomb',      'mean'); % or max
    cfg.projthresh     = ft_getopt(cfg, 'projthresh',    []);
  end

  % start with an empty structure, keep some fields
  interp = keepfields(functional, {'time', 'freq'});
  interp = copyfields(anatomical, interp, {'pos', 'tri', 'dim', 'transform', 'coordsys', 'unit'});

  % identify the inside voxels after interpolation
  interp.inside    = true(size(anatomical.pos,1),1);

  % interpolate all functional data
  for i=1:length(dat_name)
    fprintf('interpolating %s\n', dat_name{i});

    dimord = getdimord(functional, dat_name{i});
    dimtok = tokenize(dimord, '_');
    dimf   = getdimsiz(functional, dat_name{i});
    dimf(end+1:length(dimtok)) = 1; % there can be additional trailing singleton dimensions

    if prod(functional.dim)==dimf(1)
      % convert into 3-D, 4-D or 5-D array
      dimf        = [functional.dim dimf(2:end)];
      dat_array{i} = reshape(dat_array{i}, dimf);
    end

    % should be 5-D array, can have trailing singleton dimensions
    if numel(dimf)<4
      dimf(4) = 1;
    end
    if numel(dimf)<5
      dimf(5) = 1;
    end

    allav = zeros([size(anatomical.pos,1), dimf(4:end)]);
    if ~strcmp(cfg.interpmethod, 'project')
      for k=1:dimf(4)
        for m=1:dimf(5)
          fv    = double_ifnot(dat_array{i}(:,:,:,k,m)); % ensure double precision to allow sparse multiplication
          fv    = fv(functional.inside(:));
          av    = interpmat*fv;
          allav(:,k,m) = av;
        end
      end
    else
      for k=1:dimf(4)
        for m=1:dimf(5)
          fv   = double_ifnot(dat_array{i}(:,:,:,k,m));
          av   = interp_gridded(functional.transform, fv, anatomical.pos, 'dim', functional.dim, 'projmethod', 'project', 'projvec', cfg.projvec, 'projweight', cfg.projweight, 'projcomb', cfg.projcomb, 'projthresh', cfg.projthresh);
          allav(:,k,m) = av;
        end
      end
    end
    interp = setsubfield(interp, dat_name{i}, allav);
  end

elseif ~isUnstructuredFun && ~isUnstructuredAna
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % functional data defined on a volume, anatomy on a differently sampled volume
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % set default interpmethod for this situation
  cfg.interpmethod = ft_getopt(cfg, 'interpmethod', 'linear');
  if isequal(cfg.interpmethod, 'mode') && isAtlasFun
    % use a mode-based interpolation, i.e. a majority vote of the nearby
    % voxel locations.

    % first to a nearest interpolation of the voxel coordinates the other
    % way around
    tmp = anatomical;
    [tx, ty, tz] = voxelcoords(tmp.dim, tmp.transform);
    tmp.tx = tx;
    tmp.ty = ty;
    tmp.tz = tz;

    tmpcfg = [];
    tmpcfg.interpmethod = 'nearest';
    tmpcfg.parameter = {'tx' 'ty' 'tz'};
    tmpint = ft_sourceinterpolate(tmpcfg, tmp, functional);

    [ix,i1,i2] = intersect([tx(:) ty(:) tz(:)],[tmpint.tx(:) tmpint.ty(:) tmpint.tz(:)],'rows');

    interp = keepfields(anatomical, {'pos', 'dim', 'transform', 'coordsys', 'unit',  'anatomy'});
    interp.inside = false(interp.dim);
    interp.inside(i1) = true;

    for k = 1:numel(cfg.parameter)
      interp.(cfg.parameter{k}) = nan(interp.dim);
      fun = functional.(cfg.parameter{k});
      for m = 1:numel(i1)
        values = reshape(fun(tmpint.tx==tx(i1(m)) & tmpint.ty==ty(i1(m)) & tmpint.tz==tz(i1(m))),[],1);
        [M,f,c] = mode(values);
        if numel(c)==1
          interp.(cfg.parameter{k})(i1(m)) = M;
        else
          ft_warning('multiple modes per voxel, returning NaN');
          interp.(cfg.parameter{k})(i1(m)) = nan;
        end
      end
    end

  elseif isequal(cfg.interpmethod, 'mode') && ~isAtlasFun
    ft_error('the interpolation method ''mode'' is only supported for parcellations');

  else
    % start with an empty structure, keep some fields
    interp = keepfields(functional, {'time', 'freq'});
    interp = copyfields(anatomical, interp, {'pos', 'dim', 'transform', 'coordsys', 'unit', 'anatomy'});

    % convert the anatomical voxel positions into voxel indices into the functional volume
    anatomical.transform = functional.transform \ anatomical.transform;
    functional.transform = eye(4);

    [fx, fy, fz] = voxelcoords(functional.dim, functional.transform);
    [ax, ay, az] = voxelcoords(anatomical.dim, anatomical.transform);

    % estimate the subvolume of the anatomy that is spanned by the functional volume
    minfx = 1;
    minfy = 1;
    minfz = 1;
    maxfx = functional.dim(1);
    maxfy = functional.dim(2);
    maxfz = functional.dim(3);
    sel = ax(:)>=minfx & ...
      ax(:)<=maxfx & ...
      ay(:)>=minfy & ...
      ay(:)<=maxfy & ...
      az(:)>=minfz & ...
      az(:)<=maxfz;
    fprintf('selecting subvolume of %.1f%%\n', 100*sum(sel)./prod(anatomical.dim));

    if all(functional.inside(:))
      % keep all voxels marked as inside
      interp.inside = true(anatomical.dim);
    else
      % reslice and interpolate inside
      interp.inside = zeros(anatomical.dim);
      % interpolate with method nearest
      interp.inside( sel) = my_interpn(double(functional.inside), ax(sel), ay(sel), az(sel), 'nearest', cfg.feedback);
      interp.inside(~sel) = 0;
      interp.inside = logical(interp.inside);
    end

    % prepare the grid that is used in the interpolation
    fg = [fx(:) fy(:) fz(:)];
    clear fx fy fz

    % reslice and interpolate all functional volumes
    for i=1:length(dat_name)
      fprintf('reslicing and interpolating %s\n', dat_name{i});

      dimord = getdimord(functional, dat_name{i});
      dimtok = tokenize(dimord, '_');
      dimf   = getdimsiz(functional, dat_name{i});
      dimf(end+1:length(dimtok)) = 1; % there can be additional trailing singleton dimensions

      if prod(functional.dim)==dimf(1)
        % convert into 3-D, 4-D or 5-D array
        dimf = [functional.dim dimf(2:end)];
        dat_array{i} = reshape(dat_array{i}, dimf);
      end

      % should be 5-D array, can have trailing singleton dimensions
      if numel(dimf)<4
        dimf(4) = 1;
      end
      if numel(dimf)<5
        dimf(5) = 1;
      end

      av    = zeros([anatomical.dim            ]);
      allav = zeros([anatomical.dim dimf(4:end)]);
      functional.inside = functional.inside(:,:,:,1,1);

      if any(dimf(4:end)>1) && ~strcmp(cfg.feedback, 'none')
        % this is needed to prevent feedback to be displayed for every time-frequency point
        ft_warning('disabling feedback');
        cfg.feedback = 'none';
      end

      for k=1:dimf(4)
        for m=1:dimf(5)
          fv = double_ifnot(dat_array{i}(:,:,:,k,m));
          % av( sel) = my_interpn(fx, fy, fz, fv, ax(sel), ay(sel), az(sel), cfg.interpmethod, cfg.feedback);
          if islogical(dat_array{i})
            % interpolate always with method nearest
            av( sel) = my_interpn(fv, ax(sel), ay(sel), az(sel), 'nearest', cfg.feedback);
            av = logical(av);
          else
            if ~all(functional.inside(:))
              % extrapolate the outside of the functional volumes for better interpolation at the edges
              fv(~functional.inside) = griddatan(fg(functional.inside(:), :), fv(functional.inside(:)), fg(~functional.inside(:), :), 'nearest');
            end
            % interpolate functional onto anatomical grid
            av( sel) = my_interpn(fv, ax(sel), ay(sel), az(sel), cfg.interpmethod, cfg.feedback);
            av(~sel) = nan;
            av(~interp.inside) = nan;
          end
          allav(:,:,:,k,m) = av;
        end
      end
      if isfield(interp, 'freq') || isfield(interp, 'time')
        % the output should be a source representation, not a volume
        allav = reshape(allav, prod(anatomical.dim), dimf(4), dimf(5));
      end
      interp = setsubfield(interp, dat_name{i}, allav);
      % keep the description of the labels in the segmentation/parcellation
      if strcmp(cfg.interpmethod, 'nearest') && isfield(functional, [dat_name{i} 'label'])
        interp.([dat_name{i} 'label']) = functional.([dat_name{i} 'label']);
      end
    end

  end

end % computing the interpolation according to the input data

if isfield(interp, 'freq') || isfield(interp, 'time')
  % the output should be a source representation, not a volumetric representation
  if ~isfield(interp, 'pos')
    [x, y, z] = voxelcoords(interp.dim, interp.transform);
    interp.pos = [x(:) y(:) z(:)];
  end
end

if isAtlasFun
  for i=1:numel(dat_name)
    % keep the labels that describe the different tissue types
    interp = copyfields(functional, interp, [dat_name{i} 'label']);
    % replace NaNs that fall outside the labeled area with zero
    tmp = interp.(dat_name{i});
    tmp(isnan(tmp)) = 0;
    interp.(dat_name{i}) = tmp;
  end
  % remove the inside field if present
  interp = removefields(interp, 'inside');
end

if exist('interpmat', 'var')
  cfg.interpmat = interpmat;
  cfg.interpmat; % access it once to fool the cfg-tracking
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   functional anatomical
ft_postamble provenance interp
ft_postamble history    interp
ft_postamble savevar    interp

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION this function computes the location of all voxels in head
% coordinates in a memory efficient manner
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x, y, z] = voxelcoords(dim, transform)

xgrid     = 1:dim(1);
ygrid     = 1:dim(2);
zgrid     = 1:dim(3);
npix      = prod(dim(1:2));  % number of voxels in a single slice

x = zeros(dim);
y = zeros(dim);
z = zeros(dim);
X = zeros(1,npix);
Y = zeros(1,npix);
Z = zeros(1,npix);
E = ones(1,npix);
% determine the voxel locations per slice
for i=1:dim(3)
  [X(:), Y(:), Z(:)] = ndgrid(xgrid, ygrid, zgrid(i));
  tmp = transform*[X; Y; Z; E];
  x((1:npix)+(i-1)*npix) = tmp(1,:);
  y((1:npix)+(i-1)*npix) = tmp(2,:);
  z((1:npix)+(i-1)*npix) = tmp(3,:);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION for memory efficient interpolation
% the only reason for this function is that it does the interpolation in smaller chuncks
% this prevents memory problems that I often encountered here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [av] = my_interpn(fx, fy, fz, fv, ax, ay, az, interpmethod, feedback);
function [av] = my_interpn(fv, ax, ay, az, interpmethod, feedback)

num = numel(ax);            % total number of voxels
blocksize = floor(num/20);  % number of voxels to interpolate at once, split it into 20 chuncks
lastblock = 0;              % boolean flag for while loop
sel = 1:blocksize;          % selection of voxels that are interpolated, this is the first chunck
av  = zeros(size(ax));
ft_progress('init', feedback, 'interpolating');
while (1)
  ft_progress(sel(1)/num, 'interpolating %.1f%%\n', 100*sel(1)/num);
  if sel(end)>=num
    sel = sel(1):num;
    lastblock = 1;
  end
  av(sel) = interpn(fv, ax(sel), ay(sel), az(sel), interpmethod);
  if lastblock
    break
  end
  sel = sel + blocksize;
end
ft_progress('close');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION to cast array to double precision, only if needed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = double_ifnot(in)

if ~isa(in, 'double')
  out = double(in);
else
  out = in;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ result, node_num ] = sphere01_triangle_quad_icos1v ( a_xyz, ...
  b_xyz, c_xyz, factor, fun )

%*****************************************************************************80
%
%% SPHERE01_TRIANGLE_QUAD_ICOS1V: vertex rule, subdivide then project.
%
%  Discussion:
%
%    This function estimates an integral over a spherical triangle on the
%    unit sphere.
%
%    This function sets up an icosahedral grid, and subdivides each
%    edge of the icosahedron into FACTOR subedges.  These edges define a grid
%    within each triangular icosahedral face.   All of these calculations are
%    done, essentially, on the FLAT faces of the icosahedron.  Only then are
%    the triangle vertices projected to the sphere.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 April 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A_XYZ(3), B_XYZ(3), C_XYZ(3), the vertices
%    of the spherical triangle.
%
%    Input, integer FACTOR, the subdivision factor, which must
%    be at least 1.
%
%    Input, function v = fun ( x ), evaluates the integrand at the point X.
%
%    Output, real RESULT, the estimated integral.
%
%    Output, integer NODE_NUM, the number of evaluation points.
%

%
%  Destroy all row vectors.
%
  a_xyz = a_xyz(:);
  b_xyz = b_xyz(:);
  c_xyz = c_xyz(:);
%
%  Initialize the integral data.
%
  result = 0.0;
  area_total = 0.0;
  node_num = 0;
%
%  Deal with subtriangles that have same orientation as face.
%
  for f1 = 0 : factor - 1
    for f2 = 0 : factor - f1 - 1
      f3 = factor - f1 - f2;

      a2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1 + 1, f2,     f3 - 1 );
      b2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1,     f2 + 1, f3 - 1 );
      c2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1,     f2,     f3 );

      area = sphere01_triangle_vertices_to_area ( a2_xyz, b2_xyz, c2_xyz );

      node_num = node_num + 3;
      va = fun ( a2_xyz );
      vb = fun ( b2_xyz );
      vc = fun ( c2_xyz );
      result = result + area * ( va + vb + vc ) / 3.0;
      area_total = area_total + area;

      end
    end
%
%  Deal with subtriangles that have opposite orientation as face.
%
  for f3 = 0 : factor - 2
    for f2 = 1 : factor - f3 - 1
      f1 = factor - f2 - f3;

      a2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1 - 1, f2,     f3 + 1 );
      b2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1,     f2 - 1, f3 + 1 );
      c2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1,     f2,     f3 );

      area = sphere01_triangle_vertices_to_area ( a2_xyz, b2_xyz, c2_xyz );

      node_num = node_num + 3;
      va = fun ( a2_xyz );
      vb = fun ( b2_xyz );
      vc = fun ( c2_xyz );
      result = result + area * ( va + vb + vc ) / 3.0;
      area_total = area_total + area;

    end
  end

  return
end