text:
% pop_dipfit_settings() - select global settings for dipole fitting through a pop up window
%
% Usage:
%   >> OUTEEG = pop_dipfit_settings ( INEEG ); % pop up window
%   >> OUTEEG = pop_dipfit_settings ( INEEG, 'key1', 'val1', 'key2', 'val2' ... )
%
% Inputs:
%   INEEG	input dataset
%
% Optional inputs:
%   'hdmfile'  - [string] file containing a head model compatible with
%                the Fieldtrip dipolefitting() function ("vol" entry)
%   'mrifile'  - [string] file containing an anatomical MR head image.
%                The MRI must be normalized to the MNI brain. See the .mat
%                files used by the sphere and boundary element models
%                (For instance, select the sphere model and study 'EEG.dipfit').
%                If SPM2 software is installed, dipfit will be able to read
%                most MRI file formats for plotting purposes (.mnc files, etc...).
%                To plot dipoles in a subject MRI, first normalize the MRI
%                to the MNI brain using SPM2.
%   'coordformat' - ['MNI'|'Spherical'] Coordinates returned by the selected
%                head model. May be MNI coordinates or spherical coordinates
%                (For spherical coordinates, the head radius is assumed to be 85 mm.
%   'chanfile' - [string] template channel locations file. (This function will
%                check whether your channel locations file is compatible with
%                your selected head model).
%   'chansel'  - [integer vector] indices of channels to use for dipole fitting.
%                {default: all}
%   'coord_transform' - [float array] Talairach transformation matrix for
%                       aligning the dataset channel locations to the selected
%                       head model.
%   'electrodes'      - [integer array] indices of channels to include
%                       in the dipole model. {default: all}
% Outputs:
%   OUTEEG	output dataset
%
% Author: Arnaud Delorme, SCCN, La Jolla 2003-
%         Robert Oostenveld, SMI/FCDC, Nijmegen 2003

% MEG flag:
%   'gradfile' - [string] file containing gradiometer locations
%                ("gradfile" parameter in Fieldtrip dipolefitting() function)

% SMI, University Aalborg, Denmark http://www.smi.auc.dk/
% FC Donders Centre, University Nijmegen, the Netherlands http://www.fcdonders.kun.nl

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [OUTEEG, com] = pop_dipfit_settings ( EEG, varargin )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 1
   help pop_dipfit_settings;
   return;
end;

OUTEEG = EEG;
com = '';

% get the default values and filenames
dipfitdefs;

if nargin < 2

    if isstr(EEG) % setmodel
        tmpdat = get(gcf, 'userdata');
        chanfile = tmpdat.chanfile;
        tmpdat   = tmpdat.template_models;
        tmpval = get(findobj(gcf, 'tag', 'listmodels'), 'value');
        set(findobj(gcf, 'tag', 'model'), 'string', char(tmpdat(tmpval).hdmfile));
        set(findobj(gcf, 'tag', 'coord'), 'value' , fastif(strcmpi(tmpdat(tmpval).coordformat,'MNI'),2, ...
                                                     fastif(strcmpi(tmpdat(tmpval).coordformat,'CTF'),3,1)));
        set(findobj(gcf, 'tag', 'mri'  ), 'string', char(tmpdat(tmpval).mrifile));
        set(findobj(gcf, 'tag', 'meg'), 'string', char(tmpdat(tmpval).chanfile));
        set(findobj(gcf, 'tag', 'coregcheckbox'), 'value', 0);
        if tmpval < 3,
            set(findobj(gcf, 'userdata', 'editable'), 'enable', 'off');
        else,
            set(findobj(gcf, 'userdata', 'editable'), 'enable', 'on');
        end;
        if tmpval == 3,
            set(findobj(gcf, 'tag', 'headstr'), 'string', 'Subject CTF head model file (default.htm)');
            set(findobj(gcf, 'tag', 'mristr'),  'string', 'Subject MRI (coregistered with CTF head)');
            set(findobj(gcf, 'tag', 'chanstr'), 'string', 'CTF Res4 file');
            set(findobj(gcf, 'tag', 'manualcoreg'), 'enable', 'off');
            set(findobj(gcf, 'userdata', 'coreg'), 'enable', 'off');
        else,
            set(findobj(gcf, 'tag', 'headstr'), 'string', 'Head model file');
            set(findobj(gcf, 'tag', 'mristr'),  'string', 'MRI file');
            set(findobj(gcf, 'tag', 'chanstr'), 'string', 'Model template channel locations file');
            set(findobj(gcf, 'tag', 'manualcoreg'), 'enable', 'on');
            set(findobj(gcf, 'userdata', 'coreg'), 'enable', 'on');
        end;
        tmpl = tmpdat(tmpval).coord_transform;
        set(findobj(gcf, 'tag', 'coregtext'), 'string', '');
        set(findobj(gcf, 'tag', 'coregcheckbox'), 'value', 0);
        [allkeywordstrue transform] = lookupchantemplate(chanfile, tmpl);
        if allkeywordstrue,
            set(findobj(gcf, 'tag', 'coregtext'), 'string', char(vararg2str({ transform })));
            if isempty(transform)
                 set(findobj(gcf, 'tag', 'coregcheckbox'), 'value', 1);
            else set(findobj(gcf, 'tag', 'coregcheckbox'), 'value', 0);
            end;
        end;
        return;
    end;

    % detect DIPFIT1.0x structure
    % ---------------------------
    if isfield(EEG.dipfit, 'vol')
        str = [ 'Dipole information structure from DIPFIT v1.02 detected.' ...
                'Keep or erase the old dipole information including dipole locations? ' ...
                'In either case, a new dipole model can be constructed.' ];

        tmpButtonName=questdlg2( strmultiline(str, 60), 'Old DIPFIT structure', 'Keep', 'Erase', 'Keep');
        if strcmpi(tmpButtonName, 'Keep'), return; end;

    elseif isfield(EEG.dipfit, 'hdmfile')
        % detect previous DIPFIT structure
        % --------------------------------
        str = [ 'Dipole information and settings are present in the dataset. ' ...
                'Keep or erase this information?' ];
        tmpButtonName=questdlg2( strmultiline(str, 60), 'Old DIPFIT structure', 'Keep', 'Erase', 'Keep');
        if strcmpi(tmpButtonName, 'Keep'), return; end;
    end;

    % define the callbacks for the buttons
    % -------------------------------------
    cb_selectelectrodes = [ 'tmplocs = EEG.chanlocs; tmp = select_channel_list({tmplocs.label}, ' ...
                            'eval(get(findobj(gcbf, ''tag'', ''elec''), ''string'')));' ...
                            'set(findobj(gcbf, ''tag'', ''elec''), ''string'',[''[''  num2str(tmp) '']'']); clear tmplocs;' ]; % did not work
    cb_selectelectrodes = 'tmplocs = EEG.chanlocs; set(findobj(gcbf, ''tag'', ''elec''), ''string'', int2str(pop_chansel({tmplocs.labels}))); clear tmplocs;';
    cb_volmodel = [ 'tmpdat = get(gcbf, ''userdata'');' ...
                    'tmpind = get(gcbo, ''value'');' ...
                    'set(findobj(gcbf, ''tag'', ''radii''),   ''string'', num2str(tmpdat{tmpind}.r,3));' ...
                    'set(findobj(gcbf, ''tag'', ''conduct''), ''string'', num2str(tmpdat{tmpind}.c,3));' ...
                    'clear tmpdat tmpind;' ];
    cb_changeradii   = [  'tmpdat = get(gcbf, ''userdata'');' ...
                          'tmpdat.vol.r = str2num(get(gcbo, ''string''));' ...
                          'set(gcf, ''userdata'', tmpdat)' ];
    cb_changeconduct = [  'tmpdat = get(gcbf, ''userdata'');' ...
                          'tmpdat.vol.c = str2num(get(gcbo, ''string''));' ...
                          'set(gcf, ''userdata'', tmpdat)' ];
    cb_changeorigin  = [  'tmpdat = get(gcbf, ''userdata'');' ...
                          'tmpdat.vol.o = str2num(get(gcbo, ''string''));' ...
                          'set(gcf, ''userdata'', tmpdat)' ];
    % cb_fitelec = [ 'if get(gcbo, ''value''),' ...
    %                '  set(findobj(gcbf, ''tag'', ''origin''), ''enable'', ''off'');' ...
    %                'else' ...
    %                '  set(findobj(gcbf, ''tag'', ''origin''), ''enable'', ''on'');' ...
    %                'end;' ];
    valmodel     = 1;
    userdata    = [];
    if isfield(EEG.chaninfo, 'filename')
        if ~isempty(findstr(lower(EEG.chaninfo.filename), 'standard-10-5-cap385')), valmodel = 1; end;
        if ~isempty(findstr(lower(EEG.chaninfo.filename), 'standard_1005')),        valmodel = 2; end;
    end;

    geomvert = [3 1 1 1 1 1 1 1 1 1 1];

    geomhorz = {
        [1 2]
        [1]
        [1 1.3 0.5 0.5 ]
        [1 1.3 0.9 0.1 ]
        [1 1.3 0.5 0.5 ]
        [1 1.3 0.5 0.5 ]
        [1 1.3 0.5 0.5 ]
        [1 1.3 0.5 0.5 ]
        [1]
        [1]
        [1] };

    % define each individual graphical user element
    comhelp1 = [ 'warndlg2(strvcat(''The two default head models are in ''standard_BEM'' and ''standard_BESA'''',' ...
                 ''' sub-folders in the DIPFIT2 plugin folder, and may be modified there.''), ''Model type'');' ];
    comhelp3 = [ 'warndlg2(strvcat(''Any MR image normalized to the MNI brain model may be used for plotting'',' ...
                 '''(see the DIPFIT 2.0 tutorial for more information)''), ''Model type'');' ];
    comhelp2 = [ 'warndlg2(strvcat(''The template location file associated with the head model'',' ...
                 '''you are using must be entered (see tutorial).''), ''Template location file'');' ];
    commandload1 = [ '[filename, filepath] = uigetfile(''*'', ''Select a text file'');' ...
                    'if filename ~=0,' ...
                    '   set(findobj(''parent'', gcbf, ''tag'', ''model''), ''string'', [ filepath filename ]);' ...
                    'end;' ...
                    'clear filename filepath tagtest;' ];
    commandload2 = [ '[filename, filepath] = uigetfile(''*'', ''Select a text file'');' ...
                    'if filename ~=0,' ...
                    '   set(findobj(''parent'', gcbf, ''tag'', ''meg''), ''string'', [ filepath filename ]);' ...
                    'end;' ...
                    'clear filename filepath tagtest;' ];
    commandload3 = [ '[filename, filepath] = uigetfile(''*'', ''Select a text file'');' ...
                    'if filename ~=0,' ...
                    '   set(findobj(''parent'', gcbf, ''tag'', ''mri''), ''string'', [ filepath filename ]);' ...
                    'end;' ...
                    'clear filename filepath tagtest;' ];
    cb_selectcoreg = [ 'tmpmodel = get( findobj(gcbf, ''tag'', ''model''), ''string'');' ...
                       'tmploc2  = get( findobj(gcbf, ''tag'', ''meg'')  , ''string'');' ...
                       'tmploc1  = get( gcbo, ''userdata'');' ...
                       'tmptransf = get( findobj(gcbf, ''tag'', ''coregtext''), ''string'');' ...
                       '[tmp tmptransf] = coregister(tmploc1{1}, tmploc2, ''mesh'', tmpmodel,' ...
                       '                       ''transform'', str2num(tmptransf), ''chaninfo1'', tmploc1{2}, ''helpmsg'', ''on'');' ...
                       'if ~isempty(tmptransf), set( findobj(gcbf, ''tag'', ''coregtext''), ''string'', num2str(tmptransf)); end;' ...
                       'clear tmpmodel tmploc2 tmploc1 tmp tmptransf;' ];
    setmodel = [ 'pop_dipfit_settings(''setmodel'');' ];

    dipfitdefs; % contains template_model

    templatenames = { template_models.name };
    elements  = { ...
        { 'style' 'text'        'string'  [ 'Head model (click to select)' 10 '' ] } ...
        { 'style' 'listbox'     'string'  strvcat(templatenames{:}) ...
                                'callback' setmodel 'value' valmodel 'tag' 'listmodels' } { } ...
        { 'style' 'text'        'string' '________' 'tag' 'headstr' } ...
        { 'style' 'edit'        'string' '' 'tag'      'model' 'userdata' 'editable' 'enable' 'off'} ...
        { 'style' 'pushbutton'  'string' 'Browse'    'callback' commandload1       'userdata' 'editable' 'enable' 'off' } ...
        { 'style' 'pushbutton'  'string' 'Help'      'callback' comhelp1 } ...
        { 'style' 'text'        'string' 'Output coordinates' } ...
        { 'style' 'popupmenu'   'string' 'spherical (head radius 85 mm)|MNI|CTF' 'tag' 'coord' ...
          'value' 1  'userdata' 'editable' 'enable' 'off'} ...
        { 'style' 'text'        'string' 'Click to select' } { } ...
        { 'style' 'text'        'string' '________' 'tag' 'mristr' } ...
        { 'style' 'edit'        'string' '' 'tag' 'mri' } ...
        { 'style' 'pushbutton'  'string' 'Browse'       'callback' commandload3 } ...
        { 'style' 'pushbutton'  'string' 'Help'         'callback' comhelp3 } ...
        { 'style' 'text'        'string' '________', 'tag', 'chanstr' } ...
        { 'style' 'edit'        'string' '' 'tag' 'meg'  'userdata' 'editable' 'enable' 'off'} ...
        { 'style' 'pushbutton'  'string' 'Browse'       'callback' commandload2  'userdata' 'editable' 'enable' 'off'} ...
        { 'style' 'pushbutton'  'string' 'Help'         'callback' comhelp2 } ...
        { 'style' 'text'        'string' 'Co-register chan. locs. with head model' 'userdata' 'coreg' } ...
        { 'style' 'edit'        'string' '' 'tag' 'coregtext' 'userdata' 'coreg' } ...
        { 'style' 'pushbutton'  'string' 'Manual Co-Reg.' 'tag' 'manualcoreg' 'callback' cb_selectcoreg 'userdata' { EEG.chanlocs EEG.chaninfo } } ...
        { 'style' 'checkbox'    'string' 'No Co-Reg.'    'tag' 'coregcheckbox' 'value' 0  'userdata' 'coreg' } ...
        { 'style' 'text'        'string' 'Channels to omit from dipole fitting' } ...
        { 'style' 'edit'        'string' ''             'tag' 'elec' } ...
        { 'style' 'pushbutton'  'string' 'List' 'callback' cb_selectelectrodes } { } ...
        { } ...
        { 'style' 'text'        'string' 'Note: For EEG, check that the channel locations are on the surface of the head model' } ...
        { 'style' 'text'        'string' '(To do this: ''Set head radius'' to about 85 in the channel editor).' } ...
                };

    % plot GUI and protect parameters
    % -------------------------------
    userdata.template_models  = template_models;
    if isfield(EEG.chaninfo, 'filename')
         userdata.chanfile         = lower(EEG.chaninfo.filename);
    else userdata.chanfile         = '';
    end;
    optiongui = { 'geometry', geomhorz, 'uilist', elements, 'helpcom', 'pophelp(''pop_dipfit_settings'')', ...
                  'title', 'Dipole fit settings - pop_dipfit_settings()', ...
                  'userdata', userdata, 'geomvert', geomvert 'eval' 'pop_dipfit_settings(''setmodel'');' };
	[result, userdat2, strhalt, outstruct] = inputgui( 'mode', 'noclose', optiongui{:});
    if isempty(result), return; end;
    if ~isempty(get(0, 'currentfigure')) currentfig = gcf; else return; end;

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

    % decode GUI inputs
    % -----------------
    options = {};
    options = { options{:} 'hdmfile'      result{2} };
    options = { options{:} 'coordformat'  fastif(result{3} == 2, 'MNI', fastif(result{3} == 1, 'Spherical', 'CTF')) };
    options = { options{:} 'mrifile'      result{4} };
    options = { options{:} 'chanfile'     result{5} };
    if ~result{7}, options = { options{:} 'coord_transform' str2num(result{6}) }; end;
    options = { options{:} 'chansel'      setdiff(1:EEG.nbchan, str2num(result{8})) };

else
    options = varargin;
end

g = finputcheck(options, { 'hdmfile'  'string'    []         '';
                                 'mrifile'  'string'    []         '';
                                 'chanfile' 'string'    []         '';
                                 'chansel'  'integer'   []         [1:EEG.nbchan];
                                 'electrodes' 'integer'   []         [];
                                 'coord_transform' 'real' []         [];
                                 'coordformat' 'string'    { 'MNI' 'spherical' 'CTF' } 'MNI' });
if isstr(g), error(g); end;

OUTEEG = rmfield(OUTEEG, 'dipfit');
OUTEEG.dipfit.hdmfile     = g.hdmfile;
OUTEEG.dipfit.mrifile     = g.mrifile;
OUTEEG.dipfit.chanfile    = g.chanfile;
OUTEEG.dipfit.chansel     = g.chansel;
OUTEEG.dipfit.coordformat     = g.coordformat;
OUTEEG.dipfit.coord_transform = g.coord_transform;
if ~isempty(g.electrodes), OUTEEG.dipfit.chansel = g.electrodes; end;

% removing channels with no coordinates
% -------------------------------------
[tmpeloc labels Th Rd indices] = readlocs(EEG.chanlocs);
if length(indices) < length(EEG.chanlocs)
    disp('Warning: Channels removed from dipole fitting no longer have location coordinates!');
    OUTEEG.dipfit.chansel = intersect( OUTEEG.dipfit.chansel, indices);
end;

% checking electrode configuration
% --------------------------------
if 0
    disp('Checking the electrode configuration');
    tmpchan          = readlocs(OUTEEG.dipfit.chanfile);
    [tmp1 ind1 ind2] = intersect( lower({ tmpchan.labels }), lower({ OUTEEG.chanlocs.labels }));
    if isempty(tmp1)
        disp('No channel labels in common found between template and dataset channels');
        if ~isempty(findstr(OUTEEG.dipfit.hdmfile, 'BESA'))
            disp('Use the channel editor to fit a head sphere to your channel locations.');
            disp('Check for inconsistency in dipole info.');
        else
            disp('Results using standard BEM model are INACCURATE when the chan locations are not on the head surface!');
        end;
    else % common channels: performing best transformation
        TMP = OUTEEG;
        elec1 = eeglab2fieldtrip(TMP, 'elec');
        elec1 = elec1.elec;
        TMP.chanlocs = tmpchan;
        elec2 = eeglab2fieldtrip(TMP, 'elec');
        elec2 = elec2.elec;
        cfg.elec     = elec1;
        cfg.template = elec2;
        cfg.method   = 'warp';
        elec3 = electrodenormalize(cfg);

        % convert back to EEGLAB format
        OUTEEG.chanlocs = struct( 'labels', elec3.label, ...
                              'X'     , mat2cell(elec3.pnt(:,1)'), ...
                              'Y'     , mat2cell(elec3.pnt(:,2)'), ...
                              'Z'     , mat2cell(elec3.pnt(:,3)') );
        OUTEEG.chanlocs = convertlocs(OUTEEG.chanlocs, 'cart2all');
    end;

end;

com = sprintf('%s = pop_dipfit_settings( %s, %s);', inputname(1), inputname(1), vararg2str(options));

% test for wrong parameters
% -------------------------
function bool = test_wrong_parameters(hdl)

    coreg1 = get( findobj( hdl, 'tag', 'coregtext')    , 'string' );
    coreg2 = get( findobj( hdl, 'tag', 'coregcheckbox'), 'value' );
    meg    = get( findobj( hdl, 'tag', 'coord'), 'value' );

    bool = 0;
    if meg == 3, return; end;
    if coreg2 == 0 & isempty(coreg1)
         bool = 1; warndlg2(strvcat('You must co-register your channel locations', ...
                                    'with the head model (Press buttun, "Manual Co-Reg".', ...
                                    'and follow instructions); To bypass co-registration,', ...
                                    'check the checkbox " No Co-Reg".'), 'Error');
    end;

author: PatternRecognition

repo: OpenBMI


text:
function parsed_all = detect_all_by_featmap( featmap, scale, param, caffe_solver, thres )
    scale = 2.^(scale(end:-1:1)-5);
    parsed = [];
    for s = 1:numel(featmap)
        caffe_solver.reshape_as_input({{featmap{s}}});
        caffe_solver.set_input_data({{featmap{s}}});
        caffe_solver.forward_prefilled();
        parsed_t = convert_output_to_struct( param, caffe_solver, 1, thres);
        for j = 1 : length(parsed_t)
            parsed_t(j).box = round(parsed_t(j).box / scale(s));
            parsed_t(j).point = round(parsed_t(j).point / scale(s));
            parsed_t(j).box(:, end+1) = parsed_t(j).cls_score;
        end
        if s == 1
            parsed = parsed_t;
        else
            for j = 1 : length(param.gpu_id)
                parsed(j).cls_score = cat(1, parsed(j).cls_score, parsed_t(j).cls_score);
                %parsed(j).active = cat(1, parsed(j).cls_score, parsed_t.cls_score);
                parsed(j).box = cat(1, parsed(j).box, parsed_t(j).box);
                parsed(j).point = cat(1, parsed(j).point, parsed_t(j).point);
            end
        end
    end
    parsed_all = parsed;
end

author: liuyuisanai

repo: RSA-for-object-detection


text:
function [grandavg] = ft_timelockgrandaverage(cfg, varargin)

% FT_TIMELOCKGRANDAVERAGE computes ERF/ERP average and variance
% over multiple subjects or over blocks within one subject
%
% Use as
%   [grandavg] = ft_timelockgrandaverage(cfg, avg1, avg2, avg3, ...)
%
% where
%   avg1..N are the ERF/ERP averages as obtained from FT_TIMELOCKANALYSIS
%
% and cfg is a configuration structure with
%   cfg.method         = string, 'across' or 'within' (default = 'across'), see below for details
%   cfg.parameter      = string, which parameter to average (default = 'avg')
%   cfg.channel        = Nx1 cell-array with selection of channels (default = 'all'), see FT_CHANNELSELECTION for details
%   cfg.latency        = [begin end] in seconds or 'all' (default = 'all')
%   cfg.keepindividual = string, 'yes' or 'no' (default = 'no')
%   cfg.nanmean        = string, can be 'yes' or 'no' (default = 'yes')
%   cfg.normalizevar   = string, 'N' or 'N-1' (default = 'N-1')
%
% If cfg.method = 'across', a plain average is performed, i.e. the requested
% parameter in each input argument is weighted equally in the average. This is useful
% when averaging across subjects. The variance-field will contain the variance across
% the parameter of interest, and the output dof-field will contain the number of
% input arguments.
%
% If cfg.method = 'within', a weighted average is performed, i.e. the requested
% parameter in each input argument is weighted according to the degrees of freedom in
% the dof-field. This is useful when averaging within subjects across blocks, e.g.
% when each block was recorded in a separate file. The variance-field will contain
% the variance across all input observations, and the output dof-field will contain
% the total number of observations.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure. For this particular function, the input should be
% structured as a cell-array.
%
% See also FT_TIMELOCKANALYSIS, FT_TIMELOCKSTATISTICS, FT_TIMELOCKBASELINE

% Copyright (C) 2003-2006, Jens Schwarzbach
% Copyright (C) 2013, Burkhard Maess
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

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

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

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

% return immediately after distributed execution
if ~isempty(ft_getopt(cfg, 'distribute'))
  return
end

% check if the input data is valid for this function
for i=1:length(varargin)
  if isfield(varargin{i},'trial') && isfield(varargin{i},'avg') % see bug2372 (dieloz)
    varargin{i} = rmfield(varargin{i}, 'trial');
    varargin{i}.dimord = 'chan_time';
    ft_warning('not using the trials, using the single-subject average to compute the grand average');
  else
    if isfield(varargin{i},'trial') && ~isfield(varargin{i},'avg')
      ft_error('input structure %d does not contain an average, use FT_TIMELOCKANALYSIS first', i);
    end
  end
  varargin{i} = ft_checkdata(varargin{i}, 'datatype', 'timelock', 'feedback', 'no');
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'channels'}); % prevent accidental typos, see issue 1729

% set the defaults
cfg.method         = ft_getopt(cfg, 'method'        , 'across');
cfg.parameter      = ft_getopt(cfg, 'parameter'     , 'avg');
cfg.channel        = ft_getopt(cfg, 'channel'       , 'all');
cfg.latency        = ft_getopt(cfg, 'latency'       , 'all');
cfg.keepindividual = ft_getopt(cfg, 'keepindividual', 'no');
cfg.nanmean        = ft_getopt(cfg, 'nanmean'       , 'yes');
cfg.normalizevar   = ft_getopt(cfg, 'normalizevar'  , 'N-1');

if iscell(cfg.parameter)
  if numel(cfg.parameter)>1
    ft_error('only a single parameter can be specified to be averaged');
  else
    cfg.parameter = cfg.parameter{1};
  end
end

if strcmp(cfg.parameter, 'trial')
  ft_error('not supporting averaging over the repetition dimension, please use FT_TIMELOCKANALYSIS');
end

% select trials and channels of interest
orgcfg = cfg;
tmpcfg = keepfields(cfg, {'parameter', 'channel', 'tolerance', 'latency', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
[varargin{:}] = ft_selectdata(tmpcfg, varargin{:});
% restore the provenance information
[cfg, varargin{:}] = rollback_provenance(cfg, varargin{:});
% do not use the default option returned by FT_SELECTDATA, but the original one for this function
cfg.nanmean = orgcfg.nanmean;

% determine the size of the data to be averaged
datsiz = size(varargin{1}.(cfg.parameter));
dimord = getdimord(varargin{1}, cfg.parameter);
nsubj  = length(varargin);

% whether to normalize the variance with N or N-1, see VAR
normalizewithN = strcmpi(cfg.normalizevar, 'N');

% whether nans should persist in the output or be treated as missing values
if istrue(cfg.nanmean)
  mymean = @nanmean;
  myvar  = @nanvar;
  mysum  = @nansum;
else
  mymean = @mean;
  myvar  = @var;
  mysum  = @sum;
end

if istrue(cfg.keepindividual)
  fprintf('not computing average, but keeping individual %s for %d subjects\n', cfg.parameter, nsubj);

  % allocate memory to hold the data and collect it
  dat = zeros([nsubj, datsiz]);
  for s=1:nsubj
    dat(s, :, :, :) = varargin{s}.(cfg.parameter);
  end
  grandavg.individual = dat; % Nsubj x Nchan x Nsamples

else
  dat = nan([nsubj, datsiz]);
  dof = nan([nsubj, datsiz]);

  switch cfg.method
    case 'across'
      fprintf('computing average of %s across %d subjects\n', cfg.parameter, nsubj);
      for s=1:nsubj
        dat(s, :, :, :) = varargin{s}.(cfg.parameter);
      end

      % compute the mean and variance across subjects
      grandavg.avg = reshape(mymean(dat, 1),                 datsiz);
      grandavg.var = reshape(myvar(dat, normalizewithN, 1),  datsiz);
      grandavg.dof = reshape(sum(isfinite(dat), 1),          datsiz);

      if normalizewithN
        % just to be sure
        grandavg.var(grandavg.dof<=0) = NaN;
      else
        % see https://stats.stackexchange.com/questions/4068/how-should-one-define-the-sample-variance-for-scalar-input
        % the fieldtrip/external/stats/nanvar implementation behaves differently here than Mathworks VAR and NANVAR implementations
        grandavg.var(grandavg.dof<=1) = NaN;
      end

    case 'within'
      fprintf('computing average of %s within subjects over %d blocks\n', cfg.parameter, nsubj);

      for s=1:nsubj
        dat(s, :, :, :) = varargin{s}.(cfg.parameter);
        dof(s, :, :, :) = varargin{s}.dof;
      end % for nsub

      % compute the weighted mean across input arguments
      grandavg.avg = mysum(dat .* dof, 1) ./ sum(dof, 1);
      grandavg.avg = reshape(grandavg.avg, datsiz);
      grandavg.dof = sum(dof, 1);
      grandavg.dof = reshape(grandavg.dof, datsiz);

      % this follows https://en.wikipedia.org/wiki/Weighted_arithmetic_mean#Weighted_sample_variance with frequency weights
      if normalizewithN
        grandavg.var = mysum(dof .* (dat - repmat(reshape(grandavg.avg, [1 datsiz]), [nsubj 1 1])).^2, 1) ./  sum(dof, 1);
        grandavg.var = reshape(grandavg.var, datsiz);
      else
        grandavg.var = mysum(dof .* (dat - repmat(reshape(grandavg.avg, [1 datsiz]), [nsubj 1 1])).^2, 1) ./ (sum(dof, 1) - 1);
        grandavg.var = reshape(grandavg.var, datsiz);
      end

    otherwise
      ft_error('unsupported value for cfg.method')

  end % switch
end % if keepindividual

% collect the output data
grandavg = copyfields(varargin{1}, grandavg, {'time', 'freq', 'label', 'labelcmb'});

% sensor positions should only be present in the output when they are identical in all inputs
hasgrad = cellfun(@(x) isfield(x, 'grad'), varargin(:));
if all(hasgrad) % check if positions are different between subjects
  samegrad = cellfun(@(x) isequal(varargin{1}.grad, x.grad), varargin(2:end));
  if all(samegrad)
    grandavg.grad = varargin{1}.grad;
  else
    ft_warning('discarding gradiometer information because it is not identical in all inputs');
  end
end

haselec = cellfun(@(x) isfield(x, 'elec'), varargin(:));
if all(haselec) % check if positions are different between subjects
  sameelec = cellfun(@(x) isequal(varargin{1}.elec, x.elec), varargin(2:end));
  if all(sameelec)
    grandavg.elec = varargin{1}.elec;
  else
    ft_warning('discarding electrode information because it is not identical in all inputs');
  end
end

hasopto = cellfun(@(x) isfield(x, 'opto'), varargin(:));
if all(hasopto) % check if positions are different between subjects
  sameopto = cellfun(@(x) isequal(varargin{1}.opto, x.opto), varargin(2:end));
  if all(sameopto)
    grandavg.opto = varargin{1}.opto;
  else
    ft_warning('discarding optode information because it is not identical in all inputs');
  end
end

% set dimord
if strcmp(cfg.keepindividual, 'yes')
  grandavg.dimord = ['subj_', dimord];
elseif strcmp(cfg.keepindividual, 'no')
  grandavg.dimord =  dimord;
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   varargin
ft_postamble provenance grandavg
ft_postamble history    grandavg
ft_postamble savevar    grandavg

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef chebpref
%CHEBPREF   Abstract class for Chebfun preferences.
%
%   CHEBPREF is an abstract class for managing preferences of the various
%   subsystems comprising Chebfun.  On its own, it provides little more than a
%   wrapper around a MATLAB struct.  Concrete preference classes inherit from
%   CHEBPREF and add to its functionality as desired.
%
% See also CHEBFUNPREF, CHEBOPPREF.

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

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS PROPERTIES:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    properties ( Access = protected )
        % MATLAB struct to hold a list of preferences for a given subsystem.
        prefList
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Non-Static methods implemented by this class.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

        function out = subsref(pref, ind)
        %SUBSREF   Subscripted referencing for CHEBPREF.
        %   P.PROP, where P is a CHEBPREF object, returns the value of the
        %   CHEBPREF property PROP stored in P.  If PROP is not a CHEBPREF
        %   property, an error will be thrown.
        %
        %   CHEBPREF does not support any other subscripted referencing
        %   types, including '()' and '{}'.

            switch ( ind(1).type )
                case '.'
                    if ( isfield(pref.prefList, ind(1).subs) )
                        out = pref.prefList.(ind(1).subs);
                    end

                    if ( numel(ind) > 1 )
                        out = subsref(out, ind(2:end));
                    end
                otherwise
                    error('CHEBFUN:CHEBPREF:chebpref:badType', ...
                        'Invalid subscripted reference type.')
            end
        end

        function pref = subsasgn(pref, ind, val)
        %SUBSASGN   Subscripted assignment for CHEBPREF.
        %   P.PROP = VAL, where P is a CHEBPREF object, assigns the value
        %   VAL to the CHEBPREF property PROP stored in P.  If PROP is not a
        %   CHEBPREF property, an error will be thrown.
        %
        %   CHEBPREF does not support any other subscripted assignment types,
        %   including '()' and '{}'.

            switch ( ind(1).type )
                case '.'
                    if ( isfield(pref.prefList, ind(1).subs) )
                        pref.prefList = builtin('subsasgn', pref.prefList, ...
                            ind, val);
                    end
                otherwise
                    error('CHEBFUN:CHEBPREF:chebpref:badType', ...
                        'Invalid subscripted assignment type.')
            end
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Abstract methods
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Abstract = true )
        % Display information about a CHEBPREF object.
        display(pref)
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Static methods
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = true )

        function pref1 = mergePrefStructs(pref1, pref2)
        %MERGEPREFSTRUCTS   Merge preference structures.
        %   P = CHEBPREF.MERGEPREFSTRUCTS(P, Q), where P and Q are MATLAB
        %   structures, "merges" Q into P by replacing the contents of fields
        %   in P with those of identically-named fields in Q.  If Q has a field
        %   whose name does not match any of those in P, it is added to P.

            for field = fieldnames(pref2).'
                pref1.(field{1}) = pref2.(field{1});
            end
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Hidden static methods.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = true, Hidden = true  )

        function setDefaults(makeSC, manageSCDefaults, varargin)
        %SETDEFAULTS   Set default preferences.
        %   CHEBPREF.SETDEFAULTS(MAKESC, MANAGESCDEFAULTS, ...) sets the
        %   persistently stored defaults of a subclass of CHEBPREF.  MAKESC is
        %   a handle to the subclass constructor, and MANAGESCDEFAULTS is a
        %   handle to the subclass static method responsible for managing the
        %   stored defaults.
        %
        %   Users should not call CHEBPREF.SETDEFAULTS.  Call the SETDEFAULTS
        %   method of the appropriate Chebfun subsystem whose defaults need to
        %   be altered instead.

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % DEVELOPER NOTES:
        %  - This is a template method to prevent subclasses from needing to
        %    duplicate the logic required to set persistently-stored defaults
        %    As the help text says, this method is not user-facing.
        %  - This method is not likely to work if the subclass does not handle
        %    things in the exact same way that the existing CHEBFUNPREF and
        %    CHEBOPPREF subclasses do.
        %  - The function referred to by the MANAGESCDEFAULTS input needs to
        %    conform to a specific prototype.  See the documentation for
        %    CHEBFUNPREF.MANAGEDEFAULTPREFS for an example.
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            nargs = length(varargin);

            if ( nargs < 1)
                error('CHEBFUN:CHEBPREF:chebpref:notEnoughArguments', ...
                    'Not enough arguments.');
            end

            if ( nargs == 1 )
                if ( isstruct(varargin{1}) )
                    varargin{1} = makeSC(varargin{1});
                end

                if ( ischar(varargin{1}) && strcmp(varargin{1}, 'factory') )
                    manageSCDefaults('set-factory');
                elseif ( isa(varargin{1}, 'chebpref') )
                    manageSCDefaults('set', varargin{1}.prefList);
                else
                    error('CHEBFUN:CHEBPREF:chebpref:badArg', ...
                        ['When calling chebpref.setDefaults() with just ' ...
                         'one argument, that argument must be ''factory'', ' ...
                         'a subclass of CHEBPREF, or a MATLAB structure.']);
                end
            elseif ( mod(nargs, 2) == 0 )
                manageSCDefaults('set', varargin{:});
            else
                error('CHEBFUN:CHEBPREF:chebpref:unpairedArg', ...
                    'Unpaired argument in name-value pair list.');
            end
        end

    end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function D = sparseDiag(d)

% SPARSEDIAG Create a diagonal matrix that is sparse from a vector.
% FORMAT
% DESC creates a diagonal matrix that is sparse from a vector.
% ARG d : the diagonal vector from which the sparse diagonal matrix
% is formed.
% RETURN D : the sparse diagonal matrix containing the vector as
% its diagonal.
%
% SEEALSO : diag, spdiags
%
% COPYRIGHT : Neil D. Lawrence, 2005

% NDLUTIL

if length(size(d)) ~=2
  error('Input must be a vector.');
end
if size(d, 1) ~= 1 & size(d, 2) ~=1
  error('Input must be a vector.');
end

D = spdiags(d, 0, length(d), length(d));
% % Can be made more efficient.
% n = length(d);
% D = spalloc(n, n, n);
% for i = 1:n
%   D(i, i) = d(i);
% end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function h = p41_h ( n, x )

%*****************************************************************************80
%
%% P41_H evaluates the Hessian for problem 41.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real H(N,N), the N by N Hessian matrix.
%
  h = zeros ( n, n );

  h(1,1) = 1200.0 * x(1)^2 - 400.0 * x(2) + 2.0;
  h(1,2) = - 400.0 * x(1);

  h(2,1) = -400.0 * x(1);
  h(2,2) = 220.2;
  h(2,4) = 19.8;

  h(3,3) = -360.0 * x(4) + 1080.0 * x(3)^2 + 2.0;
  h(3,4) = - 360.0 * x(3);

  h(4,2) = 19.8;
  h(4,3) = - 360.0 * x(3);
  h(4,4) = 200.2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION IS RELEASED "AS IS."  THE U.S. GOVERNMENT MAKES NO WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, CONCERNING THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION, INCLUDING, WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT WILL THE U.S. GOVERNMENT BE LIABLE FOR ANY DAMAGES, INCLUDING ANY LOST PROFITS, LOST SAVINGS OR OTHER INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE, OR INABILITY TO USE, THIS SOFTWARE OR ANY ACCOMPANYING DOCUMENTATION, EVEN IF INFORMED IN ADVANCE OF THE POSSIBILITY OF SUCH DAMAGES.
%
% change log:
%   040204 tdr stub file
%
% file: get_ci_ribp_1s.m
%
function varargout = get_ci_ribp_1s(varargin)
  error('get_ci_ribp_1s is not part of current release.');
return;

author: Sable

repo: mcbench-benchmarks


text:
%kmreplicate 'Replicate a Given Matrix as Submatrices of a Larger Matrix'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros mreplicate.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'Input mmatrix data object'
% Integer: r 'Replica count along rows ', default: 1: 'replica count along rows'
% Integer: c 'Replica count along cols ', default: 1: 'replica count along columns'
% OutputFile: o 'Output', required: 'Resulting matrix output data object'
%
% Example: o = kmreplicate(i, {'i','';'r',1;'c',1;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% mreplicate - Replicate a Given Matrix as Submatrices of a Larger Matrix
%
%  DESCRIPTION
% .I mreplicate
% replicates the input matrix as a rectangular array of submatrices to form
% a larger matrix. The input maybe of any dimension NxM as long as N>1 and M>1.
% The output matrix will be (R*N)x(C*M) where R is the number of times the
% input matrix is to be replicated along the row dimension and C is the number of
% times the input matrix is to be replicated along the column dimension.
%
% The output data type is the same as the input data type.
%
%
%
%  EXAMPLES
%
% % mreplicate -i image:moon -o blah -r 2 -c 3
%
% Produces an image of the moon replicated 3 times in the horizontal direction
% and 2 times in the vertical direction.
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kmreplicate(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kmreplicate(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'r', 1;'c', 1;'o', '__output'};
maxval={0,2,2,0};
minval={0,2,2,0};
istoggle=[0,1,1,0];
was_set=istoggle * 0;
paramtype={'InputFile','Integer','Integer','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'mreplicate"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% [INPUT]
% ds = A structure representing the dataset.
% sn = A string representing the serial number of the result file.
% temp = A string representing the full path to the Excel spreadsheet used as template for the result file.
% out = A string representing the full path to the Excel spreadsheet to which the results are written, eventually replacing the previous ones.
% cvm = A string representing the method for calculating the critical values (optional, default='WB'):
%   - 'FS' for finite sample critical values;
%   - 'WB' for wild bootstrap.
% cvq = A float [0.90,0.99] representing the quantile of the critical values (optional, default=0.95).
% lag_max = An integer [0,10] representing the maximum lag order to be evaluated for the Augmented Dickey-Fuller test (optional, default=0).
% lag_sel = A string representing the lag order selection criteria for the Augmented Dickey-Fuller test (optional, default='FIX'):
%   - 'AIC' for Akaike's Information Criterion;
%   - 'BIC' for Bayesian Information Criterion;
%   - 'FIX' to use a fixed lag order;
%   - 'FPE' for Final Prediction Error;
%   - 'HQIC' for Hannan-Quinn Information Criterion.
% mbd = An integer [3,Inf) representing the minimum duration of a bubble in days (optional, default=NaN).
%   If NaN is provided, then an optimal value based on the number of observations is used.
% analyze = A boolean that indicates whether to analyse the results and display plots (optional, default=false).
%
% [OUTPUT]
% result = A structure representing the original dataset inclusive of intermediate and final calculations.
% stopped = A boolean that indicates whether the process has been stopped through user input.

function [result,stopped] = run_bubbles_detection(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('ds',@(x)validateattributes(x,{'struct'},{'nonempty'}));
        ip.addRequired('sn',@(x)validateattributes(x,{'char'},{'nonempty' 'size' [1 NaN]}));
        ip.addRequired('temp',@(x)validateattributes(x,{'char'},{'nonempty' 'size' [1 NaN]}));
        ip.addRequired('out',@(x)validateattributes(x,{'char'},{'nonempty' 'size' [1 NaN]}));
        ip.addOptional('cvm','WB',@(x)any(validatestring(x,{'FS' 'WB'})));
        ip.addOptional('cvq',0.95,@(x)validateattributes(x,{'double'},{'real' 'finite' '>=' 0.90 '<=' 0.99 'scalar'}));
        ip.addOptional('lag_max',0,@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer' '>=' 0 '<=' 10 'scalar'}));
        ip.addOptional('lag_sel','FIX',@(x)any(validatestring(x,{'AIC' 'BIC' 'FIX' 'FPE' 'HQIC'})));
        ip.addOptional('mbd',NaN,@(x)validateattributes(x,{'double'},{'real' 'scalar'}));
        ip.addOptional('analyze',false,@(x)validateattributes(x,{'logical'},{'scalar'}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    ds = validate_dataset(ipr.ds,'BubblesDetection');
    sn = ipr.sn;
    temp = validate_template(ipr.temp);
    out = validate_output(ipr.out);
    cvm = ipr.cvm;
    cvq = ipr.cvq;
    lag_max = ipr.lag_max;
    lag_sel = ipr.lag_sel;
    mbd = validate_mbd(ipr.mbd,ds.T);
    analyze = ipr.analyze;

    nargoutchk(1,2);

    [result,stopped] = run_bubbles_detection_internal(ds,sn,temp,out,cvm,cvq,lag_max,lag_sel,mbd,analyze);

end

function [result,stopped] = run_bubbles_detection_internal(ds,sn,temp,out,cvm,cvq,lag_max,lag_sel,mbd,analyze)

    result = [];
    stopped = false;
    e = [];

    ds = initialize(ds,sn,cvm,cvq,lag_max,lag_sel,mbd);
    n = ds.N;
    t = ds.T;

    bar = waitbar(0,'Initializing bubbles detection measures...','CreateCancelBtn',@(src,event)setappdata(gcbf(),'Stop', true));
    setappdata(bar,'Stop',false);
    cleanup = onCleanup(@()delete(bar));

    pause(1);
    waitbar(0,bar,'Calculating bubbles detection measures...');
    pause(1);

    try

        p = ds.Prices;

        for i = 1:n
            waitbar((i - 1) / n,bar,['Calculating bubbles detection measures for ' ds.FirmNames{i} '...']);

            if (getappdata(bar,'Stop'))
                stopped = true;
                break;
            end

            offset = min(ds.Defaults(i) - 1,t);
            p_i = p(1:offset,i);

            [bsadfs,cvs,detection,breakdown] = psy_bubbles_detection(p_i,ds.CVM,ds.CVQ,ds.LagMax,ds.LagSel,ds.MBD);
            ds.BSADFS(1:offset,i) = bsadfs;
            ds.CVS(1:offset,i) = cvs;
            ds.BUB(1:offset,i) = detection(:,1);
            ds.BMPH(1:offset,i) = detection(:,2);
            ds.BRPH(1:offset,i) = detection(:,3);
            ds.Breakdowns{i} = breakdown;

            if (getappdata(bar,'Stop'))
                stopped = true;
                break;
            end

            waitbar(i / n,bar);
        end

    catch e
    end

    if (~isempty(e))
        delete(bar);
        rethrow(e);
    end

    if (stopped)
        delete(bar);
        return;
    end

    pause(1);
    waitbar(1,bar,'Finalizing cross-sectional measures...');
    pause(1);

    try
        ds = finalize(ds);
    catch e
        delete(bar);
        rethrow(e);
    end

    pause(1);
    waitbar(1,bar,'Writing bubbles detection measures...');
    pause(1);

    try
        write_results(ds,temp,out);
        delete(bar);
    catch e
        delete(bar);
        rethrow(e);
    end

    if (analyze)
        analyze_result(ds);
    end

    result = ds;

end

%% PROCESS

function ds = initialize(ds,sn,cvm,cvq,lag_max,lag_sel,mbd)

    n = ds.N;
    t = ds.T;

    ds.Result = 'BubblesDetection';
    ds.ResultDate = now(); %#ok<TNOW1>
    ds.ResultAnalysis = @(ds)analyze_result(ds);
    ds.ResultSerial = sn;

    ds.CVM = cvm;
    ds.CVQ = cvq;
    ds.LagMax = lag_max;
    ds.LagSel = lag_sel;
    ds.MBD = mbd;

    if (strcmp(ds.LagSel,'FIX'))
        label = [' (CVM= ' ds.CVM ', CVQ=' num2str(ds.CVQ * 100) '%, L=' num2str(ds.LagMax) ')'];
    else
        label = [' (CVM= ' ds.CVM ', CVQ=' num2str(ds.CVQ * 100) '%, LM=' num2str(ds.LagMax) ', LS=' ds.LagSel ')'];
    end

    ds.LabelsMeasuresSimple = {'BUB' 'BMPH' 'BRPH' 'Breakdowns'};
    ds.LabelsMeasures = [strcat(ds.LabelsMeasuresSimple(1:3),{label}) ds.LabelsMeasuresSimple(4)];

    ds.LabelsIndicatorsSimple = {'BC' 'BCP'};
    ds.LabelsIndicators = {['BC' label] ['BCP' label]};

    ds.LabelsSheetsSimple = [ds.LabelsMeasuresSimple {'Indicators'}];
    ds.LabelsSheets = [ds.LabelsMeasures {'Indicators'}];

    ds.BSADFS = NaN(t,n);
    ds.CVS = NaN(t,n);
    ds.BUB = NaN(t,n);
    ds.BMPH = NaN(t,n);
    ds.BRPH = NaN(t,n);
    ds.Breakdowns = cell(1,n);

    ds.Indicators = NaN(t,numel(ds.LabelsIndicators));

    ds.ComparisonReferences = {'Indicators' 1 strcat({'BD-'},ds.LabelsIndicatorsSimple{1})};

end

function ds = finalize(ds)

    bc = sum(ds.BUB .* ds.Capitalizations,2,'omitnan');
    bcp = bc ./ sum(ds.Capitalizations,2,'omitnan');

    ds.Indicators = [bc bcp];

end

function write_results(ds,temp,out)

    [out_path,~,~] = fileparts(out);

    try
        if (exist(out_path,'dir') ~= 7)
            mkdir(out_path);
        end

        if (exist(out,'file') == 2)
            delete(out);
        end
    catch
        error('A system I/O error occurred while writing the results.');
    end

    copy_result = copyfile(temp,out,'f');

    if (copy_result == 0)
        error('The output file could not be created from the template file.');
    end

    lc = numel(ds.LabelsSheetsSimple);

    dates_str = cell2table(ds.DatesStr,'VariableNames',{'Date'});

    for i = 1:(lc - 2)
        sheet = ds.LabelsSheetsSimple{i};
        measure = strrep(sheet,' ','');

        tab = [dates_str array2table(ds.(measure),'VariableNames',ds.FirmNames)];
        writetable(tab,out,'FileType','spreadsheet','Sheet',sheet,'WriteRowNames',true);
    end

    breakdowns_count = cellfun(@(x)size(x,1),ds.Breakdowns);
    breakdowns_firms = arrayfun(@(x)repelem(ds.FirmNames(x),breakdowns_count(x),1),1:ds.N,'UniformOutput',false);
    breakdowns_firms = vertcat(breakdowns_firms{:});
    breakdowns_firms = cell2table(breakdowns_firms,'VariableNames',{'Firm'});

    tab = [breakdowns_firms array2table(cell2mat(ds.Breakdowns.'),'VariableNames',{'Start' 'Peak' 'End' 'Duration', 'Boom Phase' 'Burst Phase'})];
    writetable(tab,out,'FileType','spreadsheet','Sheet',ds.LabelsSheetsSimple{lc - 1},'WriteRowNames',true);

    tab = [dates_str array2table(ds.Indicators,'VariableNames',ds.LabelsIndicatorsSimple)];
    writetable(tab,out,'FileType','spreadsheet','Sheet',ds.LabelsSheetsSimple{lc},'WriteRowNames',true);

    worksheets_batch(out,ds.LabelsSheetsSimple,ds.LabelsSheets);

end

%% PLOTTING

function analyze_result(ds)

    safe_plot(@(id)plot_sequence(ds,id));
    safe_plot(@(id)plot_indicators(ds,id));

end

function plot_indicators(ds,id)

    bc = ds.Indicators(:,1);
    bcp = ds.Indicators(:,2);

    y_limits_bc = plot_limits(bc,0.1,0);
    y_limits_bcp = [0 100];

    f = figure('Name','Bubbles Detection Measures > Indicators','Units','normalized','Position',[100 100 0.85 0.85],'Tag',id);

    sub_1 = subplot(2,1,1);
    plot(sub_1,ds.DatesNum,smooth_data(bc),'Color',[0.000 0.447 0.741]);
    set(sub_1,'YLim',y_limits_bc);
    t1 = title(sub_1,ds.LabelsIndicators{1});
    set(t1,'Units','normalized');
    t1_position = get(t1,'Position');
    set(t1,'Position',[0.4783 t1_position(2) t1_position(3)]);

    sub_2 = subplot(2,1,2);
    area(sub_2,ds.DatesNum,smooth_data(bcp),'EdgeColor',[0.000 0.447 0.741],'FaceAlpha',0.5,'FaceColor',[0.749 0.862 0.933]);
    set(sub_2,'YLim',y_limits_bcp);
    t2 = title(sub_2,ds.LabelsIndicators{2});
    set(t2,'Units','normalized');
    t2_position = get(t2,'Position');
    set(t2,'Position',[0.4783 t2_position(2) t2_position(3)]);

    set([sub_1 sub_2],'XLim',[ds.DatesNum(1) ds.DatesNum(end)],'XTickLabelRotation',45);
    set([sub_1 sub_2],'XGrid','on','YGrid','on');

    if (ds.MonthlyTicks)
        date_ticks([sub_1 sub_2],'x','mm/yyyy','KeepLimits','KeepTicks');
    else
        date_ticks([sub_1 sub_2],'x','yyyy','KeepLimits');
    end

    figure_title(f,'Indicators');

    maximize_figure(f);

end

function plot_sequence(ds,id)

    n = ds.N;
    t = ds.T;
    dn = ds.DatesNum;
    mt = ds.MonthlyTicks;

    bmph = ds.BMPH;
    brph = ds.BRPH;

    data_bsadfs = ds.BSADFS;
    data_cvs = ds.CVS;
    data_ps = smooth_data(ds.Prices);
    data_ps_bm = data_ps .* bmph;
    data_ps_br = data_ps .* brph;
    data_ps_ot = NaN(t,n);

    data_bds = cell(1,n);

    for i = 1:n
        offset = min(ds.Defaults(i) - 1,t);

        bsadfs_i = data_bsadfs(1:offset,i);
        idx = find(bsadfs_i ~= 0,1,'first');
        data_bsadfs(1:idx,i) = bsadfs_i(idx + 1);
        data_bsadfs(1:offset,i) = smooth_data(data_bsadfs(1:offset,i));

        cvs_i = data_cvs(1:offset,i);
        idx = find(isnan(cvs_i),1,'last');
        data_cvs(1:idx,i) = cvs_i(idx + 1);
        data_cvs(1:offset,i) = smooth_data(data_cvs(1:offset,i));

        data_ps_ot(1:offset,i) = data_ps(1:offset,i) .* ~(bmph(1:offset,i) | bmph(1:offset,i));

        bd_i = ds.Breakdowns{i};

        bm_count = sum(bd_i(:,5));
        br_count = sum(bd_i(:,6));
        ot_count = t - (bm_count + br_count);
        data_bds{i} = [ot_count bm_count br_count];
    end

    data_ones = ones(1,n);
    data = [repmat({dn},1,n); mat2cell(data_bsadfs,t,data_ones); mat2cell(data_cvs,t,data_ones); mat2cell(data_ps,t,data_ones); mat2cell(data_ps_bm,t,data_ones); mat2cell(data_ps_br,t,data_ones); mat2cell(data_ps_ot,t,data_ones); data_bds];

    plots_title = [repmat({'Bubbles'},1,n); repmat({'Model'},1,n); repmat({'Breakdown - Percentage'},1,n); repmat({'Breakdown - Count'},1,n)];

    x_limits = [dn(1) dn(end)];

    core = struct();

    core.N = n;
    core.Data = data;
    core.Function = @(subs,data)plot_function(subs,data);

    core.OuterTitle = 'Bubbles Detection Measures > Time Series';
    core.InnerTitle = 'Time Series';
    core.SequenceTitles = ds.FirmNames;

    core.PlotsAllocation = [7 4];
    core.PlotsSpan = {[1:3 5:7 9:11 13:15] [21:23 25:27] [4 8 12] [20 24 28]};
    core.PlotsTitle = plots_title;

    core.XDates = {mt mt [] []};
    core.XGrid = {true true false false};
    core.XLabel = {[] [] [] []};
    core.XLimits = {x_limits x_limits [] []};
    core.XRotation = {45 45 [] []};
    core.XTick = {[] [] [] []};
    core.XTickLabels = {[] [] [] []};

    core.YGrid = {true true false true};
    core.YLabel = {[] [] [] []};
    core.YLimits = {[] [] [] []};
    core.YRotation = {[] [] [] []};
    core.YTick = {[] [] [] []};
    core.YTickLabels = {[] [] [] []};

    sequential_plot(core,id);

    function plot_function(subs,data)

        x = data{1};
        bsadfs = data{2};
        cvs = data{3};
        ps = data{4};
        ps_bm = data{5};
        ps_br = data{6};
        ps_ot = data{7};
        bd = data{8};

        d = find(isnan(ps),1,'first');

        if (isempty(d))
            xd = [];
        else
            xd = x(d) - 1;
        end

        sub_1 = subs(1);
        a1 = area(sub_1,x,ps_ot,'EdgeColor','none','FaceAlpha',0.5,'FaceColor',[0.749 0.862 0.933]);
        hold(sub_1,'on');
            a2 = area(sub_1,x,ps_bm,'EdgeColor','none','FaceColor',[0.200 0.627, 0.173]);
            a3 = area(sub_1,x,ps_br,'EdgeColor','none','FaceColor',[0.984 0.502 0.447]);
            plot(sub_1,x,ps,'Color',[0.000 0.000 0.000])
        hold(sub_1,'off');

        if (~isempty(xd))
            hold(sub_1,'on');
                plot(sub_1,[xd xd],get(sub_1,'YLim'),'Color',[1.000 0.400 0.400]);
            hold(sub_1,'off');
        end

        l = legend(sub_1,[a1 a2 a3],'None','Boom Phase','Burst Phase','Location','south','Orientation','horizontal');
        set(l,'Units','normalized');
        l_position = get(l,'Position');
        set(l,'Position',[l_position(1) 0.3645 l_position(3) l_position(4)]);

        sub_2 = subs(2);
        p1 = plot(sub_2,x,bsadfs,'Color',[0.000 0.447 0.741]);
        hold(sub_2,'on');
            p2 = plot(sub_2,x,cvs,'Color',[0.000 0.000 0.000],'LineStyle','--');
        hold(sub_2,'off');

        if (~isempty(xd))
            hold(sub_2,'on');
                plot(sub_2,[xd xd],get(sub_2,'YLim'),'Color',[1.000 0.400 0.400]);
            hold(sub_2,'off');
        end

        legend(sub_2,[p1 p2],'BSADF','Critical Values','Location','northwest');

        sub_3 = subs(3);
        pc1 = pie(sub_3,bd);
        hpc1 = findobj(pc1,'Type','patch');
        set(hpc1(1),'FaceAlpha',0.5,'FaceColor',[0.749 0.862 0.933]);
        set(hpc1(2),'FaceColor',[0.200 0.627, 0.173]);
        set(hpc1(3),'FaceColor',[0.984 0.502 0.447]);

        sub_4 = subs(4);
        b1 = bar(sub_4,1,bd(1));
        set(b1,'FaceAlpha',0.5,'FaceColor',[0.749 0.862 0.933]);
        hold(sub_4,'on');
            b2 = bar(sub_4,2,bd(2));
            set(b2,'FaceColor',[0.200 0.627, 0.173]);
            b3 = bar(sub_4,3,bd(3));
            set(b3,'FaceColor',[0.984 0.502 0.447]);
        hold(sub_4,'off');

    end

end

%% VALIDATION

function mbd = validate_mbd(mbd,t)

    if (~isnan(mbd))
        if (~isfinite(mbd))
            error('The value of ''mbd'' is invalid. Expected input to be finite.');
        end

        if (floor(mbd) ~= mbd)
            error('The value of ''mbd'' is invalid. Expected input to be an integer.');
        end

        b = ceil(0.2 * t);

        if ((mbd < 3) || (mbd > b))
            error(['The value of ''mbd'' is invalid. Expected input to have a value >= 5 and <= ' num2str(b) '.']);
        end
    end

end

function out = validate_output(out)

    [path,name,extension] = fileparts(out);

    if (~strcmpi(extension,'.xlsx'))
        out = fullfile(path,[name extension '.xlsx']);
    end

end

function temp = validate_template(temp)

    sheets = {'BUB' 'BMPH' 'BRPH' 'Breakdowns' 'Indicators'};
    file_sheets = validate_xls(temp,'T');

    if (~all(ismember(sheets,file_sheets)))
        error(['The template must contain the following sheets: ' sheets{1} sprintf(', %s',sheets{2:end}) '.']);
    end

    worksheets_batch(temp,sheets);

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
function [model] = PrepareNetwork4Separation(dnnFile, hasClean, useGPU)
dnn = load(dnnFile);

layer = dnn.layer;
para = dnn.para;
para.local.useFileName = 1;
% para.topology.useFileName = 1;
% para.local.seglen = 100;
% para.local.segshift = 100;
para.useGPU = useGPU;

% mixture STFT
stft_idx = ReturnLayerIdxByName(layer, 'stft');
para.out_layer_idx(1) = stft_idx(1);

% mixture log spec
log_idx = ReturnLayerIdxByName(layer, 'log');
para.out_layer_idx(2) = log_idx(1);

% enhanced log spec
hadamard_idx = ReturnLayerIdxByName(layer, 'hadamard');
if ~isempty(hadamard_idx)
    for i=1:2
        para.out_layer_idx(end+1) = hadamard_idx(i)+layer{hadamard_idx(i)}.next;
    end
    for i=1:2
        para.out_layer_idx(end+1) = hadamard_idx(i)+layer{hadamard_idx(i)}.prev(1);
    end
    para.test.mask_idx = [-1 0] + length(para.out_layer_idx);
else
    cmn_idx = ReturnLayerIdxByName(layer, 'cmn');
    if ~isempty(cmn_idx)
        for i=2:3
            para.out_layer_idx = [para.out_layer_idx cmn_idx(i)+layer{cmn_idx(i)}.prev];
        end
    else
        delta_idx = ReturnLayerIdxByName(layer, 'delta');
        para.out_layer_idx = [para.out_layer_idx delta_idx(2:3)-1];
    end
    para.test.mask_idx = [];
end

% clean STFT and log spec
input_idx = ReturnLayerIdxByName(layer, 'input');
if hasClean==0
    layer = layer(1:input_idx(end)-1);
    para.IO = RemoveIOStream(para.IO, 2);
    para.test.clean_idx = [];
    para.test.clean_stft_idx = [];
else
    para.out_layer_idx = [para.out_layer_idx log_idx(end-1:end)];
    para.test.clean_idx = [-1 0] + length(para.out_layer_idx);
    para.out_layer_idx = [para.out_layer_idx stft_idx(end-1:end)];
    para.test.clean_stft_idx = [-1 0] + length(para.out_layer_idx);
end

model.layer= layer;
model.para = para;
end

author: singaxiong

repo: SignalGraph


text:
function str = num2sci(A,precision,units)
%% FUNCTION num2sci
% Converts a number to a string in scientific format (e.g. .00325 ==> 3.25m).
% Can optionally provide precision (number of significant digits) and/or unit
% string to append.  The default precision is 4 significant digits.
%
%
% Examples:
%
% str = num2sci(165.48e-3) returns str = '165.5 m'
% str = num2sci(165.48e-3,'V') returns str = '165.5 mV'
% str = num2sci(165.48e-3,'mV') returns str = '165.5 uV'
% str = num2sci(165.48e-3,3,'V') returns str = '165 mV'
% str = num2sci([1.5 2.5],{'V','mA'}) returns str = {'1.5 V' '2.5 mA'}
%
% See also: num2str
%
%
% Written By: Jason Kaeding
%       Date: 07/10/09

%% REVISION INFORMATION:
% 11/23/09 - JAK - Added code to check units for cases when a prefix exists
%                  inside the unit string.  Examples are "kg" or "mV".  In these
%                  cases, the prefix is removed, A is recomputed based on this
%                  prefix, and then a new prefix is determined and used.
% 03/03/10 - JAK - Fixed bug in 11/23 revision which did not convert "kg"
%                  properly.

%% Check inputs
error(nargchk(1,3,nargin,'struct'));

if iscell(A)
    error('num2sci:cellInput','Input cannot be a cell array.');
end

if ischar(A)
    str = A;
    return;
end

if isempty(A)
    str = '';
    return;
end

switch nargin
    case 1
        precision = 4;
        units = '';
    case 2
        if ischar(precision) || iscellstr(precision)
            units = precision;
            precision = 4;
        else
            units = '';
        end
    case 3
        if ~ischar(units) && ~iscellstr(units)
            error('Units must be a string');
        end
end

if iscellstr(units)
    if ~all(size(units) == size(A))
        error('num2cell:badUnitsDim',...
            'Units must be the same size as input matrix.');
    end
end

%% Set up SI prefixes
prefix_vector = {'y','z','a','f','p','n','u','m','',...
    'k','M','G','T','P','E','Z','Y'};
prefix_offset = 9;

%% De-exponentiate
exponents = floor(log10(abs(A))./3);
exponents(~isfinite(A)) = 0; % +/- Inf and NaN assign to zero (no prefix)
exponents(A==0) = 0; % log10(0) = -Inf, correct this to 0 = 0*10^0
exponents = max(exponents,-8); % cap finite at 10^-24
exponents = min(exponents,8); % cap finite at 10^24

% Check for prefix in units
% Cases to check:
% 1. [prefix][A-Z][\w]* for typical units
% 2. [prefix][mgs] for SI units such as m,kg,s
% 3. [prefix][b][iyp/]?[\w/]* for bit/byte types of units (include / for bits/s)
unit_exp = 0;
if ~isempty(units)
    unit_check_results = regexp(units,...
        ['^([',[prefix_vector{:}],'])([A-Z]\w*|[mgs]|b[iyp/]?[\w/]*)$'],...
        'tokens','once');
    if ischar(units)
        unit_exp = 0;
        if ~isempty(unit_check_results)
            unit_exp = find(strcmp(unit_check_results{1},prefix_vector),1) - 9; % offset is 9
            units = unit_check_results{2}; % un-prefixed units
        end
    else
        % Check one by one
        unit_exp = zeros(size(exponents));
        for ind = 1:numel(unit_check_results)
            if ~isempty(unit_check_results{ind})
                unit_exp(ind) = find(strcmp(unit_check_results{ind}{1},...
                    prefix_vector),1) - 9; % offset is 9
                units{ind} = unit_check_results{ind}{2}; % un-prefixed units
            end
        end
    end
end

% Compute new base numbers
A = A./(10.^(3.*exponents));
exponents = exponents + unit_exp;

%% Find prefix for each element
prefix = cellfun(@(x)prefix_vector{x},num2cell(exponents+prefix_offset),...
    'UniformOutput',false);

%% Execute num2str
str = cellfun(@(x)[num2str(x,precision),' '],num2cell(A),'UniformOutput',false);

%% Build final string
% 1. Combine number with prefix & units
% 2. Check for num2str anomalies resulting from %g use (e.g. 32.4e+005), fix
% 3. Remove space if last character
str = regexprep(strcat(str,prefix,units),{'(\d+)\.?(\d+)?e[+](\d+)','\s+$'},...
    {'$1$2${repmat(''0'',1,str2double($3)-length($2))}',''});

%% De-cellify if single element
if numel(str) == 1
    str = str{1};
end


text:
% under construction...

function com = makecommand4history(p, skipfields, type, funcname)
%skipfields = {'ALLEEG', 'Saveas'};
fn     = fieldnames(p.Results);
if length(setindex)==1 && setindex(1)==1
      inputvari = 'EEG'; % Thanks to Felix Bacigalupo for this suggestion. Dic 12, 2011
else
      inputvari = inputname(1);
end
erpcom = sprintf( 'ERP = pop_averager( %s ', inputvari);
for q=1:length(fn)
      fn2com = fn{q};
      if ~ismember_bc2(fn2com, skipfields)
            fn2res = p.Results.(fn2com);
            if ~isempty(fn2res)
                  if ischar(fn2res)
                        if ~strcmpi(fn2res,'off')
                              erpcom = sprintf( '%s, ''%s'', ''%s''', erpcom, fn2com, fn2res);
                        end
                  else
                        if iscell(fn2res)
                              fn2resstr = vect2colon(cell2mat(fn2res), 'Sort','on');
                              fnformat = '{%s}';
                        else
                              fn2resstr = vect2colon(fn2res, 'Sort','on');
                              fnformat = '%s';
                        end
                        erpcom = sprintf( ['%s, ''%s'', ' fnformat], erpcom, fn2com, fn2resstr);
                  end
            end
      end
end
erpcom = sprintf( '%s );', erpcom);

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function [data] = ft_regressconfound(cfg, datain)

% FT_REGRESSCONFOUND estimates the regression weight of a set of confounds
% using a General Linear Model (GLM) and removes the estimated contribution
% from the single-trial data.
%
% Use as
%   timelock = ft_regressconfound(cfg, timelock)
% or as
%   freq     = ft_regressconfound(cfg, freq)
% or as
%   source   = ft_regressconfound(cfg, source)
%
% where timelock, freq, or, source come from FT_TIMELOCKANALYSIS,
% FT_FREQANALYSIS, or FT_SOURCEANALYSIS respectively, with keeptrials = 'yes'
%
% The cfg argument is a structure that should contain
%   cfg.confound    = matrix, [Ntrials X Nconfounds], may not contain NaNs
%
% The following configuration options are supported:
%   cfg.reject      = vector, [1 X Nconfounds], listing the confounds that
%                     are to be rejected (default = 'all')
%   cfg.normalize   = string, 'yes' or 'no', normalization to
%                     make the confounds orthogonal (default = 'yes')
%   cfg.output      = 'residual' (default), 'beta', or 'model'.
%                     If 'residual' is specified, the output is a data
%                     structure containing the residuals after regressing
%                     out the in cfg.reject listed confounds. If 'beta' or 'model'
%                     is specified, the output is a data structure containing
%                     the regression weights or the model, respectively.
%
% This method is described by Stolk et al., Online and offline tools for head
% movement compensation in MEG (Neuroimage, 2013)
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_REJECTCOMPONENT, FT_REJECTARTIFACT

% Undocumented local options:
%   cfg.ftest       = string array, {N X Nconfounds}, to F-test whether
%                     the full model explains more variance than reduced models
%                     (e.g. {'1 2'; '3 4'; '5'} where iteratively the added value of
%                     regressors 1 and 2, and then 3 and 4, etc., are tested)
%   cfg.statistics  = string, 'yes' or 'no', whether to add the statistics
%                     on the regression weights to the output (default = 'no',
%                     applies only when cfg.output = 'beta')

% Copyright (C) 2011-2017, Arjen Stolk, Robert Oostenveld, Lennart Verhagen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

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

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar datain
ft_preamble provenance datain
ft_preamble trackconfig

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

% check if the input data is valid for this function
datain = ft_checkdata(datain, 'datatype', {'timelock', 'freq', 'source'}, 'feedback', 'yes');

if isfield(cfg, 'beta') || isfield(cfg, 'model')
 ft_error('The options cfg.beta and cfg.model have been removed as of Aug 2017, please use cfg.output instead');
end

% ensure that the required options are present
cfg = ft_checkconfig(cfg, 'required', {'confound'}, 'renamed', {'Ftest','ftest'}, 'forbidden', {'beta','model'});

% specify the defaults
cfg.confound   = ft_getopt(cfg, 'confound');
cfg.reject     = ft_getopt(cfg, 'reject', 'all');
cfg.normalize  = ft_getopt(cfg, 'normalize', 'yes');
cfg.output     = ft_getopt(cfg, 'output', 'residual');
cfg.statistics = ft_getopt(cfg, 'statistics', 'no');
cfg.ftest      = ft_getopt(cfg, 'ftest');
cfg.parameter  = ft_getopt(cfg, 'parameter'); % the default is handled further down

regr = cfg.confound;
if any(isnan(regr(:)))
  ft_error('the confounds may not contain NaNs');
end
nconf     = size(regr,2);
conflist  = 1:nconf;
if strcmp(cfg.reject, 'all')
  cfg.reject = conflist(1:end); % to be removed
else
  cfg.reject = intersect(conflist, cfg.reject); % to be removed
end

fprintf('removing confound %s \n', num2str(cfg.reject));
kprs = setdiff(conflist, cfg.reject); % to be kept
fprintf('keeping confound %s \n', num2str(kprs));

% confound normalization for orthogonality
if strcmp(cfg.normalize, 'yes')
  fprintf('normalizing the confounds, except the constant \n');
  for c = 1:nconf
    AVG = mean(regr(:,c));
    STD = std(regr(:,c),0,1);
    if abs(STD/AVG)<10*eps
      fprintf('confound %s is a constant \n', num2str(c));
    else
      regr(:,c) = (regr(:,c) - AVG) / STD;
    end
    clear AVG STD;
  end
else
  fprintf('skipping normalization procedure \n');
end

switch ft_datatype(datain)
  case 'freq'
    cfg.parameter = ft_getopt(cfg, 'parameter', 'powspctrm');
  case 'timelock'
    cfg.parameter = ft_getopt(cfg, 'parameter', 'trial');
  case 'source'
    cfg.parameter = ft_getopt(cfg, 'parameter', 'pow');
end

dimord = getdimord(datain, cfg.parameter);
dimsiz = getdimsiz(datain, cfg.parameter);
dimtok = tokenize(dimord, '_');
rptdim = find(strcmp(dimtok, 'rpt'));
datdim = setdiff(1:length(dimtok), rptdim);

nrpt = dimsiz(rptdim);

dat = datain.(cfg.parameter);
if strcmp(dimtok{1}, '{pos}')
  indx = find(datain.inside);
  npos = length(indx);
  tmp = nan([nrpt npos datdim(2:end)]); % only positions inside the brain
  for i=indx'
    tmp(:,i,:,:,:) = dat{i}(:,:,:,:);
  end
  dat = tmp;
  haspermuted = false;
else
  dat = permute(dat, [rptdim datdim]);
  haspermuted = true;
end

dat = reshape(dat, nrpt, []);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GLM MODEL
%   Y = X * B + err, where Y is data, X is the model, and B are beta's
% which means
%   Best = X\Y ('matrix division', which is similar to B = inv(X)*Y)
% or when presented differently
%   Yest = X * Best
%   Yest = X * X\Y
%   Yclean = Y - Yest (the true 'clean' data is the recorded data 'Y' -
%   the data containing confounds 'Yest')
%   Yclean = Y - X * X\Y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% estimate and remove the confounds
fprintf('estimating the regression weights and removing the confounds \n');
if isempty(find(isnan(dat))) % if there are no NaNs, process all at once
  beta = regr\dat;                                                        % B = X\Y
else % otherwise process per colum set as defined by the nan distribution
  [u,i,j] = unique(~isnan(dat)','rows','first'); % find unique rows
  uniquecolumns = u'; % unique column types
  Nuniques = numel(i); % number of unique types
  beta_temp = NaN(Nuniques, nconf, size(dat,2)); % declare empty variable
  for n = 1:Nuniques % for each unique type
    rowidx = find(uniquecolumns(:,n)==1); % row indices for unique type
    colidx = find(j==n); % column indices for unique type
    if any(uniquecolumns(:,n)) % if vector contains a nonzero number
      beta_temp(n,:,colidx) = regr(rowidx,:)\dat(rowidx,colidx);         % B = X\Y
    end
  end
  beta = reshape(nansum(beta_temp,1),[nconf size(dat,2)]); % sum the betas
  clear beta_temp
end

model = regr(:, cfg.reject) * beta(cfg.reject, :);                        % model = confounds * weights = X * X\Y
Yc = dat - model;                                                         % Yclean = Y - X * X\Y

% reduced models analyses
if ~isempty(cfg.ftest)
  dfe        = nrpt - nconf;                                              % degrees of freedom
  err        = dat - regr * beta;                                         % err = Y - X * B
  tmse       = sum((err).^2)/dfe;                                         % mean squared error
  for iter = 1:numel(cfg.ftest)
    % regressors to test if they explain additional variance
    r          = str2num(cfg.ftest{iter});
    fprintf('F-testing explained additional variance of regressors %s \n', num2str(r));
    % regressors in reduced design (that is the original design)
    ri         = ~ismember(1:size(regr,2),r);
    rX         = regr(:,ri);               % reduced design
    rnr        = size(rX,2);               % number of regressors in reduced design
    % estimate reduced model betas
    rXcov      = pinv(rX'*rX);             % inverse design covariance matrix
    rb         = rXcov*rX'*dat;          	 % beta estimates using pinv
    % calculate mean squared error of reduced model
    rdfe       = size(dat,1) - size(rX,2); % degrees of freedom of the error
    rerr       = dat-rX*rb;                % residual error
    rmse       = sum(rerr'.^2,2)./rdfe;	   % mean squared error
    % F-test
    F(iter,:)  = ((rmse'-tmse)./(nconf-rnr)) ./ (tmse./(dfe-2));
    % Rik Henson defined F-test
    % F = ( ( rerr'*rerr - err'*err ) / ( nconf-rnr ) ) / ( err'*err/ ( nrpt-nconf ) );
    % convert F-value to p-value
    idx_pos    = F(iter,:) >= 0;
    idx_neg    = ~idx_pos;
    p(iter,:)     = nan(1,size(F(iter,:),2));
    p(iter,idx_pos) = (1-fcdf(F(iter,idx_pos),rnr,rdfe));
    p(iter,idx_neg) = fcdf(-F(iter,idx_neg),rnr,rdfe);
    clear rerr rmse
    % FIXME: drop in replace tcdf from the statfun/private dir
  end
  clear dfe err tmse
end

% organize the output
dataout = keepfields(datain, {'label', 'time', 'freq', 'pos', 'dim', 'transform', 'inside', 'outside', 'trialinfo', 'sampleinfo', 'dimord'});
switch cfg.output
  case 'residual'
    dataout.(cfg.parameter) = reshape(Yc, [nrpt dimsiz(datdim)]); % either powspctrm, trial, or pow
    if haspermuted
      dataout.(cfg.parameter) = ipermute(dataout.(cfg.parameter), [rptdim datdim]);
    end
    clear Yc
  case 'beta'
    dataout.beta = reshape(beta, [nconf, dimsiz(datdim)]);
    if haspermuted
      dataout.beta = ipermute(dataout.beta, [rptdim datdim]);
    end
    if strcmp(cfg.statistics, 'yes') % beta statistics
      fprintf('performing statistics on the regression weights \n');
      dfe        = nrpt - nconf;                                              % degrees of freedom
      err        = dat - regr * beta;                                         % err = Y - X * B
      mse        = sum((err).^2)/dfe;                                         % mean squared error
      covar      = diag(regr'*regr)';                                         % regressor covariance
      bvar       = repmat(mse',1,size(covar,2))./repmat(covar,size(mse,2),1); % beta variance
      tval       = (beta'./sqrt(bvar))';                                      % betas -> t-values
      prob       = (1-tcdf(tval,dfe))*2;                                      % p-values
      clear err dfe mse bvar
      % FIXME: drop in replace tcdf from the statfun/private dir
      dataout.stat = reshape(tval, [nconf dimsiz(datdim)]);
      dataout.prob = reshape(prob, [nconf dimsiz(datdim)]);
      if haspermuted
        dataout.stat = ipermute(dataout.stat, [rptdim datdim]);
        dataout.prob = ipermute(dataout.prob, [rptdim datdim]);
      end
      clear tval prob
    end
  case 'model'
    dataout.model = keepfields(datain, {'label', 'time', 'freq', 'pos', 'dim', 'transform', 'inside', 'outside', 'trialinfo', 'sampleinfo', 'dimord'});
    dataout.model.(cfg.parameter) = reshape(model, [nrpt, dimsiz(datdim)]);
    if haspermuted
      dataout.model.(cfg.parameter) = ipermute(dataout.model.(cfg.parameter), [rptdim datdim]);
    end
  otherwise
    error('output ''%s'' is not supported', cfg.output);
end

% reduced models analyses
if ~isempty(cfg.ftest)
  dataout.fvar   = reshape(F, [numel(cfg.ftest) dimsiz(datdim)]);
  dataout.pvar   = reshape(p, [numel(cfg.ftest) dimsiz(datdim)]);
  clear F p
end

% discard the gradiometer information because the weightings have been changed
if isfield(dataout, 'grad')
  ft_warning('discarding gradiometer information because the weightings have been changed');
  dataout = rmfield(dataout, 'grad');
end

% discard the electrode information because the weightings have been changed
if isfield(dataout, 'elec')
  ft_warning('discarding electrode information because the weightings have been changed');
  dataout = rmfield(dataout, 'elec');
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble trackconfig
ft_postamble previous datain

% rename the output variable to accomodate the savevar postamble
data = dataout;

ft_postamble provenance data
ft_postamble history    data
ft_postamble savevar    data

author: spm


text:
function p = prior_sqrtunif(varargin)
%PRIOR_SQRTUNIF  Uniform prior structure for the square root of the parameter
%
%  Description
%    P = PRIOR_SQRTUNIF creates uniform prior structure for the
%    square root of the parameter.
%
%  See also
%    PRIOR_*
%
% Copyright (c) 2009 Jarno Vanhatalo
% Copyright (c) 2010 Jaakko Riihimäki
% Copyright (c) 2010 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'PRIOR_SQRTUNIFORM';
  ip.addOptional('p', [], @isstruct);
  ip.parse(varargin{:});
  p=ip.Results.p;

  if isempty(p)
    init=true;
    p.type = 'Sqrt-Uniform';
  else
    if ~isfield(p,'type') && ~isequal(p.type,'Sqrt-Uniform')
      error('First argument does not seem to be a valid prior structure')
    end
    init=false;
  end

  if init
    % set functions
    p.fh.pak = @prior_sqrtunif_pak;
    p.fh.unpak = @prior_sqrtunif_unpak;
    p.fh.lp = @prior_sqrtunif_lp;
    p.fh.lpg = @prior_sqrtunif_lpg;
    p.fh.recappend = @prior_sqrtunif_recappend;
  end

end

function [w, s, h] = prior_sqrtunif_pak(p)
  w=[];
  s={};
  h=[];
end

function [p, w] = prior_sqrtunif_unpak(p, w)
  w = w;
  p = p;
end

function lp = prior_sqrtunif_lp(x, p)
  lJ  = -log(2*sqrt(x));  % log(1/(2*sqrt(x))) log(|J|) of transformation
  lp  = sum(lJ);
end

function lpg = prior_sqrtunif_lpg(x, p)
  lJg = -1./(2*x);        % gradient of log(|J|) of transformation
  lpg = lJg;
end

function rec = prior_sqrtunif_recappend(rec, ri, p)
% The parameters are not sampled in any case.
  rec = rec;
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function RHS = constantSourceTerm1D(phi)
% RHS vector for a Explicit source term
%
% SYNOPSIS:
%   RHS = constantSourceTerm1D(phi)
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% extract data from the mesh structure
Nx = phi.domain.dims(1);
G = 1:Nx+2;

% rearrange the matrix of k and build the sparse matrix for internal cells
row_index = reshape(G(2:Nx+1),Nx,1); % main diagonal (only internal cells)

% define the RHS Vector
RHS = zeros(Nx+2,1);

% assign the values of the RHS vector
RHS(row_index) = reshape(phi.value(2:Nx+1),Nx,1);

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function varargout=loadjd(filename, varargin)
%
%    data=loadjd(inputfile)
%       or
%    [data, mmap]=loadjd(inputfile, 'Param1',value1, 'Param2',value2,...)
%
%    Parse a hierarchical container data file, including JSON,
%    binary JSON (BJData/UBJSON/MessagePack) and HDF5, and output
%    the parsed data in a MATLAB/Octave data structure
%
%    author: Qianqian Fang (q.fang <at> neu.edu)
%
%    input:
%        inputfile: the input hierarchical container data file, supporting:
%                *.json,.jnii,.jdt,.jmsh,.jnirs: JSON/JData based data files, see https://neurojson.org/jdata/draft2
%                *.bjd,.bnii,.jdb,.bmsh,.bnirs: binary JData (BJData) files, see https://neurojson.org/bjdata/draft2
%                *.ubj: UBJSON-encoded files, see http://ubjson.org
%                *.msgpack: MessagePack-encoded files, see http://msgpack.org
%                *.h5,.hdf5,.snirf: HDF5 files, see https://www.hdfgroup.org/
%        options: (optional) for JSON/JData files, these are optional 'param',value pairs
%                supported by loadjson.m; for BJData/UBJSON/MessagePack files, these are
%                options supported by loadbj.m; for HDF5 files, these are options
%                supported by loadh5.m (part of EasyH5 toolbox, http://github.com/fangq/easyh5/)
%
%    output:
%        data: a structure (array) or cell (array) storing the hierarchical data
%              in the container data file
%        mmap: (optional) output storing the JSON/binary JSON memory-map table for fast
%              disk access. see help info for loadjson or loadbj for more details.
%
%    examples:
%        obj=struct('string','value','array',[1 2 3]);
%        savejd('obj', obj, 'datafile.json')
%        newobj=loadjd('datafile.json');
%
%    license:
%        BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of JSONLab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

if(nargin<1)
    error('you must provide file name');
end

if(regexpi(filename,'\.json$|\.jnii$\.jdt$\.jdat$|\.jmsh$|\.jnirs$'))
    [varargout{1:nargout}]=loadjson(filename,varargin{:});
elseif(regexpi(filename,'\.bjd$|\.bnii$\.jdb$\.jbat$|\.bmsh$|\.bnirs$'))
    [varargout{1:nargout}]=loadbj(filename,varargin{:});
elseif(regexpi(filename,'\.ubj$'))
    [varargout{1:nargout}]=loadubjson(filename,varargin{:});
elseif(regexpi(filename,'\.msgpack$'))
    [varargout{1:nargout}]=loadmsgpack(filename,varargin{:});
elseif(regexpi(filename,'\.h5$|\.hdf5$|\.snirf$'))
    if(~exist('loadh5','file'))
        error('you must first install EasyH5 from http://github.com/fangq/easyh5/');
    end
    [varargout{1:nargout}]=loadh5(filename,varargin{:});
else
    error('file suffix must be one of .json,.jnii,.jdt,.jmsh,.jnirs,.bjd,.bnii,.jdb,.bmsh,.bnirs,.ubj,.msgpack,.h5,.hdf5,.snirf');
end


text:
function p = prime_ge ( n )

%*****************************************************************************80
%
%% PRIME_GE returns the smallest prime greater than or equal to N.
%
%  Discussion:
%
%    The MATLAB version of this program is made much simpler
%    because of the availability of the IS_PRIME logical function.
%
%  Example:
%
%    N     PRIME_GE
%
%    -10    2
%      1    2
%      2    2
%      3    3
%      4    5
%      5    5
%      6    7
%      7    7
%      8   11
%      9   11
%     10   11
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 March 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number to be bounded.
%
%    Output, integer P, the smallest prime number that is greater
%    than or equal to N.
%
  p = max ( ceil ( n ), 2 );

  while ( ~ isprime ( p ) )
      p = p + 1;
  end

  return
end


text:
function mean = f_mean ( m, n )

%*****************************************************************************80
%
%% F_MEAN returns the mean of the F central PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the parameters of the PDF.
%    1 <= M,
%    1 <= N.
%    Note, however, that the mean is not defined unless 3 <= N.
%
%    Output, real MEAN, the mean of the PDF.
%
  if ( n < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'F_MEAN - Fatal error!\n' );
    fprintf ( 1, '  The mean is not defined for N < 3.\n' );
    error ( 'F_MEAN - Fatal error!' );
  end

  mean = n / ( n - 2 );

  return
end


text:
function sparse_grid_mixed_write ( dim_num, rule, alpha, beta, point_num, ...
  sparse_weight, sparse_point, file_name )

%*****************************************************************************80
%
%% SPARSE_GRID_MIXED_WRITE writes a sparse grid rule to five files.
%
%  Discussion:
%
%    The files are:
%    * the "A" file stores the ALPHA values, as a DIM_NUM x 1 list.
%    * the "B" file stores the BETA values, as a DIM_NUM x 1 list.
%    * the "R" file stores the region, as a DIM_NUM x 2 list.
%    * the "W" file stores the weights as a POINT_NUM list;
%    * the "X" file stores the abscissas as a DIM_NUM x POINT_NUM list;
%
%    The entries in the "R" file are the two corners of the DIM_NUM dimensional
%    rectangle that constitutes the integration region.  Coordinates that
%    should be infinite are set to 1.0E+30.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 March 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, int DIM_NUM, the spatial dimension.
%
%    Input, integer RULE(DIM_NUM), the rule in each dimension.
%     1, "CC",  Clenshaw Curtis, Closed Fully Nested rule.
%     2, "F2",  Fejer Type 2, Open Fully Nested rule.
%     3, "GP",  Gauss Patterson, Open Fully Nested rule.
%     4, "GL",  Gauss Legendre, Open Weakly Nested rule.
%     5, "GH",  Gauss Hermite, Open Weakly Nested rule.
%     6, "GGH", Generalized Gauss Hermite, Open Weakly Nested rule.
%     7, "LG",  Gauss Laguerre, Open Non Nested rule.
%     8, "GLG", Generalized Gauss Laguerre, Open Non Nested rule.
%     9, "GJ",  Gauss Jacobi, Open Non Nested rule.
%    10, "GW",  Golub Welsch, (presumed) Open Non Nested rule.
%    11, "CC_SE", Clenshaw Curtis Slow Exponential, Closed Fully Nested rule.
%    12, "F2_SE", Fejer Type 2 Slow Exponential, Closed Fully Nested rule.
%    13, "GP_SE", Gauss Patterson Slow Exponential, Closed Fully Nested rule.
%    14, "CC_ME", Clenshaw Curtis Moderate Exponential, Closed Fully Nested rule.
%    15, "F2_ME", Fejer Type 2 Moderate Exponential, Closed Fully Nested rule.
%    16, "GP_ME", Gauss Patterson Moderate Exponential, Closed Fully Nested rule.
%    17, "CCN", Clenshaw Curtis Nested, Linear, Closed Fully Nested rule.
%
%    Input, real ALPHA(DIM_NUM), BETA(DIM_NUM), parameters used for
%    Generalized Hermite, Generalized Laguerre, and Jacobi rules.
%
%    Input, int POINT_NUM, the number of points in the grid.
%
%    Input, real SPARSE_WEIGHT(POINT_NUM), the weights.
%
%    Input, real SPARSE_POINT(DIM_NUM,POINT_NUM), the points.
%
%    Input, string FILE_NAME, the main part of the file name.
%
  for dim = 1 : dim_num
    if ( rule(dim) == 1 )
      sparse_region(dim,1) = -1.0;
      sparse_region(dim,2) = +1.0;
    elseif ( rule(dim) == 2 )
      sparse_region(dim,1) = -1.0;
      sparse_region(dim,2) = +1.0;
    elseif ( rule(dim) == 3 )
      sparse_region(dim,1) = -1.0;
      sparse_region(dim,2) = +1.0;
    elseif ( rule(dim) == 4 )
      sparse_region(dim,1) = -1.0;
      sparse_region(dim,2) = +1.0;
    elseif ( rule(dim) == 5 )
      sparse_region(dim,1) = - r8_huge ( );
      sparse_region(dim,2) = + r8_huge ( );
    elseif ( rule(dim) == 6 )
      sparse_region(dim,1) = - r8_huge ( );
      sparse_region(dim,2) = + r8_huge ( );
    elseif ( rule(dim) == 7 )
      sparse_region(dim,1) = 0.0;
      sparse_region(dim,2) = r8_huge ( );
    elseif ( rule(dim) == 8 )
      sparse_region(dim,1) = 0.0;
      sparse_region(dim,2) = r8_huge ( );
    elseif ( rule(dim) == 9 )
      sparse_region(dim,1) = -1.0;
      sparse_region(dim,2) = +1.0;
    elseif ( rule(dim) == 10 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SPARSE_GRID_MIXED_WRITE - Fatal error!\n' );
      fprintf ( 1, '  Do not know how to specify region for rules of type 10.\n');
      error ( 'SPARSE_GRID_MIXED_WRITE - Fatal error!' );
    elseif ( rule(dim) == 11 )
      sparse_region(dim,1) = -1.0;
      sparse_region(dim,2) = +1.0;
    elseif ( rule(dim) == 12 )
      sparse_region(dim,1) = -1.0;
      sparse_region(dim,2) = +1.0;
    elseif ( rule(dim) == 13 )
      sparse_region(dim,1) = -1.0;
      sparse_region(dim,2) = +1.0;
    elseif ( rule(dim) == 14 )
      sparse_region(dim,1) = -1.0;
      sparse_region(dim,2) = +1.0;
    elseif ( rule(dim) == 15 )
      sparse_region(dim,1) = -1.0;
      sparse_region(dim,2) = +1.0;
    elseif ( rule(dim) == 16 )
      sparse_region(dim,1) = -1.0;
      sparse_region(dim,2) = +1.0;
    elseif ( rule(dim) == 17 )
      sparse_region(dim,1) = -1.0;
      sparse_region(dim,2) = +1.0;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SPARSE_GRID_MIXED_WRITE - Fatal error!\n' );
      fprintf ( 1, '  Unexpected value of RULE = %d\n', rule(dim) );
      error ( 'SPARSE_GRID_MIXED_WRITE - Fatal error!' );
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_MIXED_WRITE:\n' );

  file_name_a = strcat ( file_name, '_a.txt' );
  r8mat_write ( file_name_a, dim_num, 1, alpha );
  fprintf ( 1, '  Wrote the A file = "%s".\n', file_name_a );

  file_name_b = strcat ( file_name, '_b.txt' );
  r8mat_write ( file_name_b, dim_num, 1, beta );
  fprintf ( 1, '  Wrote the B file = "%s".\n', file_name_b );

  file_name_r = strcat ( file_name, '_r.txt' );
  r8mat_write ( file_name_r, dim_num, 2, sparse_region );
  fprintf ( 1, '  Wrote the R file = "%s".\n', file_name_r );

  file_name_w = strcat ( file_name, '_w.txt' );
  r8mat_write ( file_name_w, 1, point_num, sparse_weight );
  fprintf ( 1, '  Wrote the W file = "%s".\n', file_name_w );

  file_name_x = strcat ( file_name, '_x.txt' );
  r8mat_write ( file_name_x, dim_num, point_num, sparse_point );
  fprintf ( 1, '  Wrote the X file = "%s".\n', file_name_x );

  return
end


text:
function [cdstr, utstr] = jd2str(jdate)

% convert Julian date to string equivalent
% calendar date and universal time

% input

%  jdate = Julian date

% output

%  cdstr = calendar date string
%  utstr = universal time string

% Orbital Mechanics with Matlab

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[month, day, year] = gdate(jdate);

% serial date number

sdn = datenum(year, month, day);

% create calendar date string

cdstr = datestr(sdn, 1);

% create universal time string

utstr = datestr(day - fix(day), 'HH:MM:SS.FFF');

author: goGPS-Project