text:
classdef RunClass < TreeNodeClass

    methods

        % ----------------------------------------------------------------------------------
        function obj = RunClass(varargin)
            %
            % Syntax:
            %   obj = RunClass()
            %   obj = RunClass(filename);
            %   obj = RunClass(filename, iGroup, iSubj, iSess, iRun);
            %   obj = RunClass(run);
            %
            % Example 1:
            %   run1     = RunClass('./s1/neuro_run01.nirs',1,1,1,1);
            %   run1copy = RunClass(run1);
            %
            obj@TreeNodeClass(varargin);

            obj.type  = 'run';
            if nargin==0
                obj.name  = '';
                return;
            end
            if isa(varargin{1}, 'RunClass')
                obj.Copy(varargin{1});
                return;
            elseif isa(varargin{1}, 'FileClass')
                [~, ~, ~, obj.name] = varargin{1}.ExtractNames();
                obj.path         = varargin{1}.GetFilesPath();  % Fix wrong root path
            elseif ischar(varargin{1}) && strcmp(varargin{1},'copy')
                return;
            elseif ischar(varargin{1})
                obj.name = varargin{1};
            end
            if nargin==5
                obj.iGroup = varargin{2};
                obj.iSubj  = varargin{3};
                obj.iSess  = varargin{4};
                obj.iRun   = varargin{5};
            end

            obj.LoadAcquiredData();
            if obj.acquired.Error()
                obj = RunClass.empty();
                return;
            end
            obj.procStream = ProcStreamClass(obj.acquired);
            obj.InitTincMan();
            if isa(varargin{1}, 'FileClass')
                varargin{1}.Loaded();
            end

        end

        % ----------------------------------------------------------------------------------
        function b = Error(obj)
            if isempty(obj)
                b = -1;
                return;
            end
            b = obj.acquired.Error();
        end

        % ----------------------------------------------------------------------------------
        function err = Load(obj)
            err = 0;
            err1 = obj.LoadAcquiredData();
            err2 = obj.LoadDerivedData();
            if ~(err1==0 && err2==0)
                err = -1;
            end
        end

        % ----------------------------------------------------------------------------------
        function err = LoadDerivedData(obj)
            err = 0;
            if isempty(obj)
                return;
            end
            err = obj.procStream.Load([obj.path, obj.GetOutputFilename]);
        end

        % ----------------------------------------------------------------------------------
        function err = LoadAcquiredData(obj)
            err = -1;
            if isempty(obj)
                return;
            end

            if isempty(obj.SaveMemorySpace(obj.name))
                % Storage scheme is memory: In this case load acquisition data unconditionally.
                dataStorageScheme = 'memory';
            else
                dataStorageScheme = 'files';
            end

            if isempty(obj.acquired)
                if obj.IsNirs()
                    obj.acquired = NirsClass([obj.path, obj.name], dataStorageScheme);
                else
                    obj.acquired = SnirfClass([obj.path, obj.name], dataStorageScheme);
                end
            else
                obj.acquired.Load([obj.path, obj.name]);
            end

            if obj.acquired.Error() < 0
                obj.logger.Write( sprintf('     **** Error: "%s" failed to load - %s\n', obj.name, obj.acquired.GetErrorMsg()) );
                return;
            elseif obj.acquired.Error() > 0
                obj.logger.Write( sprintf('     **** Warning: %s in file "%s"\n', obj.acquired.GetErrorMsg(), obj.name) );
            else
                %fprintf('    Loaded file %s to run.\n', obj.name);
            end
            err = 0;
        end

        % ----------------------------------------------------------------------------------
        function FreeMemory(obj)
            if isempty(obj)
                return;
            end

            % Unload derived data
            obj.procStream.FreeMemory([obj.path, obj.GetOutputFilename()]);

            % Unload acquired data
            obj.acquired.FreeMemory([obj.path, obj.GetFilename()]);
        end

        % ----------------------------------------------------------------------------------
        function SaveAcquiredData(obj)
            if isempty(obj)
                return;
            end
            obj.procStream.input.SaveAcquiredData()
        end

        % ----------------------------------------------------------------------------------
        function b = AcquiredDataModified(obj)
            b = obj.procStream.AcquiredDataModified();
            if b
                obj.logger.Write(sprintf('Acquisition data for run %s has been modified\n', obj.name));
            end
        end

        % ----------------------------------------------------------------------------------
        % Copy processing params (procInut and procResult) from
        % N2 to N1 if N1 and N2 are same nodes
        % ----------------------------------------------------------------------------------
        function Copy(obj, obj2, conditional)
            if nargin==3 && strcmp(conditional, 'conditional')
                if obj.Mismatch(obj2)
                    return
                end
                obj.Copy@TreeNodeClass(obj2, 'conditional');
            else
                obj.Copy@TreeNodeClass(obj2);
                if isempty(obj.acquired)
                    if obj.IsNirs()
                        obj.acquired = NirsClass();
                    else
                        obj.acquired = SnirfClass();
                    end
                end
                obj.acquired.Copy(obj2.acquired);
            end
        end

        % --------------------------------------------------------------
        function CopyStims(obj, obj2)
            obj.CondNames = obj2.CondNames;
            obj.procStream.CopyStims(obj2.procStream);
        end

        % ----------------------------------------------------------------------------------
        % Subjects obj1 and obj2 are considered equivalent if their names
        % are equivalent and their sets of runs are equivalent.
        % ----------------------------------------------------------------------------------
        function B = equivalent(obj1, obj2)
            B=1;
            [p1,n1] = fileparts(obj1.name);
            [p2,n2] = fileparts(obj2.name);
            if ~strcmp([p1,'/',n1],[p2,'/',n2])
                B=0;
                return;
            end
        end

        % ----------------------------------------------------------------------------------
        function b = IsEmpty(obj)
            b = true;
            if isempty(obj)
                return;
            end
            if obj.acquired.IsEmpty()
                return;
            end
            b = false;
        end

        % ----------------------------------------------------------------------------------
        function b = IsEmptyOutput(obj)
            b = true;
            if isempty(obj)
                return;
            end
            obj.LoadDerivedData();
            if obj.procStream.IsEmptyOutput()
                return;
            end
            b = false;
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % ----------------------------------------------------------------------------------
        function b = IsNirs(obj)
            b = false;
            [~,~,ext] = fileparts(obj.name);
            if strcmp(ext,'.nirs')
                b = true;
            end
        end

        % ----------------------------------------------------------------------------------
        function varval = GetVar(obj, varname)
            varval = [];
            if isproperty(obj, varname)
                varval = eval( sprintf('obj.%s', varname) );
            end
            if isempty(varval)
                varval = obj.procStream.GetVar(varname);
            end
            if isempty(varval)
                varval = obj.acquired.GetVar(varname);
            end
        end

        % ----------------------------------------------------------------------------------
        function LoadInputVars(obj)

            % a) Find all variables needed by proc stream
            args = obj.procStream.GetInputArgs();

            % b) Find these variables in this run
            for ii = 1:length(args)
                eval( sprintf('obj.inputVars.%s = obj.GetVar(args{ii});', args{ii}) );
            end
        end

        % ----------------------------------------------------------------------------------
        function Calc(obj, options)
            if ~exist('options','var') || isempty(options)
                options = 'overwrite';
            end

            % Update call application GUI using it's generic Update function
            if ~isempty(obj.updateParentGui)
                obj.updateParentGui('DataTreeClass', [obj.iGroup, obj.iSubj, obj.iSess, obj.iRun]);
            end

            % Load acquired data
            obj.acquired.Load();

            if strcmpi(options, 'overwrite')
                % Recalculating result means deleting old results, if
                % option == 'overwrite'
                obj.procStream.output.Flush();
            end

            if obj.DEBUG
                fprintf('Calculating processing stream for group %d, subject %d, session %d, run %d\n', obj.iGroup, obj.iSubj, obj.iSess, obj.iRun);
            end

            % Find all variables needed by proc stream, find them in this run, and load them to proc stream input
            obj.LoadInputVars();

            Calc@TreeNodeClass(obj);

            if obj.DEBUG
                obj.logger.Write(sprintf('Completed processing stream for group %d, subject %d, session %d, run %d\n', obj.iGroup, obj.iSubj, obj.iSess, obj.iRun));
                obj.logger.Write('\n')
            end
        end

        % ----------------------------------------------------------------------------------
        function Print(obj, indent)
            if ~exist('indent', 'var')
                indent = 4;
            else
                indent = indent+4;
            end
            Print@TreeNodeClass(obj, indent);
        end

        % ---------------------------------------------------------------
        function PrintProcStream(obj)
            fcalls = obj.procStream.GetFuncCallChain();
            obj.logger.Write('Run processing stream:\n');
            for ii = 1:length(fcalls)
                obj.logger.Write('%s\n', fcalls{ii});
            end
        end

    end    % Public methods

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Pubic Set/Get methods for acquired data
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % ----------------------------------------------------------------------------------
        function t = GetTime(obj, iBlk)
            if nargin==1
                iBlk=1;
            end

            % Sometimes the caller is NOT the current element in which case we need to load
            % (when dataStorageScheme = 'file' mode) here explicitely.
            err = -1;
            if obj.acquired.IsEmpty()
                err = obj.acquired.LoadTime();
            end
            t = obj.acquired.GetTime(iBlk);
            if err==0
                obj.acquired.FreeMemory(obj.GetFilename());
            end
        end

        % ----------------------------------------------------------------------------------
        function t = GetTimeCombined(obj)
            % Sometimes the caller is NOT the current element in which case we need to load
            % (when dataStorageScheme = 'file' mode) here explicitely.
            err = -1;
            if obj.acquired.IsEmpty()
                err = obj.acquired.LoadTime();
            end
            t = obj.acquired.GetTimeCombined();
            if err==0
                obj.acquired.FreeMemory(obj.GetFilename());
            end
        end

        % ----------------------------------------------------------------------------------
        function t = GetAuxiliaryTime(obj)
            % Sometimes the caller is NOT the current element in which case we need to load
            % (when dataStorageScheme = 'file' mode) here explicitely. BUT this needs to be
            % optimized to NOT load the whole thing, just aux - see GetTime()/ GetTimeCombined().
            % TBD. jdubb, 08/17/2022
            err = -1;
            if obj.acquired.IsEmpty()
                obj.acquired.Load();
                err = obj.GetError();
            end
            t = obj.acquired.GetAuxiliaryTime();
            if err==0
                obj.acquired.FreeMemory(obj.GetFilename());
            end
        end

        % ----------------------------------------------------------------------------------
        function [d, t, ml] = GetRawData(obj, iBlk)
            if nargin<2
                iBlk = 1;
            end
            [d, t, ml] = obj.acquired.GetDataTimeSeries('', iBlk);
        end

        % ----------------------------------------------------------------------------------
        function [d, t, ml] = GetDataTimeSeries(obj, options, iBlk)
            if ~exist('options','var')
                options = '';
            end
            if ~exist('iBlk','var') || isempty(iBlk)
                iBlk = 1;
            end
            if isempty(options) || strcmp(options, 'reshape')
                [d, t, ml] = obj.acquired.GetDataTimeSeries(options, iBlk);
            else
                [d, t, ml] = obj.GetDataTimeSeries@TreeNodeClass(options, iBlk);
            end
        end

        % ----------------------------------------------------------------------------------
        function [iDataBlks, iCh] = GetDataBlocksIdxs(obj, iCh)
            if nargin<2
                iCh = [];
            end
            [iDataBlks, iCh] = obj.acquired.GetDataBlocksIdxs(iCh);
        end

        % ----------------------------------------------------------------------------------
        function n = GetDataBlocksNum(obj)
            n = obj.acquired.GetDataBlocksNum();
        end

        % ----------------------------------------------------------------------------------
        function SD = GetSDG(obj,option)
            % Sometimes the caller is NOT the current element in which case we need to load
            % (when dataStorageScheme = 'file' mode) here explicitely. BUT this needs to be
            % optimized to NOT load the whole thing, just SDG. TBD. jdubb, 08/17/2022
            err = -1;
            if obj.acquired.IsEmpty()
                obj.acquired.Load();
                err = obj.GetError();
            end
            SD.Lambda  = obj.acquired.GetWls();
            if exist('option','var')
                SD.SrcPos  = obj.acquired.GetSrcPos(option);
                SD.DetPos  = obj.acquired.GetDetPos(option);
            else
                SD.SrcPos  = obj.acquired.GetSrcPos();
                SD.DetPos  = obj.acquired.GetDetPos();
            end
            if err==0
                obj.acquired.FreeMemory(obj.GetFilename());
            end
        end

        % ----------------------------------------------------------------------------------
        function InitMlActMan(obj, iBlk)
            if ~exist('iBlk','var')
                iBlk = 1;
            end
            ml = obj.acquired.data.GetMeasListSrcDetPairs(iBlk);
            obj.procStream.input.SetMeasListActMan([ml, ones(size(ml, 1), 1)]);
        end

        % ----------------------------------------------------------------------------------
        function InitMlVis(obj, iBlk)
            if ~exist('iBlk','var')
                iBlk = 1;
            end
            ml = obj.acquired.data.GetMeasListSrcDetPairs(iBlk);
            obj.procStream.input.SetMeasListActMan([ml, ones(size(ml, 1), 1)]);
        end

        % ----------------------------------------------------------------------------------
        function ch = GetMeasList(obj, options, iBlk)
            if ~exist('iBlk','var') || isempty(iBlk)
                iBlk=1;
            end
            if ~exist('options','var')
                options = '';
            end

            ch = struct('MeasList',[], 'MeasListVis',[], 'MeasListActMan',[], 'MeasListActAuto',[]);

            % Sometimes the caller is NOT the current element in which case we need to load
            % (when dataStorageScheme = 'file' mode) here explicitely. BUT this needs to be
            % optimized to NOT load the whole thing, just SDG. TBD. jdubb, 08/17/2022
            err = -1;
            if obj.acquired.IsEmpty()
                obj.acquired.Load();
                err = obj.GetError();
            end

            ch.MeasList        = obj.acquired.GetMeasList(iBlk);
            ch.MeasListActMan  = obj.procStream.GetMeasListActMan(iBlk);
            ch.MeasListActAuto = obj.procStream.GetMeasListActAuto(iBlk);

            ch.MeasListActMan  = mlAct_Initialize(ch.MeasListActMan, ch.MeasList);
            ch.MeasListActAuto = mlAct_Initialize(ch.MeasListActAuto, ch.MeasList);
            if strcmp(options,'reshape')
                ch.MeasList = sortrows(ch.MeasList);
            end
            if err==0
                obj.acquired.FreeMemory(obj.GetFilename());
            end
        end

        % ----------------------------------------------------------------------------------
        function mlAct = GetActiveChannels(obj)
            % Load to memory if needed
            err = -1;
            if obj.procStream.output.IsEmpty()
                obj.Load();
                err = 0;
            end

            mlAct = obj.GetVar('mlActAuto');
            if ~isempty(mlAct)
                mlAct = mlAct{1};
            end

            % Free memory
            if err==0
                obj.FreeMemory();
            end
        end

        % ----------------------------------------------------------------------------------
        function SetStims_MatInput(obj, s, t, CondNames)
            obj.procStream.SetStims_MatInput(s, t, CondNames);
        end

        % ----------------------------------------------------------------------------------
        function ReloadStim(obj)
            % Update call application GUI using it's generic Update function
            if ~isempty(obj.updateParentGui)
                obj.updateParentGui('DataTreeClass', [obj.iGroup, obj.iSubj, obj.iSess, obj.iRun]);
            end
            if obj.DEBUG
                fprintf('group %d, subject %d, session %d, run %d\n', obj.iGroup, obj.iSubj, obj.iSess, obj.iRun);
            end
            obj.acquired.LoadStim(obj.acquired.GetFilename());
            obj.procStream.CopyStims(obj.acquired)
            pause(.5)
        end

        % ----------------------------------------------------------------------------------
        function s = GetStims(obj, t)
            % Proc stream output
            s_inp = obj.procStream.input.GetStims(t);
            s_out = obj.procStream.output.GetStims(t);

            k_inp_all  = find(s_inp~=0);
            k_out_edit = find(s_out~=0 & s_out~=1);

            % Select only those output stims which exist in the input
            b = ismember(k_out_edit, k_inp_all);

            s = s_inp;
            s(k_out_edit(b)) = s_out(k_out_edit(b));
        end

        % ----------------------------------------------------------------------------------
        function SetConditions(obj, CondNames)
            if nargin==2
                obj.procStream.SetConditions(CondNames);
            end
            obj.CondNames = unique(obj.procStream.GetConditions());
        end

        % ----------------------------------------------------------------------------------
        function CondNames = GetConditions(obj)
            CondNames = obj.procStream.GetConditions();
        end

        % ----------------------------------------------------------------------------------
        function CondNames = GetConditionsActive(obj)
            CondNames = obj.CondNames;
            t = obj.GetTime();
            s = obj.GetStims(t);
            for ii=1:size(s,2)
                if ismember(abs(1), s(:,ii))
                    CondNames{ii} = ['-- ', CondNames{ii}];
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function wls = GetWls(obj)
            wls = obj.acquired.GetWls();
        end

        % ----------------------------------------------------------------------------------
        function bbox = GetSdgBbox(obj)
            if obj.acquired.IsEmpty()
                % No need to load whole of acquired data need only probe here
                obj.acquired.LoadProbe(obj.acquired.GetFilename());
            end
            bbox = obj.acquired.GetSdgBbox();
        end

        % ----------------------------------------------------------------------------------
        function aux = GetAux(obj)
            aux = obj.acquired.GetAux();
        end

        % ----------------------------------------------------------------------------------
        function aux = GetAuxiliary(obj)
            aux = obj.acquired.GetAuxiliary();
        end

        % ----------------------------------------------------------------------------------
        function tIncAuto = GetTincAuto(obj, iBlk)
            if nargin<2
                iBlk = 1;
            end
            tIncAuto = obj.procStream.output.GetTincAuto(iBlk);
        end

        % ----------------------------------------------------------------------------------
        function tIncAutoCh = GetTincAutoCh(obj, iBlk)
            if nargin<2
                iBlk = 1;
            end
            tIncAutoCh = obj.procStream.output.GetTincAutoCh(iBlk);
        end

        % ----------------------------------------------------------------------------------
        function tIncMan = GetTincMan(obj, iBlk)
            if nargin<2
                iBlk = 1;
            end
            tIncMan = obj.procStream.input.GetTincMan(iBlk);
            if isempty(tIncMan)  % If the Tinc array is unitialized TODO find a more sensical place to do this
               obj.InitTincMan();
               tIncMan = obj.procStream.input.GetTincMan(iBlk);
            end
        end

        % ----------------------------------------------------------------------------------
        function SetTincMan(obj, idxs, iBlk, excl_incl)
            if nargin<2
                return
            end
            if nargin<4
                excl_incl = 'exclude';
            end
            tIncMan = obj.procStream.GetTincMan(iBlk);
            if strcmp(excl_incl, 'exclude')
                tIncMan(idxs) = 0;
            elseif strcmp(excl_incl, 'include')
                tIncMan(idxs) = 1;
            end
            obj.procStream.SetTincMan(tIncMan, iBlk);
        end

        % ----------------------------------------------------------------------------------
        function InitTincMan(obj)
            iBlk = 1;  % TODO implement multiple data blocks
            while 1
                t = obj.acquired.GetTime(iBlk);
                if isempty(t)
                    break
                end
                tIncMan = ones(length(t),1);
                obj.procStream.SetTincMan(tIncMan, iBlk);
                iBlk = iBlk+1;
            end
        end
    end        % Public Set/Get methods

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % All other public methods for acquired data
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % ----------------------------------------------------------------------------------
        function AddStims(obj, tPts, condition, duration, amp, more)
            if isempty(tPts)
                return;
            end
            if isempty(condition)
                return;
            end
            obj.procStream.AddStims(tPts, condition, duration, amp, more);
        end

        % ----------------------------------------------------------------------------------
        function DeleteStims(obj, tPts, condition)
            if ~exist('tPts','var') || isempty(tPts)
                return;
            end
            if ~exist('condition','var')
                condition = '';
            end
            obj.procStream.DeleteStims(tPts, condition);
        end

        % ----------------------------------------------------------------------------------
        function ToggleStims(obj, tPts, condition)
            if ~exist('tPts','var') || isempty(tPts)
                return;
            end
            if ~exist('condition','var')
                condition = '';
            end
            obj.procStream.ToggleStims(tPts, condition);
        end

        % ----------------------------------------------------------------------------------
        function MoveStims(obj, tPts, condition)
            if ~exist('tPts','var') || isempty(tPts)
                return;
            end
            if ~exist('condition','var')
                condition = '';
            end
            obj.procStream.MoveStims(tPts, condition);
        end

        % ----------------------------------------------------------------------------------
        function AddStimColumn(obj, name, initValue)
            if ~exist('name', 'var')
                return;
            end
            obj.procStream.AddStimColumn(name, initValue);
        end

        % ----------------------------------------------------------------------------------
        function DeleteStimColumn(obj, idx)
            if ~exist('idx', 'var') || idx <= 3
                return;
            end
            obj.procStream.DeleteStimColumn(idx);
        end

        % ----------------------------------------------------------------------------------
        function RenameStimColumn(obj, oldname, newname)
            if ~exist('oldname', 'var') || ~exist('newname', 'var')
                return;
            end
            obj.procStream.RenameStimColumn(oldname, newname);
        end

        % ----------------------------------------------------------------------------------
        function probe = GetProbe(obj)
            probe = obj.acquired.GetProbe();
        end

        % ----------------------------------------------------------------------------------
        function data = GetStimData(obj, icond)
            data = obj.procStream.GetStimData(icond);
        end

        % ----------------------------------------------------------------------------------
        function val = GetStimDataLabels(obj, icond)
            val = obj.procStream.GetStimDataLabels(icond);
        end

        % ----------------------------------------------------------------------------------
        function SetStimTpts(obj, icond, tpts)
            obj.procStream.SetStimTpts(icond, tpts);
        end

        % ----------------------------------------------------------------------------------
        function tpts = GetStimTpts(obj, icond)
            if ~exist('icond','var')
                icond=1;
            end
            tpts = obj.procStream.GetStimTpts(icond);
        end

        % ----------------------------------------------------------------------------------
        function SetStimDuration(obj, icond, duration, tpts)
            obj.procStream.SetStimDuration(icond, duration, tpts);
        end

        % ----------------------------------------------------------------------------------
        function duration = GetStimDuration(obj, icond)
            if ~exist('icond','var')
                icond=1;
            end
            duration = obj.procStream.GetStimDuration(icond);
        end

        % ----------------------------------------------------------------------------------
        function SetStimAmplitudes(obj, icond, amps, tpts)
            obj.procStream.SetStimAmplitudes(icond, amps, tpts);
        end

        % ----------------------------------------------------------------------------------
        function vals = GetStimAmplitudes(obj, icond)
            if ~exist('icond','var')
                icond=1;
            end
            vals = obj.procStream.GetStimAmplitudes(icond);
        end

        % ----------------------------------------------------------------------------------
        function RenameCondition(obj, oldname, newname)
            % Function to rename a condition. Important to remeber that changing the
            % condition involves 2 distinct well defined steps:
            %   a) For the current element change the name of the specified (old)
            %      condition for ONLY for ALL the acquired data elements under the
            %      currElem, be it run, subj, or group . In this step we DO NOT TOUCH
            %      the condition names of the run, subject or group .
            %   b) Rebuild condition names and tables of all the tree nodes group, subjects
            %      and runs same as if you were loading during Homer3 startup from the
            %      acquired data.
            %
            if ~exist('oldname','var') || ~ischar(oldname)
                return;
            end
            if ~exist('newname','var')  || ~ischar(newname)
                return;
            end
            newname = obj.ErrCheckNewCondName(newname);
            if obj.err ~= 0
                return;
            end
            obj.procStream.RenameCondition(oldname, newname);
        end

        % ----------------------------------------------------------------------------------
        function StimReject(obj, t, iBlk)
            obj.procStream.StimReject(t, iBlk);
        end

        % ----------------------------------------------------------------------------------
        function StimInclude(obj, t, iBlk)
            obj.procStream.StimInclude(t, iBlk);
        end

        % ----------------------------------------------------------------------------------
        function vals = GetStimValSettings(obj)
            vals = obj.procStream.input.GetStimValSettings();
        end

        % ----------------------------------------------------------------------------------
        function ExportStim(obj, options)
            global cfg
            if ~exist('options','var')
                options = '';
                if strcmpi(cfg.GetValue('Load Stim from TSV file'), 'no')
                    options = 'regenerate';
                end
            end
            SnirfFile2Tsv(obj.acquired, '', options);
        end

        % ----------------------------------------------------------------------------------
        function fname = DeleteExportStim(obj)
            fname = obj.acquired.GetStimTsvFilename();
            if ispathvalid(fname)
                try
                    fprintf('Delete  %s\n', fname)
                    delete(fname);
                catch
                end
            end
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % ----------------------------------------------------------------------------------
        function nbytes = MemoryRequired(obj, option)
            if ~exist('option','var')
                option = 'memory';
            end
            nbytes = obj.procStream.MemoryRequired();
            if strcmp(option, 'file')
                return
            end
            if isempty(obj.acquired)
                return
            end
            nbytes = nbytes + obj.acquired.MemoryRequired();
        end

        % -----------------------------------------------------------------
        function [fn_error, missing_args, prereqs] = CheckProcStreamOrder(obj)
            % Returns index of processing stream function which is missing
            % an argument, and a cell array of the missing arguments.
            % fn_error is 0 if there are no errors.

            missing_args = {};
            fn_error = 0;
            prereqs = '';

            % Processing stream begins with inputs available
            available = obj.procStream.input.GetProcInputs();
            % Inputs which are usually optional or defined elsewhere
            extras = {'iRun' 'iSubj' 'iGroup' 'mlActAuto', 'tIncAuto', 'Aaux', 'rcMap'};
            available = [available, extras];

            % For all fcalls
            for i = 1:length(obj.procStream.fcalls)

                inputs = obj.procStream.fcalls(i).GetInputs();

                % Check that each input is available
                for j = 1:length(inputs)
                    if ~any(strcmp(available, inputs{j}))
                       fn_error = obj.procStream.fcalls(i);
                       missing_args{end+1} = inputs{j}; %#ok<AGROW>
                    end
                end

                if isa(fn_error, 'FuncCallClass')
                    entry = obj.procStream.reg.GetEntryByName(fn_error.name);
                    if isfield(entry.help.sections, 'prerequisites')
                       prereqs_list = splitlines(entry.help.sections.prerequisites.str);
                       for k = 1:length(prereqs_list)
                           if ~isempty(prereqs_list{k})
                               prereqs = [prereqs, sprintf('\n'), strtrim(prereqs_list{k})];
                           end
                       end
                    end
                   return;
                end

                % Add outputs of the function to available list
                outputs = obj.procStream.fcalls(i).GetOutputs();
                for j = 1:length(outputs)
                   available{end + 1} = outputs{j};  %#ok<AGROW>
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function r = ListOutputFilenames(obj, options)
            if ~exist('options','var')
                options = '';
            end
            r = obj.GetOutputFilename(options);
            fprintf('    %s %s\n', obj.path, r);
        end

        % ----------------------------------------------------------------------------------
        function b = HaveOutput(obj)
            b1 = ~obj.procStream.output.IsEmpty();
            fname = obj.procStream.output.SetFilename([obj.path, obj.GetOutputFilename()]);
            b2 = false;
            if ispathvalid(fname)
                r = load(fname);
                b2 = ~r.output.IsEmpty();
            end
            b = b1 || b2;
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Private methods
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods (Access = public)

    end  % Private methods

end

author: BUNPC


text:
function p = polynomial_value ( m, o, c, e, n, x )

%*****************************************************************************80
%
%% POLYNOMIAL_VALUE evaluates a polynomial.
%
%  Discussion:
%
%    The polynomial is evaluated term by term, and no attempt is made to
%    use an approach such as Horner's method to speed up the process.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 December 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer O, the "order" of the polynomial.
%
%    Input, real C(O), the coefficients of the polynomial.
%
%    Input, integer E(O), the indices of the exponents of the polynomial.
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(M,N), the coordinates of the evaluation points.
%
%    Output, real P(N), the value of the polynomial at X.
%
  p = zeros ( n, 1 );

  for j = 1 : o
    f = mono_unrank_grlex ( m, e(j) );
    v = mono_value ( m, n, f, x );
    p(1:n) = p(1:n) + c(j) * v(1:n);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function x = p01_x ( n )

%*****************************************************************************80
%
%% P01_X returns the least squares solution X for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Output, real X(N,1), the least squares solution.
%
  x = [ -3.0200000; 4.4914286; -0.72857143 ];

  return
end


text:
function [ mOutputImage ] = ApplyBoxFilter( mInputImage, boxRadius, borderType, borderValue, normalizeFlag )
% ----------------------------------------------------------------------------------------------- %
% [ mFilteredImage ] = ApplyBoxFilter( mInputImage, boxBlurKernelRadius )
%   Applies Box Filter using Integral Images.
% Input:
%   - mInputImage   -   Input Image.
%                       Structure: Image Matrix (Single Channel)
%                       Type: 'Single' / 'Double'.
%                       Range: [0, 1].
%   - boxRadius     -   Box Radius.
%                       The radius of the box neighborhood for the
%                       summation process.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: {1, 2, ..., }.
% Output:
%   - mOutputImage  -   Output Image.
%                       Structure: Image Matrix (Single Channel)
%                       Type: 'Single' / 'Double'.
%                       Range: [0, 1].
% Remarks:
%   1.  References: "..."
%   2.  The running sum matches Intel IPP 'SumWindowRow' / 'SumWindowColumn'.
% TODO:
%   1.  s
%   Release Notes:
%   -   1.2.001     29/04/2016  Royi Avital
%       *   Fixed bug were the border pixels weren't calculated correctly.
%   -   1.2.000     29/04/2016  Royi Avital
%       *   Updated function input.
%   -   1.1.000     04/01/2016  Royi Avital
%       *   Using "Running Sum" instead of Integral Images / Summed Area
%           Table due to numeric issues.
%   -   1.0.000     14/03/2015  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF = 0;
ON  = 1;

BORDER_TYPE_CONSTANT    = 1;
BORDER_TYPE_CIRCULAR    = 2;
BORDER_TYPE_REPLICATE   = 3;
BORDER_TYPE_SYMMETRIC   = 4;

numRows = size(mInputImage, 1);
numCols = size(mInputImage, 2);

boxBlurKernelLength     = (2 * boxRadius) + 1;

mOutputImage    = zeros([(numRows + boxBlurKernelLength - 1), (numCols + boxBlurKernelLength - 1)]);
switch(borderType)
    case(BORDER_TYPE_CONSTANT)
        mInputImage     = padarray(mInputImage, [boxRadius, boxRadius], borderValue, 'both');
    case(BORDER_TYPE_CIRCULAR)
        mInputImage     = padarray(mInputImage, [boxRadius, boxRadius], 'circular', 'both');
    case(BORDER_TYPE_REPLICATE)
        mInputImage     = padarray(mInputImage, [boxRadius, boxRadius], 'replicate', 'both');
    case(BORDER_TYPE_SYMMETRIC)
        mInputImage     = padarray(mInputImage, [boxRadius, boxRadius], 'symmetric', 'both');
end

vRowsIdx = boxRadius + [1:numRows];
vColsIdx = boxRadius + [1:numCols];
removedIdx = boxRadius + 1;

vCurrSum = sum(mInputImage(1:boxBlurKernelLength, :), 1);
mOutputImage((boxRadius + 1), :) = vCurrSum;

for iRow = boxRadius + [2:numRows]
    vRowToAdd       = mInputImage((iRow + boxRadius) ,:);
    vRowToRemove    = mInputImage((iRow - removedIdx) ,:);
    vCurrSum        = vCurrSum + vRowToAdd - vRowToRemove;

    mOutputImage(iRow, :) = vCurrSum;
end

mInputImage = mOutputImage;

vCurrSum = sum(mInputImage(:, 1:boxBlurKernelLength), 2);
mOutputImage(:, (boxRadius + 1)) = vCurrSum;

for iCol = boxRadius + [2:numCols]
    vColToAdd       = mInputImage(:, (iCol + boxRadius));
    vColToRemove    = mInputImage(:, (iCol - removedIdx));
    vCurrSum        = vCurrSum + vColToAdd - vColToRemove;

    mOutputImage(:, iCol) = vCurrSum;
end

mOutputImage = mOutputImage(vRowsIdx, vColsIdx);

if(normalizeFlag == ON)
    mOutputImage = mOutputImage ./ (boxBlurKernelLength * boxBlurKernelLength);
end

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function planC = saveIVHMatrix(IVHNum, doseBinsV, volsHistV, planC)
%"saveIVHMatrix"
%   Store the doseBinsV and volsHistV with the binWidth shift removed.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.
%
%Usage:
%   function planC = saveIVHMatrix(IVHNum, doseBinsV, volsHistV, planC)

indexS = planC{end};

if length(doseBinsV) == 1
    doseBinsV(2) = 2*doseBinsV(1);
    volsHistV(2) = 0;
end

%Calculate width of bins, use to find middle of each bin.
binWidthsV = diff(doseBinsV);
lastBinWidth = binWidthsV(end);
binWidthsV(end+1) = lastBinWidth;

%Extract dose bin values, adding half to binwidth to get lower edge.
planC{indexS.IVH}(IVHNum).IVHMatrix(:,1) = doseBinsV - binWidthsV/2;
planC{indexS.IVH}(IVHNum).IVHMatrix(:,2) = volsHistV;
return;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% This plot a DEM map plus eq on top...
% Benoetigt gtopo302 Modifikation von gtopo30

import zmaptopo.TopoToFlag

report_this_filefun(mfilename('fullpath'));

gtopo30s([46 47],[6 7])
s4_south=20;
s3_north=50;
s2_west=30;
s1_east=50;

cegt=gtopo30s([s4_south s3_north],[s2_west s1_east]);
disp('GTOPO30 files')
disp(cegt)

%     try
%         l  = get(h1,'XLim');
%     catch
%         update(mainmap())
%         pltopo
%     end

% VON ZMAP WIEDER BENOETIGT

%     s1_east = l(2); s2_west = l(1);
%     l  = get(h1,'YLim');
%     s3_north = l(2); s4_south = l(1);
%     fac = 1;

if abs(s4_south-s3_north) > 10 | abs(s1_east-s2_west) > 10
    def = {'3'};
    ni2 = inputdlg('Decimation factor for DEM data?','Input',1,def);
    l = ni2{:};
    fac = str2double(l);
end

%     do = ['cd  ' hodi]; ; eval(do);
%     do = ['cd ' hodi fs 'dem' fs 'gtopo30']; eval(do);

%     if exist('gtopo30s2') == 2
%         fname = gtopo30s([s4_south s3_north],[ s2_west s1_east])
%         do = [' [tmap, tmapleg] = gtopo30(fname,fac,[s4_south s3_north],[ s2_west s1_east]); '];
%     else
%         do = [' [tmap, tmapleg] = gtopo30(''test'',fac,[s4_south s3_north],[ s2_west s1_east]); '];
%     end

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

% fname=gtopo30s([s4_south s3_north],[ s2_west s1_east])
% fname=upper(fname{1})
% [tmap, tmapleg] = gtopo30(fname,fac,[s4_south s3_north],[ s2_west s1_east]);
[tmap, tmapleg] = gtopo302('c:\ZMAP6\dem\gtopo30',fac,[s4_south s3_north],[s2_west s1_east]);

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

my = s4_south:1/tmapleg(1):s3_north+0.1;
mx = s2_west:1/tmapleg(1):s1_east+0.1;
vlon = mx;
vlat = my;
[m,n] = size(tmap);
toflag = TopoToFlag.five;

% entspricht pltopo mitplt='plo'

[existFlag,figNumber]=figure_exists('Topographic Map',1);

if existFlag == 0;  ac3 = 'new'; overtopo;   end
if existFlag == 1
    figure_w_normalized_uicontrolunits(to1)
    delete(gca); delete(gca);delete(gca)
end

hold on; axis off

axes('position',[0.13,  0.13, 0.65, 0.7]);
pcolor(mx(1:n),my(1:m),tmap); shading flat
demcmap(tmap);
hold on
h1topo = gca;
set(gca,'color',[ 0.341 0.776 1.000 ]')
%whitebg(gcf,[0 0 0]);

set(gca,'FontSize',12,'FontWeight','bold','TickDir','out','Ticklength',[0.02 0.02])
set(gcf,'Color','w','InvertHardcopy','off')
set(gcf,'renderer','zbuffer')

author: CelsoReyes


text:
function [T2Hot1] = T2Hot1(X,alpha)
%Hotelling's T-Squared test for one multivariate sample.
%
%   Syntax: function [T2Hot1] = T2Hot1(X,alpha)
%
%     Inputs:
%          X - multivariate data matrix.
%      alpha - significance level (default = 0.05).
%
%     Output:
%          n - sample-size.
%          p - variables.
%          T2 - Hotelling's T-Squared statistic.
%          Chi-sqr. or F - the approximation statistic test.
%          df's - degrees' of freedom of the approximation statistic test.
%          P - probability that null Ho: is true.
%
%
%    Example: For the example given by Johnson and Wichern (1992, p. 183),
%             with 20 cases (n = 20) and three variables (p = 3). We are interested
%             to test if there is any difference against the expected mean vector
%             [4 50 10] at a significance level = 0.10.
%                      --------------    --------------
%                       x1   x2   x3      x1   x2   x3
%                      --------------    --------------
%                      3.7  48.5  9.3    3.9  36.9 12.7
%                      5.7  65.1  8.0    4.5  58.8 12.3
%                      3.8  47.2 10.9    3.5  27.8  9.8
%                      3.2  53.2 12.0    4.5  40.2  8.4
%                      3.1  55.5  9.7    1.5  13.5 10.1
%                      4.6  36.1  7.9    8.5  56.4  7.1
%                      2.4  24.8 14.0    4.5  71.6  8.2
%                      7.2  33.1  7.6    6.5  52.8 10.9
%                      6.7  47.4  8.5    4.1  44.1 11.2
%                      5.4  54.1 11.3    5.5  40.9  9.4
%                      --------------    --------------
%
%             Total data matrix must be:
%              X=[3.7 48.5 9.3;5.7 65.1 8.0;3.8 47.2 10.9;3.2 53.2 12.0;3.1 55.5 9.7;
%              4.6 36.1 7.9;2.4 24.8 14.0;7.2 33.1 7.6;6.7 47.4 8.5;5.4 54.1 11.3;
%              3.9 36.9 12.7;4.5 58.8 12.3;3.5 27.8 9.8;4.5 40.2 8.4;1.5 13.5 10.1;
%              8.5 56.4 7.1;4.5 71.6 8.2;6.5 52.8 10.9;4.1 44.1 11.2;5.5 40.9 9.4];
%
%     Calling on Matlab the function:
%             T2Hot1(X)
%       Immediately it ask:
%             -Do you have an expected mean vector? (y/n):
%            For this example we must to put:
%             y  (meaning 'yes')
%            Then it ask:
%             -Give me the expected mean vector:
%            Giving the mean vector:
%             [4 50 10]
%            Otherwise (n; meaning 'no') it consider a zero mean vector.
%
%       Answer is:
% ------------------------------------------------------------------------------------
% Sample-size    Variables      T2          F           df1          df2          P
% ------------------------------------------------------------------------------------
%      20            3         9.7388     2.9045          3           17        0.0649
% ------------------------------------------------------------------------------------
% Mean vectors results significant.
%
%
%  Created by A. Trujillo-Ortiz and R. Hernandez-Walls
%             Facultad de Ciencias Marinas
%             Universidad Autonoma de Baja California
%             Apdo. Postal 453
%             Ensenada, Baja California
%             Mexico.
%             atrujo@uabc.mx
%             And the special collaboration of the post-graduate students of the 2002:2
%             Multivariate Statistics Course: Karel Castro-Morales, Alejandro Espinoza-Tenorio,
%             Andrea Guia-Ramirez.
%
%  Copyright (C) December 2002
%
%  References:
%
%  Johnson, R. A. and Wichern, D. W. (1992), Applied Multivariate Statistical Analysis.
%              3rd. ed. New-Jersey:Prentice Hall. pp. 180-181,199-200.
%

if nargin < 1,
    error('Requires at least one input argument.');
end;

if nargin < 2,
    alpha = 0.05;
end;

if (alpha <= 0 | alpha >= 1)
   fprintf('Warning: significance level must be between 0 and 1\n');
   return;
end;

[n,p]=size(X);

ask=input('Do you have an expected mean vector? (y/n): ','s');
if ask=='y'
   mu=input('Give me the expected mean vector: ');
else
   mu=zeros([1,p]);
end;

if n <= p,
   error('Warning: requires that sample-size (n) must be greater than the number of variables (p).');
   return;
else

   m=mean(X); %Mean vector from data matrix X.
   S=cov(X);  %Covariance matrix from data matrix X.
   T2=n*(m-mu)*inv(S)*(m-mu)'; %Hotelling's T-Squared statistic.

   if n >= 50 %Chi-square approximation.
      X2=T2;
      v=p; %Degrees of freedom.
      P=1-chi2cdf(X2,v); %Probability that null Ho: is true.
      disp(' ')
      fprintf('----------------------------------------------------------------------------\n');
      disp(' Sample-size    Variables      T2          Chi-sqr.         df          P')
      fprintf('----------------------------------------------------------------------------\n');
      fprintf('%8.i%13.i%15.4f%14.4f%11.i%14.4f\n\n',n,p,T2,X2,v,P);
      fprintf('----------------------------------------------------------------------------\n');
      if P >= alpha;
         disp('Mean vectors results not significant.');
      else
         disp('Mean vectors results significant.');
      end;
   else  %F approximation.
      F=(n-p)/((n-1)*p)*T2;
      v1=p;  %Numerator degrees of freedom.
      v2=n-p;  %Denominator degrees of freedom.
      P=1-fcdf(F,v1,v2);  %Probability that null Ho: is true.
      disp(' ')
      fprintf('-------------------------------------------------------------------------------------\n');
      disp(' Sample-size    Variables      T2          F           df1          df2          P')
      fprintf('-------------------------------------------------------------------------------------\n');
      fprintf('%8.i%13.i%15.4f%11.4f%9.i%14.i%14.4f\n\n',n,p,T2,F,v1,v2,P);
      fprintf('-------------------------------------------------------------------------------------\n');
      if P >= alpha;
         disp('Mean vectors results not significant.');
      else
         disp('Mean vectors results significant.');
      end;
   end;
end;

return;

author: xiuyechen

repo: FishExplorer


text:
function Map = createMap(Rob,Sen,Opt)

% CREATEMAP Create an empty Map structure.
%   Map = CREATEMAP(Rob,Sen,Lmk,Opt) creates the structure Map from the
%   information contained in Rob, Sen and Lmk, using options Opt. The
%   resulting structure is an EKF map with all empty spaces, able to host
%   all states necessary for Rob, Sen and Lmk. It contains the fields:
%       .used   flags vector to used states in the map
%       .x      state vector
%       .P      covariances matrix
%       .t      current time - used to control time updates

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.
%   Copyright 2015 Joan Sola @ IRI-UPC-CSIC.

switch lower(Opt.map.type)
    case 'ekf'
        % Help: Map has state and covariances matrix.
        % The state has:
        % - state of each robot
        % - state of each sensor being estimated
        % - state of each landmark
        % The covariances matrix has all cross-variances of the above.

        Map.type = 'ekf';

        RobStates = [Rob.state];
        SenStates = [Sen.state];

        % overall number of states needed to allocate robots, sensors and landmarks
        n = sum([RobStates.size SenStates.size Opt.map.lmkSize*Opt.map.numLmks]);

        Map.used = false(n,1);

        Map.x = zeros(n,1);
        Map.P = zeros(n,n);

    case 'graph'
        % Help: Map has nominal and manifold states.
        % The state has:
        % - state of each pose of the trajectory of each robot
        % - state of each landmark
        % The manifold state has the same, but manifold states.
        % The true state is not stored, but obtained by composing nominal
        % and manifold states.

        % TODO: sensor in map not supported yet.

        Map.type = 'graph';

        RobStates = [Rob.state];

        % overall number of states needed to allocate frames and landmarks
        nf = sum([RobStates.dsize].*Opt.map.numFrames); % number of states for frames
        nl = Opt.map.lmkDSize*Opt.map.numLmks;  % number of states for landmarks
        n = nf + nl;    % total number of states

        % Map occupancy
        Map.used = false(n,1);

        % State and manifold
        Map.x = zeros(n,1);
        Map.P = zeros(n,n);

        switch Opt.solver.decomposition

            case 'Cholesky'
                % Hessian matrix in the manifold
                Map.H  = sparse(n,n);
                Map.R  = sparse(n,n);
                Map.b  = zeros(n,1); % rhs vector.
                Map.mr = []; % range of used states in H

            case 'QR'
                m = 1000;
                Map.A  = sparse(m,n);
                Map.R  = sparse(n,n);
                Map.b  = zeros(m,1); % rhs vector.
                Map.d  = zeros(n,1); % rhs vector.
                Map.mr = []; % range of used states in A
                Map.fr = []; % range of used factors in A
                Map.pr = []; % permutation of re-ordered states in A

            case 'Schur'
                % Hessian matrix in the manifold
                Map.H    = sparse(n,n);
                Map.b    = zeros(n,1); % rhs vector.
                Map.sSff = sparse(nf,nf); % sqrt of the Schur complement
                Map.iHll = sparse(nl,nl); % Inverse of the landmarks Hessian
                Map.mr   = []; % range of used states in H

            otherwise
                error('??? Unknown graph solver ''%s''.', Opt.solver.decomposition)
        end

    otherwise

        error('??? Unknown Map type. Please use ''ekf'' or ''graph''.')

end

Map.t = 0; % Current Map's time

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
report_this_filefun(mfilename('fullpath'));
%load /nfs/alaska/home/stefan
b = a;

% Hier vielleicht dx, dy aendern ...
ve = [];dx = 0.1; dy = 0.1;

x0 = min(a.Longitude);
x1 = max(a.Longitude);
y0 = min(a.Latitude);
y1 = max(a.Latitude);

for x = x0:dx:x1
    for y = y0:dy:y1
        ve =    [ ve ; x y ];
    end
end

le = length(ve);
Y0 = zeros(le,1);

% Hier muss die richtige abminderung rein ...
for i = 1:length(b)
    di2 = deg2km((distance(ve(:,2),ve(:,1),repmat(b(i,2),le,1),repmat(b(i,1),le,1))));
    R = di2;
    r = sqrt(R.^2 + 5.57^2);
    M = b(i,6); % wenn mit felhler: + dM  (random mit standart)
    Y = -0.136 + 0.229*(M-6) - 0.778 * log10(r) ;
    Y = 10.^Y;
    c = [Y , Y0];
    mapga = max(c');
    Y0 = mapga';
end

mapga = mapga';

l1 = length(x0:dx:x1);
l2 = length(y0:dy:y1);

re = reshape(mapga,l2,l1);
rey = reshape(ve(:,2),l2,l1);
rex = reshape(ve(:,1),l2,l1);

figure
pcolor(rex,rey,re)
hold on
shading interp
overlay
colorbar


text:
function bnd = read_asa_bnd(fn)

% READ_ASA_BND reads an ASA boundary triangulation file
% converting the units of the vertices to mm

% Copyright (C) 2002, Robert Oostenveld
%
% 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$

Npnt = read_asa(fn, 'NumberPositions=', '%d');
Ntri = read_asa(fn, 'NumberPolygons=', '%d');
Unit = read_asa(fn, 'UnitPosition', '%s');

pnt = read_asa(fn, 'Positions', '%f');
if any(size(pnt)~=[Npnt,3])
  pnt_file = read_asa(fn, 'Positions', '%s');
  [path, name, ext] = fileparts(fn);
  fid = fopen_or_error(fullfile(path, pnt_file), 'rb', 'ieee-le');
  pnt = fread(fid, [3,Npnt], 'float')';
  fclose(fid);
end

tri = read_asa(fn, 'Polygons', '%f');
if any(size(tri)~=[Ntri,3])
  tri_file = read_asa(fn, 'Polygons', '%s');
  [path, name, ext] = fileparts(fn);
  fid = fopen_or_error(fullfile(path, tri_file), 'rb', 'ieee-le');
  tri = fread(fid, [3,Ntri], 'int32')';
  fclose(fid);
end

if strcmpi(Unit,'mm')
  pnt   = 1*pnt;
elseif strcmpi(Unit,'cm')
  pnt   = 100*pnt;
elseif strcmpi(Unit,'m')
  pnt   = 1000*pnt;
else
  ft_error(sprintf('Unknown unit of distance for triangulated boundary (%s)', Unit));
end

bnd.pnt = pnt;
bnd.tri = tri + 1;      % 1-offset instead of 0-offset
% bnd.tri = fliplr(bnd.tri);    % invert order of triangle vertices

author: spm


text:
%V_M2HTMLPWD - create html documentation of files in current directory
%
% M2HTML must be on the matlab path
% available from http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=4039&objectType=FILE

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: v_m2htmlpwd.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pd=pwd;
tl=pd(max([find(pd=='\') find(pd=='/')])+1:end);
cd('..');
m2html('mfiles',tl,'htmldir',[tl '\doc']);
cd(pd);
winopen(['doc\' tl '\index.html']);

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [verticeslist_org] = reorganize_verticeslist (mesh_total, A, mesh_outer, perim, verticeslist,step)

% --- Step 1: iterative reorganization of the vertices' list
% (each vertex must appear in the right order so that they can be
% subsequently linked together one by one).
% This uses mesh_vertex_nearest (D. Weber) and is based on the assumption
% that verticeslist are relatively regularly situated on a circle.

% NOTE: sometimes this function can return a small loop instead of the full
% perimeter and it's mandatory to check that with e.g. the size (number of
% vertices of the perimeter. Just starting the reorglist
% with another start_vertex usually overcomes that problem.
start_vertex=1;
reorglist=[];
while ( size (reorglist,2) ~= (size(verticeslist,2)+1) )

    if start_vertex>=size(verticeslist,2);
        step = step +1;
        [verticeslist]=SearchProjectionOnPial(mesh_total,mesh_outer,perim,step);
        start_vertex = 1;
    end

    reorglist=verticeslist(start_vertex);
    remaininglist=verticeslist;
    remaininglist(start_vertex)=[];

    % search the nearest vertex for the first vertex of the list (start_vertex)
    [nextindex,nextvalue]=mesh_vertex_nearest(mesh_total.vertices(remaininglist,:),mesh_total.vertices(verticeslist(start_vertex),:));
    clear nextvalue;

    % create the reorglist (reorganized list)
    reorglist=[reorglist remaininglist(nextindex)];

    % delete this vertex from the list of vertices (i.e. from the pool of
    % vertices to be reorganized)
    remaininglist(nextindex)=[];

    % repeat that one more time so that we are sufficiently far from the
    % start_vertex to add it to the pool of vertices to be reorganized.
    % Then continue iteratively until we close the loop by finding the start_vertex.
    [nextindex,nextvalue]=mesh_vertex_nearest(mesh_total.vertices(remaininglist,:),mesh_total.vertices(reorglist(2),:));
    clear nextvalue;
    reorglist=[reorglist remaininglist(nextindex)];
    remaininglist(nextindex)=[];

    % add the start_vertex
    remaininglist= [remaininglist verticeslist(start_vertex)];

    % continue to reorganize until start_vertex
    for z= 3: size(verticeslist,2)
        [nextindex,nextvalue]=mesh_vertex_nearest(mesh_total.vertices(remaininglist,:),mesh_total.vertices(reorglist(z),:));
        clear nextvalue;
        reorglist=[reorglist remaininglist(nextindex)];
        if remaininglist(nextindex) == verticeslist(start_vertex), break, end
            remaininglist(nextindex)=[];
    end
    start_vertex=start_vertex+1;
end

verticeslist_org = reorglist;
clear remaininglist;
clear nextindex;
clear nextvalue;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [im, scale] = readImage(imagePath)
% READIMAGE   Read and standardize image
%    [IM, SCALE] = READIMAGE(IMAGEPATH) reads the specified image file,
%    converts the result to SINGLE class, and rescales the image
%    to have a maximum height of 480 pixels, returing the corresponding
%    scaling factor SCALE.
%
%    READIMAGE(IM) where IM is already an image applies only the
%    standardization to it.

% Author: Andrea Vedaldi

% Copyright (C) 2013 Andrea Vedaldi
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

if ischar(imagePath)
  try
    im = imread(imagePath) ;
  catch
    error('Corrupted image %s', imagePath) ;
  end
else
  im = imagePath ;
end

im = im2single(im) ;

scale = 1 ;
if (size(im,1) > 480)
  scale = 480 / size(im,1) ;
  im = imresize(im, scale) ;
  im = min(max(im,0),1) ;
end

author: jianxiongxiao

repo: ProfXkit


text:
function edge_length = tetrahedron_edge_length_3d ( tetra )

%*****************************************************************************80
%
%% TETRAHEDRON_EDGE_LENGTH_3D returns edge lengths of a tetrahedron in 3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 May 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real TETRA(3,4), the tetrahedron vertices.
%
%    Output, real EDGE_LENGTH(6,1), the length of the edges.
%
  dim_num = 3;
  edge_length = zeros ( 6, 1 );

  k = 0;
  for j1 = 1 : 3
    for j2 = j1+1 : 4
      k = k + 1;
      edge_length(k) = r8vec_length ( dim_num, ...
        tetra(1:dim_num,j2) - tetra(1:dim_num,j1) );
    end
  end

  return
end


text:
% function to implement the Iterated Extended Kalman Filter (IEKF)
% Inputs:
%  OBSn - the observations (with noise)
%  xest - initial state space estimates
% Ouputs:
%  Xp - predicted states
function Xp = f_IEKF(OBSn,xest)

load avar % r1,r2, L, and T

tol = .1; % tolerance for iterations
diff = 1;
count = 0;

F = [1 T 0 0 0 0 0;
    0 1 0 0 0 0 0;
    0 0 1 T 0 0 0;
    0 0 0 1 0 0 0;
    0 0 0 0 1 0 T;
    0 0 0 0 0 1 T;
    0 0 0 0 0 0 1]; % state transition matrix

n = size(OBSn,2); % number of observations

Xp = zeros(7,n);  % make room

Pkp1 = 1e10*eye(7); %xest*xest'; %.1*ones(7,7);

%Pkp1 = xest*xest';
%Pkp1 = F*Pkp1*F';

%Pkp1 = (10*randn(7,1))*(10*randn(7,1)).';
%Pkp1 = F*Pkp1*F';

% for each observation
for i = 1:n;

    % if this is the first iteration the prior predicted estimate is xest
    % if this is not the first run the prior estimate is in Xp
    if i == 1
        xkm1 = xest;
    else
        xkm1 = Xp(:,i-1);
    end

    Pkm1 = Pkp1; % conditional covariance from last iteration

    % iterations are started with the predicted estimate from the last run
    xkn = xkm1;
    while ~(diff < tol || count > 9)

        count = count + 1;

        H = [gradest(@(x)f_h1(x),xkn); gradest(@(x)f_h2(x),xkn)];

        R = (.01*randn(2,1))*(.01*randn(2,1)).';

        Rdiag = diag(R); R = diag(Rdiag);

        K = Pkm1*H'*(H*Pkm1*H'+R)^-1;

        xkn_temp = xkm1 + K*(OBSn(:,i)-f_h(xkn)-H*(xkm1-xkn));

        diff = norm(abs(xkn_temp-xkn));

        fprintf('diff = %g \n',diff)

        xkn = xkn_temp;

    end

    H = [gradest(@(x)f_h1(x),xkn); gradest(@(x)f_h2(x),xkn)];

    Pkk = (eye(7)-K*H)*Pkm1;

    Pkp1 = F*Pkk*F';

    Xp(:,i) = F*xkn;

    clc;

    fprintf('i = %g; Count is %g \n',i,count)

    count = 0;

    diff = 1;

end

author: Sable

repo: mcbench-benchmarks


text:
function x = mods235(p)
  n=2;

  x(1)=0.1*(p(1)-1.0);
  x(2)=p(2)-p(1)*p(1);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [F,P] = spm_MH_reml(YY,X,Q,N,hE);
% Estimation of covariance components from y*y' using sampling
% FORMAT [F,P] = spm_MH_reml(YY,X,Q,N,[hE]);
%
% YY  - (m x m) sample covariance matrix Y*Y'  {Y = (m x N) data matrix}
% X   - (m x p) design matrix
% Q   - {1 x q} covariance components
% N   - number of samples
%
% hE  - prior expectation: log-normal hyper-parameterisation (with hyperpriors)
%
% F   - [-ve] free energy F = log evidence = p(Y|X,Q)
% P   - smaple of hyperparameters from thier posterioir p(h|YY,X,Q)
%--------------------------------------------------------------------------
%
% This routiens using MCMC sampling (reverible Metropolis-Hastings)
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_MH_reml.m 5033 2012-11-02 20:59:54Z karl $

% assume a single sample if not specified
%--------------------------------------------------------------------------
try
    N;
catch
    N  = 1;
end

% assume OPT = 0
%--------------------------------------------------------------------------
try
    hE;
    OPT = hE;
catch
    OPT = 0;
end

% ortho-normalise X
%--------------------------------------------------------------------------
if isempty(X)
    X = sparse(length(Q{1}),1);
else
    X = orth(full(X));
end

% remove fixed effects
%--------------------------------------------------------------------------
n     = length(Q{1});
m     = length(Q);
h     = zeros(m,1);
R     = speye(n,n) - X*X';
YY    = R*YY*R;

M.OPT = OPT;
M.Q   = Q;
M.N   = N;

% initialise and specify hyperpriors
%--------------------------------------------------------------------------
[C,h,Ph,Fr] = spm_reml(YY,X,Q,N,0,4,OPT);
if M.OPT
    M.hE  = h - 16;
    M.hP  = eye(m,m)/32;
else
    M.hE  = zeros(m,1);
    M.hP  = speye(m,m)/exp(32);
end

% sample
%--------------------------------------------------------------------------
[P,F] = spm_MH('spm_MH_reml_likelihood',h,YY,M);


text:
function [predicted_label] = rec_test(model)
% Given the full 3D shape, recognition test for a generative model.
% After the model is (generatively) trained, this function test the
% recognition ability(discriminative) of the model using a seperate test
% data.
% The classification method is free-energy.

debug = 0;
if ~isfield(model.layers{2},'uw')
    model = merge_model(model);
end
global kConv_forward2 kConv_forward_c;

num_layer = length(model.layers);

test_list = read_data_list(model.data_path, model.classnames, ...
    model.volume_size + 2 * model.pad_size, 'test', debug);

label = [];
new_list = repmat(struct('filename', '', 'label', 0), 1, 1);
now = 1;
for i = 1 : model.classes
    cnt = length(test_list{i});
    label(now:now+cnt-1,1) = i;
    new_list(now:now+cnt-1,1) = test_list{i};
    now = now + cnt;
end
n = length(label);

batch_size = 32;
batch_num = ceil(n / batch_size);
predicted_label = zeros(n, model.classes); % prediction results.
for b = 1 : batch_num
    batch_end = min(n, b * batch_size);
	batch_index = (b-1)*batch_size + 1 : batch_end;
    batch_data = read_batch(model, new_list(batch_index), false);
    this_size = size(batch_data, 1);

    % propagate/inference bottum up using recognition weight.
    for l = 2 : num_layer - 1
        if l == 2
            stride = model.layers{l}.stride;
            hidden_presigmoid = myConvolve2(kConv_forward2, batch_data, model.layers{l}.uw, stride, 'forward');
            hidden_presigmoid = bsxfun(@plus, hidden_presigmoid, permute(model.layers{l}.c, [2,3,4,5,1]));
        elseif strcmp(model.layers{l}.type, 'convolution')
            stride = model.layers{l}.stride;
            hidden_presigmoid = myConvolve(kConv_forward_c, batch_data, model.layers{l}.uw, stride, 'forward');
            hidden_presigmoid = bsxfun(@plus, hidden_presigmoid, permute(model.layers{l}.c, [2,3,4,5,1]));
        else
            batch_data = reshape(batch_data, size(batch_data,1), []);
            hidden_presigmoid = bsxfun(@plus, ...
                batch_data * model.layers{l}.uw, model.layers{l}.c);
        end
        batch_data = 1 ./ ( 1 + exp(-hidden_presigmoid) );
    end

    batch_data = reshape(batch_data, this_size, []);
    % calculate the free energy for each label hypothesis
    for c = 1 : model.classes
        try_label = zeros(this_size, model.classes);
        try_label(:, c) = 1;
        predicted_label((b-1) * batch_size+1: batch_end, c) = free_energy(model, [try_label, batch_data], num_layer);
    end

end

predicted_label = bsxfun(@rdivide, predicted_label, sum(predicted_label, 2));
[~, predicted_label] = max(predicted_label, [], 2);

acc = sum(predicted_label == label) / n;
fprintf('acc is : %f\n', acc * 100);

author: zhirongw


text:
function computeFluxConsistentReactionPresence(modelFolder,propertiesFolder,reconVersion)
% This function extracts the presence of flux consistent reactions for a
% resource of reconstructions that were refined through the semi-automatic
% refinement pipeline (1 = present in the flux consistent submodel, 0 = not
% present in the flux consistent submodel).
%
% USAGE
%   computeFluxConsistentReactionPresence(modelFolder,propertiesFolder,reconVersion)
%
% INPUTS
% modelFolder                                                                                Folder with COBRA models to be analyzed
% propertiesFolder      Folder where the retrieved reaction presences will
%                       be stored (default: current folder)
% reconVersion          Name assigned to the reconstruction resource
%
%   - AUTHOR
%   Almut Heinken, 12/2020

mkdir([propertiesFolder filesep 'ReactionMetabolitePresence'])
currentDir=pwd;
cd([propertiesFolder filesep 'ReactionMetabolitePresence'])

dInfo = dir(modelFolder);
modelList={dInfo.name};
modelList=modelList';
modelList(~(contains(modelList(:,1),{'.mat','.sbml','.xml'})),:)=[];

% check if output file already exists
if isfile(['ReactionPresence_' reconVersion '.txt'])
    reactionPresence = readInputTableForPipeline(['ReactionPresence_' reconVersion '.txt']);
    allRxns=ReactionPresence(1,2:end)';
else
    % restart from existing data if possible
    if isfile([propertiesFolder filesep 'Reactions_' reconVersion '.txt'])
        reactions = readInputTableForPipeline([propertiesFolder filesep 'Reactions_' reconVersion '.txt']);
        allRxns=reactions(:,1);
    else
        allRxns={};
        for i=1:length(modelList)
            i
            model=readCbModel([modelFolder filesep modelList{i} '.mat']);
            allRxns=unique(vertcat(allRxns,model.rxns));
        end
    end
end

% remove models that were already retrieved
modelsRenamed=strrep(modelList(:,1),'.mat','');
modelsRenamed=strrep(modelsRenamed,'.sbml','');
modelsRenamed=strrep(modelsRenamed,'.xml','');
[C,IA]=intersect(modelsRenamed,ReactionPresence(2:end,1));
modelList(IA,:)=[];

% define the intervals in which the computations will be performed
if length(modelList)>5000
    steps=2000;
elseif length(modelList)>200
    steps=200;
else
    steps=25;
end

% in case of reruns, skip if all models are already analyzed
if ~isempty(modelList)
    for i=1:steps:length(modelList)
        if length(modelList)-i>=steps-1
            endPnt=steps-1;
        else
            endPnt=length(modelList)-i;
        end

        modelsToLoad={};
        for j=i:i+endPnt
            if j <= length(modelList)
                modelsToLoad{j}=[modelFolder filesep modelList{j} '.mat'];
            end
        end

        rxnsTmp={};
        parfor j=i:i+endPnt
            model=readCbModel(modelsToLoad{j});
            [~, ~, ~, ~, model] = findFluxConsistentSubset(model);
            rxnsTmp{j}=model.rxns;
        end

        for j=i:i+endPnt
            plusonerow=size(ReactionPresence,1)+1;

            ReactionPresence{plusonerow,1}=strrep(modelList{j},'.mat','');
            for k=1:length(allRxns)
                if ~isempty(find(ismember(rxnsTmp{j},allRxns{k})))
                    ReactionPresence{plusonerow,k+1}=1;
                else
                    ReactionPresence{plusonerow,k+1}=0;
                end
            end
        end
        % export the results as a table
        writetable(cell2table(ReactionPresence),['FluxConsistent_ReactionPresence_' reconVersion],'FileType','text','WriteVariableNames',false,'Delimiter','tab');
    end
end

cd(currentDir)

end

author: opencobra

repo: cobratoolbox


text:
function [bnet, LL, engine] = learn_params_dbn_em(engine, evidence, varargin)
% LEARN_PARAMS_DBN Set the parameters in a DBN to their ML/MAP values using batch EM.
% [bnet, LLtrace, engine] = learn_params_dbn_em(engine, data, ...)
%
% data{l}{i,t} = value of node i in slice t of time-series l, or [] if hidden.
%   Suppose you have L time series, each of length T, in an O*T*L array D,
%   where O is the num of observed scalar nodes, and N is the total num nodes per slice.
%   Then you can create data as follows, where onodes is the index of the observable nodes:
%      data = cell(1,L);
%      for l=1:L
%        data{l} = cell(N, T);
%        data{l}(onodes,:) = num2cell(D(:,:,l));
%      end
% Of course it is possible for different sets of nodes to be observed in
% each slice/ sequence, and for each sequence to be a different length.
%
% LLtrace is the learning curve: the vector of log-likelihood scores at each iteration.
%
% Optional arguments [default]
%
% max_iter - specifies the maximum number of iterations [100]
% thresh - specifies the thresold for stopping EM [1e-3]
%   We stop when |f(t) - f(t-1)| / avg < threshold,
%   where avg = (|f(t)| + |f(t-1)|)/2 and f is log lik.
% verbose - display loglik at each iteration [1]
% anneal - 1 means do deterministic annealing (only for entropic priors) [0]
% anneal_rate - geometric cooling rate [0.8]
% init_temp - initial annealing temperature [10]
% final_temp - final annealing temperature [1e-3]
%

max_iter = 100;
thresh = 1e-3;
anneal = 0;
anneal_rate = 0.8;
init_temp = 10;
final_temp = 1e-3;
verbose = 1;

for i=1:2:length(varargin)
  switch varargin{i}
   case 'max_iter', max_iter = varargin{i+1};
   case 'thresh', thresh = varargin{i+1};
   case 'anneal', anneal = varargin{i+1};
   case 'anneal_rate', anneal_rate = varargin{i+1};
   case 'init_temp', init_temp = varargin{i+1};
   case 'final_temp', final_temp = varargin{i+1};
   otherwise, error(['unrecognized argument' varargin{i}])
  end
end

% take 1 EM step at each temperature value, then when temp=0, run to convergence
% When using an entropic prior, Z = 1-T, so
% T=2 => Z=-1 (max entropy)
% T=1 => Z=0 (max likelihood)
% T=0 => Z=1 (min entropy / max structure)
num_iter = 1;
LL = [];
if anneal
  temperature = init_temp;
  while temperature > final_temp
    [engine, loglik, logpost] = EM_step(engine, evidence, temperature);
    if verbose
      fprintf('EM iteration %d, loglik = %8.4f, logpost = %8.4f, temp=%8.4f\n', ...
	      num_iter, loglik, logpost, temperature);
    end
    num_iter = num_iter + 1;
    LL = [LL loglik];
    temperature = temperature * anneal_rate;
  end
  temperature = 0;
  previous_loglik = loglik;
  previous_logpost = logpost;
else
  temperature = 0;
  previous_loglik = -inf;
  previous_logpost = -inf;
end

converged = 0;
while ~converged & (num_iter <= max_iter)
  [engine, loglik, logpost] = EM_step(engine, evidence, temperature);
  if verbose
    %fprintf('EM iteration %d, loglik = %8.4f, logpost = %8.4f\n', ...
    %	    num_iter, loglik, logpost);
    fprintf('EM iteration %d, loglik = %8.4f\n', num_iter, loglik);
  end
  num_iter = num_iter + 1;
  [converged, decreased] = em_converged(loglik, previous_loglik, thresh);
  %[converged, decreased] = em_converged(logpost, previous_logpost, thresh);
  previous_loglik = loglik;
  previous_logpost = logpost;
  LL = [LL loglik];
end

bnet = bnet_from_engine(engine);

%%%%%%%%%

function [engine, loglik, logpost] = EM_step(engine, cases, temp)

bnet = bnet_from_engine(engine); % engine contains the old params that are used for the E step
ss = length(bnet.intra);
CPDs = bnet.CPD; % these are the new params that get maximized
num_CPDs = length(CPDs);

% log P(theta|D) = (log P(D|theta) + log P(theta)) - log(P(D))
% where log P(D|theta) = sum_cases log P(case|theta)
% and log P(theta) = sum_CPDs log P(CPD) - only count once even if tied!
% logpost = log P(theta,D) (un-normalized)
% This should be negative, and increase at every step.

adjustable = zeros(1,num_CPDs);
logprior = zeros(1, num_CPDs);
for e=1:num_CPDs
  adjustable(e) = adjustable_CPD(CPDs{e});
end
adj = find(adjustable);

for e=adj(:)'
  logprior(e) = log_prior(CPDs{e});
  CPDs{e} = reset_ess(CPDs{e});
end

loglik = 0;
for l=1:length(cases)
  evidence = cases{l};
  if ~iscell(evidence)
    error('training data must be a cell array of cell arrays')
  end
  [engine, ll] = enter_evidence(engine, evidence);
  assert(~isnan(ll))
  loglik = loglik + ll;
  T = size(evidence, 2);

  % We unroll ns etc because in update_ess, we refer to nodes by their unrolled number
  % so that they extract evidence from the right place.
  % (The CPD should really store its own version of ns and cnodes...)
  ns = repmat(bnet.node_sizes_slice(:), [1 T]);
  cnodes = unroll_set(bnet.cnodes_slice, ss, T);

  %hidden_bitv = repmat(bnet.hidden_bitv(1:ss), [1 T]);
  hidden_bitv = zeros(ss, T);
  hidden_bitv(isemptycell(evidence))=1;
  % hidden_bitv(i) = 1 means node i is hidden.
  % We pass this in, rather than using isemptycell(evidence(dom)), because
  % isemptycell is very slow.

  t = 1;
  for i=1:ss
    e = bnet.equiv_class(i,1);
    if adjustable(e)
      fmarg = marginal_family(engine, i, t);
      CPDs{e} = update_ess(CPDs{e}, fmarg, evidence, ns(:), cnodes(:), hidden_bitv(:));
    end
  end

  for i=1:ss
    e = bnet.equiv_class(i,2);
    if adjustable(e)
      for t=2:T
	fmarg = marginal_family(engine, i, t);
	CPDs{e} = update_ess(CPDs{e}, fmarg, evidence, ns(:), cnodes(:), hidden_bitv(:));
      end
    end
  end
end

logpost = loglik + sum(logprior(:));

for e=adj(:)'
  CPDs{e} = maximize_params(CPDs{e}, temp);
end

engine = update_engine(engine, CPDs);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [X,D,e] = sparse_solver_irls(p, A, Y, beta, terminationValue, maxIterations)
%SPARSE_SOLVER_IRLS Steered-response power map using a MVDR beamformer
%
%   This routine performs sparse recovery using the iterative reweighted
%   least squares (IRLS) algorithm, with an L_{p,2} norm, assuming
%   spatiotemporal signals Y, of M channels/sensors/microphones and T
%   temporal frames (snapshots), and a linear mixing model of the form
%   Y = A*X, where A is a MxK overcomplete basis dictionary, with K>>M.
%   The original algorithm is described in:
%
%   Chartrand, R., & Yin, W. (2008). Iteratively reweighted algorithms for
%   compressive sensing. In 2008 IEEE International Conference on Acoustics,
%   Speech and Signal Processing (ICASSP), pp. 3869-3872.
%
%   and further inspired by the microphone array processing work of
%
%   Wabnitz, A., Epain, N., McEwan, A., & Jin, C. (2011). Upscaling ambisonic
%   sound scenes using compressed sensing techniques. In 2011 IEEE Workshop
%   on Applications of Signal Processing to Audio and Acoustics (WASPAA).
%
%   Inputs:
%       p:  sparsity exponent p<=1, of the L_{p,2} norm
%       A:  overcomplete dictionary (in CS jargon) of steering vectors for
%           a dense grid of K directions, with A of size MxK.
%       Y:  measurement vectors (in CS jargon), multichannel microphone signals
%           in audio/acoustics, of size MxT
%       beta:   regularization parameter
%       terminationValue:   error threshold to stop iterations in the
%                           solver
%       maxIterations:  maximum number of iterations before stopping the
%                       solver, if terminationValue is not reached
%
%   Outputs:
%       X:  KxT sparse signal solution
%       D:  demixing matrix of size KxM
%       e:  Kx1 sparse signal powers
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SPARSE_SOLVER_IRLS.M - 13/5/2019
% Archontis Politis, archontis.politis@tuni.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<6
    maxIterations = 0;
end

nChan = size(A,1); % number of channels/sensors
nDict = size(A,2); % size of steering vector dictionary
nSnap = size(Y,2); % size of temporal snapshots (e.g. STFT frames, or time-domain samples)
% initial values for amplitudes and weights
W = eye(nDict);
minEpsilon = terminationValue;
epsilon = 1;
n = 0;

if nSnap>nChan
    % dimensionality reduction using SVD
    [U, L, ~] = svd(Y);
    L_trunc = L(:,1:nChan);
    Y2 = U*L_trunc;
else
    Y2 = Y;
end
X_prev = A\Y2; % initial LS solution
while epsilon>minEpsilon
    % regularization
    gamma = beta/(1-beta) * trace(A*W*A')/nChan;
    % gamma = beta^0.15/(1-beta^0.15) * trace(A*W*A')/nChan; % Epain & Jin
    % reconstruction matrix
    D = (W*A')/(A*W*A' + gamma*eye(nChan));
    % amplitude estimates
    X_current = D * Y2;
    % norm of difference between previous and current solution
    er_k = sqrt(sum((X_current - X_prev).^2,2));
    er = sqrt(sum(er_k.^2));
    % sparse solution signal energies
    e = sum(X_current.^2,2);
    e_max = max(e);
    epsilon = e_max/nDict; % Epain & Jin
%     if er < epsilon^(1/2)/100 % Chartrand & Yin
%         epsilon = epsilon /10;
%     end
    w = (e + epsilon).^(1-p/2);
    W = diag(w);

    n = n+1;
    disp(['Iteration ' num2str(n) ': ' num2str(epsilon)])
    if maxIterations
        if n>maxIterations, break; end
    end
end
if nSnap>nChan
    D = X_current*inv(U*L_trunc);
    X = D*Y;
else
    X = X_current;
end

author: polarch

repo: Spherical-Array-Processing


text:
%FindInInterval - Find values that fall in a given interval.
%
%  The equivalent Matlab code is trivial
%
%      i = find(values(:,1)>=interval(1)&values(:,1)<=interval(2));
%      indices = [i(1),i(end)];
%
%  but becomes extremely slow when dealing with very large lists.
%  This function can dramatically speed up things whenever one needs
%  to repeatedly find values in a long list of intervals.
%
%  USAGE
%
%    indices = FindInInterval(values,interval,from)
%
%    values         values to test, sorted in ascending order
%    interval       [start,stop] pair
%    from           optional initial index (see example below)
%
%  OUTPUT
%
%    indices        indices of the first and last values that fall
%                   in the interval
%
%  EXAMPLE
%
%    This code assumes that the variable 'interval' contains a list
%    of non-overlapping intervals sorted in ascending order,
%    i.e. interval(i+1,1) >= interval(i,2).
%
%    previous = 1;
%    for i = 1:n,
%       % Find values within this window
%       j = FindInInterval(values,interval(i,:),previous);
%       previous = j(1);
%       % ... do whatever computations here ...
%    end
%
%  SEE
%
%    See also ConsolidateIntervals, SubtractIntervals, ExcludeIntervals,
%    InIntervals, Restrict, CountInIntervals, PlotIntervals.

% Copyright (C) 2004-2010 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function x2 = dfrnt ( xx, npl )

%*****************************************************************************80
%
%% DFRNT determines the derivative of a Chebyshev series.
%
%  Discussion:
%
%    This routine computes the Chebyshev series of the derivative of a
%    function whose Chebyshev series is given.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Roger Broucke.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Roger Broucke,
%    Algorithm 446:
%    Ten Subroutines for the Manipulation of Chebyshev Series,
%    Communications of the ACM,
%    October 1973, Volume 16, Number 4, pages 254-256.
%
%  Parameters:
%
%    Input, real XX(NPL), the given Chebyshev series.
%
%    Input, integer NPL, the number of terms in the
%    Chebyshev series.
%
%    Output, real X2(NPL), the Chebyshev series for the
%    derivative.
%
  x2 = zeros ( npl, 1 );

  n = npl - 1;
  xxn = xx(npl-1);
  x2(npl-1) = 2.0 * xx(npl) * n;
  x2(npl) = 0.0;

  for k = 3 : npl
    l = npl - k + 1;
    xxl = xx(l);
    x2(l) = x2(l+2) + 2.0 * xxn * l;
    xxn = xxl;
  end

  return
end

author: johannesgerer


text:
function [dose3D, A_z3D, r3D, PB3D] = getInfluence(IM, structNum, PBWeightsV)
%Create 3D dose matrix from PB weigths (PBWeightsV) and
%a list of structures to have dose computed to (structsV).
%JOD, 14 Nov 03.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

global planC

% maskSingle3D = getUniformStr(structNum);

maskSingle3D = getSurfaceExpand(structNum,0.5,1); %%%Expandewd structure

[rowV, colV, sliceV] = find3d(maskSingle3D);

numPBs = size(IM.beamlets,2);

influenceM = sparse(prod(size(maskSingle3D)),numPBs);

s = size(maskSingle3D);

doseInStruct = zeros(s);

sampleRate = IM.beamlets(structNum,1).sampleRate;

if sampleRate ~= 1
    disp('Inflating downsampled dose distribution...')
  if rem(log2(sampleRate),1) ~= 0
    error('Sample factor must (currently) be a power of 2.')
  end
  maskSample3D = getDown3Mask(maskSingle3D, sampleRate, 1);

  tmp3D = maskSample3D .* maskSingle3D;

  [rowV,colV,sliceV] = find3d(tmp3D);

  clear tmp3D

  toInterp3D = ~maskSample3D .* maskSingle3D;

  [rInterpV,cInterpV,sInterpV] = find3d(toInterp3D);

end

dose3D = zeros(s);

for PBNum = 1 : size(IM.beamlets,2)       %Loop over beamlets

  if ~isempty(IM.beamlets(structNum,PBNum).influence)

    doseV     = double(IM.beamlets(structNum,PBNum).influence);
    indV      = IM.beamlets(structNum,PBNum).indexV;
    maxVal    = IM.beamlets(structNum,PBNum).maxInfluenceVal;
    sizeParam = IM.beamlets(structNum,PBNum).fullLength;

    doseScaledV = PBWeightsV(PBNum) * (doseV * maxVal) /(2^8 -1);

    %inflate:
    doseInflateV = zeros(1,sizeParam);
    doseInflateV(indV) = doseScaledV;

    %Then put into the dose3D matrix:
    ind2V = sub2ind(size(dose3D), rowV, colV, sliceV);

    doseInStruct(ind2V) = doseInStruct(ind2V) + doseInflateV(:);

    if strcmpi(IM.params.debug(1),'y')
      ind0V =IM.beamlets(structNum,PBNum).rIndex;
      AInflate = zeros(1,sizeParam);
      rInflate = zeros(1,sizeParam);
      AInflate(ind0V) = PBWeightsV(PBNum) * double(IM.beamlets(structNum,PBNum).A_zV);
      A_z3D(ind2V) = AInflate(:); %overwrite
      rInflate(ind0V) = PBWeightsV(PBNum) * double(IM.beamlets(structNum,PBNum).r);
      r3D(ind2V) = rInflate(:); %overwrite
      PB3D(ind2V) = PBNum * ones(length(rInflate),1);
    end

  end

end

%If sub-sampled, use 3-D interpolation to fill out dose.

if sampleRate ~= 1

  %Now do 3-D interpolation:
  sizeV = size(dose3D);

  doseInterp = doseInStruct;
  %Get rid of all points which are not sampled
  [r3V,c3V,s3V] = find3d(maskSample3D);
  whereV = sub2ind(sizeV,r3V,c3V,s3V);
  doseInterp = doseInterp(whereV);
  doseInterp = reshape(doseInterp,[sizeV(1)/2^log2(sampleRate),sizeV(2)/2^log2(sampleRate),sizeV(3)]);

  rDownV = (rInterpV + 2^log2(sampleRate) - 1)/2^log2(sampleRate);
  cDownV = (cInterpV + 2^log2(sampleRate) - 1)/2^log2(sampleRate);

  fillsV = matInterp3(rDownV,cDownV,sInterpV,doseInterp);

  ind3V = sub2ind(sizeV,rInterpV,cInterpV,sInterpV);  %index for interpolated values
  doseInStruct(ind3V) = fillsV;
  disp('Finished inflating.')
end

dose3D = dose3D + doseInStruct .* [dose3D==0];  %voxel over-adding avoided


text:
function nameList = fixNames(nameList)
% Standardizes names, ie, met, rxn and long names for them as well.
% `fixNames` makes the names lowercase, removes obtuse characters and
% whitespace and replaces them with underscores, and then removes starting
% and trailing underscores. Called by `compareCbModels`, `verifyModel`, `addSEEDInfo`.
%
% USAGE:
%
%    nameList = fixNames(nameList)
%
% INPUTS:
%    nameList:      Cell array of names.
%
% OUTPUTS:
%    nameList:      Same cell array but with fixed names.
%
% Please cite:
% `Sauls, J. T., & Buescher, J. M. (2014). Assimilating genome-scale
% metabolic reconstructions with modelBorgifier. Bioinformatics
% (Oxford, England), 30(7), 1036?8`. http://doi.org/10.1093/bioinformatics/btt747
%
% ..
%    Edit the above text to modify the response to help addMetInfo
%    Last Modified by GUIDE v2.5 06-Dec-2013 14:19:28
%    This file is published under Creative Commons BY-NC-SA.
%
%    Correspondance:
%    johntsauls@gmail.com
%
%    Developed at:
%    BRAIN Aktiengesellschaft
%    Microbial Production Technologies Unit
%    Quantitative Biology and Sequencing Platform
%    Darmstaeter Str. 34-36
%    64673 Zwingenberg, Germany
%    www.brain-biotech.de

nameList = lower(nameList) ; % Fix those names! Make everything lowercase.

% Remove wierd characters.
nameList = regexprep(nameList, ' |-(?!($|\||\[))|,|:|\(|\)', '_') ;

% Remove [ and ] that aren't part of the compartment lable.
% ie ] that are not at the end of the word.
nameList = regexprep(nameList, '\](?!($|\|))', '_') ;

% And [ that don't look ahead do a letter and a ].
nameList = regexprep(nameList, '\[(?!.\]$)', '_') ;

% Consolodate resulting underscores.
nameList = regexprep(nameList, '_____|____|___|__', '_') ;

% Delete ones at beginning and end of names.
nameList = regexprep(nameList, '^_|\|_|_\||_$|_(?=\[\w\]$)', '') ;


text:
function ini_zmap_sys()
    %    This is the  ZMAP default file used for LINUX systems.
    %    It's purpose is to modify the ZmapGlobal variables as necessary
    %    to fit the system.
    report_this_filefun();

    ZG=ZmapGlobal.Data;
    % Marker sizes
    ZG.ms6 = 3;

    %set the recursion slightly, to avoid error (specialy with ploop functions)
    set(0,'RecursionLimit',750)

    set(0,'DefaultAxesFontName','Arial');
    set(0,'DefaultTextFontName','Arial');
    set(0,'DefaultAxesTickLength',[0.01 0.01]);
    set(0,'DefaultFigurePaperPositionMode','auto');

end


text:
function [weight] = create_gc_weights(xDim, yDim, zDim, xVar, yVar, zVar)
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

    % Creates a 3D normalised distributio of size dimension^3 with a variance of var.

    xDimCentre = floor(xDim / 2) + 1;
    yDimCentre = floor(yDim / 2) + 1;
    zDimCentre = floor(zDim / 2) + 1;

    weight = zeros(xDim, yDim, zDim);

    for z = 1 : zDim
        for x = 1 : xDim
            for y = 1 : yDim
               weight(x,y,z) = 1/(xVar*sqrt(2*pi))*exp((-(x - xDimCentre) ^ 2) / (2 * xVar ^ 2)) ...
                   * 1/(yVar*sqrt(2*pi))*exp((-(y - yDimCentre) ^ 2) / (2 * yVar ^ 2)) ...
                   * 1/(zVar*sqrt(2*pi))*exp((-(z - zDimCentre) ^ 2) / (2 * zVar ^ 2));
            end
        end
    end

    % ensure that it is normalised
    total = sum(sum(sum(weight)));
    weight = weight./total;

end

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
function [mMcBdepth] = calc_McBwdepth(mCatalog, fBinning, nEvents, nOverlap, nBootSample,nMinNumberevents )
% function [mMcBdepth] = calc_McBwdepth(mCatalog, fBinning, nEvents, nOverlap, nBootSample,nMinNumberevents)
% ------------------------------------------------------------------------
% Calculate Mc and b_value with depth; Mc is EMR-medianMc determined by bootstrapping
%
% Incoming variables:
% mCatalog : Earthquake catalog
% fBinnig  : Binning interval
% nEvents  : Number of events per window
% nOverlap : Number of events that overlap in windows
% nBootSample : Number of bootstrap samples
% nMinNumberevents : Minimum number of events
%
% Outgoing variables:
% mMcBdepth(:,1) : Depth
% mMcBdepth(:,2) : EMR Median Mc
% mMcBdepth(:,3) : 16-percentile
% mMcBdepth(:,4) : 84-percentile
% mMcBdepth(:,5) : b-value
% mMcBdepth(:,6) : Standard deviation of b (Shi & Bolt, 1982)
% mMcBdepth(:,7) :a-value
%
% Author: J. Woessner
% last update: 04.04.03

% Initialze
mMcBdepth = [];

% Set fix values
fMinMag = min(mCatalog(:,6));
fMaxMag = max(mCatalog(:,6));
fMinDepth = min(mCatalog(:,7));
fMaxDepth = max(mCatalog(:,7));

% Sorting by depth
[vSortDepth,vIndiceSort] = sort(mCatalog(:,7));
mCatDep = mCatalog(vIndiceSort(:,1),:);
for nStep = 1:nEvents/nOverlap:length(mCatDep(:,6))-nEvents
    mCat = mCatDep(nStep:nStep+nEvents,:);
    % Determine Mc by bootstrapping
    [vMc, vMls, fMc_org, fStdMc_org, fSkew, vPerc, fMedianMc, fMeanMc, fStdMc, v1Sigma] = calc_BstMc(mCat,fBinning,nBootSample);
    % Mean depth
    fMeanDepth = (mCatDep(nStep,7)+mCatDep(round(nStep+nEvents),7))/2;
    % B-value determination
    % Select magnitude range to calculate b-value for EMR median Mc
    vSel = mCatDep(:,6) >= fMedianMc-0.05;
    mCat2 = mCatDep(vSel,:);
    % Check for minimum number of events
    if length(mCat(:,1)) > nMinNumberevents
        try
            [fMeanMag, fBValue, fStdDev, fAValue] =  calc_bmemag(mCat2, fBinning);
            vBvalue = [fBValue fStdDev fAValue];
        catch
            vBvalue = [NaN NaN NaN];
        end
    else
        vBvalue = [NaN NaN NaN];
    end; % END of IF
    % Result matrix
    mMcBdepth = [mMcBdepth; fMeanDepth fMedianMc v1Sigma vBvalue];
end; % End of FOR nStep