text:
function varargout = process_test_parametric2( varargin )
% PROCESS_TEST_PARAMETRIC2: Parametric two-sample tests (independent).
%
% USAGE:  OutputFiles = process_test_parametric2('Run', sProcess, sInput)
%                   p = process_test_parametric2('ComputePvalues', t, df, TestType='t', TestTail='two')

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, Dimitrios Pantazis, 2008-2019

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Parametric test: Independent';
    sProcess.Category    = 'Stat2';
    sProcess.SubGroup    = 'Test';
    sProcess.Index       = 101;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Statistics';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data',  'results',  'timefreq',  'matrix'};
    sProcess.OutputTypes = {'pdata', 'presults', 'ptimefreq', 'pmatrix'};
    sProcess.nInputs     = 2;
    sProcess.nMinFiles   = 2;

    % === GENERIC EXTRACT OPTIONS
    % Label
    sProcess.options.extract_title.Comment    = 'Select data to test:';
    sProcess.options.extract_title.Type       = 'label';
    % Options
    sProcess = process_extract_values('DefineExtractOptions', sProcess);
    % DISABLE ABSOLUTE VALUE
    sProcess.options.isabs.Value = 0;
    sProcess.options.isnorm.Value = 0;
    sProcess.options.isabs.Hidden = 1;
    sProcess.options.isnorm.Hidden = 1;

    % === EXCLUDE ZERO VALUES
    sProcess.options.iszerobad.Comment    = 'Exclude the zero values from the computation';
    sProcess.options.iszerobad.Type       = 'checkbox';
    sProcess.options.iszerobad.Value      = 1;
    sProcess.options.iszerobad.InputTypes = {'timefreq', 'matrix'};
    % === OUTPUT COMMENT
    sProcess.options.Comment.Comment = 'Comment (empty=default): ';
    sProcess.options.Comment.Type    = 'text';
    sProcess.options.Comment.Value   = '';

    % === TEST: title
    sProcess.options.test_title.Comment    = '
Test statistic:';
    sProcess.options.test_title.Type       = 'label';
    % === TEST: type
    sProcess.options.test_type.Comment = {['Student''s t-test   (equal variance)        A,B~N(m,v)
' ...
                                           't = (mean(A)-mean(B)) / (Sx * sqrt(1/nA + 1/nB))
' ...
                                           'Sx = sqrt(((nA-1)*var(A) + (nB-1)*var(B)) / (nA+nB-2))
' ...
                                           'df = nA + nB - 2'], ...
                                          ['Student''s t-test   (unequal variance)        A,B~N(m,v)
', ...
                                           't = (mean(A)-mean(B)) / sqrt(var(A)/nA + var(B)/nB)
' ...
                                           'df=(vA/nA+vB/nB)2 / ((vA/nA)2/(nA-1)+(vB/nB)2/(nB-1))'], ...
                                          ['Power F-test        A~N(0,vA), B~N(0,vB)
', ...
                                           'F = (sum(A^2)/nA) / (sum(B^2)/nB)      ~F(nA,nB)'], ...
                                          ['Power F-test (unconstrained sources)
', ...
                                           'F = (sum(Ax2+Ay2+Az2)/nA) / (sum(Bx2+By2+Bz2)/nB)
', ...
                                           'Ax,Ay,Az~N(0,vA), Bx,By,Bz~N(0,vB), F~F(3*nA,3*nB)']; ...
                                          'ttest_equal', 'ttest_unequal', 'power', 'power_unconstr'};
    sProcess.options.test_type.Type    = 'radio_label';
    sProcess.options.test_type.Value   = 'ttest_equal';
    % === TAIL FOR THE TEST STATISTIC
    sProcess.options.tail.Comment  = {'One-tailed (-)', 'Two-tailed', 'One-tailed (+)', ''; ...
                                      'one-', 'two', 'one+', ''};
    sProcess.options.tail.Type     = 'radio_linelabel';
    sProcess.options.tail.Value    = 'two';
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    % === DATA SELECTION COMMENT ===
    strData = process_extract_time('GetTimeString', sProcess);
    if isfield(sProcess.options, 'freqrange') && isfield(sProcess.options.freqrange, 'Value') && iscell(sProcess.options.freqrange.Value) && (length(sProcess.options.freqrange.Value) == 3) && (length(sProcess.options.freqrange.Value{1}) == 2)
        FreqRange = sProcess.options.freqrange.Value{1};
        if (FreqRange(1) == FreqRange(2))
            strData = [strData, ' ' num2str(FreqRange(1)) 'Hz'];
        else
            strData = [strData, ' ' num2str(FreqRange(1)) '-' num2str(FreqRange(2)) 'Hz'];
        end
    end
    if isfield(sProcess.options, 'sensortypes') && isfield(sProcess.options.sensortypes, 'Value') && ~isempty(sProcess.options.sensortypes.Value)
        strData = [strData, ' ', sProcess.options.sensortypes.Value];
    end
    if isfield(sProcess.options, 'rows') && isfield(sProcess.options.rows, 'Value') && ~isempty(sProcess.options.rows.Value)
        strData = [strData, ' ', sProcess.options.rows.Value];
    end
    if ~isempty(strData)
        strData = [' [' strData ']'];
    end

    % === ABSOLUTE VALUE ===
    % Get options
    if isfield(sProcess.options, 'isabs') && isfield(sProcess.options.isabs, 'Value') && sProcess.options.isabs.Value
        isAbsolute = 1;
        strAbs = ' abs';
    elseif isfield(sProcess.options, 'isnorm') && isfield(sProcess.options.isnorm, 'Value') && sProcess.options.isnorm.Value
        isAbsolute = 1;
        strAbs = ' norm';

    else
        isAbsolute = 0;
        strAbs = '';
    end

    % === TEST COMMENT ===
    % Get test info
    if isfield(sProcess.options, 'test_type') && isfield(sProcess.options.test_type, 'Value') && ~isempty(sProcess.options.test_type.Value)
        TestType = sProcess.options.test_type.Value;
    else
        TestType = 'ttest_unequal';
    end
    if isfield(sProcess.options, 'tail') && isfield(sProcess.options.tail, 'Value') && ~isempty(sProcess.options.tail.Value)
        TestTail = sProcess.options.tail.Value;
    else
        TestTail = [];
    end
    % Documenting test to perform
    switch (TestType)
        case {'ttest_equal', 'ttest_unequal', 'ttest_paired', 'wilcoxon_paired'}   % , 'wilcoxon'
            strHypo = '          H0:(A=B)';
            switch(TestTail)
                case 'one-',   strHypo = [strHypo, ', H1:(A<B)'];
                case 'two',    strHypo = [strHypo, ', H1:(A<>B)'];
                case 'one+',   strHypo = [strHypo, ', H1:(A>B)'];
            end
        case 'signtest'
            strHypo = '         H0:(A=B), H1:(A<>B)';
        case 'ttest_onesample'
            strHypo = '          H0:(X=0)';
            switch(TestTail)
                case 'one-',   strHypo = [strHypo, ', H1:(X<0)'];
                case 'two',    strHypo = [strHypo, ', H1:(X<>0)'];
                case 'one+',   strHypo = [strHypo, ', H1:(X>0)'];
            end
        case 'ttest_baseline'
            strHypo = '          H0:(X=Baseline)';
            switch(TestTail)
                case 'one-',   strHypo = [strHypo, ', H1:(X<Baseline)'];
                case 'two',    strHypo = [strHypo, ', H1:(X<>Baseline)'];
                case 'one+',   strHypo = [strHypo, ', H1:(X>Baseline)'];
            end
        case {'power_baseline', 'power_baseline_unconstr'}
            strHypo = '          H0:(|X|=|Baseline|)';
            switch(TestTail)
                case 'one-',   strHypo = [strHypo, ', H1:(|X|<|Baseline|)'];
                case 'two',    strHypo = [strHypo, ', H1:(|X|<>|Baseline|)'];
                case 'one+',   strHypo = [strHypo, ', H1:(|X|>|Baseline|)'];
            end
        case 'chi2_onesample'
            strHypo = '          H0:(|Zi| = 0)';
        case 'chi2_onesample_unconstr'
            strHypo = '          H0:(|Zi| = 0)';
        case {'power', 'power_unconstr'}
            strHypo = '          H0:(vA=vB)';
            switch(TestTail)
                case 'one-',   strHypo = [strHypo, ', H1:(vA<vB)'];
                case 'two',    strHypo = [strHypo, ', H1:(vA<>vB)'];
                case 'one+',   strHypo = [strHypo, ', H1:(vA>vB)'];
            end
        case {'absmean', 'absmean_param'}
            strHypo = '          H0:(|mean(A)|=|mean(B)|)';
            switch(TestTail)
                case 'one-',   strHypo = [strHypo, ', H1:(|mean(A)|<|mean(B)|)'];
                case 'two',    strHypo = [strHypo, ', H1:(|mean(A)|<>|mean(B)|)'];
                case 'one+',   strHypo = [strHypo, ', H1:(|mean(A)|>|mean(B)|)'];
            end
    end

    % No comment when forcing a one-sided test
    if ismember(TestType, {'ttest_onesample'}) && isAbsolute
        strTail = '';
    elseif strcmpi(TestType, 'chi2_onesample') || strcmpi(TestType, 'chi2_onesample_unconstr')
        strTail = '';
        strAbs = '';
    elseif strcmpi(TestType, 'signtest')
        strTail = '';
    % Comment for one-tailed tests
    elseif ismember(TestTail, {'one-','one+'})
        strTail = [' ' TestTail];
    else
        strTail = '';
    end

    % === ASSEMBLING ===
    switch (TestType)
        case 'ttest_equal',      Comment = ['t-test equal'     strTail strAbs strData strHypo];
        case 'ttest_unequal',    Comment = ['t-test unequal'   strTail strAbs strData strHypo];
        case 'ttest_paired',     Comment = ['t-test paired'    strTail strAbs strData strHypo];
        case 'ttest_onesample',  Comment = ['t-test zero'      strTail strAbs strData strHypo];
        case 'ttest_baseline',   Comment = ['t-test baseline'  strTail strAbs strData strHypo];
        case 'power_baseline',   Comment = ['power test baseline'  strTail strData strHypo];
        case 'power_baseline_unconstr',   Comment = ['power test baseline unconstr' strTail strData strHypo];
        case 'signtest',         Comment = ['signtest'         strAbs strData strHypo];
        case 'wilcoxon_paired',  Comment = ['wilcoxon paired ' strTail strAbs strData strHypo];
        % case 'wilcoxon',         Comment = ['wilcoxon'       strTail strAbs strData strHypo];
        case 'power',            Comment = ['power test'       strTail strData strHypo];
        case 'power_unconstr',   Comment = ['power test unconstr' strTail strData strHypo];
        case 'chi2_onesample',   Comment = ['Chi2-test'        strTail strData strHypo];
        case 'chi2_onesample_unconstr',  Comment = ['Chi2-test unconstr' strTail strData strHypo];
        case 'absmean',          Comment = ['absmean test'     strTail strData strHypo];
        case 'absmean_param',    Comment = ['absmean test'     strTail strData strHypo];
    end
end

%% ===== RUN =====
function sOutput = Run(sProcess, sInputsA, sInputsB) %#ok<DEFNU>
    % Initialize returned variables
    sOutput = [];

    % ===== GET OPTIONS =====
    % Get generic extract options
    OPTIONS = process_extract_values('GetExtractOptions', sProcess, sInputsA(1));
    % Exclude zero values
    if isfield(sProcess.options, 'iszerobad') && isfield(sProcess.options.iszerobad, 'Value') && ~isempty(sProcess.options.iszerobad.Value)
        OPTIONS.isZeroBad = sProcess.options.iszerobad.Value;
    else
        OPTIONS.isZeroBad = 1;
    end
    % Get test type
    OPTIONS.TestType = sProcess.options.test_type.Value;
    OPTIONS.TestTail = sProcess.options.tail.Value;
    % Invalid test/tail combinations
    if ismember(OPTIONS.TestType, {'ttest_onesample'}) && OPTIONS.isAbsolute && ismember(OPTIONS.TestTail, {'two', 'one-'})
        bst_report('Warning', sProcess, [], 'Testing |X|>0: Using a positive one-tailed test (one+) instead.');
        OPTIONS.TestTail = 'one+';
    elseif ismember(OPTIONS.TestType, {'chi2_onesample', 'chi2_onesample_unconstr'}) && ismember(OPTIONS.TestTail, {'two', 'one-'})
        bst_report('Warning', sProcess, [], 'Testing |X|>0: Using a positive one-tailed test (one+) instead.');
        OPTIONS.TestTail = 'one+';
    end
    % Time-frequency: Warning if processing power
    isTfPower = false;
    if strcmpi(sInputsA(1).FileType, 'timefreq')
        TfMat = in_bst_timefreq(sInputsA(1).FileName, 0, 'Measure');
        if isequal(TfMat.Measure, 'power')
            isTfPower = true;
        end
    end
    % Get average function
    switch (OPTIONS.TestType)
        case {'ttest_equal', 'ttest_unequal', 'ttest_onesample', 'ttest_paired', 'ttest_baseline', 'absmean', 'absmean_param'}
            if isTfPower
                bst_report('Warning', sProcess, [], ['You are testing power values, while a more standard analysis is to test the magnitude (ie. sqrt(power)).' 10 ...
                    'Option #1: Recompute the time-frequency maps using the option "Measure: Magnitude".' 10 ...
                    'Option #2: Run the process "Extract > Measure from complex values", with option "Magntiude".']);
                isTfPower = false;
            end
            if OPTIONS.isAbsolute
                AvgFunction = 'norm';
                isAvgVariance = 1;
            else
                AvgFunction = 'mean';
                isAvgVariance = 1;
            end
        case {'power', 'power_baseline', 'power_baseline_unconstr', 'power_unconstr', 'chi2_onesample', 'chi2_onesample_unconstr'}
            AvgFunction = 'rms';
            isAvgVariance = 0;
    end
    isAvgWeighted = 0;
    % Baseline: Only for test against baseline
    if isfield(sProcess.options, 'baseline') && isfield(sProcess.options.baseline, 'Value') && iscell(sProcess.options.baseline.Value) && ~isempty(sProcess.options.baseline.Value) && ~isempty(sProcess.options.baseline.Value{1})
        Baseline = sProcess.options.baseline.Value{1};
    else
        Baseline = [];
    end
    % Unconstrained chi2 test: force the computation of norm
    if ismember(OPTIONS.TestType, {'chi2_onesample_unconstr', 'power_baseline_unconstr'})
        OPTIONS.isAbsolute = 1;
    end

    % ===== CHECK INPUT FILES =====
    % Make sure that file type is indentical for both sets
    if ~isempty(sInputsA) && ~isempty(sInputsB) && ~strcmpi(sInputsA(1).FileType, sInputsB(1).FileType)
        bst_report('Error', sProcess, [], 'Cannot process inputs from different types.');
        return;
    end
    % Check the number of files in input
    if (length(sInputsA) < 2) && ~strcmpi(OPTIONS.TestType, 'ttest_baseline')
        bst_report('Error', sProcess, [], 'Not enough files in input.');
        return;
    end
    % Load time vector from the first file: if same as input, discard input
    TimeVector = in_bst(sInputsA(1).FileName, 'Time');
    if ~isempty(OPTIONS.TimeWindow) && (abs(TimeVector(1) - OPTIONS.TimeWindow(1)) < 1e-4) && (abs(TimeVector(end) - OPTIONS.TimeWindow(2)) < 1e-4)
        OPTIONS.TimeWindow = [];
    end
    % Load freq range from the first file: if same as input, discard input
    if ~isempty(OPTIONS.FreqRange)
        % Load Freqs field from the input file
        TfMat = in_bst_timefreq(sInputsA(1).FileName, 0, 'Freqs');
        if iscell(TfMat.Freqs)
            BandBounds = process_tf_bands('GetBounds', TfMat.Freqs);
            FreqList = unique(BandBounds(:));
        else
            FreqList = TfMat.Freqs;
        end
        if (abs(OPTIONS.FreqRange(1) - FreqList(1)) < 1e-4) && (abs(OPTIONS.FreqRange(2) - FreqList(end)) < 1e-4)
            OPTIONS.FreqRange = [];
        end
    end

    % ===== INPUT DATA =====
    % If there is nothing special done with the files: files can be handled directly by bst_avg_files
    if isempty(OPTIONS.TimeWindow) && isempty(OPTIONS.ScoutSel) && isempty(OPTIONS.SensorTypes) && isempty(OPTIONS.Rows) && isempty(OPTIONS.FreqRange) && ~OPTIONS.isAvgTime && ~OPTIONS.isAvgRow && ~OPTIONS.isAvgFreq && ~isTfPower
        InputSetA = {sInputsA.FileName};
        if ~isempty(sInputsB)
            InputSetB = {sInputsB.FileName};
        else
            InputSetB = [];
        end
        OutputType = sInputsA(1).FileType;
    % Else: Call process "Extract values" first
    else
        % Do not concatenate the output
        OPTIONS.Dim = 0;
        % Call extraction process: FilesA
        [InputSetA, OutputType] = process_extract_values('Extract', sProcess, sInputsA, OPTIONS);
        if isempty(InputSetA)
            return;
        end
        % Read FilesB
        if ~isempty(sInputsB)
            InputSetB = process_extract_values('Extract', sProcess, sInputsB, OPTIONS);
            if isempty(InputSetB)
                return;
            end
        else
            InputSetB = [];
        end
        % Adjust time-frequency already in 'power', for power stats.
        if isTfPower
            bst_report('Info', sProcess, [], 'Data is already power values, adapting power test (not squaring again).');
            for iIn = 1:numel(InputSetA)
                [InputSetA{iIn}.TF, isError] = process_tf_measure('Compute', InputSetA{iIn}.TF, InputSetA{iIn}.Measure, 'magnitude', true);
                if isError
                    bst_report('Error', sProcess, sInputsA(1), ['Error converting time-frequency measure ' InputSetA{iIn}.Measure 'to magnitude.']);
                end
            end
            for iIn = 1:numel(InputSetB)
                [InputSetB{iIn}.TF, isError] = process_tf_measure('Compute', InputSetB{iIn}.TF, InputSetB{iIn}.Measure, 'magnitude', true);
                if isError
                    bst_report('Error', sProcess, sInputsA(1), ['Error converting time-frequency measure ' InputSetB{iIn}.Measure 'to magnitude.']);
                end
            end
        end
    end

    % === COMPUTE TEST ===
    % Branch between dependent(=paired) and independent tests
    switch (OPTIONS.TestType)

        % ===== INDEPENDENT TESTS =====
        case {'ttest_equal', 'ttest_unequal', 'absmean', 'absmean_param', 'power', 'power_unconstr'}
            % Compute mean and var for both files sets
            [StatA, MessagesA] = bst_avg_files(InputSetA, [], AvgFunction, isAvgVariance, isAvgWeighted, OPTIONS.isMatchRows, OPTIONS.isZeroBad);
            [StatB, MessagesB] = bst_avg_files(InputSetB, [], AvgFunction, isAvgVariance, isAvgWeighted, OPTIONS.isMatchRows, OPTIONS.isZeroBad);
            % Add messages to report
            if ~isempty(MessagesA)
                if isempty(StatA)
                    bst_report('Error', sProcess, sInputsA, MessagesA);
                    return;
                else
                    bst_report('Warning', sProcess, sInputsA, MessagesA);
                end
            end
            if ~isempty(MessagesB)
                if isempty(StatB)
                    bst_report('Error', sProcess, sInputsB, MessagesB);
                    return;
                else
                    bst_report('Warning', sProcess, sInputsB, MessagesB);
                end
            end
            if ~isequal(size(StatA.mean), size(StatB.mean))
                bst_report('Error', sProcess, [], 'Files A and B do not have the same number of signals or time samples.');
                return;
            end
            % Detect if the source model is unconstrained
            isUnconstrained = panel_scout('isUnconstrained', StatA);
            % Do not allow unconstrained sources without a norm
            if isUnconstrained && ~OPTIONS.isAbsolute
                bst_report('Error', sProcess, [], ['Cannot run this test on unconstrained sources:' 10 'you must compute the norm of the three orientations first.']);
                return;
            end
            % Bad channels: For recordings, keep only the channels that are good in BOTH A and B sets
            if strcmpi(sInputsA(1).FileType, 'data') && ~isempty(StatA.ChannelFlag) && ~isempty(StatB.ChannelFlag)
                ChannelFlag = StatA.ChannelFlag;
                ChannelFlag(StatB.ChannelFlag == -1) = -1;
                isGood = (ChannelFlag == 1);
            else  % case {'results', 'timefreq', 'matrix'}
                ChannelFlag = [];
                isGood = true(size(StatA.mean, 1), 1);
                isGood((StatA.nGoodSamples < 2) | (StatB.nGoodSamples < 2)) = 0;
            end

            % === COMPUTE TEST ===
            % Display progress bar
            bst_progress('start', 'Processes', 'Computing test...');
            % Get average results
            mA = StatA.mean(isGood,:,:);
            mB = StatB.mean(isGood,:,:);
            nA = repmat(StatA.nGoodSamples(isGood,:,:), [1, size(mA,2), size(mA,3)]);
            nB = repmat(StatB.nGoodSamples(isGood,:,:), [1, size(mB,2), size(mB,3)]);
            % Get variance (if needed)
            if isAvgVariance
                vA = StatA.var(isGood,:,:);
                vB = StatB.var(isGood,:,:);
                % Remove null variances
                iNull = find((vA == 0) | (vB == 0));
                vA(iNull) = eps;
                vB(iNull) = eps;
            else
                iNull = [];
            end

            % Compute test statistic
            switch (OPTIONS.TestType)
                % === T-TEST: EQUAL VARIANCE ===
                case 'ttest_equal'
                    df   = nA + nB - 2 ;
                    pvar = ((nA-1).*vA + (nB-1).*vB) ./ df;
                    tmap = (mA-mB) ./ sqrt(pvar .* (1./nA + 1./nB));
                    % Calculate p-values from t-values
                    pmap = ComputePvalues(tmap, df, 't', OPTIONS.TestTail);
                    % Units: t
                    DisplayUnits = 't';

                % === T-TEST: UNEQUAL VARIANCE ===
                case 'ttest_unequal'
                    df = (vA./nA + vB./nB).^2 ./ ...
                         ((vA./nA).^2./(nA-1) + (vB./nB).^2./(nB-1));
                    tmap = (mA-mB) ./ sqrt(vA./nA + vB./nB);
                    % Calculate p-values from t-values
                    pmap = ComputePvalues(tmap, df, 't', OPTIONS.TestTail);
                    % Units: t
                    DisplayUnits = 't';

                % ===== POWER TEST (A/B) =====
                case {'power', 'power_unconstr'}
                    % If you have xi, n normal random variables with zero mean and unit variance N(0,1), then:
                    % X1 = sum_i(xi^2) is chi-square random variable with n degrees of freedom
                    % https://en.wikipedia.org/wiki/Chi-squared_distribution
                    %
                    % If X1 and X2 are chi-square random variables with nA and nB degrees of freedom, then
                    % F = (X1/nA) / (X2/nB) is F-distributed with nA numerator degrees of freedom and nB denominator degrees of freedom.
                    % https://en.wikipedia.org/wiki/F-distribution
                    %
                    % Use case here: If A~N(0,vA) and B~N(0,vB)
                    % Then F = (sum(A^2)/nA) / (sum(B^2)/nB)  ~ F(nA,nB)
                    %
                    % Can be used for two things:
                    % 1) Testing for variance difference:   H0:(A~N(0,vA), B~N(0,vB), vA=vB)
                    %    => Samples must be zero-mean (mA=0, mB=0) so probably normalized before testing
                    % 2) Testing for power difference:      H0:(A~N(0,1) and B~(0,1))
                    %    => Samples must be normalized before testing

                    % The output of bst_avg_files is the RMS of A and B:  mA = sqrt(sum(A^2)/nA)
                    % F statistic = (sum(A^2)/nA) / (sum(B^2)/nB)
                    %             = (mA^2) / (mB^2)
                    tmap = mA.^2 ./ mB.^2;
                    % Degrees of freedom
                    if strcmpi(OPTIONS.TestType, 'power_unconstr')
                        df = {3*nA,3*nB};
                    else
                        df = {nA,nB};
                    end
                    % Calculate p-values from t-values
                    pmap = ComputePvalues(tmap, df, 'F', OPTIONS.TestTail);
                    % Units: F
                    DisplayUnits = 'F';

                % === ABSOLUTE MEAN TEST ===
                case 'absmean_param'
                    % EXPLANATIONS:
                    %   Assume the individual samples xi follow a normal distribution with mean "m" and variance "s^2":  X ~ N(m,s^2)
                    %   Then mean(x) is also normal with mean m and variance sm^2 = s^2/N:    mean(X) ~ N(m,s^2/N)
                    %
                    %   When we apply abs(mean(x)), we are folding this normal distribution to make it positive.
                    %   Details are discussed here: https://en.wikipedia.org/wiki/Folded_normal_distribution
                    %   If y=|x|, with x~N(m,sm^2), then y has a new distribution with mean my and variance sy^2:
                    %      my = sm*sqrt(2/pi)*exp(-m^2/(2*sm^2)) - m*erf(-m/sqrt(2*sm^2))
                    %         = s/sqrt(N)*sqrt(2/pi)*exp(-m^2/(2*s^2/N)) - m*erf(-m/sqrt(2*s^2/N))
                    %      sy^2 = m^2 + sm^2 - my^2
                    %           = m^2 + s^2/N - my^2
                    %
                    % RESTRICTIONS
                    %   - A and B are normally distributed (same as t-test assumptions)
                    %   - Cannot be applied if an absolute has been applied already, we need the original values

                    % Test to check that there was no abs already applied, we need the relative values
                    if all(mA(:) > 0) && all(mB(:) > 0)
                        bst_report('Error', sProcess, [], ['This test is designed for values that are positive and negative.' 10 'It cannot be applied to values for which we have already discarded the sign.' 10 'If all your measures you are testing are always strictly positive, then use a Student t-test.']);
                        return;
                    end

                    % Mean of: abs(mean(A))-abs(mean(B))
                    % mAabs = sA/sqrt(N)*sqrt(2/pi)*exp(-mA^2/(2*sA^2/nA)) - mA*erf(-mA/sqrt(2*sA^2/nA))
                    %       = sqrt(vA./nA.*(2/pi)) .* exp(-mA.^2/(2.*vA./nA)) - mA*erf(-mA./sqrt(2.*vA./nA))
                    mAabs = sqrt(vA./nA.*(2/pi)) .* exp(-mA.^2./(2.*vA./nA)) - mA.*erf(-mA./sqrt(2.*vA./nA));
                    mBabs = sqrt(vB./nB.*(2/pi)) .* exp(-mB.^2./(2.*vB./nB)) - mB.*erf(-mB./sqrt(2.*vB./nB));
                    mAB = mAabs - mBabs;
                    % Variance of: abs(mean(A))-abs(mean(B))
                    vAabs = mA.^2 + vA./nA - mAabs.^2;
                    vBabs = mB.^2 + vB./nB - mBabs.^2;
                    sdAB = sqrt(vAabs + vBabs);
                    S = (abs(mA) - abs(mB) - mAB) ./ sdAB;  %S should be zero mean, unit variance under the null hypothesis

                    % [H,P] = ztest(S,0,1);   m = 0; sigma = 1;
                    % zval = (S - m) ./ (sigma ./ sqrt(length(S)));
                    tmap = S .* sqrt(length(S));
                    % Two-tailed test
                    pmap = 2 * (1/2 * erfc(-1 * -abs(tmap) / sqrt(2)));     % 2 * normcdf(-abs(zval),0,1);
                    % No need to recompute the values on the fly
                    df = [];
                    % Units: z
                    DisplayUnits = 'z';

                otherwise
                    error('Not supported yet');
            end
            % Remove values with null variances
            if ~isempty(iNull)
                tmap(iNull) = 0;
                pmap(iNull) = 1;
            end

        % ===== PAIRED/ONE-SAMPLE TESTS =====
        case {'ttest_paired', 'ttest_onesample', 'ttest_baseline', 'chi2_onesample', 'chi2_onesample_unconstr', 'power_baseline', 'power_baseline_unconstr'}
            % Number of samples must be equal
            if (length(sInputsA) ~= length(sInputsB)) && ismember(OPTIONS.TestType, {'ttest_paired'})
                bst_report('Error', sProcess, [], 'For a paired test, the number of files must be the same in the two groups.');
                return;
            end
            % Compute the mean and variance of (samples A - samples B)
            [StatA, MessagesA] = bst_avg_files(InputSetA, InputSetB, AvgFunction, isAvgVariance, isAvgWeighted, OPTIONS.isMatchRows, OPTIONS.isZeroBad);
            % Add messages to report
            if ~isempty(MessagesA)
                if isempty(StatA)
                    bst_report('Error', sProcess, [], MessagesA);
                    return;
                else
                    bst_report('Warning', sProcess, [], MessagesA);
                end
            end
            % Display progress bar
            bst_progress('start', 'Processes', 'Computing test...');
            % Bad channels and other properties
            switch lower(sInputsA(1).FileType)
                case {'data', 'pdata'}
                    ChannelFlag = StatA.ChannelFlag;
                    isGood = (ChannelFlag == 1);
                case {'results', 'timefreq', 'matrix', 'presults', 'ptimefreq', 'pmatrix'}
                    ChannelFlag = [];
                    isGood = true(size(StatA.mean, 1), 1);
                    isGood(StatA.nGoodSamples < 2) = -1;
            end

            % === COMPUTE TEST ===
            % Display progress bar
            bst_progress('start', 'Processes', 'Computing test...');
            % Get results
            mean_diff = StatA.mean(isGood,:,:);
            nA = repmat(StatA.nGoodSamples(isGood,:,:), [1, size(mean_diff,2), size(mean_diff,3)]);
            nB = [];
            % Get variance (if needed)
            if isAvgVariance
                std_diff = sqrt(StatA.var(isGood,:,:));
                % Remove null variances
                iNull = find(std_diff == 0);
                std_diff(iNull) = eps;
            else
                iNull = [];
            end

            % Get pre-stimulus baseline (for tests vs baseline)
            if ismember(OPTIONS.TestType, {'ttest_baseline', 'power_baseline', 'power_baseline_unconstr'})
                if ~isempty(Baseline)
                    % Get baseline bounds
                    iBaseline = bst_closest(Baseline, StatA.Time);
                    if (iBaseline(1) == iBaseline(2))
                        bst_report('Error', sProcess, [], 'The baseline must be included in the time window on which you run the test.');
                        return;
                    end
                    iBaseline = iBaseline(1):iBaseline(2);
                else
                    bst_report('Warning', sProcess, [], 'Baseline is not defined, using the entire time definition.');
                    iBaseline = 1:length(TimeVector);
                end
            end

            % Compute test statistic
            switch (OPTIONS.TestType)
                case {'ttest_paired', 'ttest_onesample'}
                    % Compute t-test
                    tmap = mean_diff ./ std_diff .* sqrt(nA);
                    df = nA - 1;
                    % Test if the statistics make sense
                    if all(tmap(:) == 0)
                        bst_report('Error', sProcess, [], 'The T-statistics is zero for all the tests.');
                        return;
                    end
                    % Calculate p-values from t-values
                    pmap = ComputePvalues(tmap, df, 't', OPTIONS.TestTail);
                    % Units: t
                    DisplayUnits = 't';

                case {'chi2_onesample', 'chi2_onesample_unconstr'}
                    % https://en.wikipedia.org/wiki/Chi-squared_distribution
                    % =>  If Zi~N(0,1) i=1..n  =>  Q=sum(Zi^2) ~ Chi2(n)
                    % Variable "mean_diff" contains RMS(data)=sqrt(sum(data.^2)/n)
                    % =>  If data is ~N(0,1)   =>  (mean_diff^2 * n) ~ Chi2(n)
                    tmap = mean_diff .^ 2 .* nA;
                    % Number of degrees of freedom
                    if strcmpi(OPTIONS.TestType, 'chi2_onesample_unconstr')
                        df = 3 * nA;
                    else
                        df = nA;
                    end
                    % Calculate p-values from F-values
                    pmap = ComputePvalues(tmap, df, 'chi2', OPTIONS.TestTail);
                    % Units: t
                    DisplayUnits = 'chi2';

                case 'ttest_baseline'
                    % TEST:  Y = mean_trials(X)
                    %        t = (Y - mean_time(Y(baseline)) / std_time(Y(baseline)))
                    % Compute variance over baseline (pre-stim interval)
                    meanBaseline = mean(mean_diff(:,iBaseline,:), 2);
                    stdBaseline = std(mean_diff(:,iBaseline,:), 0, 2);
                    % Remove null variances
                    iNull = find(stdBaseline == 0);
                    stdBaseline(iNull) = eps;
                    % Compute t-statistics (formula from wikipedia)
                    tmap = bst_bsxfun(@minus, mean_diff, meanBaseline);
                    tmap = bst_bsxfun(@rdivide, tmap, stdBaseline);
                    df = repmat(length(iBaseline) - 1, size(tmap));
                    % Calculate p-values from t-values
                    pmap = ComputePvalues(tmap, df, 't', OPTIONS.TestTail);
                    % Units: t
                    DisplayUnits = 't';

                case {'power_baseline', 'power_baseline_unconstr'}
                    % TEST:  Y = sum_trials(X^2)
                    %        F = Y / mean_time(Y(baseline))     F~F(Ntrials,Ntrials)
                    %
                    % The output of bst_avg_files is the RMS of X:  data = sqrt(sum_trials(X^2)/Ntrials)
                    % => Y = data^2 * Ntrials
                    % => F = data^2 / mean(data^2(baseline))

                    % Square the RMS
                    data = mean_diff .^ 2;
                    % Compute mean over baseline
                    meanBaseline = mean(data(:,iBaseline,:), 2);
                    % Remove null denominators
                    iNull = find(meanBaseline == 0);
                    meanBaseline(iNull) = eps;
                    % Compute F statistic
                    tmap = bst_bsxfun(@rdivide, data, meanBaseline);
                    % Degrees of freedom
                    if strcmpi(OPTIONS.TestType, 'power_baseline_unconstr')
                        df = {3*nA,3*nA};
                    else
                        df = {nA,nA};
                    end
                    % Calculate p-values from F-values
                    pmap = ComputePvalues(tmap, df, 'F', OPTIONS.TestTail);
                    % No need to recompute the values on the fly
                    df = [];
                    % Units: t
                    DisplayUnits = 'F';

                otherwise
                    error('Not supported yet');
            end
            % Remove values with null variances
            if ~isempty(iNull)
                tmap(iNull) = 0;
                pmap(iNull) = 1;
            end
    end

    % Return full matrices
    if all(isGood)
        tmap_full = tmap;
        pmap_full = pmap;
        df_full   = df;
        nA_full   = nA;
        nB_full   = nB;
    else
        tmap_full = zeros(size(StatA.mean));
        tmap_full(isGood,:,:) = tmap;
        if ~isempty(df)
            df_full = zeros(size(StatA.mean));
            df_full(isGood,:,:) = df;
        else
            df_full = [];
        end
        if ~isempty(pmap)
            pmap_full = ones(size(StatA.mean));
            pmap_full(isGood,:,:) = pmap;
        else
            pmap_full = [];
        end
        if ~isempty(nA)
            nA_full = zeros(1,size(StatA.mean,1));
            nA_full(isGood) = nA(1:size(nA,1));
        else
            nA_full = [];
        end
        if ~isempty(nB)
            nB_full = zeros(1,size(StatA.mean,1));
            nB_full(isGood) = nB(1:size(nB,1));
        else
            nB_full = [];
        end
    end

    % === CONVERT BACK MATRIX => DATA ===
    % If processing recordings with only some sensor types selected
    if strcmpi(sInputsA(1).FileType, 'data') && strcmpi(OutputType, 'matrix') && ~isempty(OPTIONS.SensorTypes) && ~OPTIONS.isAvgTime && ~OPTIONS.isAvgRow && ~OPTIONS.isAvgFreq
        % Get the list of selected sensors
        dataTypes = strtrim(str_split(OPTIONS.SensorTypes, ',;'));
        % If only major data types were selected: save results in "data" format
        if ~isempty(dataTypes) && all(ismember(dataTypes, {'MEG','EEG','MEG MAG''MEG GRAD','MEG GRAD2','MEG GRAD3','SEEG','ECOG','NIRS'}))
            % Load channel file
            ChannelMat = in_bst_channel(sInputsA(1).ChannelFile);
            % Find channel names in the output row names
            [tmp,iChan,iRow] = intersect({ChannelMat.Channel.Name}, StatA.RowNames);
            % Convert output data matrices
            tmap_tmp = zeros(length(ChannelMat.Channel), size(tmap_full,2), size(tmap_full,3));
            tmap_tmp(iChan,:,:) = tmap_full(iRow,:,:);
            tmap_full = tmap_tmp;
            if ~isempty(pmap_full)
                pmap_tmp = zeros(size(tmap_tmp));
                pmap_tmp(iChan,:,:) = pmap_full(iRow,:,:);
                pmap_full = pmap_tmp;
            end
            if ~isempty(df_full)
                df_tmp = zeros(size(tmap_tmp));
                df_tmp(iChan,:,:) = df_full(iRow,:,:);
                df_full = df_tmp;
            end
            % New channel flag
            tmpChannelFlag = -1 .* ones(length(ChannelMat.Channel), 1);
            if ~isempty(ChannelFlag) && (length(ChannelFlag) == length(iChan))
                tmpChannelFlag(iChan) = ChannelFlag(iRow);
            else
                tmpChannelFlag(iChan) = 1;
            end
            ChannelFlag = tmpChannelFlag;
            % Convert Stat structure
            OutputType = 'data';
            StatA.RowNames = [];
        end
    end

    % === OUTPUT STRUCTURE ===
    % Initialize output structure
    sOutput = db_template('statmat');
    sOutput.pmap         = pmap_full;
    sOutput.tmap         = tmap_full;
    sOutput.df           = df_full;
    sOutput.Correction   = 'no';
    sOutput.Type         = OutputType;
    sOutput.ChannelFlag  = ChannelFlag;
    sOutput.Time         = StatA.Time;
    sOutput.ColormapType = 'stat2';
    sOutput.DisplayUnits = DisplayUnits;
    sOutput.nComponents  = StatA.nComponents;
    sOutput.GridAtlas    = StatA.GridAtlas;
    sOutput.Freqs        = StatA.Freqs;
    sOutput.TFmask       = StatA.TFmask;
    % Row names
    if isfield(StatA, 'RowNames') && ~isempty(StatA.RowNames)
        if strcmpi(OutputType, 'matrix')
            sOutput.Description = StatA.RowNames;
        elseif strcmpi(OutputType, 'timefreq')
            sOutput.RowNames = StatA.RowNames;
        end
    end
    % Save options
    sOutput.Options = OPTIONS;
    % Save the number of good samples used for both sets
    sOutput.Options.nGoodSamplesA = nA_full;
    sOutput.Options.nGoodSamplesB = nB_full;
end

%% ===== COMPUTE P-VALUES ====
function p = ComputePvalues(t, df, TestDistrib, TestTail)
    % Default: two-tailed tests
    if (nargin < 4) || isempty(TestTail)
        TestTail = 'two';
    end
    % Default: F-distribution
    if (nargin < 3) || isempty(TestDistrib)
        TestDistrib = 'f';
    end
    % Nothing to test
    if strcmpi(TestTail, 'no')
        p = zeros(size(t));
        return;
    end

    % Different distributions
    switch lower(TestDistrib)
        % === T-TEST ===
        case 't'
            % Calculate p-values from t-values
            switch (TestTail)
                case 'one-'
                    % Inferior one-tailed t-test:   p = tcdf(t, df);
                    % Equivalent without the statistics toolbox (FieldTrip formula)
                    p = 0.5 .* ( 1 + sign(t) .* betainc( t.^2 ./ (df + t.^2), 0.5, 0.5.*df ) );
                case 'two'
                    % Two-tailed t-test:     p = 2 * (1 - tcdf(abs(t),df));
                    % Equivalent without the statistics toolbox
                    p = betainc( df ./ (df + t .^ 2), df./2, 0.5);
                    % FieldTrip equivalent: p2 = 1 - betainc( t.^2 ./ (df + t.^2), 0.5, 0.5.*df );
                case 'one+'
                    % Superior one-tailed t-test:    p = 1 - tcdf(t, df);
                    % Equivalent without the statistics toolbox (FieldTrip formula)
                    p = 0.5 .* ( 1 - sign(t) .* betainc( t.^2 ./ (df + t.^2), 0.5, 0.5.*df ) );
            end

        % === F-TEST ===
        case 'f'
            v1 = df{1};
            v2 = df{2};
            % Evaluate for which values we can compute something
            k = ((t > 0) & ~isinf(t) & (v1 > 0) & (v2 > 0));
            % Initialize returned p-values
            p = ones(size(t));
            % Calculate p-values from F-values
            switch (TestTail)
                case 'one-'
                    % Inferior one-tailed F-test
                    % p = fcdf(t, v1, v2);
                    p(k) = 1 - betainc(v2(k)./(v2(k) + v1(k).*t(k)), v2(k)./2, v1(k)./2);
                case 'two'
                    % Two tailed F-test
                    % p = 2*min(fcdf(F,df1,df2),fpval(F,df1,df2))
                    p(k) = 2 * min(...
                            1 - betainc(v2(k)./(v2(k) + v1(k).*t(k)), v2(k)./2, v1(k)./2), ...
                            1 - betainc(v1(k)./(v1(k) + v2(k)./t(k)), v1(k)./2, v2(k)./2));
                case 'one+'
                    % Superior one-tailed F-test
                    % p = fpval(t, v1, v2);
                    %   = fcdf(1/t, v2, v1);
                    p(k) = 1 - betainc(v1(k)./(v1(k) + v2(k)./t(k)), v1(k)./2, v2(k)./2);
            end

        % === CHI2-TEST ===
        case 'chi2'
            % Calculate p-values from Chi2-values
            %   chi2cdf(x,n) = gammainc(t/2, n/2)
            switch (TestTail)
                case 'one-'
                    % Inferior one-tailed Chi2-test:    p = gammainc(t./2, df./2);
                    error('Not relevant.');
                case 'two'
                    % Two-tailed Chi2-test
                    error('Not relevant.');
                case 'one+'
                    % Superior one-tailed Chi2-test:    p = 1 - gammainc(t./2, df./2);
                    p = 1 - gammainc(t./2, df./2);
            end
    end
end

author: brainstorm-tools


text:
function [ title, key, totcrd, ptrcrd, indcrd, valcrd, rhscrd, mxtype, ...
  nrow, ncol, nnzero, neltvl, ptrfmt, indfmt, valfmt, rhsfmt, rhstyp, ...
  nrhs, nrhsix, colptr, rowind, values, rhsval, rhsptr, rhsind, guess, ...
  exact ] = hb_file_read ( input_unit )

%*****************************************************************************80
%
%% HB_FILE_READ reads an HB file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Iain Duff, Roger Grimes, John Lewis,
%    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
%    October 1992.
%
%  Parameters:
%
%    Input, integer INPUT_UNIT, the unit from data is read.
%
%    Output, character ( len = 72 ) TITLE, a title for the matrix.
%
%    Output, character ( len = 8 ) KEY, an identifier for the matrix.
%
%    Output, integer TOTCRD, the total number of lines of data.
%
%    Output, integer PTRCRD, the number of input lines for pointers.
%
%    Output, integer INDCRD, the number of input lines for row indices.
%
%    Output, integer VALCRD, the number of input lines for numerical values.
%
%    Output, integer RHSCRD, the number of input lines for right hand sides.
%
%    Output, character ( len = 3 ) MXTYPE, the matrix type.
%    First character is R for Real, C for complex, P for pattern only.
%    Second character is S for symmetric, U for unsymmetric, H for
%      Hermitian, Z for skew symmetric, R for rectangular.
%    Third character is A for assembled and E for unassembled
%      finite element matrices.
%
%    Output, integer NROW, the number of rows or variables.
%
%    Output, integer NCOL, the number of columns or elements.
%
%    Output, integer NNZERO.  In the case of assembled sparse matrices,
%    this is the number of nonzeroes.  In the case of unassembled finite
%    element matrices, in which the right hand side vectors are also
%    stored as unassembled finite element vectors, this is the total
%    number of entries in a single unassembled right hand side vector.
%
%    Output, integer NELTVL, the number of finite element matrix entries,
%    set to 0 in the case of assembled matrices.
%
%    Output, character ( len = 16 ) PTRFMT, the format for reading pointers.
%
%    Output, character ( len = 16 ) INDFMT, the format for reading indices.
%
%    Output, character ( len = 20 ) VALFMT, the format for reading values.
%
%    Output, character ( len = 20 ) RHSFMT, the format for reading values
%    of the right hand side.
%
%    Output, character ( len = 3 ) RHSTYP, the right hand side type.
%    First character is F for full storage or M for same as matrix.
%    Second character is G if starting "guess" vectors are supplied.
%    Third character is X if exact solution vectors are supplied.
%    Ignored if NRHS = 0.
%
%    Output, integer NRHS, the number of right hand sides.
%
%    Output, integer NRHSIX, the number of row indices (set to 0
%    in the case of unassembled matrices.)  Ignored if NRHS = 0.
%
%    Output, integer COLPTR(NCOL+1), COLPTR(I) points to the location of
%    the first entry of column I in the sparse matrix structure.
%
%    Output, integer ROWIND(NNZERO) or ROWIND(NELTVL), the row index of
%    each item.
%
%    Output, real VALUES(NNZERO) or VALUES(NELTVL), the nonzero values
%    of the matrix.
%
%    If RHSTYP(1:1) == 'F':
%
%      Output, integer RHSPTR(*), is not used.
%
%      Output, integer RHSIND(*), is not used.
%
%      Output, real RHSVAL(NROW,NRHS), contains NRHS dense right hand
%      side vectors.
%
%    If RHSTYP(1:1) = 'M' and MXTYPE(3:3) = 'A':
%
%      Output, integer RHSPTR(NRHS+1), RHSPTR(I) points to the location of
%      the first entry of right hand side I in the sparse right hand
%      side vector.
%
%      Output, integer RHSIND(NRHSIX), indicates, for each entry of
%      RHSVAL, the corresponding row index.
%
%      Output, real RHSVAL(NRHSIX), contains the value of the right hand
%      side entries.
%
%    If RHSTYP(1:1) = 'M' and MXTYPE(3:3) = 'E':
%
%      Output, integer RHSPTR(*), is not used.
%
%      Output, integer RHSIND(*), is not used.
%
%      Output, real RHSVAL(NNZERO,NRHS), contains NRHS unassembled
%      finite element vector right hand sides.
%
%    Output, real GUESS(NROW,NRHS), the starting guess vectors.
%
%    Output, real EXACT(NROW,NRHS), the exact solution vectors.
%

%
%  Read the header block.
%
  [ title, key, totcrd, ptrcrd, indcrd, valcrd, rhscrd, mxtype, ...
    nrow, ncol, nnzero, neltvl, ptrfmt, indfmt, valfmt, rhsfmt, rhstyp, ...
    nrhs, nrhsix ] = hb_header_read ( input_unit );
%
%  Read the matrix structure.
%
  [ colptr, rowind ] = hb_structure_read ( input_unit, ncol, mxtype, ...
    nnzero, neltvl, ptrcrd, ptrfmt, indcrd, indfmt );
%
%  Read the matrix values.
%
  values = hb_values_read ( input_unit, valcrd, mxtype, nnzero, neltvl, ...
    valfmt );
%
%  Read the right hand sides.
%
  [ rhsval, rhsptr, rhsind ] = hb_rhs_read ( input_unit, nrow, nnzero, ...
    nrhs, nrhsix, rhscrd, ptrfmt, indfmt, rhsfmt, mxtype, rhstyp );
%
%  Read the starting guesses.
%
  guess = hb_guess_read ( input_unit, nrow, nrhs, rhscrd, rhsfmt, rhstyp );
%
%  Read the exact solutions.
%
  exact = hb_exact_read ( input_unit, nrow, nrhs, rhscrd, rhsfmt, rhstyp );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function T = tn_vector()

    T.expectation = @TN_vector_expectation;
    function [exp] = TN_vector_expectation(mus, taus)

        sigmas = 1.0 ./ sqrt(taus);
        x = - mus ./ sigmas;
        %lambdax = norm.pdf(x)/(0.5*erfc(x/sqrt(2)));
        lambdax = pdf('Normal',x,0,1) ./ (0.5 * erfc(x/sqrt(2)));
        exp = mus + sigmas .* lambdax;

        for i = 1 : length(exp)
            v = exp(i);
            mu = mus(i);
            tau = taus(i);
            sigma = sigmas(i);

            if mu < -30 * sigma
                v = 1/(abs(mu)*tau);
            else
                % do nothing;
            end

            if isnan(v)
                v = 0;
            end

            exp(i) = v;
        end

        exp(isinf(exp)) = 0;

    end

    T.variance = @TN_vector_variance;
    function [var] = TN_vector_variance(mus, taus)

        sigmas = 1.0 ./ sqrt(taus);
        x = - mus ./ sigmas;
        %lambdax = norm.pdf(x)/(0.5*erfc(x/math.sqrt(2)));
        lambdax = pdf('Normal',x,0,1) ./ (0.5 * erfc(x/sqrt(2)));
        deltax = lambdax .* (lambdax-x);
        var = sigmas.^2 .* ( 1 - deltax );

        for i = 1 : length(var)
            v = var(i);
            mu = mus(i);
            tau = taus(i);
            sigma = sigmas(i);

            if mu < -30 * sigma
                v = (1/(abs(mu)*tau)).^2;
            else
                % do nothing;
            end

            if isnan(v)
                v = 0;
            end

            var(i) = v;
        end

        var(isinf(var)) = 0;
    end
end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function kl = cross_entropy(p, q, symmetric)
% CROSS_ENTROPY Compute the Kullback-Leibler divergence between two discrete prob. distributions
% kl = cross_entropy(p, q, symmetric)
%
% If symmetric = 1, we compute the symmetric version. Default: symmetric = 0;

tiny = exp(-700);
if nargin < 3, symmetric = 0; end
p = p(:);
q = q(:);
if symmetric
  kl  = (sum(p .* log((p+tiny)./(q+tiny))) + sum(q .* log((q+tiny)./(p+tiny))))/2;
else
  kl  = sum(p .* log((p+tiny)./(q+tiny)));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [net, info] = train_DCFNet(varargin)
%Training DCFNet
training_env();
opts.networkType = 12;
opts.lossType = 1;
opts.inputSize = 125;
opts.padding = 2; % padding size for crop
opts.gpus = 1;
[opts] = vl_argparse(opts, varargin) ;

opts.expDir = fullfile('..', 'data', sprintf('DCFNet-net-%d-%d-%1.1f', opts.networkType, opts.inputSize, opts.padding)) ;
opts.imdbPath = fullfile('..', 'data', sprintf('imdb-vid2015-DCFNet-%d-%1.1f', opts.inputSize, opts.padding), 'imdb.mat');

trainOpts.momentum = 0.9;
trainOpts.weightDecay = 0.0005;
trainOpts.numEpochs = 50;
trainOpts.learningRate = logspace(-2, -5, trainOpts.numEpochs) ; % from SiameseFC
trainOpts.batchSize = 32;
trainOpts.gpus = [opts.gpus]; %only support single gpu
opts.train = trainOpts;

if ~isfield(opts.train, 'gpus')
    opts.train.gpus = [];
elseif numel(opts.train.gpus) ~=0
    gpuDevice(opts.train.gpus);
end

% --------------------------------------------------------------------
%                                                 Prepare net and data
% --------------------------------------------------------------------
net = init_DCFNet('networkType', opts.networkType,...
    'inputSize', opts.inputSize, 'padding', opts.padding);

if exist(opts.imdbPath, 'file')
    imdb = load(opts.imdbPath) ;
else
    imdb = getImdbDCFNet('output_size', opts.inputSize,...
        'padding', opts.padding) ;
    if ~exist(fileparts(opts.imdbPath), 'dir')
        mkdir(fileparts(opts.imdbPath));
    end
    save(opts.imdbPath, '-v7.3', '-struct', 'imdb') ;
end

% --------------------------------------------------------------------
%                                                                Train
% --------------------------------------------------------------------
[net, info] = cnn_train_dag(net, imdb, getBatch(opts), ...
    'expDir', opts.expDir, opts.train, 'val', find(imdb.images.set == 2)) ;

transition_net('expDir', opts.expDir);
transition_net_multi('expDir', opts.expDir);
end

% --------------------------------------------------------------------
function fn = getBatch(opts)
% --------------------------------------------------------------------
bopts = struct('numGpus', numel(opts.train.gpus), 'batchSize', opts.train.batchSize, 'maxStep', 10) ;
fn = @(x,y) getDagNNBatch(bopts,x,y) ;
end

% --------------------------------------------------------------------
function inputs = getDagNNBatch(opts, imdb, batch)
% --------------------------------------------------------------------
rand_next = randi([1, opts.maxStep], size(batch));
rand_next = min(rand_next, imdb.images.up_index(batch));

if opts.numGpus > 0
    target = gpuArray(single(imdb.images.images(:,:,:,batch)));
    search = gpuArray(single(imdb.images.images(:,:,:,batch+rand_next)));
else
    target = single(imdb.images.images(:,:,:,batch));
    search = single(imdb.images.images(:,:,:,batch+rand_next));
end
target = bsxfun(@minus, target, imdb.images.data_mean);
search = bsxfun(@minus, search, imdb.images.data_mean);
inputs = {'target', target, 'search', search} ;
end

author: foolwood

repo: DCFNet

save_path: github-repos/MATLAB/foolwood-DCFNet


text:
function geometry_test016 ( )

%*****************************************************************************80
%
%% GEOMETRY_TEST016 tests CIRCLE_IMP_POINTS_2D, POLYGON_AREA_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
  center(1:2,1) =  [ 5.0; -2.0 ];
  r = 2.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'GEOMETRY_TEST016\n' );
  fprintf ( 1, '  CIRCLE_IMP_POINTS_2D gets points on a circle;\n' );
  fprintf ( 1, '  POLYGON_AREA_2D finds the area of a polygon.\n' );

  circle_imp_print_2d ( r, center, '  The implicit circle:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The area = %f\n', pi * r * r );

  n = 8;

  p = circle_imp_points_2d ( r, center, n );

  r8mat_print ( 2, n, p, '  Sample results:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  For any N, the sampled points define a polygon\n' );
  fprintf ( 1, '  whose area approximates the circle area.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  N      Area\n' );
  fprintf ( 1, '\n' );

  for n = 3 : 24

    p = circle_imp_points_2d ( r, center, n );
    result = polygon_area_2d ( n, p );
    fprintf ( 1, '  %6d  %12f\n', n, result );

  end

  return
end


text:
function k = cmpndKernDiagCompute(kern, x)

% CMPNDKERNDIAGCOMPUTE Compute diagonal of CMPND kernel.
% FORMAT
% DESC computes the diagonal of the kernel matrix for the compound kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : a vector containing the diagonal of the kernel matrix
% computed at the given points.
%
% SEEALSO : cmpndKernParamInit, kernDiagCompute, kernCreate, cmpndKernCompute
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% KERN

i = 1;
if ~isempty(kern.comp{i}.index)
  % only part of the data is involved with the kernel.
  k  = kernDiagCompute(kern.comp{i}, x(:, kern.comp{i}.index));
else
  % all the data is involved with the kernel.
  k  = kernDiagCompute(kern.comp{i}, x);
end
for i = 2:length(kern.comp)
  if ~isempty(kern.comp{i}.index)
    % only part of the data is involved with the kernel.
    k  = k + kernDiagCompute(kern.comp{i}, x(:, kern.comp{i}.index));
  else
    % all the data is involved with the kernel.
    k  = k + kernDiagCompute(kern.comp{i}, x);
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function x = slopeinit(xs,X)
%SLOPEINIT    Initialization of slope expansion
%
%  x = slopeinit(xs,X)
%
%This statement with interval quantities xs and X followed by
%  a function evaluation
%
%  y = f(x)
%
%using arithmetic operations and standard functions (overloaded by
%slope operators) implies 3 assertions:
%
%  1)  f(xs) in y.c
%  2)  f(X)  in y.r
%  3)  For all xp in X there exists s in y.s with  f(xp) = f(xs) + s*(xp-xs) .
%
%The expansion point xs and expansion interval X must be real quantities,
%  non-interval input data is converted into interval data to ensure
%  correctness of 1..3).
%
%The call
%
%  slopeinit
%
%without input parameters initializes the slope expansion to xs = X = [].
%
%The slope of arrays is stored in the 'next dimension'. So y.s is
%  3-dimensional for a slope row vector y. Since Matlab does not
%  support multi-dimensional sparse arrays, y.s is not accessible
%  if y has more than one column.
%
%As a simple example of slopes
%
%  xs = [ -3 ; 7.15 ] ;  X = intval(' [ -3.1,-2.9 ]  [ 7,7.1 ] ') ;
%  u = slopeinit( xs , X );
%
%initializes the slope package to have two dependent variables. The
%expansion point is xs=[-3;7.15] with expansion interval X.
%
%For a function f from R^2 to R, y=f(u) computed by overloaded slope
%operators has access y.c, y.r and y.s with the properties listed above.
%Internally, slopes are computed and stored according to
%
%  S.M. Rump: Expansion and Estimation of the Range of Nonlinear Functions,
%    Math. Comp. 65(216), p. 1503-1512 (1996).
%
%In contrast to the usual definition of slopes (see Neumaier's book),
%
%  - rather than one n-dimensional slope, n one-dimensional slopes
%      are calculated treating the function
%        f(X_1,...,X_k-1,z,xs_k+1,...,xs_n)
%      as a function in one variable z.
%  - intersections for the two definitions of multiplication and division
%      sharpens results (note this is not possible for n-dimensional slopes)
%  - intersections of the computed range y.r with y.c+y.s(X-xs) sharpens
%      results
%  - the sharp formulas for slopes of convex and concave functions are
%      used (for details, see the paper).
%
%Slope are useful for verified inclusion of clustered or multiple zeros.
%Sometimes, the range estimation is better than naive interval arithmetic
%or gradient expansion.
%
%For a simple plot of one-dimensional slopes, see slopeplot.
%

% written  12/06/98     S.M. Rump
% modified 09/28/01     S.M. Rump  matrices and multi-dimensional arrays
% modified 04/06/04     S.M. Rump  handling of derivatives of sparse arrays
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 08/26/12     S.M. Rump  global variables removed
% modified 10/03/12     S.M. Rump  SlopeSparseArrayDeriv removed
%

  if nargin==0
    INTLAB_SLOPE.NUMVAR = 0;
    INTLAB_SLOPE.Xxs = [];
    setappdata(0,'INTLAB_SLOPE',INTLAB_SLOPE);
    return
  end

  if ~isreal(xs) | ~isreal(X)
    error('Complex numbers not allowed in slope expansion')
  end

  if ( ~isa(xs,'double') & ~isa(xs,'intval') ) | ...
     ( ~isa(X,'double') & ~isa(X,'intval') )
    error('invalid initialization of slopes: input must be double or intval')
  end

  if ~isequal(size(xs),size(X))
    error('invalid initialization of slopes: dimensions do not match')
  end

  % initialize INTLAB constants
  INTLAB_SLOPE.NUMVAR = prod(size(xs));
  dummy.xs = intval(xs);
  dummy.X = intval(X);
  INTLAB_SLOPE.Xxs = dummy.X - dummy.xs;
  INTLAB_SLOPE.Xxs = INTLAB_SLOPE.Xxs(:);
  setappdata(0,'INTLAB_SLOPE',INTLAB_SLOPE);
  x = slope( dummy, 'slopeinit' );

repo: OpenMAS


text:
function s = char(q,eps,varargin)
% quaternion to char

if nargin>1 && ~isnumeric(eps)
  varargin = [{eps},varargin];
  eps = 1;
elseif nargin == 1
  eps = 1;
end

if length(q) == 1
  [alpha,beta,gamma] = Euler(q);
  if check_option(varargin,'nodegree')
    degchar = '';
  else
    degchar = mtexdegchar;
  end

  s = ['(' xnum2str([alpha,beta,gamma]./degree,...
    'precision',eps,'delimiter',[degchar ',']) degchar ')'];

else
  s = ['Rotations: ',size2str(q)];
end

author: mtex-toolbox

repo: mtex

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


text:
function A = tt_cp_vec_to_fac(x,Z)
%TT_CP_VEC_TO_FAC Converts a vector to a cell array of factor matrices.
%
%   A = TT_CP_VEC_TO_FAC(X,Z) converts the vector X into a cell array
%   of factor matrices consistent with the size of the tensor Z.
%
%   See also TT_FAC_TO_VEC, TT_CP_FUN, TT_CP_OPT.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Set-up
P = length(x);
N = ndims(Z);
sz = size(Z);

%% Determine R
R = P / sum(sz);

%% Create A
A = cell(N,1);
for n = 1:N
    idx1 = sum(sz(1:n-1))*R + 1;
    idx2 = sum(sz(1:n))*R;
    A{n} = reshape(x(idx1:idx2),sz(n),R);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [m] = esvm_mine_train_iteration(m, training_function)
%% ONE ITERATION OF: Mine negatives until cache is full and update the current
% classifier using training_function (do_svm, do_rank, ...). m must
% contain the field m.train_set, which indicates the current
% training set of negative images
% Returns the updated model (where m.mining_queue is updated mining_queue)
%
% Copyright (C) 2011-12 by Tomasz Malisiewicz
% All rights reserved.
%
% This file is part of the Exemplar-SVM library and is made
% available under the terms of the MIT license (see COPYING file).
% Project homepage: https://github.com/quantombone/exemplarsvm

% Start wtrace (trace of learned classifier parameters across
% iterations) with first round classifier, if not present already
if ~isfield(m.model,'wtrace')
  m.model.wtrace{1} = m.model.w;
  m.model.btrace{1} = m.model.b;
end

if length(m.mining_queue) == 0
  fprintf(1,' ---Null mining queue, not mining!\n');
  return;
end

%If the skip is enabled, we just update the model
if m.mining_params.train_skip_mining == 0
  [hn, m.mining_queue, mining_stats] = ...
      esvm_mine_negatives({m}, m.mining_queue, m.train_set, ...
                     m.mining_params);

  m = add_new_detections(m, cat(2,hn.xs{1}{:}), cat(1,hn.bbs{1}{: ...
                   }));
else
  mining_stats.num_visited = 0;
  fprintf(1,'WARNING: train_skip_mining==0, just updating model\n');
end

m = update_the_model(m, mining_stats, training_function);

if isfield(m,'dataset_params') && m.dataset_params.display == 1
  dump_figures(m);
end

function [m] = update_the_model(m, mining_stats, training_function)
%% UPDATE the current SVM, keep max number of svs, and show the results

if ~isfield(m,'mining_stats')
  m.mining_stats{1} = mining_stats;
else
  m.mining_stats{end+1} = mining_stats;
end

m = training_function(m);

% Append new w to trace
m.model.wtrace{end+1} = m.model.w;
m.model.btrace{end+1} = m.model.b;

% if (m.mining_params.dfun == 1)
%   r = m.model.w(:)'*bsxfun(@minus,m.model.svxs,m.model.x(:,1)).^2 - ...
%       m.model.b;
% else
%   r = m.model.w(:)'*m.model.svxs - m.model.b;
% end
% m.model.svbbs(:,end) = r;

function dump_figures(m)

% figure(1)
% clf
% show_cool_os(m)

% if (mining_params.dump_images == 1) || ...
%       (mining_params.dump_last_image == 1 && ...
%        m.iteration == mining_params.train_max_mine_iterations)
%   set(gcf,'PaperPosition',[0 0 10 3]);
%   print(gcf,sprintf('%s/%s.%d_iter=%05d.png', ...
%                     mining_params.final_directory,m.curid,...
%                     m.objectid,m.iteration),'-dpng');
% end

figure(2)
clf
Isv1 = esvm_show_det_stack(m,7);

imagesc(Isv1)
axis image
axis off
iter = length(m.model.wtrace)-1;
title(sprintf('Ex %s.%d.%s SVM-iter=%03d',m.curid,m.objectid,m.cls,iter))
drawnow
snapnow

if (m.mining_params.dump_images == 1) || ...
      (m.mining_params.dump_last_image == 1 && ...
       m.iteration == m.mining_params.train_max_mine_iterations)

  imwrite(Isv1,sprintf('%s/%s.%d_iter_I=%05d.png', ...
                    m.mining_params.final_directory, m.curid,...
                    m.objectid, m.iteration), 'png');
end

function m = add_new_detections(m, xs, bbs)
% Add current detections (xs,bbs) to the model struct (m)
% making sure we prune away duplicates, and then sort by score
%
% Tomasz Malisiewicz (tomasz@cmu.edu)

%First iteration might not have support vector information stored
if ~isfield(m.model, 'svxs') || isempty(m.model.svxs)
  m.model.svxs = [];
  m.model.svbbs = [];
end

m.model.svxs = cat(2,m.model.svxs,xs);
m.model.svbbs = cat(1,m.model.svbbs,bbs);

%Create a unique string identifier for each of the supports
names = cell(size(m.model.svbbs,1),1);
for i = 1:length(names)
  bb = m.model.svbbs(i,:);
  names{i} = sprintf('%d.%.3f.%d.%d.%d',bb(11),bb(8), ...
                             bb(9),bb(10),bb(7));
end

[unames,subset,j] = unique(names);
m.model.svbbs = m.model.svbbs(subset,:);
m.model.svxs = m.model.svxs(:,subset);

[aa,bb] = sort(m.model.w(:)'*m.model.svxs,'descend');
m.model.svbbs = m.model.svbbs(bb,:);
m.model.svxs = m.model.svxs(:,bb);

author: quantombone

repo: exemplarsvm


text:
function mprob = convBaron(prob,opts)
%CONVBARON Convert OPTI problem to BARON Nonlinear problem
%
%   mprob = convBaron(prob,opts)

%   Copyright (C) 2013 Jonathan Currie (IPL)

%Ensure all args passed
if(nargin < 2)
    error('You must supply both the problem + options');
end
if(~isstruct(prob) || ~isstruct(opts))
    error('Both prob + opts must be structures!');
end

%Get Warning Level
if(strcmpi(opts.warnings,'all'))
    warn = 2;
elseif(strcmpi(opts.warnings,'critical'))
    warn = 1;
else
    warn = 0;
end

%Assign Common Args
mprob.A = prob.A;
mprob.rl = prob.rl;
mprob.ru = prob.ru;
mprob.lb = prob.lb;
mprob.ub = prob.ub;
if(isfield(prob,'int') && isfield(prob.int,'str'))
    mprob.xtype = prob.int.str;
end

%Check if the user wants to use BARON LP/QP/QCQP solver
switch(lower(prob.type))
    case {'lp','milp','bilp','qp','qcqp','miqp','miqcqp'}
        if(warn > 1)
            optiwarn('OPTI:QCQP2NLP','Quadratic Program converted to (MI)NLP to suit BARON interface');
        end
        [prob,opts] = QCQP2NLP(prob,opts);
    case {'uno','nlp','minlp'}
        %nothing
    otherwise
        error('Problem type ''%s'' is not supported by BARON',prob.type);
end

%Assign remaining args
mprob.fun = prob.fun;
mprob.nlcon = prob.nlcon;
mprob.cl = prob.cl;
mprob.cu = prob.cu;
mprob.x0 = prob.x0;

% Set BARON options (optiset options override common ones)
%Build baronset struct, adding any passed options (ignored by baronset if not recognised)
options = baronset(opts.solverOpts);
%Add OPTISET options
options.epsa = opts.tolafun;
options.epsr = opts.tolrfun;
options.inttol = opts.tolint;
options.maxiter = opts.maxiter;
options.maxtime = opts.maxtime;
options.prlevel = dispLevel(opts.display);
%Save Options Structure
mprob.options = options;

function  print_level = dispLevel(lev)
%Return IPOPT compatible display level
switch(lower(lev))
    case'off'
        print_level = 0;
    case 'iter'
        print_level = 1;
    case 'final'
        print_level = 1;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [alpha_fin,b_fin] = rsp_max(K,Y,C,T,S_min),
% Recursive Stabilization Procedure
bias_eps=10^(-7);
m= size(K,1);
k = ceil(log(m/S_min)/log(T/(T-1)));
current_position = [0,0];
Node_done = zeros(k,T);
Current_Path = zeros(k*T,m);
not_finished = 1;
count = 0;
Path_from_root = zeros(k,1);
%%
%% Computation of the kernel
%%
while not_finished,

    cp1 = current_position(1);
    cp2 = current_position(2);

    if (cp1~=k),
        active = find(Node_done(cp1+1,1:T)==0);
    else
        active = [];
    end;%if (cp1~=k)

    if ~isempty(active),

        % The part of the execution tree is not completely explored
        if cp1==0,
            tmp=[1:m];
        else
            tmp = find(Current_Path(cp2 + (cp1-1)*T,:)==1);
        end; % if cp1==0,

        % Build the training set for the node
        tab_tmp = zeros(1,m);
        train_ind_tmp = tmp([1:floor((active(1)-1)*length(tmp)/T),floor(active(1)*length(tmp)/T):length(tmp)]);
        %% If all the points are from the same class then do not go further
        is_plus = find(Y(train_ind_tmp)==1);

        if (length(is_plus)==0|(length(is_plus)==length(train_ind_tmp))),
            Node_done(cp1+1,1:T)=1;
            Current_Path([1:T] + (cp1)*T,:) = zeros(T,m);
        else
            tab_tmp(train_ind_tmp)=1;
            current_position(1) = cp1+1;
            Path_from_root(current_position(1)) = current_position(2);
            current_position(2) = active(1);
            Current_Path(current_position(2) + (current_position(1)-1)*T,:) = tab_tmp;
        end;
    else
        if (cp1==0),
            tmp = [1:m];
        else
            tmp = find(Current_Path(current_position(2)+(current_position(1)-1)*T,:)==1);
        end;
        Ylearn = Y(tmp,:);
        Clearn = C*m/length(tmp);
        Klearn = K(tmp,tmp);

        if current_position(1)==k,
            %keyboard;
            [alpha] = quadsolve(Klearn,-ones(length(tmp),1),Ylearn',0,Clearn);
            alpha_bias=alpha;
        else
            A = Current_Path(current_position(1)*T+1:current_position(1)*T+T,tmp);
            for t=1:T,
                g(t) = A(t,:)*(Klearn*A(t,:)');
            end;
            % scale the values of gamma_i (related to the way the optimization is done)
            if C==Inf,
                scc = 1;
            else
                scc = Clearn;
            end;
            H = [[Klearn, Klearn*A'/scc];[A*Klearn/scc,A*Klearn*A'/scc^2]];
            c = -[ones(length(tmp),1);0.5*g'/scc];
            A = [[Ylearn',zeros(1,T)];[zeros(1,length(tmp)),ones(1,T)]];
            %keyboard;
            btmp = [0;scc];
            [alpha_tmp,y] = quadsolve(H,c,A,btmp,Clearn);
            alpha_bias=alpha_tmp;
           % keyboard;
            alpha = alpha_tmp(1:length(tmp)) + Current_Path(current_position(1)*T+1:current_position(1)*T+T,tmp)'*alpha_tmp(length(tmp)+1:length(alpha_tmp))/scc;
%            checka = -alpha'*Klearn*Current_Path(current_position(1)*T+1:current_position(1)*T+T,tmp)' + 0.5*g;
            %keyboard;
        end; % if current_pos...==k-1

        if current_position ~=0,
            Current_Path(current_position(2)+(current_position(1)-1)*T,tmp) = alpha';
            Node_done(current_position(1)+1,1:T)=0;
            Node_done(current_position(1),current_position(2)) = 1;
            count=count+1;
            current_position(2) = Path_from_root(current_position(1));
            current_position(1) = current_position(1)-1;
        else,
            not_finished=0;
            alpha_fin = alpha;
            b_fin = -y(1);
        end; % if current_pos...~=0

    end;% if ~isempty(active)
end; %while not_finished
% End function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [gx, gy, gz] = gaussPtetra(X1, X2, X3, X4, ng)

%% USAGE: generate Gaussian nodes on tetra mesh
%
% INPUTS:
% p --- nt-by-3 vector
% t --- nt-by-4 vector
% ng --- number of Gaussian nodes in an triangle
%            possible values = 1, 4, 5, 11, 15.
%
% OUTPUTS:
% gx --- nt-by-ng vector: stores x coordinates of Gaussian nodes
% gy --- nt-by-ng vector: stores y coordinates of Gaussian nodes
% gz --- nt-by-ng vector: stores z coordinates of Gaussian nodes
%
% Last Modified: 07/02/2020 by Xu Zhang
%%

G = zeros(size(X1,1),3*ng); i = 1;

if ng == 1 % accurate up to polynomial degree p = 1

    w1 = 1/4; % point at centroid
    G(:,[i,i+ng,i+2*ng]) = w1*(X1+X2+X3+X4);

elseif ng == 4 % accurate up to polynomial degree p = 2

    w1 = 0.585410196624969; w2 = 0.138196601125011; % median line
    G(:,[i,i+ng,i+2*ng]) = w1*X1 + w2*(X2+X3+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*X2 + w2*(X1+X3+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*X3 + w2*(X1+X2+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*X4 + w2*(X1+X2+X3);

elseif ng == 5 % accurate up to polynomial degree p = 3

    w1 = 1/4; % point at centroid
    G(:,[i,i+ng,i+2*ng]) = w1*(X1+X2+X3+X4); i = i+1;

    w1 = 1/2; w2 = 1/6;% median line
    G(:,[i,i+ng,i+2*ng]) = w1*X1 + w2*(X2+X3+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*X2 + w2*(X1+X3+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*X3 + w2*(X1+X2+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*X4 + w2*(X1+X2+X3);

elseif ng == 11 % accurate up to polynomial degree p = 4

    w1 = 1/4; % point at centroid
    G(:,[i,i+ng,i+2*ng]) = w1*(X1+X2+X3+X4); i = i+1;

    w1 = 0.785714285714286; w2 = 0.071428571428571; % median line
    G(:,[i,i+ng,i+2*ng]) = w1*X1 + w2*(X2+X3+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*X2 + w2*(X1+X3+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*X3 + w2*(X1+X2+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*X4 + w2*(X1+X2+X3); i = i+1;

    w1 = 0.399403576166799; w2 = 0.100596423833201; % line opposite sides
    G(:,[i,i+ng,i+2*ng]) = w1*(X1+X2) + w2*(X3+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*(X1+X3) + w2*(X2+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*(X1+X4) + w2*(X2+X3); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*(X2+X3) + w2*(X1+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*(X2+X4) + w2*(X1+X3); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*(X3+X4) + w2*(X1+X2);

elseif ng == 15 % accurate up to polynomial degree p = 5

    w1 = 1/4; % point at centroid
    G(:,[i,i+ng,i+2*ng]) = w1*(X1+X2+X3+X4); i = i+1;

    w1 = 0; w2 = 1/3; % median line
    G(:,[i,i+ng,i+2*ng]) = w1*X1 + w2*(X2+X3+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*X2 + w2*(X1+X3+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*X3 + w2*(X1+X2+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*X4 + w2*(X1+X2+X3); i = i+1;

    w1 = 0.727272727272727; w2 = 0.090909090909091; % median line
    G(:,[i,i+ng,i+2*ng]) = w1*X1 + w2*(X2+X3+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*X2 + w2*(X1+X3+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*X3 + w2*(X1+X2+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*X4 + w2*(X1+X2+X3); i = i+1;

    w1 = 0.066550153573664; w2 = 0.433449846426336; % line opposite sides
    G(:,[i,i+ng,i+2*ng]) = w1*(X1+X2) + w2*(X3+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*(X1+X3) + w2*(X2+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*(X1+X4) + w2*(X2+X3); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*(X2+X3) + w2*(X1+X4); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*(X2+X4) + w2*(X1+X3); i = i+1;
    G(:,[i,i+ng,i+2*ng]) = w1*(X3+X4) + w2*(X1+X2);

end
gx = G(:,1:ng); gy = G(:,ng+1:2*ng); gz = G(:,2*ng+1:3*ng);

repo: ifem


text:
function D = LRSDL_updateD(Y, D, D0, Y_range, D_range, opts, X, X0, pars)
    if nargin == 0
        addpath('../utils');
        addpath('../sparse_coding');
        tic
        d = 30;
        N = 7;
        k = 5;
        k0 = 5;
        C = 3 ;
        Y = normc(rand(d,N*C));
        D = normc(rand(d,k*C));
        D0 = normc(rand(d, k0));
        Y_range = N* (0:C);
        D_range = k* (0:C);
        opts.lambda1 = 0.01;
        opts.lambda2 = 0.002;
        opts.lambda3 = 0.1;
        opts.max_iter = 250;
        opts.show = true;
        X = randn(size(D,2), size(Y,2));
        X0 = randn(size(D0, 2), size(Y, 2));
    end

    Z = Y - D0*X0;
    W = Z - D*X;

    %% ========= update D1 ==============================

    X1 = X(1:D_range(2), :);
    X11 = X1(:, 1: Y_range(2));
    D1 = D(:, 1: D_range(2));
    W = W + D1*X1;
    E = W*X1' + get_block_col(Z, 1, Y_range)*X11';
    F = 2*X1*X1';
    D_prev = updateD_EF(D1, E, F, pars.max_iter);
    D(:, 1: D_range(2)) = D_prev;

    %% ========= update Di, i > 1 ==============================
    for i = 2: numel(Y_range) - 1
        D_cur = get_block_col(D, i, D_range);
        D_com = [-D_prev, D_cur];
        X_com = X(D_range(i-1) + 1: D_range(i+1), :);
        W = W + D_com*X_com; %W = W - D_{i-1}X^{i-1} + D_iX^i
        Xi = X(D_range(i) + 1: D_range(i+1), :);
        Xii = Xi(:, Y_range(i) + 1: Y_range(i+1));
        E = W*Xi' + get_block_col(Z, i, Y_range)*Xii';
        F = 2*Xi*Xi';
        D_prev = updateD_EF(D_cur, E, F, pars.max_iter);
        D(:, D_range(i) + 1: D_range(i+1)) = D_prev;
    end
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
%                           solveCoordinateDescent.m
%
%  Implementation of the Coordinate Descent algorithm proposed in the paper.
%
%% I/O
%  Inputs:
%     A:    m x n matrix or a function handle to a method that
%           returns A*x.
%     At:   The adjoint (transpose) of 'A'. If 'A' is a function handle, 'At'
%           must be provided.
%     b0:   m x 1 real,non-negative vector consists of all the measurements.
%     x0:   n x 1 vector. It is the initial guess of the unknown signal x.
%     opts: A struct consists of the options for the algorithm. For details,
%           see header in solvePhaseRetrieval.m or the User Guide.
%
%     Note: When a function handle is used, the
%     value of 'At' (a function handle for the adjoint of 'A') must be
%     supplied.
%
%  Outputs:
%     sol:  n x 1 vector. It is the estimated signal.
%     outs: A struct consists of the convergence info. For details,
%           see header in solvePhaseRetrieval.m or the User Guide.
%
%
%  See the script 'testCoordinateDescent.m' for an example of proper usage
%  of this function.
%
%% Notations
%  Notations mainly follow the paper's notation.
%
%% Algorithm Description
%  CD is an iterative procedure that successively minimizes the objective
%  function along coordinate directions. A single unknown is solved at each
%  iteration while all other variables are kept fixed. As a result, only
%  minimization of a univariate quartic polynomial is needed which is
%  easily achieved by finding the closed-form roots of a cubic equation.
%
%  Specifically, the method has the following steps: It keeps running until
%  the normalized gradient becomes smaller than the tolerance (1) At each
%  iteration, use the selected rule to choose an index i. (2) Minimize the
%  objective f with respect to the ith variable while keeping
%      all other 2n-1 (both real and imaginary parts are variables so there
%      are 2n in total) variables fixed by solving the cubic equation to
%      get alpha that minimize the objective along the ith variable.
%  (3) update the estimate along the ith variable by alpha.
%
%% References
%  Paper Title:   Coordinate Descent Algorithms for Phase Retrieval
%  Place:         Chapter II.B
%  Authors:       Wen-Jun Zeng, H. C. So
%  arXiv Address: https://arxiv.org/abs/1706.03474
%
% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

%% -----------------------------START-----------------------------------

function [sol, outs] = solveCoordinateDescent(A, At, b0, x0, opts)

    validateOptions(opts); % Check the validity of algorithm-specific options

    % Initialization
    m = length(b0);
    n = length(x0);
    sol = x0;
    Ax = A(sol);

    % Initialize values potentially computed at each round.
    currentTime = [];
    currentResid = [];
    currentReconError = [];
    currentMeasurementError = [];

    % Initialize vectors for recording convergence information
    [solveTimes,measurementErrors,reconErrors,residuals] = initializeContainers(opts);

    maxDiff = -inf;

    C = zeros(m, 3);

    % Used to compute gradient of objective function (used by greedy index
    % choice rule)
    f = @(z) (abs(z).^2 - b0.^2) .* z;

    startTime = tic; % Begin timer
    for iter = 1 : opts.maxIters
        switch lower(opts.indexChoice)
            case 'random'
                if opts.isComplex
                    index = randi(2*n);
                else
                    index = randi(n);
                end
            case 'cyclic'
                if opts.isComplex
                    index = mod(iter-1, 2*n) + 1;
                else
                    index = mod(iter-1, n) + 1;
                end
            case 'greedy'
                grad = At(f(A(sol)));

                if opts.isComplex
                    grad_bar = [real(grad); imag(grad)];
                else
                    grad_bar = grad;
                end
                [~, index] = max(abs(grad_bar));
        end

        vals = conj(A(double(1:n == mod(index-1, n)+1)'));
        for j = 1 : m
            if index > n
                C(j, 2) = 2*imag(Ax(j) * vals(j));
            else
                C(j, 2) = 2*real(Ax(j) * vals(j));
            end
            C(j, 3) = abs(vals(j))^2;
            C(j, 1) = abs(Ax(j))^2;
        end

        d_4 = sum(C(:, 3).^2);
        d_3 = sum(2 * C(:, 3) .* C(:, 2));
        d_2 = sum(C(:, 2).^2 + 2 * C(:, 3) .* (C(:, 1) - b0.^2));
        d_1 = sum(2 * C(:, 2) .* (C(:, 1) - b0.^2));

        % Desired alpha is a root of polynomial
        alphas = roots([4*d_4 3*d_3 2*d_2 d_1])';
        % Select only the real roots
        alphas = alphas(imag(alphas) == 0);
        % Function of alpha to be minimized
        g = @(x) d_4*x.^4 + d_3*x.^3 + d_2*x.^2 + d_1*x;
        % Find index of best alpha
        [~, idx] = min(g(alphas));
        alpha = alphas(idx);

        % Update x
        if (index > n)
            a_j = A(double(1:n == index - n)');
            sol(index - n) = sol(index - n) + 1i * alpha;
            Ax = Ax + (1i * alpha) * a_j;
        else
            a_j = A(double(1:n == index)');
            sol(index) = sol(index) + alpha;
            Ax = Ax + alpha * a_j;
        end

        diff = abs(alpha);
        maxDiff = max(diff, maxDiff);

        % Record convergence information and check stopping condition
        % If xt is provided, reconstruction error will be computed and used for stopping
        % condition. Otherwise, residual will be computed and used for stopping
        % condition.
        if ~isempty(opts.xt)
            x = sol;
            xt = opts.xt;
            %  Compute optimal rotation
            alpha = (x(:)'*xt(:))/(x(:)'*x(:));
            x = alpha*x;
            currentReconError = norm(x-xt)/norm(xt);
            if opts.recordReconErrors
                reconErrors(iter) = currentReconError;
            end
        end

        if isempty(opts.xt) | opts.recordResiduals
            currentResid = diff / max(maxDiff, 1.0e-30);
        end

        if opts.recordResiduals
            residuals(iter) = currentResid;
        end

        currentTime = toc(startTime);                % Record elapsed time so far
        if opts.recordTimes
            solveTimes(iter) = currentTime;
        end
        if opts.recordMeasurementErrors
            currentMeasurementError = norm(abs(A(sol)) - b0) / norm(b0);
            measurementErrors(iter) = currentMeasurementError;
        end

        % Display verbose output if specified
        if opts.verbose == 2
            displayVerboseOutput(iter, currentTime, currentResid, currentReconError, currentMeasurementError);
        end

        %  Test stopping criteria.
        if stopNow(opts, currentTime, currentResid, currentReconError)
            break;
        end

    end

    % Create output according to the options chosen by user
    outs = generateOutputs(opts, iter, solveTimes, measurementErrors, reconErrors, residuals);

    % Display verbose output if specified
    if opts.verbose == 1
        displayVerboseOutput(iter, currentTime, currentResid, currentReconError, currentMeasurementError);
    end
end

% Validate algorithm-specific options
function validateOptions(opts)
    validIndexChoices = {'cyclic', 'random', 'greedy'};
    checkIfInList('indexChoice',opts.indexChoice,validIndexChoices);
end

author: tomgoldstein

repo: phasepack-matlab


text:
classdef DateTime
    % This class enclapsulates date time value but behaves in a very
    % similar way as typical array of timestamps in datenum format
  %======================================================================
    %{
		Copyright (c) 2011
		This program is a result of a joined cooperation of Energocentrum
		PLUS, s.r.o. and Czech Technical University (CTU) in Prague.
        The program is maintained by Energocentrum PLUS, s.r.o. and
        licensed under the terms of MIT license. Full text of the license
        is included in the program release.

        Author(s):
		Jiri Cigler, Dept. of Control Engineering, CTU Prague & Automatic Control Laboratory, ETH Zurich
		Jan  Siroky, Energocentrum PLUS s.r.o.

        Implementation and Revisions:

        Auth  Date        Description of change
        ----  ---------   -------------------------------------------------
        jc    01-Mar-11   First implementation
        jc    30-Sep-11   Added function colon
        jc    07-Jan-12   Added functions addtodate,datevec,weekday
    %}
    %======================================================================

    % TODO

%ate and Time Operations
%addtodateModify date number by fieldcalendarCalendar for specified monthclockCurrent time as date vectorcputimeElapsed CPU timedateCurrent date stringdatenumConvert date and time to serial date numberdatestrConvert date and time to string formatdatevecConvert date and time to vector of componentseomdayLast day of monthetimeTime elapsed between date vectorsnowCurrent date and timeweekdayDay of wee

    properties
        serialDate
    end
    methods
        function this = DateTime(varargin)
            if numel(varargin)==1 && isa(varargin{1},'java.util.Date')
                    df = java.text.SimpleDateFormat( 'yyyy-MM-dd HH:mm:ss' );
                    tz = java.util.TimeZone.getTimeZone ('UTC');
                    df.setTimeZone( tz );
                    this.serialDate=datenum(char(df.format(varargin{1})));
                    %this.serialDate = datenum(char(varargin{1}.toString)) - varargin{1}.getTimezoneOffset/60/24;
                    %disp ( [ char(varargin{1}.toGMTString), '---' char(varargin{1}.toString ), '---', char(varargin{1}.toLocaleString )]);
                    %if (varargin{1}.getTimezoneOffset)~=-120
                    %    disp(1);
                    %end
            else
                this.serialDate=datenum(varargin{:});
            end
        end
        function this = plus(this,val)
            o =@plus;
            this = doFun(this,o,val);
        end
        function this = minus(this,val)
            o =@minus;
            this = doFun(this,o,val);
        end
        function this = times(this,val)
            o =@times;
            this = doFun(this,o,val);
        end

        function this = mtimes(this,val)
            o =@mtimes;
            this = doFun(this,o,val);
        end

        function this = mrdivide(this,val)
            o =@mrdivide;
            this = doFun(this,o,val);
        end

        function this = rdivide(this,val)
            o =@rdivide;
            this = doFun(this,o,val);
        end

        function this = horzcat(this,varargin)
            %this.serialDate = [this.serialDate, n.serialDate];
            for i=1:numel(varargin)
                this.serialDate = [this.serialDate, varargin{i}.serialDate];
            end
        end

        function out = colon(this,step,to)
            vect = [double(this):double(step):double(to)]';

            out =DateTime(vect);
        end

        function this = vertcat(this,varargin)
            for i=1:numel(varargin)
                this.serialDate = [this.serialDate; varargin{i}.serialDate];
            end
        end

        function this = ctranspose(this)
            this.serialDate = this.serialDate';
        end

        function this = transpose(this)
            this.serialDate = this.serialDate';
        end
        function  disp(this)
            disp([this.serialDate])
        end
        function out = double(this)
            out = this.serialDate;
        end
        function out = length(this)
            out = length(this.serialDate);
        end

        function out = size(this,varargin)
            out = size(this.serialDate,varargin{:});
        end

        function out = numel(this)
            out = numel(this.serialDate);
        end
        function out = isreal(this)
            out = isreal(this.serialDate);
        end
        function out = isnan(this)
            out = isnan(this.serialDate);
        end
        function out = isfinite(this)
            out = isfinite(this.serialDate);
        end

        function out = le(this,B)
            if isa(B,'DateTime')
                out = le(this.serialDate,B.serialDate);
            else
                out = le(this.serialDate,B);
            end
        end

        function out = lt(this,B)
            fun=@lt;
            if isa(B,'DateTime')
                out = fun(this.serialDate,B.serialDate);
            else
                out = fun(this.serialDate,B);
            end
        end
        function out = gt(this,B)
            fun=@gt;
            if isa(B,'DateTime')
                out = fun(this.serialDate,B.serialDate);
            else
                out = fun(this.serialDate,B);
            end
        end
        function out = eq(this,B)
            fun=@eq;
            if isa(B,'DateTime')
                out = fun(this.serialDate,B.serialDate);
            else
                out = fun(this.serialDate,B);
            end
        end
        function out = diff(this)
            out = diff(this.serialDate);
        end

        function out = norm(this,varargin)
            out = norm(this.serialDate,varargin{:});
        end

        function [this k] = sort(this,varargin)
            [this.serialDate k] = sort(this.serialDate,varargin{:});
        end

        function this = subsref(this,S)
            if isa(S.subs{1},'DateTime')
                S.subs{1}=double(S.subs{1});
            end

            this.serialDate =  subsref(this.serialDate,S);
        end

        function idx = subsindex(this)
            idx = double(this)-1;
        end

        function endidx = end(this,k,n)
            if size(this.serialDate,1)==1 || size(this.serialDate,2)==1
                endidx=numel(this.serialDate);
            else
                endidx = size(this.serialDate,k);
            end
        end

        function this = subsasgn(this, S, B)
            if not(isa(B,'DateTime'))
                B=DateTime(B);
            end

            this.serialDate =subsasgn(this.serialDate, S, B);
        end

        function res = bsxfun(fun,A,B)
            res = fun(A,B);
        end

        function out =superiorfloat (x,y,xi)
            if isa(x,'DateTime') && isa(xi,'DateTime')
                out = superiorfloat(x.serialDate,y,xi.serialDate);
            elseif isa(x,'DateTime') && not(isa(xi,'DateTime'))
                out = superiorfloat(x.serialDate,y,xi);
            elseif not(isa(x,'DateTime')) && isa(xi,'DateTime')
                out = superiorfloat(x,y,xi.serialDate);
            else
                out = superiorfloat(x,y,xi);
            end
        end

        function this = floor(this)
            this.serialDate = floor(this.serialDate);
        end
        function this = max(this,varargin)
            this.serialDate = max(this.serialDate,varargin{:});
        end
        function this = min(this,varargin)
            this.serialDate = min(this.serialDate,varargin{:});
        end
        function out = datestr(this,varargin)
            out = datestr(this.serialDate,varargin{:});
        end

        function out = addtodate(this,varargin)
            out = addtodate(this.serialDate,varargin{:});
        end
        function varargout= datevec(this,varargin)
            nout = nargout;
            if nout <=1
                varargout{1} = datevec(this.serialDate,varargin{:});
            elseif nout ==2
                [varargout{1} varargout{2}] = datevec(this.serialDate,varargin{:});

            elseif nout ==3
                [varargout{1} varargout{2} varargout{3}] = datevec(this.serialDate,varargin{:});

            elseif nout ==4
                [varargout{1} varargout{2} varargout{3} varargout{4}] = datevec(this.serialDate,varargin{:});

            elseif nout ==5
                [varargout{1} varargout{2} varargout{3} varargout{4} varargout{5} ] = datevec(this.serialDate,varargin{:});

            elseif nout ==6
                [varargout{1} varargout{2} varargout{3} varargout{4} varargout{5} varargout{6} ] = datevec(this.serialDate,varargin{:});
            else
                error('Unknown function call');
            end
        end

    end

    methods (Access = private)
        function this = doFun (this,o, val)
            if isa(val,'DateTime') && isa(this,'DateTime')
                this.serialDate=o(this.serialDate, val.serialDate);
            elseif isa(val,'DateTime') && not(isa(this,'DateTime'))
                val.serialDate=o(this, val.serialDate);
                this = val;
            elseif not(isa(val,'DateTime')) && (isa(this,'DateTime'))
                this.serialDate=o(this.serialDate, val);
            else
                this.serialDate=DateTime(o(this, val));
            end
        end

    end

end

author: ziyuw

repo: rembo

save_path: github-repos/MATLAB/ziyuw-rembo