text:
% POP_STUDY  - create a new STUDY set structure defining a group of related EEG datasets.
%                The STUDY set also contains information about each of the datasets: the
%                subject code, subject group, experimental condition, and session. This can
%                be provided interactively in a pop-up window or be automatically filled
%                in by the function. Defaults: Assume a different subject for each
%                dataset and only one condition; leave subject group and session fields
%                empty. Additional STUDY information about the STUDY name, task and
%                miscellaneous notes can also be saved in the STUDY structure.
% Usage:
%  >> [ STUDY ALLEEG ] = pop_study([],[], 'gui', 'on'); % create new study interactively
%  >> [ STUDY ALLEEG ] = pop_study(STUDY, ALLEEG, 'gui', 'on'); % edit study interactively
%  >> [ STUDY ALLEEG ] = pop_study(STUDY, ALLEEG, 'key', 'val', ...); % edit study
%
% Optional Inputs:
%   STUDY                - existing study structure.
%   ALLEEG               - vector of EEG dataset structures to be included in the STUDY.
%
% Optional Inputs:
%   All "'key', 'val'" inputs of STD_EDITSET may be used.
%
% Outputs:
%   STUDY                - new STUDY set comprising some or all of the datasets in
%                          ALLEEG, plus other information about the experiments.
%   ALLEEG               - an updated ALLEEG structure including the STUDY datasets.
%
% Graphic interface buttons:
%  "STUDY set name"      - [edit box] name for the STUDY structure {default: ''}
%  "STUDY set task name" - [edit box] name for the task performed by the subject {default: ''}
%  "STUDY set notes"     - [edit box] notes about the experiment, the datasets, the STUDY,
%                          or anything else to store with the rest of the STUDY information
%                          {default: ''}
%  "dataset filename/browse" - Load dataset from specified filename. It may be
%                          an EEGLAB dataset or any file supported by File-IO (if
%                          it is installed from the plugin manager).
%  "subject"             - [edit box] subject code associated with the dataset. If no
%                          subject code is provided, each dataset will assumed to be from
%                          a different subject {default: 'S1', 'S2', ..., 'Sn'}
%  "session"             - [edit box] dataset session. If no session information is
%                          provided, all datasets that belong to one subject are assumed to
%                          have been recorded within one session {default: []}
%  "run"                 - [edit box] dataset run. If no run information is
%                          provided, all datasets that belong to one subject are assumed to
%                          have been recorded within one run {default: []}
%  "condition"           - [edit box] dataset condition. If no condition code is provided,
%                          all datasets are assumed to be from the same condition {default:[]}
%  "group"               - [edit box] the subject group the dataset belongs to. If no group
%                          is provided, all subjects and datasets are assumed to belong to
%                          the same group. {default: []}
%  "Save this STUDY set to disk file" - [check box] If checked, save the new STUDY set
%                          structure to disk. If no filename is provided, a window will
%                          pop up to ask for it.
%
% See also: std_editset, POP_LOADSTUDY, POP_PRECLUST, POP_CLUST
%
% Authors: Arnaud Delorme, Hilit Serby, Scott Makeig, SCCN, INC, UCSD, July 22, 2005

% Copyright (C) Arnaud Delorme, SCCN, INC, UCSD, July 22, 2005, arno@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% Coding notes: Useful information on functions and global variables used.

function [STUDY, ALLEEG, com]  = pop_study(STUDY, ALLEEG, varargin)

com = '';

if nargin < 1
    help pop_study;
    return;
end

% type of call (gui, script or internal)
% --------------------------------------
mode = 'internal_command';
if ~ischar(STUDY) %initial settings
    mode = 'script';
    if nargin > 2
        for index = 1:length(varargin)
            if ischar(varargin{index})
                if strcmpi(varargin{index}, 'gui')
                    mode = 'gui';
                end
            end
        end
    end
end

if isempty(STUDY)
    newstudy       = 1;
    STUDY.name     = '';
    STUDY.task     = '';
    STUDY.notes    = '';
    STUDY.filename = '';
    STUDY.cluster  = [];
    STUDY.history = 'STUDY = [];';
else
    newstudy       = 0;
end

if strcmpi(mode, 'script') % script mode
    [STUDY, ALLEEG] = std_editset(STUDY, ALLEEG, varargin{:});
    return;
elseif strcmpi(mode, 'gui') % GUI mode

    % check events
    fieldList = { 'condition' 'group' 'subject' 'session' 'run' };
    if ~isempty(ALLEEG) && isstruct(ALLEEG) && isfield(ALLEEG, 'event')
        for iField = 1:length(fieldList)
            if any(cellfun(@(x)isfield(x, fieldList{iField}), {ALLEEG.event}))
                fprintf(2, 'Rename field "%s" in datasets'' event structure as it conflicts with the STUDY field bearing the same name\n', fieldList{iField});
            end
        end
    end

    % show warning if necessary
    % -------------------------
    if isreal(ALLEEG)
        if ALLEEG == 0
            res = questdlg2( strvcat('Datasets currently loaded will be removed from EEGLAB memory.', ...
                                     'Are you sure you want to continue?'), ...
                                     'Discard loaded EEGLAB datasets?', 'Cancel', 'Yes', 'Yes');
            if strcmpi(res, 'cancel'), return; end
        end
        ALLEEG = [];
        alleegEmpty = 1;
    else
        alleegEmpty = 0;
    end

    % set initial datasetinfo
    % -----------------------
    if isfield(STUDY, 'datasetinfo')
        datasetinfo = STUDY.datasetinfo;
        different = 0;
        for k = 1:length(ALLEEG)
            if ~strcmpi(datasetinfo(k).filename, ALLEEG(k).filename), different = 1; break; end
            if ~strcmpi(datasetinfo(k).subject,   ALLEEG(k).subject),   different = 1; break; end
            if ~strcmpi(datasetinfo(k).condition, ALLEEG(k).condition), different = 1; break; end
            if ~strcmpi(char(datasetinfo(k).group), char(ALLEEG(k).group)),     different = 1; break; end
            if ~isequal(datasetinfo(k).session, ALLEEG(k).session),             different = 1; break; end
            if ~isequal(datasetinfo(k).run, ALLEEG(k).run),                     different = 1; break; end
        end
        if different
            info = 'from_STUDY_different_from_ALLEEG';
        else
            info = 'from_STUDY';
        end
        if ~isfield(datasetinfo, 'comps');
            datasetinfo(1).comps = [];
        end
    else
        info = 'from_ALLEEG';
        if length(ALLEEG) > 0
            datasetinfo(length(ALLEEG)).filename  = [];
            datasetinfo(length(ALLEEG)).filepath  = [];
            datasetinfo(length(ALLEEG)).subject   = [];
            datasetinfo(length(ALLEEG)).session   = [];
            datasetinfo(length(ALLEEG)).run       = [];
            datasetinfo(length(ALLEEG)).condition = [];
            datasetinfo(length(ALLEEG)).group     = [];
            for k = 1:length(ALLEEG)
                datasetinfo(k).filename  = ALLEEG(k).filename;
                datasetinfo(k).filepath  = ALLEEG(k).filepath;
                datasetinfo(k).subject   = ALLEEG(k).subject;
                datasetinfo(k).session   = ALLEEG(k).session;
                datasetinfo(k).run       = ALLEEG(k).run;
                datasetinfo(k).condition = ALLEEG(k).condition;
                datasetinfo(k).group     = ALLEEG(k).group;
            end
            if ~isfield(datasetinfo, 'comps');
                datasetinfo(1).comps = [];
            end
        else
            datasetinfo = [];
        end
    end

    nextpage    = 'pop_study(''nextpage'', gcbf);';
    prevpage    = 'pop_study(''prevpage'', gcbf);';
    delset      = 'pop_study(''clear'',    gcbf, get(gcbo, ''userdata''));';
    loadset     = 'pop_study(''load'',     gcbf, get(guiind, ''userdata''), get(guiind, ''string'')); clear guiind;';
    loadsetedit = [ 'guiind = gcbo;' loadset ];
    subcom      = 'pop_study(''subject''  , gcbf, get(gcbo, ''userdata''), get(gcbo, ''string''));';
    sescom      = 'pop_study(''session''  , gcbf, get(gcbo, ''userdata''), get(gcbo, ''string''));';
    runcom      = 'pop_study(''run''      , gcbf, get(gcbo, ''userdata''), get(gcbo, ''string''));';
    condcom     = 'pop_study(''condition'', gcbf, get(gcbo, ''userdata''), get(gcbo, ''string''));';
    grpcom      = 'pop_study(''group''    , gcbf, get(gcbo, ''userdata''), get(gcbo, ''string''));';
    compcom     = 'pop_study(''component'', gcbf, get(gcbo, ''userdata''), get(gcbo, ''string''));';
    cb_del      = 'pop_study(''delclust'' , gcbf, ''showwarning'');';
    cb_dipole   = 'pop_study(''dipselect'', gcbf, ''showwarning'');';

	browsestudy = [ '[filename, filepath] = uiputfile2(''*.study'', ''Use existing STUDY set to import dataset information -- pop_study()''); ' ...
                      'set(findobj(''parent'', gcbf, ''tag'', ''usestudy_file''), ''string'', [filepath filename]);' ];
	saveSTUDY = [ 'set(findobj(''parent'', gcbf, ''userdata'', ''save''), ''enable'', fastif(get(gcbo, ''value'')==1, ''on'', ''off''));' ];
	browsesave = [ '[filename, filepath] = uiputfile2(''*.study'', ''Save STUDY with .study extension -- pop_clust()''); ' ...
                   'if filename ~= 0,' ...
                   '   set(findobj(''parent'', gcbf, ''tag'', ''studyfile''), ''string'', [filepath filename]);' ...
                   'end;' ...
                   'clear filename filepath;' ];

    texthead = fastif(newstudy, 'Create a new STUDY set', 'Edit STUDY set information - remember to save changes');
	guispec = { ...
        {'style' 'text' 'string' texthead 'FontWeight' 'Bold' 'HorizontalAlignment' 'center'} ...
        {} {'style' 'text' 'string' 'STUDY set name:' } { 'style' 'edit' 'string' char(STUDY.name) 'tag' 'study_name' } ...
        {} {'style' 'text' 'string' 'STUDY set task name:' } { 'style' 'edit' 'string' char(STUDY.task) 'tag' 'study_task' } ...
        {} {'style' 'text' 'string' 'STUDY set notes:' } { 'style' 'edit' 'string' char(STUDY.notes) 'tag' 'study_notes' } {}...
        {} ...
        {'style' 'text' 'string' 'dataset filename' 'userdata' 'addt'} {'style' 'text' 'string' 'browse' 'userdata' 'addt'} ...
        {'style' 'text' 'string' 'subject'    'userdata' 'addt'} ...
        {'style' 'text' 'string' 'session'    'userdata' 'addt'} ...
        {'style' 'text' 'string' 'run'        'userdata' 'addt'} ...
        {'style' 'text' 'string' 'condition'  'userdata' 'addt'} ...
        {'style' 'text' 'string' 'group'      'userdata' 'addt'} ...
        {'style' 'pushbutton' 'string' 'Select by r.v.' 'userdata' 'addt' 'callback' cb_dipole } ...
        {} };
    geomUnit  = [0.2 1 3.5];
    geomUnit2 = [0.2 1.05 0.35 0.4 0.35 0.25 0.6 0.4 0.6 0.3];
	guigeom = { [1] geomUnit geomUnit geomUnit [1] geomUnit2 };

    % create edit boxes
    % -----------------
    for index = 1:10
        guigeom = { guigeom{:} geomUnit2 };
        select_com = ['[inputname, inputpath] = uigetfile2(''*.set;*.SET'', ''Choose dataset to add to STUDY -- pop_study()'');'...
                      'if inputname ~= 0,' ...
                      '   guiind = findobj(''parent'', gcbf, ''tag'', ''set' int2str(index) ''');' ...
                      '   set( guiind,''string'', [inputpath inputname]);' ...
                          loadset ...
                      'end; clear inputname inputpath;'];
        numstr = int2str(index);
        guispec = { guispec{:}, ...
        {'style' 'text'       'string' numstr  'tag' [ 'num'   int2str(index) ] 'userdata' index }, ...
        {'style' 'edit'       'string' ''      'tag' [ 'set'   int2str(index) ] 'userdata' index 'callback' loadsetedit}, ...
        {'style' 'pushbutton' 'string' '...'   'tag' [ 'brw'   int2str(index) ] 'userdata' index 'Callback' select_com}, ...
        {'style' 'edit'       'string' ''      'tag' [ 'sub'   int2str(index) ] 'userdata' index 'Callback' subcom}, ...
        {'style' 'edit'       'string' ''      'tag' [ 'sess'  int2str(index) ] 'userdata' index 'Callback' sescom}, ...
        {'style' 'edit'       'string' ''      'tag' [ 'run'   int2str(index) ] 'userdata' index 'Callback' runcom}, ...
        {'style' 'edit'       'string' ''      'tag' [ 'cond'  int2str(index) ] 'userdata' index 'Callback' condcom}, ...
        {'style' 'edit'       'string' ''      'tag' [ 'group' int2str(index) ] 'userdata' index 'Callback' grpcom}, ...
        {'style' 'pushbutton' 'string' 'All comp.'   'tag' [ 'comps' int2str(index) ] 'userdata' index 'Callback' compcom}, ...
        {'style' 'pushbutton' 'string' 'CLear' 'tag' [ 'clear' int2str(index) ] 'userdata' index 'callback' delset} };
    end

    guispec = { guispec{:}, ...
                {'style' 'text'       'string'  'Important note: Removed datasets will not be saved before being deleted from EEGLAB memory' }, ...
                {}, ...
                {'style' 'pushbutton' 'string'  '<'      'Callback' prevpage 'userdata' 'addt'}, ...
                {'style' 'text'       'string'  'Page 1' 'tag' 'page' 'horizontalalignment' 'center' }, ...
                {'style' 'pushbutton' 'string'  '>'      'Callback' nextpage 'userdata' 'addt'}, {}, ...
                {}, ...
                {'style' 'checkbox'   'value'   0        'tag' 'delclust' 'callback' cb_del }, ...
                {'style' 'text'       'string'  'Delete cluster information (to allow loading new datasets, set new components for clustering, etc.)' } };
	guigeom = { guigeom{:} [1] [1 0.2 0.3 0.2 1] [1] [0.14 3] };

    fig_arg{1} = ALLEEG;      % datasets
    fig_arg{2} = datasetinfo; % datasetinfo
    fig_arg{3} = 1;           % page
    fig_arg{4} = {};          % all commands
    fig_arg{5} = (length(STUDY.cluster) > 1); % are cluster present
    fig_arg{6} = STUDY;       % are cluster present
    fig_arg{7} = [];          % tags

    % generate GUI
    % ------------
    optiongui = { 'geometry', guigeom, ...
                  'uilist'  , guispec, ...
                  'helpcom' , 'pophelp(''pop_study'')', ...
                  'title'   , 'Create a new STUDY set -- pop_study()', ...
                  'userdata', fig_arg };
	[~, ~, ~, ~, instruct, alltags] = inputgui( 'mode', 'plot', optiongui{:}); % plot only to get tags

    % save tags
    currentfig = gcf;
    fig_arg{7} = alltags;
    set(currentfig, 'userdata', fig_arg);

    % update GUI and wait for user input
    pop_study('delclust', currentfig);
    pop_study('redraw', currentfig);
    [result, userdat2, ~, outstruct] = inputgui( 'mode', currentfig, optiongui{:});
    if isempty(result), return; end

    while test_wrong_parameters(currentfig)
    	[result, userdat2, ~, outstruct] = inputgui( 'mode', currentfig, optiongui{:});
        if isempty(result), return; end
    end
    close(currentfig);

    % convert GUI selection to options
    % --------------------------------
    allcom = simplifycom(userdat2{4});
    options = {};
    if ~strcmpi(result{1}, STUDY.name ), options = { options{:} 'name'        result{1} }; end
    if ~strcmpi(result{2}, STUDY.task ), options = { options{:} 'task'        result{2} }; end
    if ~strcmpi(result{3}, STUDY.notes), options = { options{:} 'notes'       result{3} }; end
    if ~isempty(allcom),                 options = { options{:} 'commands'    allcom    }; end
%     if isnumeric(outstruct(1).studyfile)
%         if outstruct(1).studyfile == 1,  options = { options{:} 'resave'      'on' }; end
%     else
%         if ~isempty(outstruct(1).studyfile), options = { options{:} 'filename' outstruct(1).studyfile }; end
%     end
    options = { options{:} 'updatedat' 'on' };

    if outstruct(1).delclust == 1
        options = { options{:} 'rmclust' 'on' };
    else
        options = { options{:} 'rmclust' 'off' };
    end

    % ---
    if ~isequal(outstruct, instruct) && (outstruct(1).delclust ~= 1) % notice that isequal is sensitive to fields order. isequaln isn't backward compatible
        strfields = fieldnames(outstruct);
        for i = 1:length(strfields)
            strdiff(i) = strcmp(outstruct.(strfields{i}),instruct.(strfields{i}));
        end
        % If the information of the STUDY differ, then remove information from clusters
        % Ignoring ('STUDY set name','STUDY set task','STUDY set notes')
        if any(~strdiff(3:end-2))
            options{find(strcmp(options,'rmclust'))+1} = 'on';
        end
    end
    % ---

    % check channel labels
    % --------------------
    TMPALLEEG   = userdat2{1};
    if isfield(TMPALLEEG, 'chanlocs')
        allchans = { TMPALLEEG.chanlocs };
        if any(cellfun('isempty', allchans))
            txt = strvcat('Some datasets do not have channel labels. Do you wish to generate', ...
                          'channel labels automatically for all datasets ("1" for channel 1,', ...
                          '"2" for channel 2, ...). Datasets will be overwritten on disk.', ...
                          'If you abort, the STUDY will not be created.');
            res = questdlg2(txt, 'Dataset format problem', 'Yes', 'No, abort', 'Yes');
            if strcmpi(res, 'yes'), options = { options{:} 'addchannellabels' 'on' 'savedat' 'on'};
            else return;
            end
        end
    end

    % run command and create history
    % ------------------------------
    if alleegEmpty
        com = sprintf( '[STUDY ALLEEG] = std_editset( STUDY, [], %s );\n[STUDY ALLEEG] = std_checkset(STUDY, ALLEEG);', vararg2str(options) );
    else
        com = sprintf( '[STUDY ALLEEG] = std_editset( STUDY, ALLEEG, %s );\n[STUDY ALLEEG] = std_checkset(STUDY, ALLEEG);', vararg2str(options) );
    end
    [STUDY, ALLEEG] = std_editset(STUDY, ALLEEG, options{:});

    % check sessions
    % --------------
    allSessions = [ STUDY.datasetinfo.session ];
    if ~isempty(allSessions)
        if ~ismember(1, allSessions)
            fprintf(2, 'Sessions are usually numbered starting at 1; Do not use session to store custom information\n');
        end
    end

else % internal command

    com = STUDY;
    hdl = ALLEEG; %figure handle

    % userdata info
    % -------------
    userdat     = get(hdl, 'userdata');
    ALLEEG      = userdat{1};
    datasetinfo = userdat{2};
    page        = userdat{3};
    allcom      = userdat{4};
    clusterpresent = userdat{5};
    STUDY       = userdat{6};
    tags        = userdat{7};

    switch  com
        case 'subject'
            guiindex  = varargin{1};
            realindex = guiindex+(page-1)*10;

            datasetinfo(realindex).subject   = varargin{2};
            if get(findobj(hdl, 'tag', 'copy_to_dataset'), 'value')
                ALLEEG(realindex).subject    = varargin{2};
            end
            allcom = { allcom{:}, { 'index' realindex 'subject' varargin{2} } };
            userdat{1} = ALLEEG;
            userdat{2} = datasetinfo;
            userdat{4} = allcom;
            set(hdl, 'userdata', userdat);

        case 'run'
            guiindex  = varargin{1};
            realindex = guiindex+(page-1)*10;
            datasetinfo(realindex).run = str2num(varargin{2});
            if get(findobj(hdl, 'tag', 'copy_to_dataset'), 'value')
                ALLEEG(realindex).run  = str2num(varargin{2});
            end
            allcom = { allcom{:}, { 'index' realindex 'run' str2num(varargin{2}) } };
            userdat{1} = ALLEEG;
            userdat{2} = datasetinfo;
            userdat{4} = allcom;
            set(hdl, 'userdata', userdat);

        case 'session'
            guiindex  = varargin{1};
            realindex = guiindex+(page-1)*10;
            datasetinfo(realindex).session   = str2num(varargin{2});
            if get(findobj(hdl, 'tag', 'copy_to_dataset'), 'value')
                ALLEEG(realindex).session    = str2num(varargin{2});
            end
            allcom = { allcom{:}, { 'index' realindex 'session' str2num(varargin{2}) } };
            userdat{1} = ALLEEG;
            userdat{2} = datasetinfo;
            userdat{4} = allcom;
            set(hdl, 'userdata', userdat);

        case 'group'
            guiindex  = varargin{1};
            realindex = guiindex+(page-1)*10;
            datasetinfo(realindex).group   = varargin{2};
            if get(findobj(hdl, 'tag', 'copy_to_dataset'), 'value')
                ALLEEG(realindex).group    = varargin{2};
            end
            allcom = { allcom{:}, { 'index' realindex 'group' varargin{2} } };
            userdat{1} = ALLEEG;
            userdat{2} = datasetinfo;
            userdat{4} = allcom;
            set(hdl, 'userdata', userdat);

        case 'condition'
            guiindex  = varargin{1};
            realindex = guiindex+(page-1)*10;
            datasetinfo(realindex).condition   = varargin{2};
            if get(findobj(hdl, 'tag', 'copy_to_dataset'), 'value')
                ALLEEG(realindex).condition     = varargin{2};
            end
            allcom = { allcom{:}, { 'index' realindex 'condition' varargin{2} } };
            userdat{1} = ALLEEG;
            userdat{2} = datasetinfo;
            userdat{4} = allcom;
            set(hdl, 'userdata', userdat);

        case 'dipselect'
            STUDY.datasetinfo = datasetinfo;

            res = inputdlg2_with_checkbox( { strvcat('Enter maximum residual (topo map - dipole proj.) var. (in %)', ...
                                       'NOTE: This will delete any existing component clusters!') }, ...
                                       'pop_study():  Pre-select components', 1, { '15' },'pop_study' );

            if isempty(res), return; end
            if res{2} == 1
                STUDY = std_editset(STUDY, ALLEEG, 'commands', { 'inbrain', 'on', 'dipselect' str2num(res{1})/100 'return' });
                allcom = { allcom{:}, { 'inbrain', 'on', 'dipselect' str2num(res{1})/100 } };
            else
                STUDY = std_editset(STUDY, ALLEEG, 'commands', { 'inbrain', 'off','dipselect' str2num(res{1})/100 'return' });
                allcom = { allcom{:}, { 'inbrain', 'off', 'dipselect' str2num(res{1})/100 } };
            end

            datasetinfo   = STUDY.datasetinfo;

            userdat{2} = datasetinfo;
            userdat{4} = allcom;
            set(hdl, 'userdata', userdat);
            set(findobj(hdl, 'tag', 'delclust'), 'value', 1);
            pop_study('delclust', hdl);
            pop_study('redraw', hdl);

        case 'component'
            guiindex  = varargin{1};
            realindex = guiindex+(page-1)*10;

            for index = 1:size(ALLEEG(realindex).icaweights,1)
                complist{index} = [ 'IC ' int2str(index) ];
            end
            [tmps,tmpv] = listdlg2('PromptString', 'Select components', 'SelectionMode', ...
                                    'multiple', 'ListString', strvcat(complist), 'initialvalue', datasetinfo(realindex).comps);
            if tmpv ~= 0 % no cancel

                % find other subjects with the same session
                % -----------------------------------------
                for index = 1:length(datasetinfo)
                    if realindex == index || (strcmpi(datasetinfo(index).subject, datasetinfo(realindex).subject) && ...
                                ~isempty(datasetinfo(index).subject) && ...
                                isequal( datasetinfo(index).session, datasetinfo(realindex).session ) && ...
                                isequal( datasetinfo(index).run, datasetinfo(realindex).run ) )
                        datasetinfo(index).comps = tmps;
                        allcom = { allcom{:}, { 'index' index 'comps' tmps } };
                        set(findobj('tag', [ 'comps' int2str(index) ]), ...
                            'string', formatbut(tmps), 'horizontalalignment', 'left');
                    end
                end

            end
            userdat{2} = datasetinfo;
            userdat{4} = allcom;
            set(hdl, 'userdata', userdat);
            pop_study('redraw', hdl);

        case 'clear'
            guiindex  = varargin{1};
            realindex = guiindex+(page-1)*10;

            datasetinfo(realindex).filename  = '';
            datasetinfo(realindex).filepath  = '';
            datasetinfo(realindex).subject   = '';
            datasetinfo(realindex).session   = [];
            datasetinfo(realindex).run       = [];
            datasetinfo(realindex).condition = '';
            datasetinfo(realindex).group     = '';
            datasetinfo(realindex).comps     = [];

            allcom = { allcom{:}, { 'remove' realindex } };
            userdat{1} = ALLEEG;
            userdat{2} = datasetinfo;
            userdat{4} = allcom;
            set(hdl, 'userdata', userdat);
            pop_study('redraw', hdl);

        case 'nextpage'
            userdat{3} = page+1;
            set(hdl, 'userdata', userdat);
            pop_study('redraw', hdl);

        case 'prevpage'
            userdat{3} = max(1,page-1);
            set(hdl, 'userdata', userdat);
            pop_study('redraw', hdl);

        case 'load'
            guiindex  = varargin{1};
            filename  = varargin{2};
            realindex = guiindex+(page-1)*10;

            % load dataset
            % ------------
            TMPEEG = pop_loadset('filename', filename, 'loadmode', 'info');
            ALLEEG = eeg_store(ALLEEG, eeg_checkset(TMPEEG), realindex);

            % update datasetinfo structure
            % ----------------------------
            datasetinfo(realindex).filename  = ALLEEG(realindex).filename;
            datasetinfo(realindex).filepath  = ALLEEG(realindex).filepath;
            datasetinfo(realindex).subject   = ALLEEG(realindex).subject;
            datasetinfo(realindex).session   = ALLEEG(realindex).session;
            datasetinfo(realindex).run       = ALLEEG(realindex).run;
            datasetinfo(realindex).condition = ALLEEG(realindex).condition;
            datasetinfo(realindex).group     = ALLEEG(realindex).group;
            datasetinfo(realindex).comps     = [];

            allcom = { allcom{:}, { 'index' realindex 'load' filename } };
            userdat{1} = ALLEEG;
            userdat{2} = datasetinfo;
            userdat{4} = allcom;
            set(hdl, 'userdata', userdat);
            pop_study('redraw', hdl);

         case 'delclust'
            if clusterpresent
                if ~get(findobj(hdl, 'tag', 'delclust'), 'value')
                    for k = 1:10
                        set(tags.(['set'   num2str(k)]), 'style', 'text');
                        set(tags.(['comps' num2str(k)]), 'enable', 'off');
                        set(tags.(['sess'  num2str(k)]), 'enable', 'off');
                        set(tags.(['brw'   num2str(k)]), 'enable', 'off');
                    end
                else
                    for k = 1:10
                        set(tags.(['set'   num2str(k)]), 'style', 'edit');
                        set(tags.(['comps' num2str(k)]), 'enable', 'on');
                        set(tags.(['sess'  num2str(k)]), 'enable', 'on');
                        set(tags.(['brw'   num2str(k)]), 'enable', 'on');
                    end
                end
            else
                set(findobj(hdl, 'tag', 'delclust'), 'value', 0)
                if nargin > 2
                    warndlg2('No cluster present');
                end
            end

        case 'redraw'
            for k = 1:10
                kk = k+(page-1)*10; % real index
                if kk > length(datasetinfo)
                    set(tags.(['num' num2str(k)]), 'string', int2str(kk));
                    set(tags.(['set' num2str(k)]), 'string', '');
                    set(tags.(['sub' num2str(k)]), 'string','');
                    set(tags.(['sess' num2str(k)]), 'string','');
                    set(tags.(['run'  num2str(k)]), 'string','');
                    set(tags.(['cond' num2str(k)]), 'string','');
                    set(tags.(['comps' num2str(k)]), 'string','');
                    set(tags.(['group' num2str(k)]), 'string','');
                else
                    set(tags.(['num' num2str(k)]), 'string', int2str(kk));
                    set(tags.(['set' num2str(k)]), 'string', fullfile(char(datasetinfo(kk).filepath), char(datasetinfo(kk).filename)));
                    set(tags.(['sub' num2str(k)]), 'string', char(datasetinfo(kk).subject));
                    set(tags.(['sess' num2str(k)]), 'string', int2str(datasetinfo(kk).session));
                    set(tags.(['run' num2str(k)]),  'string', int2str(datasetinfo(kk).run));
                    set(tags.(['cond' num2str(k)]), 'string', char(datasetinfo(kk).condition));
                    set(tags.(['comps' num2str(k)]), 'string', formatbut(datasetinfo(kk).comps));
                    set(tags.(['group' num2str(k)]), 'string', char(datasetinfo(kk).group));
                end
            end
            if page<10
                 pagestr =  [ ' Page ' int2str(page) ];
            else pagestr =  [ 'Page ' int2str(page) ];
            end
            set(tags.('page'), 'string', pagestr );
    end
end

% remove empty elements in allcom
% -------------------------------
function allcom = simplifycom(allcom);

    for index = length(allcom)-1:-1:1
        if strcmpi(allcom{index}{1}, 'index') && strcmpi(allcom{index+1}{1}, 'index')
            if allcom{index}{2} == allcom{index+1}{2} % same dataset index
                allcom{index}(end+1:end+length(allcom{index+1})-2) = allcom{index+1}(3:end);
                allcom(index+1) = [];
            end
        end
    end

% test for wrong parameters
% -------------------------
function bool = test_wrong_parameters(hdl)
    userdat     = get(hdl, 'userdata');
    datasetinfo = userdat{2};
    datastrinfo = userdat{1};

    bool = 0;
    for index = 1:length(datasetinfo)
        if ~isempty(datasetinfo(index).filename)
            if isempty(datasetinfo(index).subject) && bool == 0
                bool = 1; warndlg2('All datasets must have a subject name or code', 'Error');
            end
        end
    end

    nonempty     = cellfun('isempty', { datasetinfo.filename });
    anysession   = any(~cellfun('isempty', { datasetinfo(nonempty).session }));
    allsession   = all(~cellfun('isempty', { datasetinfo(nonempty).session }));
    anyrun       = any(~cellfun('isempty', { datastrinfo(nonempty).run}));
    allrun       = all(~cellfun('isempty', { datastrinfo(nonempty).run}));
    anycondition = any(~cellfun('isempty', { datasetinfo(nonempty).condition }));
    allcondition = all(~cellfun('isempty', { datasetinfo(nonempty).condition }));
    anygroup     = any(~cellfun('isempty', { datasetinfo(nonempty).group }));
    allgroup     = all(~cellfun('isempty', { datasetinfo(nonempty).group }));
    anydipfit    = any(~cellfun('isempty', { datastrinfo(nonempty).dipfit}));
    alldipfit    = all(~cellfun('isempty', { datastrinfo(nonempty).dipfit}));

    if anygroup && ~allgroup
         bool = 1; warndlg2('If one dataset has a group label, they must all have one', 'Error');
    end
    if anycondition && ~allcondition
         bool = 1; warndlg2('If one dataset has a condition label, they must all have one', 'Error');
    end
    if anysession && ~allsession
         bool = 1; warndlg2('If one dataset has a session index, they must all have one', 'Error');
    end
    if anyrun && ~allrun
         bool = 1; warndlg2('If one dataset has a session index, they must all have one', 'Error');
    end
    if anydipfit && ~alldipfit
         bool = 1; warndlg2('Dipole''s data across datasets is not uniform');
    end
function strbut = formatbut(complist)
    if isempty(complist)
        strbut = 'All comp.';
    else
        if length(complist) > 3, strbut = [ 'Comp.: ' int2str(complist(1:2)) ' ...' ];
        else                     strbut = [ 'Comp.: ' int2str(complist) ];
        end
    end

%---------------------- helper functions -------------------------------------

function [result] = inputdlg2_with_checkbox(Prompt,Title,LineNo,DefAns,funcname);

if nargin < 4
   help inputdlg2;
   return;
end;
if nargin < 5
	funcname = '';
end

if length(Prompt) ~= length(DefAns)
	error('inputdlg2: prompt and default answer cell array must have the same size');
end

geometry = {};
listgui = {};

% determine if vertical or horizontal
% -----------------------------------
geomvert = [];
for index = 1:length(Prompt)
	geomvert = [geomvert size(Prompt{index},1) 1];  % default is vertical geometry
end
if all(geomvert == 1) && length(Prompt) > 1
	geomvert = []; % horizontal
end

for index = 1:length(Prompt)
	if ~isempty(geomvert) % vertical
		geometry = { geometry{:} [ 1] [1 ]};
	else
		geometry = { geometry{:} [ 1 0.6 ]};
	end
	listgui = { listgui{:} { 'Style', 'text', 'string', Prompt{index}}  ...
				{ 'Style', 'edit', 'string', DefAns{index} } { 'Style', 'checkbox', 'string','Keep only in-brain dipoles (requires Fieldtrip extension).','value',1 }  };
end
geometry = [1 1 1];geomvert = [2 1 1];
result = inputgui(geometry, listgui, ['pophelp(''' funcname ''');'], Title, [], 'normal', geomvert);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [ Iframe, Annotation, I, bbox_orig, bbox, OrigPoints ] = getImageFrame( options, idx, annotation_struct )
%GETIMAGEFRAME return normalized frame and its annotation
%   Given image from db and detected frame of face on it this function
%   expands detected frame by options.bw_margin and than crops image to
%   this frame and rescales it to size of options.bw (base window)
%
% INPUT
%
%   options     ...     structure containing parameters - bw (base_window),
%                       bw_margin (base window margin), components (marix with bboxes of components),
%                       image_path (path to directory containing jpg files of faces)
%   idx         ...     index of image in db
%   image       ...     cell image with filename and  detected frame
%
% OUTPUT
%
%   Iframe      ...     detected frame normalized to bw size after
%                       expansion of bw_margin in both dimensions
%   Annotation  ...     structure with annotation in similar format as
%                       original XML annotation
%   I           ...     original Image (rgb)
%   bbox_orig   ...     original bbox of detected face
%   bbox        ...     extended bbox by options.bw_margin
%   OrigPoints  ...     original landmarks (in unnormalized image)
%
% 02-08-10 Michal Uricar
% 20-10-10 Michal Uricar, normalized image frame get by getNormalizedFrame.m
% 12-07-11 Michal Uricar, corners dataset
% 21-03-12 Michal Uricar, LFW annotation in one file

    M = 10;
    fname = strtok(annotation_struct.names{idx}, '.');
    bbox_orig = double(annotation_struct.bbox(idx, :));
    filename = [options.image_path annotation_struct.names{idx}];
    I = rgb2gray(imread(filename));

    % get normalized image frame
    try
        [Iframe, bbox] = getNormalizedFrame(I, bbox_orig, options);
    catch
        Iframe = [];
        Annotation = [];
        bbox = [];
        OrigPoints = [];
        return;
    end;

    % compute coordinates of landmarks in normalized frame
    OrigPoints = [annotation_struct.eye_r(idx, :)' annotation_struct.eye_l(idx, :)' ,...
        annotation_struct.canthus_rr(idx, :)' annotation_struct.canthus_rl(idx, :)' ,...
        annotation_struct.canthus_lr(idx, :)' annotation_struct.canthus_ll(idx, :)' ,...
        annotation_struct.mouth(idx, :)' annotation_struct.mouth_corner_r(idx, :)'  ,...
        annotation_struct.mouth_corner_l(idx, :)' annotation_struct.nose(idx, :)'];
    O = [bbox(1); bbox(2)];
    scalefac = [ options.bw(1)/(bbox(3)-bbox(1)+1); options.bw(2)/(bbox(4)-bbox(2)+1) ];

    if (size(OrigPoints, 2) ~= M)
        error('??? Count of landmarks does not equal M. Corrupted xml annotation?');
        Annotation = [];
        fprintf('Bad annotation found.\n');
        return;
    end;

    % save annotation
    Points = zeros(2, M);
    for j = 1 : M
        Points(:, j) = (OrigPoints(:, j) - O) .* scalefac;
    end;

    if (numel(Points(Points < 0)) > 0)
        Annotation = [];
        fprintf('Bad annotation found.\n');
        return;
    end;

    [h w d] = size(Iframe);

    Annotation.P = round(Points);
    Annotation.image.filename = [fname '.jpg'];
    Annotation.image.width = w;
    Annotation.image.height = h;
end

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
% signed velocity in the direction of the closest fly, according to type
function [data,units] = compute_veltoward(trx,n,type)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);

for i1 = 1:nflies,
  fly1 = flies(i1);

  % fly closest to fly1 according to type
  closestfly = trx(fly1).(['closestfly_',type]);

  % velocity of fly1
  dx1 = diff(trx(fly1).x_mm,1,2);
  dy1 = diff(trx(fly1).y_mm,1,2);
  x_mm1 = trx(fly1).x_mm;
  y_mm1 = trx(fly1).y_mm;

  % loop over all flies
  for i2 = 1:nflies,

    fly2 = flies(i2);
    if i1 == i2, continue; end

    % frames where this fly is closest
    idx = find(closestfly(1:end-1) == fly2);

    % don't use the last frame of fly2
    off = trx(fly1).firstframe - trx(fly2).firstframe;
    idx(idx+off == trx(fly2).nframes) = [];

    if isempty(idx), continue; end

    % unit vector in direction of fly2 from fly1
    off = trx(fly1).firstframe - trx(fly2).firstframe;
    dx2 = trx(fly2).x_mm(off+idx)-x_mm1(idx);
    dy2 = trx(fly2).y_mm(off+idx)-y_mm1(idx);
    dz2 = sqrt(dx2.^2 + dy2.^2);
    dx2 = dx2 ./ dz2;
    dy2 = dy2 ./ dz2;
    dx2(dz2==0) = 0;
    dy2(dz2==0) = 0;

    % project velocity of fly1 onto this vector
    data{i1}(idx) = dx1(idx).*dx2 + dy1(idx).*dy2;

  end
end

units = parseunits('mm/s');

author: kristinbranson

repo: JAABA


text:
clc;
close all;
clearvars;
rng('default');
n = 10;
% Let's get a symmetric positive definite matrix
A = gallery('gcdmat', n);
% Compute its Cholesky decomposition
L = chol(A, 'lower')
b = randn(n, 1);
% solve the problem
x = L \b
x = spx.la.tris.forward_row(L, b)

author: indigits

repo: sparse-plex


text:
function [ node_x, element_node ] = gmsh_data_read ( gmsh_filename, ...
  node_dim, node_num, element_order, element_num )

%*****************************************************************************80
%
%% GMSH_DATA_READ reads data from a GMSH file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string GMSH_FILENAME, the GMSH filename.
%
%    Input, integer NODE_DIM, the spatial dimension.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer ELEMENT_ORDER, the order of the elements.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Output, real NODE_X(NODE_DIM,NODE_NUM), the node coordinates.
%
%    Output, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM),
%    the nodes that make up each element.
%
  node_x = zeros ( node_dim, node_num );
  element_node = zeros ( element_order, element_num );
%
%  Open the file.
%
  input = fopen ( gmsh_filename, 'rt' );

  if ( input < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GMSH_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the input file "%s".\n', gmsh_filename );
    error ( 'GMSH_DATA_READ - Error!' );
    return
  end

  level = 0;

  while ( 1 )

    text = fgetl ( input );

    if ( text == -1 )
      break
    end

    if ( level == 0 )
      if ( s_begin ( text(1:6), '$Nodes' ) )
        level = 1;
        j = 0;
      end
    elseif ( level == 1 )
      [ value, count ] = sscanf ( text, '%d' );
      level = 2;
    elseif ( level == 2 )
      if ( s_begin ( text(1:9), '$EndNodes' ) )
        break
      else
        j = j + 1;
        [ value, count ] = sscanf ( text, '%d  %f  %f  %f' );
        indx = value(1);
        node_x(1,j) = value(2);
        if ( 2 < count )
          node_x(2,j) = value(3);
          if ( 3 < count )
            node_x(3,j) = value(4);
          end
        end
      end
    end

  end
%
%  Now read element information.
%
  level = 0;

  while ( 1 )

    text = fgetl ( input );

    if ( text == -1 )
      fprintf ( 'ran out\n' );
      break
    end

    if ( level == 0 )
      if ( s_begin ( text(1:9), '$Elements' ) )
        level = 1;
        j = 0;
      end
    elseif ( level == 1 )
      [ value, count ] = sscanf ( text, '%d' );
      level = 2;
    elseif ( level == 2 )
      if ( s_begin ( text(1:12), '$EndElements' ) )
        break
      else
        j = j + 1;
        [ value, count ] = sscanf ( text, '%d' );
        for i = 1 : element_order
          element_node(i,j) = value(5+i);
        end
      end
    end

  end

  fclose ( input );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function xCommon=findLongestCommonSubsequence(x,y)
%%FINDLONGESTCOMMONSUBSEQUENCE Given two vectors of numbers, characters,
%              etc., find the longest sequence of elements in both of them.
%              The elements must appear in order need NOT BE CONSECUTIVE.
%
%INPUTS:      x,y  Two vectors.
%
%OUTPUTS: xCommon  The longest common subsequence shared by the vectors.
%                  The elements in the subsequence need not be consecutive.
%                  For example, the longest common subsequence of "dunkin"
%                  and "doughnuts" is "dun".
%
%The algorithm is taken from Chapter 15.4 of [1], but modified so that the
%sequence comes out in order and to deal with Matlab addressing things from
%1 instead of 0.
%
%REFERENCES:
%[1] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein,
%    Introduction to Algorithms, 2nd ed. Cambridge, MA: The MIT Press,
%    2001.
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

m=length(x);
n=length(y);

c=zeros(m+1,n+1);
for i=m:-1:1
    for j=n:-1:1
        if(x(i)==y(j))
            c(i,j)=c(i+1,j+1)+1;
        else
            c(i,j)=max(c(i,j+1),c(i+1,j));
        end
    end
end

i=1;
j=1;
%Setting xCommon to x rather than alocating with zeros assured that xCommon
%will have the same type and shape as x. This is good when x is a character
%string.
xCommon=x;
numFound=0;
while(i<=m&&j<=n)
    if(x(i)==y(j))
       numFound=numFound+1;
       xCommon(numFound)=x(i);
       i=i+1;
       j=j+1;
    elseif(c(i+1,j)>=c(i,j+1))
        i=i+1;
    else
        j=j+1;
    end
end

%Get rid of extra elements in xCommon.
xCommon=xCommon(1:numFound);
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
 function ob = Gblur(mask, varargin)
%function ob = Gblur(mask, options)
%|
%| Construct Gblur object for image restoration.
%|
%| See Gblur_test.m for example usage.
%|
%| in
%|	mask	size(image)	logical array of object support.
%|
%| options
%|	'chat'		verbose printing of debug messages
%|	'psf'		point spread function (aka impulse response)
%|	'type'		type of blur:
%|				'conv,same'	usual case (default)
%|				'conv,per'	(periodic end conditions)
%|				'conv,valid'	(realistic extended input)
%|				'fft,same'	(periodic end conditions)
%|				todo: allow replicated end conditions!
%|				'imfilter,same'	(requires image toolbox)
%|				'imfilter,circ'	(requires image toolbox)
%|				'imfilter,mirror' (requires image toolbox)
%|					(do not use - adjoint fails)
%|				'sparse'	todo: return sparse matrix
%|	'imfilter_options'	options to 'imfilter' (if used)
%|
%| out
%|	ob [nd np]	np = sum(mask(:)), so it is already "masked"
%|			nd = prod(size(mask)) for 'conv,same' type
%|
%| Copyright 2005-4-22, Jeff Fessler, University of Michigan

if nargin == 1 && streq(mask, 'test'), Gblur_test, return, end
if nargin < 1, ir_usage, end

arg.mask = mask;

% option defaults
arg.chat = 0;
arg.psf = 1; % identity matrix
arg.type = 'conv,same';
arg.class = 'fatrix2';
arg.imfilter_options = {};

% options specified by name/value pairs
arg = vararg_pair(arg, varargin);

if ndims(arg.psf) ~= ndims(mask), error 'psf dim mismatch', end
if any(size(arg.psf) > size(mask)), error 'psf too large', end

if streq(arg.type, 'imfilter,circ') && exist('imfilter') ~= 2
	warn 'no imfilter (no image processing toolbox) so using fft'
	arg.type = 'fft,same';
end

if streq(arg.type, 'imfilter,mirror') && exist('imfilter') ~= 2
	fail 'no imfilter (no image processing toolbox)'
end

if streq(arg.type, 'imfilter,same') && exist('imfilter') ~= 2
	warn 'no imfilter (no image processing toolbox) so using slower convn'
	arg.type = 'conv,same';
end

arg.psf = Gblur_mask_psf_odd(arg.psf);

psf_flip = conj(flipdims(arg.psf, 'odd', 1));
does_many = false;

arg.odim = size(mask);
if ndims(mask) == 2 && size(mask,2) == 1 % trick for 1D
	arg.odim = size(mask,1);
end
arg.idim = arg.odim;

switch arg.type

case 'conv,same'
	arg.psf_flip = psf_flip;
	handle_forw = @(arg,x) convn(x, arg.psf, 'same');
	handle_back = @(arg,y) fatrix2_maskit(arg.mask, ...
			convn(y, arg.psf_flip, 'same'));
	does_many = true;

case 'conv,per'
	arg.psf_flip = psf_flip;
	handle_forw = @(arg,x) ir_conv(x, arg.psf, 'per', true);
	handle_back = @(arg,y) fatrix2_maskit(arg.mask, ...
			ir_conv(y, arg.psf_flip, 'per', true));

case 'conv,valid'
	arg.psf_flip = psf_flip;
	handle_forw = @(arg,x) convn(x, arg.psf, 'valid');
	handle_back = @(arg,y) convn(y, arg.psf_flip);
	arg.odim = arg.idim - size(arg.psf) + 1;
	does_many = true;

case 'fft,same'
	[arg, handle_forw, handle_back] = Gblur_setup_fft_same(arg);

case 'imfilter,circ'
	arg.psf_flip = psf_flip;
	handle_forw = @(arg,x) imfilter(x, arg.psf, 'conv', 'same', ...
		'circular', arg.imfilter_options{:});
	handle_back = @(arg,y) fatrix2_maskit(arg.mask, ...
			imfilter(y, arg.psf_flip, 'conv', 'same', ...
				'circular', arg.imfilter_options{:}));

case 'imfilter,mirror'
	arg.psf_flip = psf_flip;
	handle_forw = @(arg,x) imfilter(x, arg.psf, 'conv', 'same', ...
		'symmetric', arg.imfilter_options{:});
	handle_back = @(arg,y) fatrix2_maskit(arg.mask, ...
			imfilter(y, arg.psf_flip, 'conv', 'same', ...
				'symmetric', arg.imfilter_options{:}));

case 'imfilter,same'
	arg.psf_flip = psf_flip;
	handle_forw = @(arg,x) imfilter(x, arg.psf, 'conv', 'same', ...
		arg.imfilter_options{:});
	handle_back = @(arg,y) fatrix2_maskit(arg.mask, ...
			imfilter(y, arg.psf_flip, 'conv', 'same', ...
				arg.imfilter_options{:}));

case 'sparse'
	ob = Gblur_sparse(arg.psf, arg.idim);
	ob = ob(:,arg.mask(:));
	return

otherwise
	error 'unknown blur type'
end

switch arg.class

case 'Fatrix'

	arg.nd = prod(arg.odim);
	arg.np = sum(mask(:));
	dim = [arg.nd arg.np]; % trick: make it masked by default!

	ob = Fatrix(dim, arg, 'caller', 'Gblur', ...
		'abs', @Gblur_abs, ...
		'forw', @Gblur_forw_Fatrix, 'back', @Gblur_back_Fatrix);

case 'fatrix2'

	ob = fatrix2('idim', arg.idim, 'odim', arg.odim, ...
		'arg', arg, 'does_many', does_many, ...
		'abs', @Gblur_abs, 'forw', handle_forw, 'back', handle_back);

otherwise
	fail('bug')
end

% Gblur_mask_psf_odd(psf)
% having an odd-sized psf facilites adjoint
function psf = Gblur_mask_psf_odd(psf)
nd = ndims(psf);
for id=1:nd
	sz = size(psf, id);
	if ~mod(sz, 2) % even
		warn('size(psf, %d) = %d is even; appending 0', id, sz)
		tmp = size(psf);
		tmp(id) = 1;
		z = zeros(tmp);
		psf = cat(id, psf, z);
	end
end

% Gblur_sparse()
% a_ij = b[ n(i) - n(j), m(i) - m(j) ]
% n(i) = (i-1) mod N
% m(i) = floor((i-1) / N)
function sp = Gblur_sparse(psf, idim)
if numel(idim) ~= 2
	fail 'not done'
end
fail 'todo: not done'
%sp = sparse(i, j, s, length(ii), size(ob.arg.G, 2));

% Gblur_abs(): |A| for abs(A)
function ob = Gblur_abs(ob)
ob.arg.psf = abs(ob.arg.psf);
if isfield(ob.arg, 'psf_flip')
	ob.arg.psf_flip = abs(ob.arg.psf_flip);
end
if isfield(ob.arg, 'psf_fft') % fixed 2012-04-08
	ob.arg.psf_fft = Gblur_setup_fft_same_fft(ob.arg.psf, ob.mask);
end

% Gblur_setup_fft_same()
function [arg, handle_forw, handle_back] = Gblur_setup_fft_same(arg)

arg.psf_fft = Gblur_setup_fft_same_fft(arg.psf, arg.mask);

handle_forw = @(arg,x) ifftn(fftn(x) .* arg.psf_fft);
handle_back = @(arg,y) fatrix2_maskit(arg.mask, ...
			ifftn(fftn(y) .* conj(arg.psf_fft)));

% Gblur_setup_fft_same_fft()
function psf_fft = Gblur_setup_fft_same_fft(psf, mask)

% put psf in center of array
tmp = zeros(size(mask));
switch ndims(psf)
case 2
	[n1 n2] = size(mask);
	[p1 p2] = size(psf);
%	h1 = (size(psf,1)-1)/2;
%	h2 = (size(psf,2)-1)/2;
%	i1 = floor(n1/2) + 1 + [-h1:h1]
%	i2 = floor(n2/2) + 1 + [-h2:h2]
	i1 = floor(n1/2) + 1 + (1:p1) - floor((p1+1)/2);
	i2 = floor(n2/2) + 1 + (1:p2) - floor((p2+1)/2);
	tmp(i1,i2) = psf;
case 3
	if any(~rem(size(psf),2)), error 'psf size must be odd', end
	h1 = (size(psf,1)-1)/2;
	h2 = (size(psf,2)-1)/2;
	h3 = (size(psf,3)-1)/2;
	i1 = floor(size(mask,1)/2) + 1 + [-h1:h1];
	i2 = floor(size(mask,2)/2) + 1 + [-h2:h2];
	i3 = floor(size(mask,3)/2) + 1 + [-h3:h3];
	tmp(i1,i2,i3) = psf;
otherwise
	error 'only 2d & 3d done'
end
psf_fft = fftn(ifftshift(tmp));

% Gblur_forw_Fatrix(): y = A * x
function y = Gblur_forw_Fatrix(arg, x)

[x ei] = embed_in(x, arg.mask, arg.np);

switch arg.type
case 'conv,same'
	y = convn(x, arg.psf, 'same');
case 'conv,per'
	y = ir_conv(x, arg.psf, 'per', true);
case 'fft,same'
	if ndims(x) > ndims(arg.mask);
		y = fft_conv_multi(x, arg.psf_fft);
	else
		y = ifftn(fftn(x) .* arg.psf_fft);
	end
case 'imfilter,circ'
	y = ir_imfilter_many(x, arg.psf, 'conv', 'same', ...
		'circular', arg.imfilter_options{:});
case 'imfilter,same'
	y = ir_imfilter_many(x, arg.psf, 'conv', 'same', ...
		arg.imfilter_options{:});
otherwise
	error 'bug'
end

y = ei.shape(y);

% Gblur_back_Fatrix(): x = A' * y
function x = Gblur_back_Fatrix(arg, y)

[y eo] = embed_out(y, arg.odim);

switch arg.type
case 'conv,same'
	x = convn(y, arg.psf_flip, 'same');
case 'conv,per'
	x = ir_conv(y, arg.psf_flip, 'per', true);
case 'fft,same'
	if ndims(y) > ndims(arg.mask);
		x = fft_conv_multi(y, conj(arg.psf_fft));
	else
		x = ifftn(fftn(y) .* conj(arg.psf_fft));
	end
case 'imfilter,circ'
	x = ir_imfilter_many(y, arg.psf_flip, 'conv', 'same', ...
		'circular', arg.imfilter_options{:});
case 'imfilter,same'
	x = ir_imfilter_many(y, arg.psf_flip, 'conv', 'same', ...
		arg.imfilter_options{:});
otherwise
	error 'bug'
end

x = eo.shape(x, arg.mask, arg.np);

% fft_conv_multi()
function y = fft_conv_multi(x, H)
dim = size(x);
y = zeros(size(x));
y = [];
for ii=1:dim(end)
	tmp = ifftn(fftn(stackpick(x, ii)) .* H);
	y = cat(length(dim), y, tmp); % slow: keeps enlarging y.
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [ KRlm ] = AMICO_RotateKernel( K, AUX, idx_IN, idx_OUT, isIsotropic )

	if ( isIsotropic == false )
		% fit SH and rotate kernel to 181*181 directions
		KRlm = zeros( size(AUX.fit,1), 181, 181, 'single' );
		for ox = 1:181
		for oy = 1:181
			Ylm_rot = AUX.Ylm_rot{ ox, oy };
			for s = 1 : numel(idx_IN)
				Klm = AUX.fit * K( idx_IN{s} );		% fit SH of shell to rotate
				Rlm = zeros( size(Klm) );
				idx = 1;
				for l = 0 : 2 : AUX.lmax
					const = sqrt(4.0*pi/(2.0*l+1.0)) * Klm( (l*l + l + 2.0)/2.0 );
					for m = -l : l
						Rlm(idx) = const * Ylm_rot(idx);
						idx = idx+1;
					end
				end
				KRlm( idx_OUT{s}, ox, oy ) = single( Rlm );
			end
		end
		end
	else
		% simply fit SH
		KRlm = zeros( size(AUX.fit,1), 1, 1, 'single' );
		Ylm_rot = AUX.Ylm_rot{ 1, 1 };
		for s = 1:numel(idx_IN)
			KRlm( idx_OUT{s}, 1, 1 ) = single( AUX.fit * K( idx_IN{s} ) );
		end
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
%% The class fibre
%
%%
%

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
classdef PTKPoint
    % PTKPoint. A class for storing points
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        CoordX
        CoordY
        CoordZ
    end

    methods
        function obj = PTKPoint(c_x, c_y, c_z)
            if nargin > 0
                obj.CoordX = c_x;
                obj.CoordY = c_y;
                obj.CoordZ = c_z;
            end
        end
    end

    methods (Static)
        function diff = Difference(point_1, point_2)
            diff = PTKPoint(point_1.CoordX - point_2.CoordX, point_1.CoordY - point_2.CoordY, point_1.CoordZ - point_2.CoordZ);
        end

        function mag = Magnitude(point)
            mag = norm([point.CoordX, point.CoordY, point.CoordZ]);
        end

        function point = SetCoordinate(point, dimension, value)
            switch dimension
                case PTKImageOrientation.Coronal
                    point.CoordY = value;
                case PTKImageOrientation.Sagittal
                    point.CoordX = value;
                case PTKImageOrientation.Axial
                    point.CoordZ = value;
                otherwise
                    error('Unknown dimension');
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function Yq = QuantLloyd (Nlev, FPDF, QSym)
%  Iterate to find the output levels for a minimum mean square error
%  quantizer.
%
% This subroutine searches for a set of quantizer output levels which are
% the conditional means (centroids) of the quantizer decision regions. The
% decision boundaries are assumed to lie mid-way between output levels.
%
% The procedure is based on a Lloyd-Max iteration.
% (1) Given the start of an interval and an output level, find the end of
%     the interval, such that the output level is the centroid of the
%     probability density function in that interval.
% (2) Use the output level and upper interval edge to find the output
%     level in the next interval (the interval edge must lie midway
%     between output levels).
% (3) With the start and output level of the next interval
%     determined, repeat step (1) for this interval.
% The iteration continues by modifying the initial output level based
% on whether the centroid property of the last interval was satisfied.
%    S. P. Lloyd, "Least Squares Quantization in PCM", IEEE Trans. Inform.
%      Theory, vol. 28, no. 2, pp. 129-137, March 1982.
%    J. Max, "Quantizing for Minimum Distortion", IEEE Trans. Inform.
%      Theory, vol. 6, no. l, pp. 7-12, March 1960.
%
% Nlev - Number of quantizer output levels. For symmetric quantizers
%   this is the number of levels above the mean.
% FPDF - Cell array of function pointers {Farea, Fmean, Fvar}
% QSym - Symmetry flag (optional, default 0)
%   0 - Quantizer not constrained to be symmetric
%   1 - Quantizer is symmetric with an odd number of levels. The middle
%       level is fixed at the mean. This routine finds the Nlev output
%       levels above the mean.
%   2 - Quantizer is symmetric with an even number of levels. The middle
%       decision level is at the mean. This routine finds the Nlev output
%       levels above the mean.
%
% Yq - Nlev output levels in ascending order

% There are 3 symmetry cases to consider for symmetry
% QSym == 0: No symmetry. The first interval is defined by a lower
%   boundary XL = -Inf and an output level (centroid) Yc. The centroid
%   position is iterated.
% QSym == 1: Symmetrical quantizer with fixed output level at the mean.
%   If we place another output level Yc, the decision level XL lies
%   midway between the mean and Yc. When we iterate Yc, the decision level
%   is also be readjusted.
% QSym == 2: Symmetric quantizer with a decision boundary XL at the mean.
%   The output level Yc is iterated.

if (nargin < 3)
  QSym = 0;
end

Fmean = FPDF{2};
Fvar = FPDF{3};

% Parameters
MaxIter = 100;
TolR = 1e-5;

% Set the optimization parameters
Xmean = feval(Fmean, -Inf, Inf);
sd = feval(Fvar, -Inf, Inf) - Xmean^2;

% Convergence criterion
Tol = TolR * sd;

if (QSym == 0)                  % No symmetry
  XL = -Inf;
  Yc = Xmean - 4 * sd;
  Xstep = sd;
elseif (QSym == 1)              % Symmetric, odd number of coefficients
  Xstep = 2 * sd / Nlev;
  Yc = Xmean + Xstep;
  XL = 0.5 * (Xmean + Yc);
elseif (QSym == 2)              % Symmetric, even number of coefficients
  Xstep = 2 * sd / Nlev;
  XL = Xmean;
  Yc = Xmean + Xstep;
end

% Initialization
Ytrial = Yc;
Ybase = Yc;
Xstep = 0.5 * Xstep;
FL = false;         % True if an upper bound has been found
FU = false;         % True if a lower bound has been found

for (Iter = 1:MaxIter)

% Find a set of output levels satisfying the necessary conditions
% for a minimum mean square error quantizer
% XU is the largest quantizer decision level
%  Yc = Ytrial;
  if (QSym == 1)
    XL = 0.5 * (Xmean + Ytrial);
  end

  [Yq, XU] = QuantLevel(Nlev, FPDF, XL, Ytrial);

  if (isnan(XU) || Yq(end) == XU)
    FU = true;
  else
    FL = true;
    Ybase = Ytrial;
  end

% Check for convergence, adjust the step size
  if (FL && FU)
    if (Xstep < Tol)
      break
    end
    Xstep = 0.5 * Xstep;
    Ytrial = Ybase + Xstep;
  elseif (FL)     % Need to extend the search upward
    Xstep = 2 * Xstep;
    Ytrial = Ybase + Xstep;
  else            % Need to back up the initial value
    Ybase = Ytrial;
    Xstep = 2 * Xstep;
    if (~isinf(XL))
      Xstep = min(Xstep, 0.5*(Ybase - XL)); % Don't step below XL
    end
    Ytrial = Ybase - Xstep;
  end

end

if (Iter >= MaxIter)
  error('QuantLloyd: Failed to converge');
end
if (~FU || ~FL)
  error('QuantLloyd: Feasible solution not found');
end
fprintf('QuantLloyd: Converged, %d iterations\n', Iter);

return

% ----- -----
function [Yq, XU] = QuantLevel(Nlev, FPDF, XL, Yc)
% Find a set of quantizer output levels, given the lower boundary and
% centroid of the first interval.
%
% Given the lower boundary (decision level) for an initial interval and the
% output level (centroid) of that interval, this routine first finds the
% upper decision boundary for that interval. These values are telescoped to
% give the lower boundary and output level for the second interval. This
% process continues for each interval.
%
% The last decision level is returned by this routine. If this value is
% finite, the last interval does not fully encompass the tail of the
% probability density function. If the last decision level is NaN, the last
% output level is not the centroid of the last region extenting to infinity.
%
% Nlev - Number of quantizer output levels. For symmetric quantizers this
%   is the number of levels above the mean.
% FPDF - Cell array of function handles {Farea, Fmean, Fvar}
% XL - Lower boundary of the first interval
% Yc - Centroid of the first interval
%
% Yq - Nlev quantizer output levels in ascending order. Trailing values may
%   be NaN if it is not possible to have intervals with these output levels
%   as the centroids of the corresponding intervals. The quantizer decision
%   levels lie midway between output levels.
% XU - Largest decision level (greater than or equal to Yq(Nlev)) or
%   NaN.

% If XU is finite, then the levels should be adjusted upward (increase
% Yc). If XU is NaN, the last output level is not the centroid of the last
% interval, and the levels should be adjusted downward. Another case occurs
% if an entire interval is zero probability. Then the upper decision level
% falls on the output level for that interval. Test for XU == Yq(Nlev).

Yq = zeros(1, Nlev);
for (i = 1:Nlev)

% Given the lower decision level and the output level for an interval,
% find the upper decision level
  if (isnan(XL))
    XU = XL;
    Yc = XL;
  else
    XU = QuantInterval(FPDF, XL, Yc);
  end

% Given the interval limits just found, telescope to find the
% output level to be used in the next iteration
  Yq(i) = Yc;
  XL = XU;
  Yc = XU + (XU - Yc);

end

return

% ----- ------
function Xb = QuantInterval (FPDF, Xa, Yc)
% Find the upper edge of an interval which has a given centroid.
%
% This routine solves for upper limit of the integral
%    Xb
%   Int (x-Yc) p(x) dx = 0.
%    Xa
%
% The routine is designed to allow for Yc to be less than Xa, in which case
% Xb <= Yc, or for Yc to be greater than Xa, in which case Xb >= Yc.
%
% FPDF - Cell array of function handles {Farea, Fmean, Fvar}
% Xa - Lower boundary of the interval
% Yc - Centroid of the interval
%
% Xb - Returned value representing the upper boundary of the interval. This
%   value is set to NaN if there is no solution for Xb on the opposite side
%   of Yc from Xa.

% Parameters
XstepR = 1.2;    % Initial relative step size
TolF = 1e-5;     % Function amplitude relative tolerance
TolX = 1e-5;     % Position relative tolerance
MaxIter = 100;
EpsdF = 0.01;    % Choose between bisection or linear interpolation
EpsdX = 0.05;    % For linear interpolation, constrain the relative
                 % step size to EpsdX <= dX <= 1-EpsdX.

% Searching for a solution of the equation F(Xb) = 0. Consider the
% case that Yc > Xa. The function F(Xb) is zero at Xb = Xa. It becomes
% negative with increasing Xb (since p(x) is positive). It takes on its
% most negative value at Xb = Yc. It then decreases as Xb increases. It
% is the second zero crossing we seek.

% The search procedure has as its stopping criteria:
%  a) abs(F(Xb) < TolF*abs(F(YC)).
%  b) The interval of uncertainty is less than TolX*abs(Xb)
%  c) The probability from the present trial point to Inf or -Inf (as
%     appropriate) is zero. In this case Xb is set to NaN.
%  d) The maximum number of iterations is exceeded.
%  e) F(Yc)=0. This indicates that the probability density function is
%     zero in the interval (Xa,Yc).

if (Yc == Xa)
  Xb = Yc;
  return
end

Farea = FPDF{1};
Fmean = FPDF{2};
Fvar = FPDF{3};

% Evaluate the function at Yc to check if the probability
% is zero, and to determine the stopping criterion Feps
Fm = feval(Fmean, Xa, Yc) - Yc*feval(Farea, Xa, Yc);  % Should be negative
if (Fm >= 0)
  if (Fm > 0)
    error('QuantInterval: Error, function value at centroid positive');
  end
  Xb = Yc;          % Fm == 0;
  return
end

% Set up the boundaries of the search and the step size
FL = Fm;            % Value at lower boundary (initially negative)
XL = Yc;            % Lower boundary of search
if (~isinf(Xa))

  Xb = Yc + XstepR * (Yc - Xa);  % Initial trial upper boundary

else

% Find the standard deviation of the distribution
  Xmean = feval(Fmean, -Inf, Inf);
  sd = sqrt(feval(Fvar, -Inf, Inf) - Xmean^2);
  Xb = Yc + sign(Yc-Xa) * sd;   % Initial test upper boundary
end

FR = -1;            % Same sign as FL to indicate no zero found
XR = Xb;
Feps = TolF * abs(Fm); % Tolerance on integral value
Xstep = 0.5 * (Xb - Yc);

% Search loop
for (Iter = 1:MaxIter)

  Fp = Fm;          % Previous Fm
  Fm = feval(Fmean, Xa, Xb) - Yc*feval(Farea, Xa, Xb);

% Update the end-points of the search
if (Fm >= 0)
  FR = Fm;
  XR = Xb;
else
  FL = Fm;
  XL = Xb;
end

% ----- ------
  if (FR >= 0)
% Straddling a root

% Check for convergence in the function value
% Check for convergence of the position
    if (abs(Fm) <= Feps) || ...
       (~isinf(XL) && ~isinf(XR) && ...
        abs(XR-XL) <= TolX*max(abs(XR), abs(XL)))
      return
    end

% The root location is sought by cautious linear interpolation; however
% if successive function values are nearly equal, bisection is used.
    if (abs(Fp-Fm) > EpsdF * abs(FL-FR))
      dX = FL / (FL-FR);    % Estimated zero crossing position
      dX = max(min(dX, 1-EpsdX), EpsdX); % Avoid region close to XL or XR
    else
      dX = 0.5;             % Bisection
    end
    Xb = XL + dX*(XR-XL);

% ------ ------
  else
% Not straddling a root

% Right boundary undefined
% Check for successive identical returned values
    if (Yc > Xa)
      if (Fm == Fp && feval(Farea, Xb, Inf) <= 0)
        Xb = NaN;
        return
      end
    else
      if (Fm == Fp && feval(Farea, -Inf, Xb) <= 0)
        Xb = NaN;
        return
      end
    end

% Increase the step size
    Xb = Xb + Xstep;
    Xstep = 2 * Xstep;

  end

end

error('QuantInterval: Failed to converge');

author: Sable

repo: mcbench-benchmarks


text:
function tests = test_ft_preproc_resample

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_preproc_resample

if nargout
  % assume that this is called by RUNTESTS
  tests = functiontests(localfunctions);
else
  % assume that this is called from the command line
  fn = localfunctions;
  for i=1:numel(fn)
    feval(fn{i});
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function testOptions(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nchan   = 8;
nsample = 1000;
dat     = randn(nchan, nsample) + 1;

result = {};
result{end+1} = ft_preproc_resample(dat, 1000, 100, 'resample');
result{end+1} = ft_preproc_resample(dat, 1000, 250, 'resample');
result{end+1} = ft_preproc_resample(dat, 1000, 500, 'resample');
result{end+1} = ft_preproc_resample(dat, 1000, 100, 'decimate');
result{end+1} = ft_preproc_resample(dat, 1000, 250, 'decimate');
result{end+1} = ft_preproc_resample(dat, 1000, 500, 'decimate');
result{end+1} = ft_preproc_resample(dat, 1000, 100, 'downsample');
result{end+1} = ft_preproc_resample(dat, 1000, 250, 'downsample');
result{end+1} = ft_preproc_resample(dat, 1000, 500, 'downsample');
result{end+1} = ft_preproc_resample(dat, 1000, 100, 'fft');
result{end+1} = ft_preproc_resample(dat, 1000, 250, 'fft');
result{end+1} = ft_preproc_resample(dat, 1000, 500, 'fft');

% all iterations were done with (slightly) different options, hence the results should not be equal
for i=1:numel(result)
  for j=(i+1):numel(result)
    assert(~isequal(result{i}, result{j}), 'the results %d and %d should not be equal', i, j);
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [out] = saturation_10(p1,p2,p3,S,In)
%saturation_10

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Saturation excess flow from a store with different degrees
%               of saturation (min-max exponential variant)
% Constraints:  -
% @(Inputs):    p1   - maximum contributing fraction area [-]
%               p2   - minimum contributing fraction area [-]
%               p3   - exponentia scaling parameter [-]
%               S    - current storage [mm]
%               In   - incoming flux [mm/d]

out = min(p1,p2+p2.*exp(p3.*S)).*In;

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function [struct_irf_record D_record gamma_record]=irfres_relmagnitude_panel(beta_gibbs,sigma_gibbs,It,Bu,IRFperiods,n,m,p,k,signrestable,signresperiods, relmagrestable, relmagresperiods)

% function [struct_irf_record D_record gamma_record Qdraw Qsuccess]=bear.irfres(beta_gibbs,sigma_gibbs,It,Bu,IRFperiods,n,m,p,k,signrestable,signresperiods)
% runs the gibbs sampler to obtain draws from the posterior distribution of IRFs, orthogonalised with a sign restriction setting
% inputs:  - matrix 'beta_gibbs': record of the gibbs sampler draws for the beta vector
%          - matrix 'sigma_gibbs': record of the gibbs sampler draws for the sigma matrix (vectorised)
%          - integer 'It': total number of iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - integer 'Bu': number of burn-in iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - integer 'IRFperiods': number of periods for IRFs
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'm': number of exogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'p': number of lags included in the model (defined p 7 of technical guide)
%          - integer 'k': number of coefficients to estimate for each equation in the BVAR model (defined p 7 of technical guide)
%          - cell 'signrestable': table recording the sign restriction input from the user
%          - cell 'signresperiods': table containing the periods corresponding to each restriction
% outputs: - cell 'struct_irf_record': record of the gibbs sampler draws for the orthogonalised IRFs
%          - matrix 'D_record': record of the gibbs sampler draws for the structural matrix D
%          - matrix 'gamma_record': record of the gibbs sampler draws for the structural disturbances variance-covariance matrix gamma
%          - integer 'Qdraw': total number of draws of the Q matrix
%          - integer 'Qsuccess': number of successful draws of the Q matrix

% preliminary tasks
% create first the cell that will store the results from the simulations
struct_irf_record=cell(n,n);
% storage cell
storage1=cell(It-Bu,1);
storage2=cell(It-Bu,1);

% now identify all the periods concerned with restrictions
% first expand the non-empty entries in signresperiods since they are only expressed in intervals: transform into list
% for instance, translate [1 4] into [1 2 3 4]; I don't think this can done without a loop
temp=cell2mat(signresperiods(~cellfun(@isempty,signresperiods)));
periods=[];
for ii=1:size(temp,1)
periods=[periods temp(ii,1):temp(ii,2)];
end
% suppress duplicates and sort
periods=sort(unique(periods))';
% count the total number of restriction periods (required for IRF matrix)
nperiods=size(periods,1);

% now identify all the periods concerned with relative magnitude
% restrictions
% first expand the non-empty entries in magresperiods since they are only expressed in intervals: transform into list
% for instance, translate [1 4] into [1 2 3 4];
temp=cell2mat(relmagresperiods(~cellfun(@isempty,relmagresperiods)));
mperiods=[];
for ii=1:size(temp,1)
mperiods=[mperiods temp(ii,1):temp(ii,2)];
end
% suppress duplicates and sort
mperiods=sort(unique(mperiods))';
% count the total number of restriction periods (required for IRF matrix)
rmperiods=size(periods,1);

% Identify the restriction matrices
% create five cells, corresponding to the three possible restrictions:
% one cell for sign restrictions, three cells for magnitude restrictions, one cell for zero restrictions
Scell=cell(1,n);
Mcell=cell(1,n);
Mlcell=cell(1,n);
Mucell=cell(1,n);
Zcell=cell(1,n);

% Check if value and periods restrictions correspond to each other
if sum(sum(~cellfun(@isempty,signresperiods) == ~cellfun(@isempty,signrestable))) == n^2
    % All cells with sign restrictions also specify the horizon over which
    % these are applied
else
    disp('Warning: Value restrictions do not correspond to period restrictions one to one')
    pause(1)
end

% loop over rows and columns of the period matrix
for ii=1:n
   for jj=1:n
      % if entry (ii,jj) of the period matrix and of the value matrix is not empty...
      if ~isempty(signresperiods{ii,jj}) && ~isempty(signrestable{ii,jj})
      % ... then there is a restriction over one (or several) periods
      % loop overt those periods
         for kk=signresperiods{ii,jj}(1,1):signresperiods{ii,jj}(1,2)
         % identify the position of the considered period within the list of all periods (required to build the matrix)
         position=find(periods==kk);
         % now create the restriction matrix: this will depend on the type of restriction
            % if it is a positive sign restriction...
            if strcmp(signrestable{ii,jj},'+')
            % ... then input a 1 entry in the corresponding S matrix
            Scell{1,jj}=[Scell{1,jj};zeros(1,n*nperiods)];
            Scell{1,jj}(end,(position-1)*n+ii)=1;
            % if it is a negative sign restriction...
            elseif strcmp(signrestable{ii,jj},'-')
            % ... then input a -1 entry in the corresponding S matrix
            Scell{1,jj}=[Scell{1,jj};zeros(1,n*nperiods)];
            Scell{1,jj}(end,(position-1)*n+ii)=-1;
            % if it is a zero restriction...
            elseif strcmp(signrestable{ii,jj},'0')
            % ... then input a 1 entry in the corresponding Z matrix
            Zcell{1,jj}=[Zcell{1,jj};zeros(1,n*nperiods)];
            Zcell{1,jj}(end,(position-1)*n+ii)=1;
            % else, a non-empty entry being neither a sign nor a zero restriction has to be a magnitude restriction
            else
            % fill the corresponding M matrices:
            % input a 1 in M
            Mcell{1,jj}=[Mcell{1,jj};zeros(1,n*nperiods)];
            Mcell{1,jj}(end,(position-1)*n+ii)=1;
            % input the lower value of the interval in Ml
            temp=str2num(signrestable{ii,jj});
            Mlcell{1,jj}=[Mlcell{1,jj};temp(1,1)];
            % input the upper value of the interval in Mu
            Mucell{1,jj}=[Mucell{1,jj};temp(1,2)];
            end
         end
      end
   end
end

%create matrix entry for relative magnitude restrictions (on impact)
[r clm] = find(~cellfun('isempty',relmagrestable));
%2. Indentify which entry corresponds to the positive magnitude
%restriction (which shock is supposed to have a larger impact on which
%variable)
num_magres=length(r)/2; %number of relative magnitude restrictions
IndextempL=double.empty;
kk=1; %number of the restrictions
IndextempS=double.empty;
kk=1; %%number of restriction

rowsS = [];
columnsS = [];
for jj=1:num_magres %%loop over number of magnitude restrictions
strtemp = strcat('S',num2str(jj)); %%find entry in the table corresponding to the Stronger than restriction
Stronger = strcmp(relmagrestable, strtemp);
[rowS columnS] = find(Stronger==1);
rowsS = [rowsS rowS];
columnsS = [columnsS columnS];
end

rowsW = [];
columnsW = [];
for jj=1:num_magres
strtemp = strcat('W',num2str(jj));
Weaker = strcmp(relmagrestable, strtemp);
[rowW columnW] = find(Weaker==1);
rowsW = [rowsW rowW];
columnsW = [columnsW columnW];
end

% now check what kind of restrictions apply among sign, zero and magnitude restrictions
% check for sign restrictions: if there are any, at least one entry in the cell Scell is non-empty
if sum(~cellfun(@isempty,Scell))~=0
signres=1;
else
signres=0;
end
% similarly check for zero restrictions
if sum(~cellfun(@isempty,Zcell))~=0
zerores=1;
else
zerores=0;
end
% and finally, check for magnitude restrictions
if sum(~cellfun(@isempty,Mcell))~=0
magnres=1;
else
magnres=0;
end
if length(columnsS)~=0
relmagnres=1;
else
relmagnres=0;
end

% initiate Gibbs algorithm
not_successful = 0;
hbar = bear.parfor_progressbar(It-Bu,'Progress of Sign Restriction Draws');  %create the progress bar
for ii=1:It-Bu
% initiate the variable 'success'; this variable will be used to check whether the restrictions are satisfied
% if there are only zero restrictions, they will be satisfied by construction, and 'success' will simply be ignored
success=0;

% how the algorithm will be conducted will depend on the types of restrictions implemented

   % if there are only zero restrictions, the algorithm is simple as no checking is required: the conditions are satisfied by construction
   if zerores==1 && signres==0 && magnres==0 && relmagnres==0
   % draw beta and sigma
   beta=beta_gibbs(:,ii);
   sigma=reshape(sigma_gibbs(:,ii),n,n);
   hsigma=chol(bear.nspd(sigma),'lower');
   % obtain orthogonalised IRFs
   [irfmatrix ortirfmatrix]=bear.irfsim(beta,hsigma,n,m,p,k,max(IRFperiods,max(periods)));
   % generate the stacked IRF matrix
   stackedirfmat=[];
      for jj=1:numel(periods)
      stackedirfmat=[stackedirfmat;ortirfmatrix(:,:,periods(jj,1)+1)];
      end
   % draw an entire random matrix Q satisfying the zero restrictions
   [Q]=bear.qzerores(n,Zcell,stackedirfmat);
   % there is no need to verify the restrictions: there are satisfied by construction

   % if there are sign/magnitude restrictions, possibly associated with zero restrictions
   else
   % the algorithm becomes a bit more complicated as conditions now need to be checked
   % to maintain efficiency, the algorithm proceeds recursively shock by shock, and stops as soon as a condition on the considered shock fails
      % repeat algorithm for the iteration as long as not all conditions are satisfied
      while success==0
      not_successful = not_successful+1;
      % switch 'success' to 1; it will be turned back to zero if at any time Q is detected as a candidate not satisfying the restrictions
      success=1;
      % draw randomly the vector of VAR coefficients: draw a random index
      index=floor(rand*(It-Bu))+1;
      % then draw a random set of beta and sigma corresponding to this index (this is done to make it possible to draw, if required, an infinite number of values from the gibbs sampler record, with equal probability on each value)
      beta=beta_gibbs(:,index);
      sigma=reshape(sigma_gibbs(:,index),n,n);
      hsigma=chol(bear.nspd(sigma),'lower');
      % obtain orthogonalised IRFs
      [irfmatrix ortirfmatrix]=bear.irfsim(beta,hsigma,n,m,p,k,max(IRFperiods,max(periods)));
      % generate the stacked IRF matrix
      stackedirfmat=[];
         for jj=1:numel(periods)
         stackedirfmat=[stackedirfmat;ortirfmatrix(:,:,periods(jj,1)+1)];
         end
      % initiate Qj
      Qj=[];
      % now start looping over the shocks and checking sequentially whether conditions on these shocks hold
      % stop as soon as one restriction fails
         jj=1;
         while success==1 && jj<=n
         % build column j of the random matrix Q
         [qj]=bear.qrandj(n,Zcell{1,jj},stackedirfmat,Qj);
         % obtain the candidate column fj
         fj=stackedirfmat*qj;
         % check restrictions: first sign restrictions
         [success qj]=bear.checksignres(Scell{1,jj},qj,fj);
         % if 'success' is still equal to 1, also check for magnitude restrictions
            if success==1
            [success]=bear.checkmagres(Mcell{1,jj},Mlcell{1,jj},Mucell{1,jj},fj);
            end
         % also, if 'success' is still equal to 1, update Qj by concatenating qj
            if success==1
            Qj=[Qj qj];
            end

%once all n columns are build and fullfill the sign restrictions, check
%relative magnitudes
            jj=jj+1;
           if size(Qj,2)==n && success==1  && relmagnres==1
             %disp('I reached magnitude restrictions')
             D=hsigma*Qj;
             [~, ortirfmatrixmagnitude]=bear.irfsim(beta,D,n,m,p,k,max(IRFperiods,max(mperiods)));
             % generate the stacked IRF matrix
              stackedirfmatmagn=[];
              for kk=1:numel(mperiods)
              stackedirfmatmagn=[stackedirfmatmagn;ortirfmatrixmagnitude(:,:,mperiods(kk,1)+1)];
              end
              [success]=bear.checkrelmag(stackedirfmatmagn,columnsS, columnsW, rowsS, rowsW, n, mperiods);
              if success ==1
              %disp('Sign and Magnitude Restrictions fullfilled')
              else
              disregarded=ortirfmatrixmagnitude;
              end
            end
        end
      % repeat this loop until a succesful draw is obtained
   end
% with succesful Qj at hand, eventually set Q as Qj
Q=Qj;
end

  % store
   for jj=1:IRFperiods
   storage1{ii,1}(:,:,jj)=ortirfmatrix(:,:,jj)*Q;
   end
   storage2{ii,1}=hsigma*Q;
   hbar.iterate(1);   % update progress by one iteration
end
close(hbar);   %close progress bar

% reorganise storage
% loop over iterations
for ii=1:It-Bu
   % loop over IRF periods
   for jj=1:IRFperiods
      % loop over variables
      for kk=1:n
         % loop over shocks
         for ll=1:n
         struct_irf_record{kk,ll}(ii,jj)=storage1{ii,1}(kk,ll,jj);
         end
      end
   end
D_record(:,ii)=storage2{ii,1}(:);
gamma_record(:,ii)=bear.vec(eye(n));
end

fprintf('Accepted Draws in Percent of Total Number of Draws: %f', 100*(It-Bu)/(not_successful + It-Bu))

author: european-central-bank

repo: BEAR-toolbox


text:
function [P,R,S] = lagrangepoly(X,Y,XX)
%LAGRANGEPOLY  Lagrange interpolation polynomial fitting a set of points
%   [P,R,S] = LAGRANGEPOLY(X,Y)  where X and Y are row vectors
%   defining a set of N points uses Lagrange's method to find
%   the N-1th order polynomial in X that passes through these
%   points.  P returns the N coefficients defining the polynomial,
%   in the same order as used by POLY and POLYVAL (highest order first).
%   Then, polyval(P,X) = Y.  R returns the x-coordinates of the N-1
%   extrema of the resulting polynomial (roots of its derivative),
%   and S returns the y-values  at those extrema.
%
%   YY = LAGRANGEPOLY(X,Y,XX) returns the values of the polynomial
%   sampled at the points specified in XX -- the same as
%   YY = POLYVAL(LAGRANGEPOLY(X,Y)).
%
%   Example:
%   To find the 4th-degree polynomial that oscillates between
%   1 and 0 across 5 points around zero, then plot the interpolation
%   on a denser grid inbetween:
%     X = -2:2;  Y = [1 0 1 0 1];
%     P = lagrangepoly(X,Y);
%     xx = -2.5:.01:2.5;
%     plot(xx,polyval(P,xx),X,Y,'or');
%     grid;
%   Or simply:
%     plot(xx,lagrangepoly(X,Y,xx));
%
%   Note: if you are just looking for a smooth curve passing through
%   a set of points, you can get a better fit with SPLINE, which
%   fits piecewise polynomials rather than a single polynomial.
%
%   See also: POLY, POLYVAL, SPLINE

% 2006-11-20 Dan Ellis dpwe@ee.columbia.edu
% $Header: $

%  For more info on Lagrange interpolation, see Mathworld:
%  http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html

% Make sure that X and Y are row vectors
if size(X,1) > 1;  X = X'; end
if size(Y,1) > 1;  Y = Y'; end
if size(X,1) > 1 || size(Y,1) > 1 || size(X,2) ~= size(Y,2)
  error('both inputs must be equal-length vectors')
end

N = length(X);

pvals = zeros(N,N);

% Calculate the polynomial weights for each order
for i = 1:N
  % the polynomial whose roots are all the values of X except this one
  pp = poly(X( (1:N) ~= i));
  % scale so its value is exactly 1 at this X point (and zero
  % at others, of course)
  pvals(i,:) = pp ./ polyval(pp, X(i));
end

% Each row gives the polynomial that is 1 at the corresponding X
% point and zero everywhere else, so weighting each row by the
% desired row and summing (in this case the polycoeffs) gives
% the final polynomial
P = Y*pvals;

if nargin==3
  % output is YY corresponding to input XX
  YY = polyval(P,XX);
  % assign to output
  P = YY;
end

if nargout > 1
  % Extra return arguments are values where dy/dx is zero
  % Solve for x s.t. dy/dx is zero i.e. roots of derivative polynomial
  % derivative of polynomial P scales each power by its power, downshifts
  R = roots( ((N-1):-1:1) .* P(1:(N-1)) );
  if nargout > 2
    % calculate the actual values at the points of zero derivative
    S = polyval(P,R);
  end
end


text:
function M = getmassmatvec(elem2edge,area,Dlambda,elemType,K)
%% GETMASSMATVEC Get the mass matrix of vector finite element space
%
% M = GETMASSMATVEC(elem2edge,area,Dlambda,elemType,K) get mass matrix of
% the finite element space specified by elemType. The coefficient K is
% piecewise constant.
%
% The elemType can be:
%
% - "RT0": The lowest order Raviart-Thomas element
% - "BDM1": The lowest order Brezzi-Douglas-Michel element
% - "BDM1B": The BDM element enriched by the curl of cubic bubble function
% - "ND0": The lowest order Nedelec element
% - "RT1": Linear-Quadratic Rvairart-Thomas element
%
% Note that RT0 and ND0 share the same mass matrix.
%
% RT0, BMD1 are created by Ming Wang at July, 2012. Improved by Long Chen.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('K','var'), K = []; end
NE = double(max(elem2edge(:)));
NT = size(elem2edge,1);
DiDj = zeros(NT,3,3);
for i = 1:3
    for j = i:3
        if isempty(K)
            DiDj(:,i,j) = dot(Dlambda(:,:,i),Dlambda(:,:,j),2);
        else
            switch size(K,2)
                case 1 % scalar K
                    DiDj(:,i,j) = dot(Dlambda(:,:,i),Dlambda(:,:,j),2)./K;
                case 2 % diagonal matrix
                    DiDj(:,i,j) =  Dlambda(:,1,i).*Dlambda(:,1,j)./K(:,1) ...
                                 + Dlambda(:,2,i).*Dlambda(:,2,j)./K(:,2);
                case 3 % K is 2 by 2 SPD matrix
                    detK = K(:,1).*K(:,2) - K(:,3).^2;
                    DiDj(:,i,j) = (Dlambda(:,1,i).*Dlambda(:,1,j).*K(:,2) ...
                                 + Dlambda(:,2,i).*Dlambda(:,2,j).*K(:,1) ...
                                 - Dlambda(:,1,i).*Dlambda(:,2,j).*K(:,3) ...
                                 - Dlambda(:,2,i).*Dlambda(:,1,j).*K(:,3))./detK;

            end
        end
        DiDj(:,j,i) = DiDj(:,i,j);
    end
end
localEdge = [2 3; 1 3; 1 2]; % ascend ordering

%% RT0 and ND0
if strcmp(elemType,'RT0') || strcmp(elemType,'ND0')
    M = sparse(NE,NE);
    for i = 1:3
        for j = i:3
            % local to global index map and its sign
            ii = double(elem2edge(:,i));
            jj = double(elem2edge(:,j));
            i1 = localEdge(i,1); i2 = localEdge(i,2); % [i1,i2] is the edge opposite to vertex i.
            j1 = localEdge(j,1); j2 = localEdge(j,2);
            % computation of mass matrix --- (phi_i, phi_j)
            Mij = 1/12*area.*( (1+(i1==j1))*DiDj(:,i2,j2) ...
                             - (1+(i1==j2))*DiDj(:,i2,j1) ...
                             - (1+(i2==j1))*DiDj(:,i1,j2) ...
                             + (1+(i2==j2))*DiDj(:,i1,j1));
            if (j==i)
                M = M + sparse(ii,jj,Mij,NE,NE);
            else
                M = M + sparse([ii;jj],[jj;ii],[Mij; Mij],NE,NE);
            end
        end
    end
end

%% BDM1
if strcmp(elemType,'BDM1') || strcmp(elemType,'BDM1B') || strcmp(elemType,'RT1')
    M = sparse(2*NE,2*NE);
    for i = 1:3
        for j = i:3
            % local to global index map and its sign
            ii = double(elem2edge(:,i));
            jj = double(elem2edge(:,j));
            i1 = localEdge(i,1); i2 = localEdge(i,2); % [i1,i2] is the edge opposite to vertex i.
            j1 = localEdge(j,1); j2 = localEdge(j,2);
            % computation of mass matrix, note that (rot u, rot v) = (grad u, grad v)
            % (phi_i, phi_j)
            Mij = 1/12*area.*( (1+(i1==j1))*DiDj(:,i2,j2) ...
                             - (1+(i1==j2))*DiDj(:,i2,j1) ...
                             - (1+(i2==j1))*DiDj(:,i1,j2) ...
                             + (1+(i2==j2))*DiDj(:,i1,j1));
            if (j==i)
                M = M + sparse(ii,jj,Mij,2*NE,2*NE);
            else
                M = M + sparse([ii;jj],[jj;ii],[Mij; Mij],2*NE,2*NE);
            end
            % (psi_i,psi_j)
            Mij = 1/12*area.*( (1+(i1==j1))*DiDj(:,i2,j2) ...
                             + (1+(i1==j2))*DiDj(:,i2,j1) ...
                             + (1+(i2==j1))*DiDj(:,i1,j2) ...
                             + (1+(i2==j2))*DiDj(:,i1,j1));
            if (j==i)
                M = M + sparse(ii+NE,jj+NE,Mij,2*NE,2*NE);
            else
                M = M + sparse([ii;jj]+NE,[jj;ii]+NE,[Mij; Mij],2*NE,2*NE);
            end
        end
    end
    for i = 1:3
        for j = 1:3
            % local to global index map and its sign
            ii = double(elem2edge(:,i));
            jj = double(elem2edge(:,j));
            i1 = localEdge(i,1); i2 = localEdge(i,2);
            j1 = localEdge(j,1); j2 = localEdge(j,2);
            % (psi_i,phi_j)
            Mij = 1/12*area.*( (1+(i1==j1))*DiDj(:,i2,j2) ...
                             - (1+(i1==j2))*DiDj(:,i2,j1) ...
                             + (1+(i2==j1))*DiDj(:,i1,j2) ...
                             - (1+(i2==j2))*DiDj(:,i1,j1));
            M = M + sparse([ii+NE;jj],[jj;ii+NE],[Mij; Mij],2*NE,2*NE);
        end
    end
end

%% BDM1B
if strcmp(elemType,'BDM1B')
    newM = sparse(2*NE+NT,2*NE+NT);
    newM(1:2*NE,1:2*NE) = M;
    M = newM;
    % (phi_i, bubble)
    for i = 1:3
        ii = double(elem2edge(:,i));
        jj = double((1:NT)');
        i1 = localEdge(i,1); i2 = localEdge(i,2); i3 = 6-i1-i2;
        Mij = 9/10*area.*(DiDj(:,i2,i3)...
                         +DiDj(:,i2,i2)...
                         -DiDj(:,i1,i3)...
                         -DiDj(:,i1,i1));
        M = M + sparse([ii;jj+2*NE],[jj+2*NE;ii],[Mij;Mij],2*NE+NT,2*NE+NT);
    end
    clear ii jj i1 i2 i3 Mij;
    % (psi_i, bubble)
    for i = 1:3
        ii = double(elem2edge(:,i));
        jj = double((1:NT)');
        i1 = localEdge(i,1); i2 = localEdge(i,2); i3 = 6-i1-i2;
        Mij = 9/10*area.*(DiDj(:,i2,i3)...
                        + DiDj(:,i1,i3)...
                        + DiDj(:,i2,i2)...
                        + DiDj(:,i1,i1)...
                        + DiDj(:,i1,i2));
        M = M + sparse([ii+NE;jj+2*NE],[jj+2*NE;ii+NE],[Mij;Mij],2*NE+NT,2*NE+NT);
    end
    clear ii jj i1 i2 i3 Mij;
    % (bubble, bubble)
    ii = 1+2*NE:NT+2*NE;
    Mij = 81/10*area.*(DiDj(:,1,1)...
                    + DiDj(:,1,2)...
                    + DiDj(:,1,3)...
                    + DiDj(:,2,2)...
                    + DiDj(:,2,3)...
                    + DiDj(:,3,3));
    if ~isempty(K)
        Mij = Mij./K;
    end
    M = M + sparse(ii,ii,Mij,2*NE+NT,2*NE+NT);
end

%% RT1
elem2dof = [elem2edge elem2edge+NE (1:NT)'+2*NE (1:NT)'+2*NE+NT];
locBasesIdx = [1 2 0; 1 3 0; 2 3 0; ... % phi and psi
               2 1 3; 3 1 2];           % chi
if strcmp(elemType,'RT1')
    newM = sparse(2*NE+2*NT,2*NE+2*NT);
    newM(1:2*NE,1:2*NE) = M;
    M = newM;
    % (phi_i, chi_j)
    for i = 1:3
        for j = 1:2
            ii = double(elem2dof(:,i));
            jj = double(elem2dof(:,j+6));
            i1 = locBasesIdx(i,1); i2 = locBasesIdx(i,2); %i3 = locBasesIdx(i,3);
            j1 = locBasesIdx(j,1); j2 = locBasesIdx(j,2); j3 = locBasesIdx(j,3);
            Mij = intlambda([j1,i1,j2],2)*DiDj(:,i2,j3) ...
                 -intlambda([j1,i1,j3],2)*DiDj(:,i2,j2) ...
                 -intlambda([j1,i2,j2],2)*DiDj(:,i1,j3) ...
                 +intlambda([j1,i2,j3],2)*DiDj(:,i1,j2);
            M = M + sparse([ii;jj],[jj;ii],[Mij;Mij],2*NE+2*NT,2*NE+2*NT);
        end
    end
    clear ii jj i1 i2 i3 Mij;
    % (psi_i, chi_j)
    for i = 1:3
        for j = 1:2
            ii = double(elem2dof(:,i+3));
            jj = double(elem2dof(:,j+6));
            i1 = locBasesIdx(i,1); i2 = locBasesIdx(i,2); %i3 = locBasesIdx(i,3);
            j1 = locBasesIdx(j,1); j2 = locBasesIdx(j,2); j3 = locBasesIdx(j,3);
            Mij = intlambda([j1,i1,j2],2)*DiDj(:,i2,j3) ...
                 -intlambda([j1,i1,j3],2)*DiDj(:,i2,j2) ...
                 +intlambda([j1,i2,j2],2)*DiDj(:,i1,j3) ...
                 -intlambda([j1,i2,j3],2)*DiDj(:,i1,j2);
            M = M + sparse([ii;jj],[jj;ii],[Mij;Mij],2*NE+2*NT,2*NE+2*NT);
        end
    end
    clear ii jj i1 i2 i3 Mij;
    for i = 1:2
        for j = 1:2
            % (chi_j,chi_i)
            ii = double(elem2dof(:,i+3));
            jj = double(elem2dof(:,j+6));
            i1 = locBasesIdx(i,1); i2 = locBasesIdx(i,2); i3 = locBasesIdx(i,3);
            j1 = locBasesIdx(j,1); j2 = locBasesIdx(j,2); j3 = locBasesIdx(j,3);
            Mij = intlambda([i1,j1,i2,j2],2)*DiDj(:,i3,j3) ...
                 -intlambda([i1,j1,i2,j3],2)*DiDj(:,i3,j2) ...
                 -intlambda([i1,j1,i3,j2],2)*DiDj(:,i2,j3) ...
                 +intlambda([i1,j1,i3,j3],2)*DiDj(:,i2,j2);
            M = M + sparse([ii;jj],[jj;ii],[Mij;Mij],2*NE+2*NT,2*NE+2*NT);
        end
    end
end

repo: ifem


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('Ex0702: Average Objects')

SpharmMatDir = '../';
CodeDir = [SpharmMatDir 'code'];
addpath(CodeDir);

DataFolder = 'data';

%% (1) Input data directory
InDataDir = [SpharmMatDir DataFolder '/Ex0702/hip05_prm'];

%% (2) List input file names
inFiles = dir([InDataDir '/*_prm.mat']); inNames={};
for i=1:length(inFiles)
    inNames{end+1} = [InDataDir '/' inFiles(i).name];
end

%% (3) Display input objects (optional)
    %Available values for Space- 'object';'param';'both'
dispConfs.Space = 'object';
    % Available values for Mesh- 'orig';'quad32';'quad64';'quad128';'quad256';'quad512';'icosa1';'icosa2'; ...
    %    'icosa3';'icosa4';'icosa5';'icosa6'
dispConfs.Mesh = 'orig';
    % Available values for Shape- 'solid';'mesh';'both'
dispConfs.Shade = 'both';
    % Available values for Overlay- 'none';'adc_paramap'
dispConfs.Overlay = 'none';
    % Available values for Export- 'screen';'png';'both'
dispConfs.Export = 'png';
dispConfs.Degree = [];
dispConfs.Template = '';

SpharmMatUtilDisplayObjs(dispConfs, inNames, CodeDir);

%% (4) Average objects
confs.OutputName='atlas.mat';
confs.OutDirectory = [SpharmMatDir DataFolder '/Ex0702/hip06_reg'];

if ~exist(confs.OutDirectory,'dir')
    mkdir(confs.OutDirectory);
end

SpharmMatUtilAverageObjs(confs, inNames);

%% (5) Display output objects (optional)
    %Available values for Space- 'object';'param';'both'
dispConfs.Space = 'object';
    % Available values for Mesh- 'orig';'quad32';'quad64';'quad128';'quad256';'quad512';'icosa1';'icosa2'; ...
    %    'icosa3';'icosa4';'icosa5';'icosa6'
dispConfs.Mesh = 'icosa4';
    % Available values for Shape- 'solid';'mesh';'both'
dispConfs.Shade = 'both';
    % Available values for Overlay- 'none';'adc_paramap'
dispConfs.Overlay = 'adc_paramap';
    % Available values for Export- 'screen';'png';'both'
dispConfs.Export = 'png';
dispConfs.Degree = [];
dispConfs.Template = '';

inFiles = dir([confs.OutDirectory '/*_als.mat']); inNames={};
for i=1:length(inFiles)
    inNames{end+1} = [confs.OutDirectory '/' inFiles(i).name];
end

SpharmMatUtilDisplayObjs(dispConfs, inNames, CodeDir);

inNames={}; inNames{1} = [confs.OutDirectory '/atlas.mat'];

SpharmMatUtilDisplayObjs(dispConfs, inNames, CodeDir);

rmpath(CodeDir);

clear;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus