text:
function [accuracy,Mdl] = GiveMeCfn(XTrain,yTrain,XTest,yTest,cfnParams,beVerbose)
% GiveMeCfn    Returns classification results from training a classifier on
%               training/test data
%
%---INPUTS:
% XTrain -- training data matrix
% yTrain -- training data labels
% XTest -- testing data matrix
% yTest -- testing data labels
% cfnParams -- parameters for classification
% beVerbose -- [default: false] whether to give text output of progress
%
%---OUTPUTS:
% accuracy -- IF cfnParams.computePerFold = true:
%               a 2x1 cell of accuracy results in each fold for training data (element 1)
%               and test data (element 2).
%             IF cfnParams.computePerFold = false:
%               A real number giving the average accuracy in the test folds.
% Mdl -- the trained classification model.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

%-------------------------------------------------------------------------------
%% Check Inputs:
%-------------------------------------------------------------------------------
% Ensure y is a column vector
if size(yTrain,1) < size(yTrain,2)
    yTrain = yTrain';
end
if nargin < 3
    XTest = [];
end
if nargin < 4
    yTest = [];
end
if nargin < 6
    beVerbose = false;
end

%-------------------------------------------------------------------------------
%% Define the classification model
%------------------------------------------------------------------------------
if strcmp(cfnParams.whatClassifier,'fast-linear')
    %--------------------------------------------------------------------------
    % Special case -- the `classify` function is faster than others
    %--------------------------------------------------------------------------
    if cfnParams.numFolds > 0
        if cfnParams.doReweight
            Mdl = fitcdiscr(XTrain,yTrain,'Prior','empirical',...
            'DiscrimType','linear','Weights',InverseProbWeight(yTrain),...
            'KFold',cfnParams.numFolds);
        else
            Mdl = fitcdiscr(XTrain,yTrain,'Prior','empirical',...
            'DiscrimType','linear','KFold',cfnParams.numFolds);
        end
    else
        if nargout == 1
            % It's faster to use classify if we only want the accuracy
            accuracy = BF_LossFunction(yTest,classify(yTest,XTrain,yTrain),...
                            cfnParams.whatLoss,cfnParams.classLabels);
            return;
        else
            Mdl = fitcdiscr(XTrain,yTrain,'Prior','empirical');
        end
    end
else
    if cfnParams.numClasses==2
        %----------------------------------------------------------------------
        % Special case: a two-class model
        %----------------------------------------------------------------------
        switch cfnParams.whatClassifier
        case 'knn'
            if beVerbose
                fprintf(1,'Using three neighbors for knn\n');
            end
            if cfnParams.numFolds > 0
                Mdl = fitcknn(XTrain,yTrain,'NumNeighbors',3,'KFold',cfnParams.numFolds);
            else
                Mdl = fitcknn(XTrain,yTrain,'NumNeighbors',3);
            end
        case 'tree'
            if cfnParams.numFolds > 0
                Mdl = fitctree(XTrain,yTrain,'KFold',cfnParams.numFolds);
            else
                Mdl = fitctree(XTrain,yTrain);
            end
        case {'linear','linclass'}
            if cfnParams.numFolds > 0
                Mdl = fitcdiscr(XTrain,yTrain,'FillCoeffs','off','SaveMemory','on','KFold',cfnParams.numFolds);
            else
                Mdl = fitcdiscr(XTrain,yTrain,'FillCoeffs','off','SaveMemory','on');
            end
        case 'logistic'
            if cfnParams.doReweight
                if cfnParams.numFolds > 0
                    Mdl = fitclinear(XTrain,yTrain,'Learner','logistic','Regularization','ridge',...
                                'Lambda','auto','Weights',InverseProbWeight(yTrain),'KFold',cfnParams.numFolds);
                else
                    Mdl = fitclinear(XTrain,yTrain,'Learner','logistic','Regularization','ridge',...
                                'Lambda','auto','Weights',InverseProbWeight(yTrain));
                end
            else
                if cfnParams.numFolds > 0
                    Mdl = fitclinear(XTrain,yTrain,'Learner','logistic','Regularization','ridge',...
                                'Lambda','auto','KFold',cfnParams.numFolds);
                else
                    Mdl = fitclinear(XTrain,yTrain,'Learner','logistic','Regularization','ridge',...
                                'Lambda','auto');
                end
            end
        case {'svm-linear','svm-linear-highdim'}
            % This one uses the fitclinear function (fast; suited for high-dimensional datasets)
            % lambdaDefault = 1/size(XTrain,1);
            % Weight observations by inverse class probability:
            if cfnParams.doReweight
                if cfnParams.numFolds > 0
                    Mdl = fitclinear(XTrain,yTrain,'Learner','svm','Lambda','auto','Regularization','ridge',...
                                'Weights',InverseProbWeight(yTrain),'KFold',cfnParams.numFolds);
                else
                    Mdl = fitclinear(XTrain,yTrain,'Learner','svm','Lambda','auto','Regularization','ridge',...
                                'Weights',InverseProbWeight(yTrain));
                end
            else
                if cfnParams.numFolds > 0
                    Mdl = fitclinear(XTrain,yTrain,'Learner','svm','Lambda','auto',...
                                        'Regularization','ridge','KFold',cfnParams.numFolds);
                else
                    Mdl = fitclinear(XTrain,yTrain,'Learner','svm','Lambda','auto','Regularization','ridge');
                end
            end
        case {'svm','svm-linear-lowdim'}
            % For low-dim datasets (e.g., using a reduced feature set), the 'traditional' linear SVM (fitcsvm)
            % can be preferable to the above (fitclinear)
            if cfnParams.doReweight
                if cfnParams.numFolds > 0
                    Mdl = fitcsvm(XTrain,yTrain,'KernelFunction','linear','Weights',...
                                    InverseProbWeight(yTrain),'KFold',cfnParams.numFolds);
                else
                    Mdl = fitcsvm(XTrain,yTrain,'KernelFunction','linear','Weights',InverseProbWeight(yTrain));
                end
            else
                if cfnParams.numFolds > 0
                    Mdl = fitcsvm(XTrain,yTrain,'KernelFunction','linear','KFold',cfnParams.numFolds);
                else
                    Mdl = fitcsvm(XTrain,yTrain,'KernelFunction','linear');
                end
            end
        case 'svm-rbf'
            % Weight observations by inverse class probability:
            if cfnParams.doReweight
                if cfnParams.numFolds > 0
                    Mdl = fitcsvm(XTrain,yTrain,'KernelFunction','rbf','Weights',...
                                    InverseProbWeight(yTrain),'KFold',cfnParams.numFolds);
                else
                    Mdl = fitcsvm(XTrain,yTrain,'KernelFunction','rbf','Weights',...
                                                    InverseProbWeight(yTrain));
                end
            else
                if cfnParams.numFolds > 0
                    Mdl = fitcsvm(XTrain,yTrain,'KernelFunction','rbf','KFold',cfnParams.numFolds);
                else
                    Mdl = fitcsvm(XTrain,yTrain,'KernelFunction','rbf');
                end
            end
        case 'diaglinear'
            if cfnParams.numFolds > 0
                Mdl = fitcnb(XTrain,yTrain,'KFold',cfnParams.numFolds);
            else
                Mdl = fitcnb(XTrain,yTrain);
            end
        otherwise
            error('Unknown classifier: ''%s''',cfnParams.whatClassifier);
        end
    else
        %----------------------------------------------------------------------
        % Define and fit a classification model involving 3 or more classes
        %----------------------------------------------------------------------
        % Define the model:
        switch cfnParams.whatClassifier
        case 'knn'
            t = templateKNN('NumNeighbors',3,'Standardize',true);
            if beVerbose, fprintf(1,'Using knn classifier\n'); end
        case 'tree'
            t = templateTree('Surrogate','off');
            if beVerbose, fprintf(1,'Using a tree classifier\n'); end
        case {'linear','linclass'}
            t = templateDiscriminant('DiscrimType','linear','FillCoeffs','off','SaveMemory','on');
            if beVerbose, fprintf(1,'Using a linear classifier\n'); end
        case 'diaglinear'
            % Naive Bayes classifier:
            t = templateNaiveBayes('DistributionNames','normal');
            if beVerbose, fprintf(1,'Using a naive Bayes classifier\n'); end
        case {'svm','svm-linear','svm-linear-lowdim'}
            t = templateSVM('KernelFunction','linear');
            if beVerbose, fprintf(1,'Using a linear svm classifier\n'); end
        case 'svm-rbf'
            t = templateSVM('KernelFunction','rbf');
            if beVerbose, fprintf(1,'Using a rbf svm classifier\n'); end
        case 'logistic'
            error('Logistic regression is not applicable to problems with >2 classes')
        otherwise
            error('Unknown classifier: ''%s''',cfnParams.whatClassifier);
        end

        % Fit the model:
        if ismember(cfnParams.whatClassifier,{'svm-linear','svm-rbf','linear','linclass','diaglinear'}) && cfnParams.doReweight
            % Reweight to give equal weight to each class (in case of class imbalance)
            if cfnParams.numFolds > 0
                Mdl = fitcecoc(XTrain,yTrain,'Learners',t,'Weights',InverseProbWeight(yTrain),'KFold',cfnParams.numFolds);
            else
                Mdl = fitcecoc(XTrain,yTrain,'Learners',t,'Weights',InverseProbWeight(yTrain));
            end
        else
            % Don't reweight:
            if cfnParams.numFolds > 0
                Mdl = fitcecoc(XTrain,yTrain,'Learners',t,'KFold',cfnParams.numFolds);
            else
                Mdl = fitcecoc(XTrain,yTrain,'Learners',t);
            end
        end
    end
end

%-------------------------------------------------------------------------------
%% Evaluate performance on test data (in-sample, or from CV folds)
%-------------------------------------------------------------------------------
% Predict the test data:
if cfnParams.numFolds == 0
    if isempty(XTest) || isempty(yTest)
        % No need to compute this if you're just after the model
        accuracy = [];
        return
    else
        yPredict = predict(Mdl,XTest);
        accuracy = BF_LossFunction(yTest,yPredict,cfnParams.whatLoss,cfnParams.classLabels);
    end
else
    % Test data is mixed through the training data provided using k-fold cross validation
    % Output is the accuracy/loss measure for each fold
    if cfnParams.computePerFold
        % Compute separately for each test fold:
        yPredict = kfoldPredict(Mdl);
        accuracyTestFolds = arrayfun(@(x) BF_LossFunction(yTrain(Mdl.Partition.test(x)),...
                            yPredict(Mdl.Partition.test(x)),cfnParams.whatLoss,...
                            cfnParams.classLabels),1:cfnParams.numFolds);

        % Compute for each training fold:
        accuracyTrainFolds = zeros(cfnParams.numFolds,1);
        for k = 1:cfnParams.numFolds
            XTrain_k = XTrain(Mdl.Partition.training(k),:);
            yTrain_k = yTrain(Mdl.Partition.training(k));
            yPredict_k = predict(Mdl.Trained{k},XTrain_k);
            accuracyTrainFolds(k) = BF_LossFunction(yTrain_k,yPredict_k,cfnParams.whatLoss,...
                            cfnParams.classLabels);
        end

        % Combine train/test accuracies:
        accuracy = {accuracyTrainFolds,accuracyTestFolds};

        % Warning for training fold over-fitting:
        if ~cfnParams.suppressWarning && all(accuracyTrainFolds==100)
            warning('100% in-sample accuracy alert (just FYI: no need for any action!)')
            % Could suggest regularization or reducing the number of features
        end
    else
        yPredict = kfoldPredict(Mdl);
        % Aggregate all outputs from folds and compute accuracy once from that aggregate:
        accuracy = BF_LossFunction(yTrain,yPredict,cfnParams.whatLoss,cfnParams.classLabels);

        % % Above is equivalent to computing the accuracy in each test fold, then computing the mean
        % % across folds:
        % accuracyTestFolds = arrayfun(@(x) BF_LossFunction(yTrain(Mdl.Partition.test(x)),...
        %                 yPredict(Mdl.Partition.test(x)),cfnParams.whatLoss,...
        %                 cfnParams.classLabels),1:cfnParams.numFolds);
        % accuracy = mean(accuracyTestFolds);
    end
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function anal = er_chopTSeries2(tSeries,trials,params,varargin);
% anal = er_chopTSeries2(tSeries,trials,params,[options]);
%
% View-independent verson of er_chopTSeries.
% chops up entered tSeries according to the assigned parfiles.
%
% Inputs are tSeries (matrix where rows are time points, cols are
% different tSeries, say from diff't rois/voxels), and trials
% (a struct obtained from er_concatParfiles, containing design matrix
% information).
%
% returns an analysis struct with the following fields:
%
%   wholeTc: vector of whole time course, concatenated across scans,
%            for selected scans, voxels
%   allTcs:  3D matrix of time courses for every trial. The rows
%            are different time points, the columns are different
%            trials, and the slices (z dim) are different conditions.
%            As with other trials, the data is taken from the specified
%            time window, incl. prestim frames.
%   meanTcs: matrix of mean time courses for each condition. Rows
%            are different time points, columns are different conds.
%   sems:    corresponding standard errors of the mean for meanTcs.
%   timeWindow: vector specifying the time in seconds, relative to the
%               trial onset, from which each trial / mean trial time
%               course is taken. [default is -4:16].
%   peakPeriod: time, in seconds, where the HRF is expected to peak.
%               This will be used for t-test and amplitude results below.
%               [default is 8:14].
%   bslPeriod:  time, in seconds, to consider as a baseline period.
%               [default is 2:6].
%   amps:       Estimated amplitudes for each trial, taken
%               as the mean amplitude during the peak period
%               minus the mean amplitude during the baseline period.
%               This is a 2D matrix where rows are different trials
%               and columns are different conditions.
%   relamps:    relative fMRI amplitudes (e.g. dot products with the
%               mean time course) for each trial, in the same format
%               as amps. These should be less sensitive to the choice
%               of peak and baseline periods than amps, but may not be
%               accurate if different conditions have fundamentally
%               different response shapes (e.g., if there's a baseline
%               period where everything is decreasing and all other
%               conditions are increasing.)
%   Hs:         1 x nConds binary vector reflecting whether
%               each condition had a mean significant activation
%               during peak relative to baseline periods
%               (one-sided t-test alpha = 0.05 by default but
%               may be entered as an optional argument).
%   ps:         Corresponding p-values for the t-tests for Hs.
%   SNR:        Mean signal to noise ratio for all trials.
%               (mean signal for peakPeriod - bslPeriod)/(stdev bslPeriod)
%   SNRdb:      Expression of SNR as decibels: 20 * log10(SNR).
%
% Many params can be entered as optional arguments. In these cases,
% call them with the form ...'arg name',[value],.... The fields
% that can be entered are: timeWindow, peakPeriod, bslPeriod, alpha.
%
% Further options:
%   barebones:          do only a minimal analysis, extracting
%                       mean time courses, amplitudes, and SEMs.
%                       This is useful for across-voxel analyses.
%
%   normBsl,[1 or 0]:   if 1, will align trials during the baseline
%                       period [default is 1].
%   alpha,[val]:        alpha value for t-tests [default 0.05].
%   onsetDelta,[val]:   automatically shift the onsets in parfiles
%                       relative to the time course by this amount
%                       (e.g. to compensate for HRF rise time).
%   'mrvWaitbar':          put up a mrvWaitbar instead of showing progress in
%                       the command line.
%   'findPeaks':        when calculating response amplitudes, figure
%                       out the peak amplitude separately for each
%                       condition, taking the peak and surrounding 2 points.
%   [params struct]:    input a struct containing all of these
%                       params. See er_getParams.
%
% 06/17/04 ras: wrote it.
% 07/28/04 ras: clarified annotation.
% 01/25/05 ras: can now input params struct.
% 05/25/05 ras: fixed calculation of relative amplitudes.
if ieNotDefined('params')
    params = er_defaultParams;
end

%%%%% params/defaults %%%%%
barebones = 1;                  % if 0, do full analysis; if 1, do minimal analysis
normBsl = params.normBsl;       % flag to zero baseline or not
alpha = params.alpha;           % threshold for significant activations
bslPeriod = params.bslPeriod;   % period to use as baseline in t-tests, in seconds
peakPeriod = params.peakPeriod; % period to look for peaks in t-tests, in seconds
timeWindow = params.timeWindow; % seconds relative to trial onset to take for each trial
onsetDelta = params.onsetDelta; % # secs to shift onsets in parfiles, relative to time course
snrConds = params.snrConds;     % For calculating SNR, which conditions to use (if empty, use all)
waitbarFlag = 0;        % flag to show a graphical mrvWaitbar to show load progress
findPeaksFlag = 0;      % when computing amps, find peak period separately for each cond
TR = trials.TR;

%%%%% parse the options %%%%%
varargin = unNestCell(varargin);
for i = 1:length(varargin)
    if isstruct(varargin{i})
        % assume it's a params struct
        names = fieldnames(varargin{i});
        for j = 1:length(names)
            cmd = sprintf('%s = varargin{i}.%s;',names{j},names{j});
            eval(cmd);
        end

    elseif ischar(varargin{i})
        switch lower(varargin{i})
        case 'barebones', barebones = 1;
        case 'normbsl', normBsl = varargin{i+1};
        case 'alpha', alpha = varargin{i+1};
        case 'peakperiod', peakPeriod = varargin{i+1};
        case 'bslperiod', bslPeriod = varargin{i+1};
        case 'timewindow', timeWindow = varargin{i+1};
        case 'onsetdelta', onsetDelta = varargin{i+1};
        case 'snrconds', snrConds = varargin{i+1};
        case 'mrvWaitbar', waitbarFlag = 1;
        case 'findpeaks', findPeaksFlag = 1;
        otherwise, % ignore
        end
    end
end

%%%%% account for format of time series
% in this code compared to elsewhere (dumb, I know, but there's a reason)
if size(tSeries,2)==1 & size(tSeries,1) > 1 % column vector
    wholeTc = tSeries'; % flip the col vector around
elseif (size(tSeries,1)>1) & (size(tSeries,2)>1)
    % what the heck, let's go recursively on the columns
    for i = 1:size(tSeries,2)
        anal(i) = er_chopTSeries2(tSeries(:,i),trials,params,varargin);
    end
    return
else
    wholeTc = tSeries;
end
clear tSeries;

% account for onset shift
if mod(onsetDelta,TR) ~= 0
    % ensure we shift by an integer # of frames
    onsetDelta = TR * round(onsetDelta/TR);
end
trials.onsetSecs = trials.onsetSecs + onsetDelta;
trials.onsetFrames = trials.onsetFrames + onsetDelta/TR;

%%%%% get nConds from trials struct
condNums = unique(trials.cond(trials.cond >= 0));
nConds = length(condNums);

%%%%% get a set of label names, if they were specified in the parfiles
for i = 1:nConds
    ind = find(trials.cond==condNums(i));
    labels{i} = trials.label{ind(1)};
end

%%%%% convert params expressed in secs into frames
% first, find the min and max time in seconds expressed in time window:
% with this change, we ignore the intervening time-window values (e.g.,
% timeWindow=[-4 20] and timeWindow=[-4:2:20] produce the same result, so
% we're more flexible when it comes to non-integer TRs)
t1 = min(timeWindow);  t2 = max(timeWindow);

% second, re-express t1 and t2 in terms of TRs / MR frames
f1 = fix(t1 / TR);  f2 = fix(t2 / TR);

% now have a 'frame window' specifying the time window in integer MR frames
frameWindow = f1:f2;

% figure out auxilliary values: how many frames before the stimulus starts?
prestim = max( -f1, 0 ); % if f1 > 0, there is no prestim (only zuul :)

% which frames reflect the peak and baseline periods?
pk1 = fix( min(peakPeriod) / TR );
pk2 = fix( max(peakPeriod) / TR );
peakFrames = find(ismember(frameWindow, pk1:pk2));

bs1 = fix( min(bslPeriod) / TR );
bs2 = fix( max(bslPeriod) / TR );
bslFrames = find(ismember(frameWindow, bs1:bs2));

%%%%% remove trials at the very end of a scan, without
%%%%% enough data to fill the time window
cutOff = find(trials.onsetFrames+frameWindow(end) > length(wholeTc)+1);
if ~isempty(cutOff)
    keep = setdiff(1:length(trials.cond),cutOff);
    trials.cond = trials.cond(keep);
    trials.onsetFrames = trials.onsetFrames(keep);
    trials.onsetSecs = trials.onsetSecs(keep);
end

%%%%% build allTcs matrix of time points x trials x  conditions
%%%%% take (frameWindow) secs from each trial
allTcs = [];

for i = 1:nConds
   cond = condNums(i);
   ind = find(trials.cond==cond);
   for j = 1:length(ind)
       tstart = max(trials.onsetFrames(ind(j)),1);
       tend = min([tstart+frameWindow(end),length(wholeTc)]);
       rng = tstart:tend;

       % add prestim
       if tstart < prestim+1
           % for 1st trial, no baseline available -- set to 0
           allTcs(:,j,i) = [zeros(1,prestim) wholeTc(rng)]';
       else
           % augment the range by previous [prestim] frames
           fullrng = rng(1)-prestim:rng(end);
           allTcs(1:length(fullrng),j,i) = wholeTc(fullrng)';
       end

       % remove baseline estimate, if selected
       if normBsl
           % estimate DC offset by prestim baseline vals
           DC = nanmean(allTcs(bslFrames,j,i));
           allTcs(:,j,i) = allTcs(:,j,i) - DC;
       end
   end
end

%%%%% find 'empty' trials, set to NaNs
% (Empty trials will result if some conditions have more
% trials than others -- in the conditions w/ fewer trials,
% the allTcs matrix will be padded with 0s to keep it a cube).
for y = 1:size(allTcs,2)
    for z = 1:size(allTcs,3)
        if all(allTcs(:,y,z)==0)
            allTcs(:,y,z) = NaN;
        end
    end
end

%%%%% get mean time courses, sems for each condition
meanTcs = zeros(length(frameWindow),nConds);
sems = zeros(length(frameWindow),nConds);
maxNTrials = size(allTcs,2);

for i = 1:nConds
    nTrials = size(allTcs,2) - sum(any(isnan(allTcs(:,:,i))));
    if maxNTrials > 1
        meanTcs(:,i) = nanmean(allTcs(:,:,i)')';
        sems(:,i) = nanstd(allTcs(:,:,i)')' ./ sqrt(nTrials);
    else
        meanTcs(:,i) = allTcs(:,:,i);
    end
end

%%%%% calc amplitudes, do t-tests of post-baseline v. baseline
for i = 1:nConds
    if findPeaksFlag==1
        % find the peak separately for
        % each condition
        maxVal = max(meanTcs(2:end-1,i));
        maxT = find(meanTcs(:,i)==maxVal);
        peak = allTcs(maxT-1:maxT+1,:,i);
        bsl = allTcs(bslFrames,:,i);
    else
        bsl = allTcs(bslFrames,:,i);
        peak = allTcs(peakFrames,:,i);
    end
    amps(:,i) = (mean(peak) - mean(bsl))';

    if ~barebones
        [Hs(i) ps(i)] = ttest2(bsl(:),peak(:),alpha,-1);
    end
end

%%%%% compute Signal-to-Noise Ratio
if ~barebones
    if isempty(snrConds)
        snrConds = find(condNums > 0);
    else
        % index into condNumbers (e.g. so you can select 0)
        snrConds = find(ismember(condNums,snrConds));
    end

	allBsl = meanTcs(bslFrames,snrConds);
	allPk = meanTcs(peakFrames,snrConds);
	SNR = abs(mean(allPk(:)) - mean(allBsl(:))) / std(allBsl(:));
end

%%%%% compute relamps
% the resulting matrix will be of size
% nTrials x nConds
if ~barebones
%	relamps = er_relamps(allTcs);
end

%%%%% assign everything to the output struct
anal.wholeTc = double(wholeTc);
anal.allTcs = allTcs;
anal.meanTcs = double(meanTcs);
anal.sems = sems;
anal.labels = labels;
anal.timeWindow = timeWindow(mod(timeWindow,TR)==0);
anal.peakPeriod = peakPeriod(mod(peakPeriod,TR)==0);
anal.bslPeriod = bslPeriod(mod(bslPeriod,TR)==0);
anal.condNums = condNums;
anal.amps = amps;

if ~barebones
	anal.Hs = zeros(1, nConds); %Hs;
	anal.ps = ones(1, nConds); %ps;
	% anal.relamps = relamps;
	anal.SNR = SNR;
	anal.SNRdb = 20 * log10(SNR);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% PURPOSE  : 	Restore messy event code values.
%                 pop_setcodebit sets the bit(s) at position(s) "bitindex" in each EEG.event(i).type value to 0 (off).
%                 "bitindex" must contain number(s) between 1 and 16.
%
% FORMAT   :
%
% EEG = pop_setcodebit(EEG, bitindex, newvalue)
%
% or
%
% pop_rt2text(ERPLAB, varargin)
%
%
% INPUTS   :
%
%    EEG           - continuous dataset
%    bitindex      - bit position(s). From 1-16
%    newvalue      - 1 or 0 (zero).
%
%
% OUTPUTS :
%
%    EEG           - updated continuous dataset%
%
%
% Example:
%
% In Biosemi system you have a 16-bit word for sending your event codes. However, you generally only use event codes from 1 to 255
% (8-bit numbers). In this case, the upper byte (bits 9 to 16) should be silent, and these bits should be zero.
% Unfortunately, sometimes this does not happen and you get different and larger event codes.
% You may find some discussions and proposed solutions about this subject in many blogs on the internet.
% However, as a general way to deal with this issue, setcodebit.m will help you to control every single bit from this 16-bit word.
% Hence, you will be able to set each bit, or a group of them, to either "0" or "1", using a single line command.
% For example, if you want to assure that your event codes keep values from 1 to 255, you must set the upper byte
% (bits 9 to 16) to zero (cleaning any spuriously activated bit), using a command like the following:
%
%
% EEG = pop_setcodebit(EEG, 9:16, 0);
%
%
% Note, EEG = pop_setcodebit(EEG, 9:16);  will work as well, since "0" is the default value for "newvalue".
%
%
% See also setcodebitGUI.m setcodebit.m
%
% *** This function is part of ERPLAB Toolbox ***
% Author: Javier Lopez-Calderon & Johanna Kreither
% Center for Mind and Brain
% University of California, Davis,
% Davis, CA
% 2011
%
% Thanks to Eric Foo for helping with testing and the help section.

function [EEG, com ] = pop_setcodebit(EEG, bitindex, newvalue, varargin)
com = '';
if nargin<1
        help setcodebit
        return
end
if isobject(EEG) % eegobj
        whenEEGisanObject % calls a script for showing an error window
        return
end
if ~isempty(EEG.epoch)
        msgboxText = 'pop_setcodebit.m only works for continuous dataset.';
        title = 'ERPLAB: setcodebitGUI few inputs';
        errorfound(msgboxText, title);
        return
end
if ischar(EEG.event(1).type)
        msgboxText = 'Your event codes are not numeric.';
        title = 'ERPLAB: pop_setcodebit few inputs';
        errorfound(msgboxText, title);
        return
end
if nargin==1
        % call GUI
        answer = setcodebitGUI;
        if isempty(answer)
                disp('User selected Cancel')
                return
        end

        todo =  answer{1};
        if strcmpi(todo, 'all')
                bitindex = 1:16;
                newvalue = 0;
        elseif strcmpi(todo, 'lower')
                bitindex = 1:8;
                newvalue = 0;
        elseif strcmpi(todo, 'upper')
                bitindex = 9:16;
                newvalue = 0;
        else
                if ~isempty(answer{2})
                        bitindex = answer{2};
                        newvalue = 0;
                        [EEG, com ] = pop_setcodebit(EEG, bitindex, newvalue, 'History', 'off');
                end
                if ~isempty(answer{3})
                        bitindex = answer{3};
                        newvalue = 1;
                        [EEG, com2 ] = pop_setcodebit(EEG, bitindex, newvalue, 'History', 'off');
                        com = sprintf('%s\n%s', com, com2);
                end
                return
        end

        %
        % Somersault
        %
        [EEG, com ] = pop_setcodebit(EEG, bitindex, newvalue, 'History', 'gui');
        return
end
p = inputParser;
p.FunctionName  = mfilename;
p.CaseSensitive = false;
p.addRequired('EEG');
p.addRequired('bitindex', @isnumeric);
p.addRequired('newvalue', @isnumeric);
p.addParamValue('History', 'script', @ischar);             % history from scripting
p.parse(EEG, bitindex, newvalue, varargin{:});

if strcmpi(p.Results.History,'implicit')
        shist = 3; % implicit
elseif strcmpi(p.Results.History,'script')
        shist = 2; % script
elseif strcmpi(p.Results.History,'gui')
        shist = 1; % gui
else
        shist = 0; % off
end
if isempty(bitindex)
        msgboxText = 'Error: you must specify one bit index, at least.';
        title = 'ERPLAB: pop_setcodebit few inputs';
        errorfound(msgboxText, title);
        return
else
        if isnumeric(bitindex)
                if min(bitindex)<1 || max(bitindex)>16
                        msgboxText = 'Error: bit index must be a positive integer between 1 and 16.';
                        title = 'ERPLAB: setcodebitGUI few inputs';
                        errorfound(msgboxText, title);
                        return
                end
        else
                msgboxText = 'Error: bit index must be numeric.';
                title = 'ERPLAB: pop_setcodebit few inputs';
                errorfound(msgboxText, title);
                return
        end
        bitindex = unique_bc2(bitindex);
end
if ischar(newvalue)
        msgboxText = 'Error: new value must be numeric.';
        title = 'ERPLAB: setcodebitGUI few inputs';
        errorfound(msgboxText, title);
        return
else
        if length(newvalue)~=1
                msgboxText = 'Error: new value must be a single value, either 0 or 1.';
                title = 'ERPLAB: pop_setcodebit few inputs';
                errorfound(msgboxText, title);
                return
        end
        if newvalue~=0 && newvalue~=1
                msgboxText = 'Error: new value must be a single value, either 0 or 1.';
                title = 'ERPLAB: pop_setcodebit few inputs';
                errorfound(msgboxText, title);
                return
        end
end

EEG = setcodebit(EEG, bitindex, newvalue);

% History
bitindexstr = vect2colon(bitindex);
newvalue    = vect2colon(newvalue);
com = sprintf('%s = pop_setcodebit( %s, %s, %s);', inputname(1), inputname(1), bitindexstr, newvalue);

% get history from script. EEG
switch shist
        case 1 % from GUI
                com = sprintf('%s %% GUI: %s', com, datestr(now));
                %fprintf('%%Equivalent command:\n%s\n\n', com);
                displayEquiComERP(com);
        case 2 % from script
                EEG = erphistory(EEG, [], com, 1);
        case 3
                % implicit
                % EEG = erphistory(EEG, [], com, 1);
                % fprintf('%%Equivalent command:\n%s\n\n', com);
        otherwise %off or none
                com = '';
                return
end

%
% Completion statement
%
msg2end
return

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function test_failed = test_libltfat_rtdgtreal(varargin)
test_failed = 0;
return;
fprintf(' ===============  %s ================ \n',upper(mfilename));

definput.flags.complexity={'double','single'};
[flags]=ltfatarghelper({},definput,varargin);
dataPtr = [flags.complexity, 'Ptr'];

[~,~,enuminfo]=libltfatprotofile;
LTFAT_FIRWIN = enuminfo.LTFAT_FIRWIN;
rdgt_phasetype = enuminfo.rtdgt_phasetype;

glarr =     [500,  512, 1024, 90];
Marr =      [1000, 512, 2048, 101];
aarr =      [100 , 256,  256, 40];
Warr =      [10  ,   2,    1,  3];

f = greasy;

a = 16;
gl = 60;
gdl = gl;
M = 64;
L = dgtlength(numel(f),a,M);
M2 = floor(M/2) + 1;
g = zeros(gl,1);
gd = zeros(gdl,1);
gPtr = libpointer(dataPtr,g);
gdPtr = libpointer(dataPtr,gd);
funname = makelibraryname('firwin',flags.complexity,0);
calllib('libltfat',funname,LTFAT_FIRWIN.LTFAT_HANN,gl,gPtr);

funname = makelibraryname('gabdual_painless',flags.complexity,0);
prd=calllib('libltfat',funname,gPtr,gl,a,M,gdPtr);
if prd
    warning('This is not painless frame');
    g2Ptr = libpointer(dataPtr,fir2long(gPtr.Value,L));
    gdPtr = libpointer(dataPtr,fir2long(gd,L));
    funname = makelibraryname('gabdual_long',flags.complexity,0);
    calllib('libltfat',funname,g2Ptr,L,1,a,M,gdPtr);
    gdl = L;
end

cframes = signal2frames(f,gl,a,M);
% cframes2 = bsxfun(@times,cframes,M*gPtr.Value.*gdPtr.Value);
% f2 = frames2signal(cframes2,gl,a);
% norm(f-postpad(f2,numel(f)))
N = size(cframes,2);

%ctrue = fftreal(ifftshift(bsxfun(@times,cframes,fftshift(gPtr.Value)),1));
ctrue = dgtreal(f,gPtr.Value,a,M,'timeinv');
p = libpointer();
pinv = libpointer();

funname = makelibraryname('rtdgtreal_init',flags.complexity,0);
calllib('libltfat',funname,gPtr,gl,M,rdgt_phasetype.LTFAT_RTDGTPHASE_ZERO,p);
funname = makelibraryname('rtidgtreal_init',flags.complexity,0);
calllib('libltfat',funname,gdPtr,gdl,M,rdgt_phasetype.LTFAT_RTDGTPHASE_ZERO,pinv);

c = zeros(2*M2,N);
f2 = zeros(gdl,N);
fPtr = libpointer(dataPtr,cframes);
cPtr = libpointer(dataPtr,c);
f2Ptr = libpointer(dataPtr,f2);

calllib('libltfat','ltfat_rtdgtreal_execute_d',p,fPtr,N,cPtr);
norm(ctrue - interleaved2complex(cPtr.Value),'fro')
calllib('libltfat','ltfat_rtidgtreal_execute_d',pinv,cPtr,N,f2Ptr);

frec = frames2signal(f2Ptr.Value,gdl,a);
figure(1);plot([f-postpad(frec,numel(f))]);

fprintf('Coef error. %d \n', norm(ctrue-interleaved2complex(cPtr.Value),'fro'));
fprintf('Rec error. %d \n', norm(f-postpad(frec,numel(f))));

figure(2);plotdgtreal(abs(ctrue)-abs(interleaved2complex(cPtr.Value)),a,M,'linabs');
%imagesc(abs(interleaved2complex(cPtr.Value))./abs(ctrue))

calllib('libltfat','ltfat_rtdgtreal_done_d',p);
calllib('libltfat','ltfat_rtidgtreal_done_d',pinv);

%norm(long2fir(gdtrue,gl) - gdPtr.Value)

function cframes = signal2frames(f,gl,a,M)

[Ls,~] = size(f);
L=dgtlength(Ls,a,M);
N = L/a;
f = postpad(f,L);

cframes = zeros(gl,N,1);

idxrange = -floor(gl/2):ceil(gl/2)-1;

for n=0:N-1
    idx = mod(n*a + idxrange,L) + 1;
    cframes(:,n+1) = f(idx);
end

function f = frames2signal(cframes,gl,a)

N = size(cframes,2);
L = N*a;
f = zeros(L,1);

idxrange = -floor(gl/2):ceil(gl/2)-1;
for n=0:N-1
    idx = mod(n*a + idxrange,L) + 1;
    f(idx) = f(idx) + cframes(:,n+1);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function kl = kldivGaussian(mean1, cov1, mean2, cov2)

% KLDIVGAUSSIAN Give the KL divergence between two Gaussians.
% FORMAT
% DESC returns the Kullback-Leibler divergence between two
% Gaussians with given means and covariances.
% ARG mean1 : mean of the first Gaussian.
% ARG cov1 : covariance of the first Gaussian.
% ARG mean2 : mean of the second Gaussian.
% ARG cov2 : covariance of the second Gaussian.
%
% SEEALSO : logdet, pdinv
%
% COPYRIGHT : Neil D. Lawrence, 2005

% NDLUTIL

[invCov2, U] = pdinv(cov2);
logDet2 = logdet(cov2, U);
logDet1 = logdet(cov1);
N = size(cov1, 1);
meanDiff = mean1 - mean2;
if size(meanDiff, 1) == 1
  meanDiff = meanDiff';
end
kl = -0.5*(logDet1 - logDet2 - trace(cov1*invCov2) ...
           + N - meanDiff'*invCov2*meanDiff);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%  FOLD Apply a binary function left-associatively to a vector
%    FOLD(F, V, INITIALVALUE) iteratively computes
%    RES = F(V(1), V(2))
%    RES = F(RES, V(3))
%    ....
%    RES = F(RES, V(END))
%    and, finally, returns RES. If V is empty, INITIALVALUE is returned.
%    For non-empty V, the third argument INITIALVALUE may be left out.
%
%    Example:
%    syms x
%    fold(@or, x == 1:3) gives x == 1 | x==2 | x==3
%    fold(@and, [], true) gives true
%    fold(@intersect, {[1, 3], [2, 3, 29], [3, 1, 17]}) gives 3
%
%    See also PROD, SUM.
%
%    Reference page in Doc Center
%       doc fold
%
%    Other functions named fold
%
%       symfun/fold    ts/fold
%

author: jmaih


text:
function inds = meshBoundaryEdgeIndices(varargin)
%MESHBOUNDARYEDGEINDICES Indices of boundary edges of a mesh.
%
%   INDS = meshBoundaryVertexIndices(V, F)
%   INDS = meshBoundaryVertexIndices(V, E, F)
%
%   Example
%     % create centered icosahedron
%     [v, f] = createIcosahedron;
%     v(:,3) = v(:,3) - mean(v(:,3));
%     % convert to simili-sphere
%     [v2, f2] = subdivideMesh(v, f, 3);
%     v3 = normalizeVector3d(v2);
%     % clip with plane
%     plane = createPlane([0 0 0], [-1 -2 3]);
%     [vc, fc] = clipMeshVertices(v3, f2, plane, 'shape', 'plane');
%     figure; drawMesh(vc, fc); axis equal; view(3);
%     % draw boundary vertices
%     ec = meshEdges(vc, fc);
%     inds = meshBoundaryEdgeIndices(vc, ec, fc);
%     edges = [vc(ec(inds, 1), :) vc(ec(inds, 2), :)];
%     hold on; drawEdge3d(edges, 'linewidth', 2, 'color', 'b');
%
%   See also
%     meshes3d, meshBoundary, meshBoundaryVertexIndices, meshEdgeFaces

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2019-05-01,    using Matlab 8.6.0.267246 (R2015b)
% Copyright 2019 INRA - Cepia Software Platform.

[vertices, edges, faces] = parseMeshData(varargin{:});

% Compute edge-vertex map if not specified
if isempty(edges)
    edges = meshEdges(vertices, faces);
end

% compute edges to faces map
edgeFaces = meshEdgeFaces(vertices, edges, faces);

inds = find(sum(edgeFaces == 0, 2) > 0);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function jdate = julian (month, day, year)

% Julian date

% Input

%  month = calendar month [1 - 12]
%  day   = calendar day [1 - 31]
%  year  = calendar year [yyyy]

% Output

%  jdate = Julian date

% special notes

%  (1) calendar year must include all digits

%  (2) will report October 5, 1582 to October 14, 1582
%      as invalid calendar dates and stop

% Orbital Mechanics with Matlab

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

y = year;
m = month;
b = 0;
c = 0;

if (m <= 2)
   y = y - 1;
   m = m + 12;
end

if (y < 0)
   c = -.75;
end

% check for valid calendar date

if (year < 1582)
   % null
elseif (year > 1582)
   a = fix(y / 100);
   b = 2 - a + floor(a / 4);
elseif (month < 10)
   % null
elseif (month > 10)
   a = fix(y / 100);
   b = 2 - a + floor(a / 4);
elseif (day <= 4)
   % null
elseif (day > 14)
   a = fix(y / 100);
   b = 2 - a + floor(a / 4);
else
   clc; home;

   fprintf('\n\n  this is an invalid calendar date!!\n');

   keycheck;

   return;
end

jd = fix(365.25 * y + c) + fix(30.6001 * (m + 1));

jdate = jd + day + b + 1720994.5;

author: goGPS-Project


text:
function knee_angle(nlp, bounds)
    % constraints for step length and step width

    domain = nlp.Plant;
    x = domain.States.x;

    % knee angle
    knee = [x('qBRight') - x('qARight')
        x('qBLeft') - x('qALeft')];
    knee_fun = SymFunction(['kneeAngles_',domain.Name], knee, {x});
    addNodeConstraint(nlp, knee_fun, {'x'}, 'all', ...
        bounds.constrBounds.knee.lb, ...
        bounds.constrBounds.knee.ub,'Linear');

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function value = r4_modp ( x, y )

%*****************************************************************************80
%
%% R4_MODP returns the nonnegative remainder of R4 division.
%
%  Discussion:
%
%    If
%      REM = R4_MODP ( X, Y )
%      RMULT = ( X - REM ) / Y
%    then
%      X = Y * RMULT + REM
%    where REM is always nonnegative.
%
%    The MOD function computes a result with the same sign as the
%    quantity being divided.  Thus, suppose you had an angle A,
%    and you wanted to ensure that it was between 0 and 360.
%    Then mod(A,360.0) would do, if A was positive, but if A
%    was negative, your result would be between -360 and 0.
%
%    On the other hand, R4_MODP(A,360.0) is between 0 and 360, always.
%
%  Example:
%
%        I         J     MOD R4_MODP  R4_MODP Factorization
%
%      107        50       7       7    107 =  2 *  50 + 7
%      107       -50       7       7    107 = -2 * -50 + 7
%     -107        50      -7      43   -107 = -3 *  50 + 43
%     -107       -50      -7      43   -107 =  3 * -50 + 43
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the number to be divided.
%
%    Input, real Y, the number that divides X.
%
%    Output, real VALUE, the nonnegative remainder when X is divided by Y.
%
  if ( y == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_MODP - Fatal error!\n' );
    fprintf ( 1, '  R4_MODP ( X, Y ) called with Y = %f\n', y );
    error ( 'R4_MODP - Fatal error!' );
  end

  value = mod ( x, y );

  if ( value < 0.0 )
    value = value + abs ( y );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function E = evaluate_weighted_energy(uv, im, occ, bfhsz, mfsz, sigma_i)

% edge region: weighted median filtering, the weights are determined by
%  spatial distance, intensity distance, occlusion state
% smooth region:

%   Author: Deqing Sun, Department of Computer Science, Brown University
%   Contact: dqsun@cs.brown.edu
%   $Date: 2009$
%   $Revision $
%
% Copyright 2009-2010, Brown University, Providence, RI. USA
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

sigma_x = 7;   %  spatial distance (7)

dilate_sz = 5*[1 1];  % dilation window size for flow edge region [5 5]

sz = size(im);
sz = sz(1:2);

if nargin < 3
    occ = ones(sz);
end;

if nargin < 4
    bfhsz = 10; % half window size
end;

%mfsz = [7 7]; % for test

if nargin < 5
    uvo = uv;
else
%     uvo(:,:,1) = medfilt2(uv(:,:,1), mfsz, 'symmetric');
%     uvo(:,:,2) = medfilt2(uv(:,:,2), mfsz, 'symmetric');
end;

if nargin < 6
    sigma_i = 5; %3.5; % intensity distance (10) 5 better than 10 and 20
end;

% % % WMF first then MF
% % uvo = uv;

e1 = edge(uv(:,:,1), 'sobel');
e2 = edge(uv(:,:,2), 'sobel');
e  = e1|e2;
mask = imdilate(e, ones(dilate_sz) );

% nearest 4 neighbors smaller than [3 3]
% tmp = [ 0 1 0; 1 1 1; 0 1 0];
% mask = imdilate(e, tmp );

% mask = e; % no dilation

% Update non-edge regions first
% uv(repmat(mask, [1 1 2])) = uvo(repmat(mask, [1 1 2]));

% Uncomment below to apply WMF to all regions
% mask = ones(size(im));

% Select boundary regions
[indx_row, indx_col] = find(mask ==1); %

% [H W] = size(im);
pad_u  = padarray(uv(:,:,1), bfhsz*[1 1], 'symmetric', 'both');
pad_v  = padarray(uv(:,:,2), bfhsz*[1 1], 'symmetric', 'both');
pad_im = padarray(im, bfhsz*[1 1], 'symmetric', 'both');
pad_occ= padarray(occ, bfhsz*[1 1], 'symmetric', 'both');

% fprintf('%d\n', length(indx_row));

% Divide into several groups for memory reasons ~70,000 causes out of memory

Indx_Row = indx_row;
Indx_Col = indx_col;
N        = length(Indx_Row); % number of elements to process
n        = 4e4;              % number of elements per batch
nB       = ceil(N/n);

E = 0;

for ib = 1:nB;
    istart = (ib-1)*n + 1;
    iend   = min(ib*n, N);
    indx_row = Indx_Row(istart:iend);
    indx_col = Indx_Col(istart:iend);

    neighbors_u = zeros((bfhsz*2+1)^2, length(indx_row));
    neighbors_v = neighbors_u;
    weights     = neighbors_u;

    for i = 1:length(indx_row)

        % crop window
        r1 = indx_row(i);
        r2 = indx_row(i) + 2*bfhsz;
        c1 = indx_col(i);
        c2 = indx_col(i) + 2*bfhsz;

        rc = indx_row(i) + bfhsz; % row  center
        cc = indx_col(i) + bfhsz; % column center
        ic = pad_im(rc, cc, :);       % intensity of the center pixel

        tmp_u = pad_u(r1:r2, c1:c2);
        tmp_v = pad_v(r1:r2, c1:c2);
        tmp_i = pad_im(r1:r2, c1:c2, :);

        % spatial weight
        [C R] = meshgrid(c1:c2, r1:r2);
        w = exp( -((C-cc).^2+(R-rc).^2)/2/sigma_x^2 );
        % Uncomment below: no spatial weight for test
        % w = ones(size(w));

        % intensity weight; comment below to disable the term
        w = w.* mean( exp(- (tmp_i - repmat(ic, [size(tmp_i,1) size(tmp_i, 2) 1])).^2/2/sigma_i^2), 3);

        % occluded weight;  comment below to disable the term
        w = w.*pad_occ(r1:r2, c1:c2);

        % Normalize
        w = w/sum(w(:));

%         neighbors_u(:,i) = tmp_u(:);
%         neighbors_v(:,i) = tmp_v(:);
%         weights(:,i)     = w(:);

        E = E + sum( abs(tmp_u(:) - pad_u(rc, cc)).* w(:));
        E = E + sum( abs(tmp_v(:) - pad_v(rc, cc)).* w(:));

    end;

end;

% for 5x5 equal weight

% Select non-boundary regions
[indx_row, indx_col] = find(mask ~=1);
bfhsz = floor(mfsz(1)/2);

pad_u  = padarray(uv(:,:,1), bfhsz*[1 1], 'symmetric', 'both');
pad_v  = padarray(uv(:,:,2), bfhsz*[1 1], 'symmetric', 'both');
pad_im = padarray(im, bfhsz*[1 1], 'symmetric', 'both');
pad_occ= padarray(occ, bfhsz*[1 1], 'symmetric', 'both');

% fprintf('%d\n', length(indx_row));

% Divide into several groups for memory reasons ~70,000 causes out of memory

Indx_Row = indx_row;
Indx_Col = indx_col;
N        = length(Indx_Row); % number of elements to process
n        = 4e4;              % number of elements per batch
nB       = ceil(N/n);

for ib = 1:nB;
    istart = (ib-1)*n + 1;
    iend   = min(ib*n, N);
    indx_row = Indx_Row(istart:iend);
    indx_col = Indx_Col(istart:iend);

    neighbors_u = zeros((bfhsz*2+1)^2, length(indx_row));
    neighbors_v = neighbors_u;
    weights     = neighbors_u;

    for i = 1:length(indx_row)

        % crop window
        r1 = indx_row(i);
        r2 = indx_row(i) + 2*bfhsz;
        c1 = indx_col(i);
        c2 = indx_col(i) + 2*bfhsz;

        rc = indx_row(i) + bfhsz; % row  center
        cc = indx_col(i) + bfhsz; % column center
        ic = pad_im(rc, cc, :);       % intensity of the center pixel

        tmp_u = pad_u(r1:r2, c1:c2);
        tmp_v = pad_v(r1:r2, c1:c2);
        tmp_i = pad_im(r1:r2, c1:c2, :);

        % spatial weight
        [C R] = meshgrid(c1:c2, r1:r2);
        w = exp( -((C-cc).^2+(R-rc).^2)/2/sigma_x^2 );
        % Uncomment below: no spatial weight for test
        w = ones(size(w));

%         % intensity weight; comment below to disable the term
%         w = w.* mean( exp(- (tmp_i - repmat(ic, [size(tmp_i,1) size(tmp_i, 2) 1])).^2/2/sigma_i^2), 3);
%
%         % occluded weight;  comment below to disable the term
%         w = w.*pad_occ(r1:r2, c1:c2);

        % Normalize
        w = w/sum(w(:));

%         neighbors_u(:,i) = tmp_u(:);
%         neighbors_v(:,i) = tmp_v(:);
%         weights(:,i)     = w(:);

        E = E + sum( abs(tmp_u(:) - pad_u(rc, cc)).* w(:));
        E = E + sum( abs(tmp_v(:) - pad_v(rc, cc)).* w(:));

    end;

end;

author: kristinbranson

repo: JAABA


text:
function [sumdetector,sumphi,sumspectr,peaks,bgpeaks] = proceed_spectra(spec,bg,range,peakpositions)
% procede Dubna spectra
%
% Input
%  spec  - spectra
%  bg    - background
%  range -
%  peakpositions - double
%
% Output
%  sumdetector - sum over all detectors (range x 72)
%  sumphi      - sum over all phi (range x 19)
%  sumspectr   - sum over all spectra (72 x 19)
%  peaks       - peak sums (peaks x 72 x 19)
%

coeff = fclencurt(19+18,0,1);
sumdetector = reshape(reshape(spec,[],size(spec,3))*coeff(1:19),size(spec,1),size(spec,2));
sumphi = sum(spec,2);
spec = spec - bg;
sumspectr = squeeze(sum(bg(range,:,:),1));
sumspectr = sumspectr ./ mean(sumspectr(:));
peaks = zeros(length(peakpositions),size(spec,2),size(spec,3));
for p = 1:size(peaks,1)
  peaks(p,:) = sum(spec(peakpositions(p,1):peakpositions(p,2),:),1);
end

bgpeaks = zeros(length(peakpositions),size(spec,2),size(spec,3));
for p = 1:size(peaks,1)
  bgpeaks(p,:) = sum(bg(peakpositions(p,1):peakpositions(p,2),:),1);
end

author: mtex-toolbox

repo: mtex

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


text:
function imden=den2d_block(x,qmf,thd,L,sigma)

if (exist('x')~=1), 	error('Provide input image'); 		end;
if (exist('qmf')~=1), 	qmf = MakeONFilter('Symmlet',6); 	end;
if (exist('thd')~=1), 	thd = 4.5;				end;

[nx ny] = size(x);
if nx~=ny
	disp('Matrix must be square');
	return;
end

n 	= nx;
J	= nextpow2(n);
if (exist('L')~=1), 	L   = floor(sqrt(2*J)); 		end; % Size of block.
scale 	= floor(log2(L));   % Coarsest scale L#2^Jc matches the size of the block.

% wc=FWT2_PO(x,scale,qmf);
wc = x;

for j = 0:J-1-scale
 HH = wc(2^(J-j-1)+1 : 2^(J-j), 2^(J-j-1)+1 : 2^(J-j));
 HH = block_partition(HH,L);
 %HHmask = sqrt(sum(HH.^2)) >= thd*L*sigma;
 HHmask = max(1-thd*log(n)*sigma^2./sum(HH.^2),0);
 HH = HH .* repmat(HHmask,L*L,1);
 HH = invblock_partition(HH,2^(J-j-1),L);
 wc(2^(J-j-1)+1 : 2^(J-j), 2^(J-j-1)+1 : 2^(J-j)) = HH;

 HL = wc(2^(J-j-1)+1 : 2^(J-j), 1 : 2^(J-j-1));
 HL = block_partition(HL,L);
 %HLmask = sqrt(sum(HL.^2)) >= thd*L*sigma;
 HLmask = max(1-thd*log(n)*sigma^2./sum(HL.^2),0);
 HL = HL .* repmat(HLmask,L*L,1);
 HL = invblock_partition(HL,2^(J-j-1),L);
 wc(2^(J-j-1)+1 : 2^(J-j), 1 : 2^(J-j-1)) = HL;

 LH = wc(1 : 2^(J-j-1), 2^(J-j-1)+1 : 2^(J-j));
 LH = block_partition(LH,L);
 %LHmask = sqrt(sum(LH.^2)) >= thd*L*sigma;
 LHmask = max(1-thd*log(n)*sigma^2./sum(LH.^2),0);
 LH = LH .* repmat(LHmask,L*L,1);
 LH = invblock_partition(LH,2^(J-j-1),L);
 wc(1 : 2^(J-j-1), 2^(J-j-1)+1 : 2^(J-j)) = LH;
end

imden = IWT2_PO(wc,scale,qmf);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = block_partition(x,L)

n 	= length(x);
nblocks = floor(n/L);

out = [];
for i=0:nblocks-1
  out = [out x(i*L+1:(i+1)*L,:)];
end

out = reshape(out(:),L^2,nblocks^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = invblock_partition(x,n,L)

nblocks = floor(n/L);

buf = reshape(x,L,L*nblocks^2);

out = [];
for i=0:nblocks-1
  out = [out;buf(:,i*L*nblocks+1:(i+1)*L*nblocks)];
end

author: gpeyre

repo: matlab-toolboxes


text:
function theta = deg2rad(omega)

% DEG2RAD Transform degrees to radians.
% FORMAT
% DESC onverts degrees to radians.
% ARG omega : angle in degrees.
% ARG theta : angle in radians.
%
% COPYRIGHT : Neil D. Lawrence, 2005

% NDLUTIL

theta = omega/180*pi;


text:
function average = spm_cfg_eeg_average
% Configuration file for M/EEG epoching
%__________________________________________________________________________
% Copyright (C) 2008-2016 Wellcome Trust Centre for Neuroimaging

% Stefan Kiebel
% $Id: spm_cfg_eeg_average.m 6926 2016-11-09 22:13:19Z guillaume $

rev = '$Rev: 6926 $';
D = cfg_files;
D.tag = 'D';
D.name = 'File Name';
D.filter = 'mat';
D.num = [1 1];
D.help = {'Select the M/EEG mat file.'};

standard = cfg_const;
standard.tag = 'standard';
standard.name = 'Standard';
standard.val  = {false};
standard.help = {''};

ks = cfg_entry;
ks.tag = 'ks';
ks.name = 'Offset of the weighting function';
ks.strtype = 'r';
ks.val = {3};
ks.num = [1 1];
ks.help = {'Parameter determining the how far the values should be from the median, '...
    'to be considered outliers (the larger, the farther).'};

bycondition = cfg_menu;
bycondition.tag = 'bycondition';
bycondition.name = 'Compute weights by condition';
bycondition.help = {'Compute weights for each condition separately or for all conditions together.'};
bycondition.labels = {'Yes', 'No'};
bycondition.values = {true, false};
bycondition.val = {false};

savew = cfg_menu;
savew.tag = 'savew';
savew.name = 'Save weights';
savew.help = {'Save weights in a separate dataset for quality control.'};
savew.labels = {'Yes', 'No'};
savew.values = {true, false};
savew.val = {false};

removebad = cfg_menu;
removebad.tag = 'removebad';
removebad.name = 'Remove bad data';
removebad.help = {'Replace data marked as bad by NaNs before averaging',...
    'Warning: beware if there is no good data, NaNs may stay in the average.'};
removebad.labels = {'Yes', 'No'};
removebad.values = {true, false};
removebad.val = {false};

robust = cfg_branch;
robust.tag = 'robust';
robust.name = 'Robust';
robust.val = {ks, bycondition, savew, removebad};
robust.help = {''};

userobust = cfg_choice;
userobust.tag = 'userobust';
userobust.name = 'Averaging type';
userobust.help = {'choose between using standard and robust averaging.'};
userobust.values = {standard, robust};
userobust.val = {standard};

plv = cfg_menu;
plv.tag = 'plv';
plv.name = 'Compute phase-locking value';
plv.help = {'Compute phase-locking value rather than average the phase',...
    'This option is only relevant for TF-phase datasets.'};
plv.labels = {'Yes', 'No'};
plv.values = {true, false};
plv.val = {false};

prefix         = cfg_entry;
prefix.tag     = 'prefix';
prefix.name    = 'Filename Prefix';
prefix.help    = {'Specify the string to be prepended to the filenames of the averaged dataset. Default prefix is ''m''.'};
prefix.strtype = 's';
prefix.num     = [1 Inf];
prefix.val     = {'m'};

average = cfg_exbranch;
average.tag = 'average';
average.name = 'Averaging';
average.val = {D, userobust, plv, prefix};
average.help = {'Average epoched EEG/MEG data.'};
average.prog = @eeg_average;
average.vout = @vout_eeg_average;
average.modality = {'EEG'};

function out = eeg_average(job)
% construct the S struct
S.D = job.D{1};
if isfield(job.userobust, 'robust')
    S.robust = job.userobust.robust;
else
    S.robust = false;
end

S.circularise = job.plv;
S.prefix = job.prefix;

out.D = spm_eeg_average(S);
out.Dfname = {fullfile(out.D.path, out.D.fname)};

function dep = vout_eeg_average(job)
% Output is always in field "D", no matter how job is structured
dep = cfg_dep;
dep.sname = 'Average Data';
% reference field "D" from output
dep.src_output = substruct('.','D');
% this can be entered into any evaluated input
dep.tgt_spec   = cfg_findspec({{'strtype','e'}});

dep(2) = cfg_dep;
dep(2).sname = 'Averaged Datafile';
% reference field "Dfname" from output
dep(2).src_output = substruct('.','Dfname');
% this can be entered into any file selector
dep(2).tgt_spec   = cfg_findspec({{'filter','mat'}});

author: spm


text:
% POP_ICATHRESH - main menu for choosing threshold for component
%                      rejection in EEGLAB.
%
% Usage:
%   >> [OUTEEG rej] = pop_icathresh(INEEG, threshval, rejmethod,
%                                      rejvalue, interact);
%
% Inputs:
%   INEEG      - input dataset
%   threshval  - values of thresholds for each of the 3 statistical
%                measures. Default is [] and the program uses the value
%                in the dataset.
%   rejmethod  - either 'percent', 'dataset' or 'current'. 'percent'
%                will reject a given percentage of components with
%                the highest value in one or several statistical
%                measure. 'dataset' will use an other dataset for
%                calibration. 'current' will use the current dataset
%                for calibration. Default is 'current'.
%   rejvalue   - percentage if rejmethod is 'percent', dataset number
%                if rejmethod is 'dataset' (no input if 'current'). If it
%                is a percentage, '25' for instance means that the 25%
%                independent components with the highest values of one
%                statistical measure are  rejected. Note that, one can
%                also enter one percentage per statistical value (such as
%                25 20 30).
%   interact   - interactive windows or just rejection
%
% Inputs:
%   OUTEEG     - output dataset with updated thresholds
%   rej        - new rejection array
%
% Graphic interface:
%   The graphic interface is divided into 3 parts. On the top, the
% experimenter chooses the basis for the component rejection (see
% rejmethod input). On the middle panel, the user can tune thresholds
% manually and plot the distribution of statistical values. Each time
% that setting of one of these two top panels is modified, the curve
% on the third panel are redrawn.
%   The third panel is composed of 3 graphs, one for each statistical
% measure. The blue curve on each graph indicates the accuracy of the
% measure to detect artifactual components and non-artifactual
% components of a dataset. Note that 'artifactual components are
% defined with the rejection methods fields on the top (if you use
% the percentage methods, these artifactual components may not actually
% be artifactual). For accurate rejection, one should first reject
% artifactual component manually and then set up the threshold to
% approximate this manual rejection. The green curve indicate the
% rejection when all measure are considered. Points on the blue and
% on the green curves indicate the position of the current threshold on
% the parametric curve. Not that for the green cumulative curve, this
% point is at the same location for all graphs.
%   How to tune the threshold? To tune the threshold, one should try to
% maximize the detection of non-artifactual components (because it is
% preferable to miss some artifactual components than to classify as
% artifactual non-artifact components).
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: EEGLAB

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

%PROBLEM: when the % is set and we change manually the threshold,
% the software comes back to the percentage rejection

% 01-25-02 reformated help & license -ad

function [EEG, rej, com] = pop_icathresh( EEG, threshval, rejmethod, rejvalue, interact);

com = [];
rej = EEG.reject.gcompreject;
if nargin < 1
	help pop_icathresh;
	return;
end

if nargin < 2
	threshval = [];
end
if nargin < 3
	rejmethod = 'current';
end
if nargin < 4
	rejvalue = 25;
end
if nargin < 5
	interact = 1;
end
if ~isempty(threshval)
	EEG.reject.threshentropy = threshval(1);
	EEG.reject.threshkurtact = threshval(2);
	EEG.reject.threshkurtdist = threshval(3);
end

tagmenu = 'pop_icathresh';

if ~isempty( findobj('tag', tagmenu))
	error('cannot open two identical windows, close the first one first');
end

% the current rejection will be stored in userdata of the figure
% --------------------------------------------------------------
gcf = figure('visible', 'off', 'numbertitle', 'off', 'name', 'Choose thresholds', 'tag', tagmenu, 'userdata', rej);
pos = get(gca, 'position');
q = [pos(1) pos(2) 0 0];
s = [pos(3) pos(4) pos(3) pos(4)]./100;
axis off;

% definition of callbacks
% -----------------------

cb_cancel  = [  'userdat = get(gcbo, ''userdata'');' ... % restore thresholds
				'EEG.reject.threshentropy  = userdat{1};' ...
				'EEG.reject.threshkurtact  = userdat{2};' ...
				'EEG.reject.threshkurtdist = userdat{3};' ...
				'clear userdat; close(gcbf);' ];

drawgraphs  = [ 'if isempty( gcbf ), fig = gcf; else fig = gcbf; end;' ... % determine the figure handler
                'if get( findobj(''parent'', fig, ''tag'', ''Ipercent''), ''value'') == 1, ' ... % test if percentage
				'   perc = str2num(get( findobj(''parent'', fig, ''tag'', ''Ipercenttext''), ''String''));' ...
				'   if length(perc < 2), perc = [perc perc perc]; end;' ...
				'   perc = round((100-perc) * length( EEG.stats.compenta) / 100);' ... % convert the percentage
				'   method = zeros( size( EEG.stats.compenta ) );' ...
				'   [tmprej tmpindex] = sort(EEG.stats.compenta(:));' ...
				'   method( tmpindex(perc(1)+1:end) ) = 1;' ...
				'   set( findobj(''parent'', fig, ''tag'', ''entstring''), ''string'', num2str(tmprej(perc(1)+1)));' ... % update threshold
				'   [tmprej tmpindex] = sort(EEG.stats.compkurta(:));' ...
				'   method( tmpindex(perc(2)+1:end) ) = 1;' ...
				'   set( findobj(''parent'', fig, ''tag'', ''kurtstring''), ''string'', num2str(tmprej(perc(1)+1)));' ... % update threshold
				'   [tmprej tmpindex] = sort(EEG.stats.compkurtdist(:));' ...
				'   method( tmpindex(perc(3)+1:end) ) = 1;' ...
				'   set( findobj(''parent'', fig, ''tag'', ''kurtdiststring''), ''string'', num2str(tmprej(perc(1)+1)));' ... % update threshold
				'   allvalues = [EEG.stats.compenta(:) EEG.stats.compkurta(:) EEG.stats.compkurtdist(:) ];' ...
				'   clear perc tmprej tmpindex;' ...
                'end;' ...
				'if get( findobj(''parent'', fig, ''tag'', ''Iother''), ''value'') == 1, ' ... % test if other dataset
				'   di = str2num(get( findobj(''parent'', fig, ''tag'', ''Iothertext''), ''String''));' ...
				'   if isempty( di ), clear fig; return; end;' ...
				'   method    = ALLEEG(di).reject.gcompreject'';' ...
				'   allvalues = [ALLEEG(di).stats.compenta(:) ALLEEG(di).stats.compkurta(:) ALLEEG(di).stats.compkurtdist(:) ];' ...
				'   clear di;' ...
                'end;' ...
				'if get( findobj(''parent'', fig, ''tag'', ''Icurrent''), ''value'') == 1, ' ... % test if current dataset
				'   method = EEG.reject.gcompreject'';' ...
				'   allvalues = [EEG.stats.compenta(:) EEG.stats.compkurta(:) EEG.stats.compkurtdist(:) ];' ...
                'end;' ...
				'axes( findobj( ''parent'', fig, ''tag'', ''graphent'')); cla;' ...
                '[tmp1 tmp2 tmp3] = drawproba( method, allvalues, [0 EEG.reject.threshkurtact EEG.reject.threshkurtdist ], ' ...
				'   { ''&'' ''|'' }, EEG.reject.threshentropy, ''Entropy of activity'',''Artifact detection (%)'', ''Non-artifact detection (%)'');' ...
				'set(gca, ''tag'', ''graphent'');' ...
				'axes( findobj( ''parent'', fig, ''tag'', ''graphkurt'')); cla;' ...
                '[tmp1 tmp2 tmp3] = drawproba( method, allvalues, [EEG.reject.threshentropy 0 EEG.reject.threshkurtdist ], ' ...
				'   { ''&'' ''|'' }, EEG.reject.threshkurtact, ''Kurtosis of activity'', ''Artifact detection (%)'', '''');' ...
				'set(gca, ''tag'', ''graphkurt'');' ...
				'axes( findobj( ''parent'', fig, ''tag'', ''graphkurtdist'')); cla;' ...
                '[tmp1 tmp2 tmp3] = drawproba( method, allvalues, [EEG.reject.threshentropy EEG.reject.threshkurtact 0 ], ' ...
				'   { ''&'' ''|'' }, EEG.reject.threshkurtdist, ''Kurtosis of topography'', ''Artifact detection (%)'', '''');' ...
				'set(gca, ''tag'', ''graphkurtdist'');' ...
				'clear method allvalues fig;' ...
				];

cb_percent = [  'if isempty( gcbf ), fig = gcf; else fig = gcbf; end;' ... % determine the figure handler
                'set( findobj(''parent'', fig, ''tag'', ''Ipercent''), ''value'', 1);'...
				'set( findobj(''parent'', fig, ''tag'', ''Iother''), ''value'', 0);' ...
				'set( findobj(''parent'', fig, ''tag'', ''Icurrent''), ''value'', 0);' ...
				'set( findobj(''parent'', fig, ''tag'', ''Ipercenttext''), ''enable'', ''on'');' ...
				'set( findobj(''parent'', fig, ''tag'', ''Iothertext''), ''enable'', ''off''); clear fig;' drawgraphs ];

cb_other   = [  'if isempty( gcbf ), fig = gcf; else fig = gcbf; end;' ... % determine the figure handler
                'set( findobj(''parent'', fig, ''tag'', ''Iother''), ''value'', 1);'...
				'set( findobj(''parent'', fig, ''tag'', ''Ipercent''), ''value'', 0);' ...
				'set( findobj(''parent'', fig, ''tag'', ''Icurrent''), ''value'', 0);' ...
				'set( findobj(''parent'', fig, ''tag'', ''Ipercenttext''), ''enable'', ''off'');' ...
				'set( findobj(''parent'', fig, ''tag'', ''Iothertext''), ''enable'', ''on''); clear fig;' drawgraphs ];

cb_current = [  'if isempty( gcbf ), fig = gcf; else fig = gcbf; end;' ... % determine the figure handler
                'set( findobj(''parent'', fig, ''tag'', ''Icurrent''), ''value'', 1);'...
				'set( findobj(''parent'', fig, ''tag'', ''Iother''), ''value'', 0);' ...
				'set( findobj(''parent'', fig, ''tag'', ''Ipercent''), ''value'', 0);' ...
				'set( findobj(''parent'', fig, ''tag'', ''Ipercenttext''), ''enable'', ''off'');' ...
				'set( findobj(''parent'', fig, ''tag'', ''Iothertext''), ''enable'', ''off''); clear fig;' drawgraphs ];

cb_entropy  = [ 'EEG.reject.threshentropy = str2num(get(gcbo, ''string''));' ...
				drawgraphs ];
cb_kurtact  = [ 'EEG.reject.threshkurtact = str2num(get(gcbo, ''string''));' ...
				drawgraphs ];
cb_kurtdist  = [ 'EEG.reject.threshkurtdist = str2num(get(gcbo, ''string''));' ...
				drawgraphs ];

if interact
	cb_calrej = [ 'ButtonName=questdlg2( ''This will erase previous projections'', ''Confirmation'', ''CANCEL'', ''OK'', ''OK'');' ]
else
	cb_calrej = [ 'ButtonName= ''OK''' ];
end
cb_calrej = [ cb_calrej ...
    		  'switch ButtonName,' ...
       	 	  '   case ''OK'',' ...
			  '    rej1 = find(EEG.stats.compenta > EEG.reject.threshentropy);' ...
			  '    rej2 = find(EEG.stats.compkurta > EEG.reject.threshkurtact);' ...
			  '    rej3 = find(EEG.stats.compkurtdist > EEG.reject.threshkurtdist);' ...
			  '    EEG.reject.gcompreject = (rej1 & rej2) | rej3;' ...
			  '    clear rej1 rej2 rej3;' ...
			  'end; clear ButtonName;' ];

% default value for rejection methods
% -----------------------------------
rejvalother = '';
rejvalpercent = '25';
switch rejmethod
	case 'percent', rejvalpercent = num2str( rejvalue);
	case 'dataset', rejvalother = num2str( rejvalue);
end

% -----------------------------------------------------
allh = supergui(gcf, { 	[1] ...
					[2 1] [2 1] [2 1] ...
					[1] ...
					[1] ...
					[1.5 0.5 1] [1.5 0.5 1] [1.5 0.5 1] ...
					[1] [1] [1] [1] [1] [1] [1] [1] [1] [1] ...
					[1 1 1 1 1] ...
					}, [], ...
	{ 'Style', 'text', 'string', 'Calibration method', 'FontSize', 13, 'fontweight', 'bold' }, ...
	...
	{ 'style', 'checkbox', 'String', '% of artifactual components (can also put one % per rejection)', 'tag', 'Ipercent', 'value', 1, 'callback', cb_percent}, ...
	{ 'style', 'edit', 'String', rejvalpercent, 'tag', 'Ipercenttext', 'callback', drawgraphs }, ...
	...
	{ 'style', 'checkbox', 'String', 'Specific dataset (enter dataset number)', 'tag', 'Iother', 'value', 0, 'callback', cb_other}, ...
	{ 'style', 'edit', 'String', rejvalother, 'tag', 'Iothertext', 'enable', 'off', 'callback', drawgraphs }, ...
	...
	{ 'style', 'checkbox', 'String', 'Current dataset', 'tag', 'Icurrent', 'value', 0, 'callback', cb_current}, ...
	{ }, ...
	...
	{ }, ...
	...
	{ 'Style', 'text', 'string', 'Threshold values', 'FontSize', 13, 'fontweight', 'bold' }, ...
	...
	{ 'style', 'text', 'String', 'Entropy of activity threshold' }, ...
	{ 'style', 'edit', 'String', num2str(EEG.reject.threshentropy), 'tag', 'entstring', 'callback',  cb_entropy }, ...
	{ 'style', 'pushbutton', 'String', 'Histogram', 'callback', 'figure; hist(EEG.stats.compenta,20); title(''Entropy of activity'');' }, ...
	...
	{ 'style', 'text', 'String', 'Kurtosis of activity threshold' }, ...
	{ 'style', 'edit', 'String', num2str(EEG.reject.threshkurtact), 'tag', 'kurtstring', 'callback', cb_kurtact }, ...
	{ 'style', 'pushbutton', 'String', 'Histogram', 'callback', 'figure; hist(EEG.stats.compkurta,200); title(''Kurtosis of activity'');' },  ...
	...
	{ 'style', 'text', 'String', 'Kurtosis of topography threshold' }, ...
	{ 'style', 'edit', 'String', num2str(EEG.reject.threshkurtdist), 'tag', 'kurtdiststring', 'callback', cb_kurtdist }, ...
	{ 'style', 'pushbutton', 'String', 'Histogram', 'callback', 'figure; hist(EEG.stats.compkurtdist,20); title(''Kurtosis of topography'');' }, ...
	...
	{ }, { }, { }, { }, { }, { }, { }, { }, { }, { }, ...
	...
	{ 'style', 'pushbutton', 'String', 'Cancel', 'callback', cb_cancel, 'userdata', { EEG.reject.threshentropy EEG.reject.threshkurtact EEG.reject.threshkurtdist }}, ...
	{ 'style', 'pushbutton', 'String', 'Auto thresh'  , 'callback', '', 'enable', 'off' }, ...
	{ 'style', 'pushbutton', 'String', 'Help'  , 'callback', 'pophelp(''pop_icathresh'');' }, ...
	{ 'style', 'pushbutton', 'String', 'Accept thresh.' , 'callback', 'close(gcbf);' }, ...
	{ 'style', 'pushbutton', 'String', 'Calc. rejection' , 'callback', 'close(gcbf);' } ...
	...
	);

h = axes('position', [0  15 30 30].*s+q, 'tag', 'graphent');
h = axes('position', [35 15 30 30].*s+q, 'tag', 'graphkurt');
h = axes('position', [70 15 30 30].*s+q, 'tag', 'graphkurtdist');

rejmethod
switch rejmethod
	case 'dataset', eval([ 'gcbf = [];' cb_other]);
	case 'current', eval([ 'gcbf = [];' cb_current]);
	otherwise, eval([ 'gcbf = [];' cb_percent]);
end

return;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
% Test file for @deltafun/chebcoeffs.m.

function pass = test_chebcoeffs(pref)

if (nargin < 1)
    pref = chebfunpref();
end

a = -4; b = 4;

f = fun.constructor(@(x) sin(x), struct('domain', [a, b]));
mag = .9*(a + (b-a)*rand(3,3));
loc = .9*(a + (b-a)*rand(1,3));

df = deltafun(f, struct('deltaMag', mag, 'deltaLoc', loc));

pass(1) = all( chebcoeffs(f) == chebcoeffs(df) );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [  ] = gsp_plot_sgram( G,A,param )
%GSP_PLOT_SGRAM Plot graph spectrogram
%   Usage:  gsp_plot_sgram( G,A );
%           gsp_plot_sgram( G,A,param );
%
%   Input parameters:
%         G     : Graph
%         A     : Graph windowed Fourrier transform
%         param : Structure of optional parameter
%   Output parameters:
%         none
%
%   *param* is a structure of optional parameter with
%
%   * *param.colorbar*: Use the colorbar (default 1)
%
%   Example:::
%
%           N = 15;
%           G = gsp_ring(2*N);
%           G = gsp_compute_fourier_basis(G);
%           x = [0:N,(N-1):-1:1]';
%           s = 3;
%           g = exp(-(x-1).^2/s^2);
%           f = gsp_localize(G,g,N);
%           c = gsp_gwft(G,f,g);
%           gsp_plot_sgram(G,c);
%
%   See also: gsp_plot_signal, gsp_plot_graph, gsp_plot_signal_spectral
%

% Author: Nathanael Perraudin
% Date  : 09.12.2013
% testing: test_plotting

% Optional parameter handling
if nargin<3
    param=struct;
end

if ~isfield(param, 'colorbar'), param.colorbar = 1; end;

imagesc(1:size(A,2), 0:size(A,1)-1,abs( A));

% Hack to overpass a matlab bug with latex interpretex
latex = get(gca,'DefaultTextInterpreter');
set(gca,'DefaultTextInterpreter','Tex');

xlabel('Nodes');
ylabel('Freqencies');

set(gca,'DefaultTextInterpreter',latex);

if param.colorbar
    colorbar
end

end

repo: gspbox


text:
function [MOIndex] = lipschitz(u, y, maxlag, model, fig)
% A method to determine the lag space, based on Lipschitz quotients
%
%% Syntax
% [MOIndex] = lipschitz(u, y, maxlag)
%
%% Description
%  Given a set of corresponding inputs and outputs the function calculates
%  so called Lipschitz number for each combination of m and l where m
%  represents the number of delayed outputs, and l the number of delayed
%  inputs for the case of dynamic system: y(t) = f(y(t-1),...,y(t-l),
%  u(t-1),..., u(t-m)).
%  To small m and l result in a large Lipschitz nuber, while to large lag
%  spaces do not have greater effect on the Lipschitz nuber. In order to
%  determine the proper lag space we have to look for the knee point, where
%  Lipscitz number stops decreasing.
%
% Input:
% * u       ... the system input(column vector)
% * y       ... the system output(column vector)
% * maxlag  ... the max lag space to investigate
% * model   ... optional - 'arx' if we only want to investigate m = l case
%
% Output:
% * MOIndex ... the maxlag by maxlag matrix containing calculated Lipscitz
%               numbers (Model order index) for each combination of m and l

%% Signature
% Written by Tomaz Sustar
% Based on the algorithm by Xiangdong He and Haruhiko Asada

if(nargin<4), model='unknown'; end;

NN = length(y);                                         % number of samples
MOIndex = zeros(maxlag);                    % matrix of lipschitz's indexes

for m=1:maxlag,                                 % number of delayed outputs
  % m,
  for l=0:maxlag,                                 % number of delayed inputs

    % if we are investigating arx model srtucture only, we calculate MOIndex only when m = l
    if(strcmp(model, 'arx') && l ~= m), continue; end;

    lag = max(l,m);       % the greater from m, l
    % Because of the lag we can construct only NN - lag input output pairs
    N    = NN-lag;    % number of input - output pairs.
    p    = floor(0.02*N); % number of Lipschitz quotients used to determine model order index

    [input target] = construct([m l], u, y); % construct regressors and target

    % calculation of Lipschitz quotients

    Q = zeros(N);      % initialize Q matrix for storing Lipschitz qotients

    for i=1:N-1,
      % for each input/output pair calculate the their Lipschitz quotients all
      % further inputs/outputs pairs. In this way all possible Lipschitz
      % quotients q(i,j) are calculated.

      Q(i,i+1:N)=(target(i)-target(i+1:N)).^2 ./ ...
       sum((repmat(input(i,:), N-i, 1)-input(i+1:N,:)).^2, 2);

    end

    Q_max = Q(Q~=0);                                         % remove zeros
    Q_max = (-sort(-Q_max(:)));        % sort qoutients in descending order
    Q_max = sqrt(Q_max(1:p));                  % take p - largest quotients

    n = m+l;
    MOIndex(m,l+1)=prod(sqrt(n)*Q_max)^(1/p); % calculates order index and stores it to the matrix

  end % end for l

end % end for m

% draw some figures

if(~strcmp(model, 'arx'))

  figure('Name', 'Model order index vs. lag space')
  surf(1:maxlag, 0:maxlag,MOIndex');
  view([-600 40]);
  set(gca, 'Zscale','log');
  set(gca, 'XTick', 1:maxlag)
  set(gca, 'XTick', 1:maxlag)
  xlabel('l - number of past outputs')
  ylabel('m - number of past inputs')
  zlabel('Model Order Index')
end

if(nargin > 4)
  figure(fig);
else
  figure('Name', 'Model order index vs. lag space - arx case')
end
semilogy(diag(MOIndex));
xlabel('m = l - number of past inputs and outputs');
ylabel('Model order index');
set(gca, 'XTick', 1:maxlag);
grid on;

author: Dynamic-Systems-and-GP

repo: GPdyn


text:
function [M_c,C_c_nu_c,g_c,Jc_c,dJc_nu_c,nu_c] = fromFloatingToCentroidalDynamics(M, h, g, Jc, dJc_nu, nu, T, dT)
%FROMFLOATINGTOCENTROIDALDYNAMICS  converts dynamic equation parameters to the
%                       corresponding values in centroidal frame of reference.
%
% [M_c,C_c_nu_c,g_c,Jc_c,dJc_nu_c,nu_c] = FROMFLOATINGTOCENTROIDALDYNAMICS
% (M, h, g, Jc, dJc_nu, nu, T, dT) takes as an input the robot dynamics.
%  The output is the robot dynamics in centroidal frame of reference.
%

% ------------Initialization----------------
ndof   = size(g,1)-6;
invT   = eye(ndof+6)/T;
invTt  = eye(ndof+6)/(T');

%% Control terms conversion
% mass matrix
M_c            = invTt*M*invT;
M_c(1:6,7:end) = zeros(6,ndof);
M_c(7:end,1:6) = zeros(ndof,6);
M_c(1:3,1:3)   = M(1,1)*eye(3);
M_c(1:3,4:6)   = zeros(3);
M_c(4:6,1:3)   = zeros(3);
Mb             = M(1:6,1:6);
Mbj            = M(1:6,7:end);

nu_c           = T*nu;
gravAcc        = norm(invTt*g)/M(1,1);
e3             = zeros(ndof+6,1);
e3(3)          = 1;
g_c            = M(1,1)*gravAcc*e3;

%coriolis terms
C_nu           = h - g;
C_nu_j         = C_nu(7:end);
C_nu_b         = C_nu(1:6);
C_c_nu_c_dT    = invTt*C_nu - M_c*dT*nu;
C_c_nu_c       = [ zeros(3,1);
                  C_c_nu_c_dT(4:6);
                  C_nu_j-(Mbj')*(Mb\C_nu_b)];

%new dT*nu computation for Jacobian
dT_nu          = M_c\(C_nu-C_c_nu_c); % M_c\(inv(T)*C_nu-C_c_nu_c);
Jc_c           = Jc*invT;
dJc_nu_c       = dJc_nu - Jc*invT*dT_nu;

end

author: robotology

repo: whole-body-controllers


text:
function engine = kalman_inf_engine(bnet)
% KALMAN_INF_ENGINE Inference engine for Linear-Gaussian state-space models.
% engine = kalman_inf_engine(bnet)
%
% 'onodes' specifies which nodes are observed; these must be leaves.
% The remaining nodes are all hidden. All nodes must have linear-Gaussian CPDs.
% The hidden nodes must be persistent, i.e., they must have children in
% the next time slice. In addition, they may not have any children within the current slice,
% except to the observed leaves. In other words, the topology must be isomorphic to a standard LDS.
%
% There are many derivations of the filtering and smoothing equations for Linear Dynamical
% Systems in the literature. I particularly like the following
% - "From HMMs to LDSs", T. Minka, MIT Tech Report, (no date), available from
%    ftp://vismod.www.media.mit.edu/pub/tpminka/papers/minka-lds-tut.ps.gz

[engine.trans_mat, engine.trans_cov, engine.obs_mat, engine.obs_cov, engine.init_state, engine.init_cov] = ...
    dbn_to_lds(bnet);

% This is where we will store the results between enter_evidence and marginal_nodes
engine.one_slice_marginal = [];
engine.two_slice_marginal = [];

engine = class(engine, 'kalman_inf_engine', inf_engine(bnet));

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function meshdemo_2d ( )

%*****************************************************************************80
%
%% MESHDEMO_2D displays some 2D examples of the use of DISTMESH.
%
%  Licensing:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Modified:
%
%    27 October 2011
%
%  Reference:
%
%    Per-Olof Persson, Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Local parameters:
%
%    Local, integer ITERATION_MAX, the maximum number of iterations that
%    DISTMESH should take.  (The program might take fewer iterations if it
%    detects convergence.)
%
%    Local, real H, the desired initial spacing.  If a uniform mesh density is used,
%    this is the approximate spacing throughout the region.
%
%    Local, pointer FH, an inline formula, or the name of an M file, which calculates
%    the mesh density function.
%
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'MESHDEMO_2D\n' );
  fprintf ( 1, '  Demonstrations of the use of DISTMESH.\n' );
%
%  Problem 1, the circle, with spacings H = 0.4, 0.2, 0.1.
%
  iteration_max = 200;
  h = 0.40;
  p01_demo ( iteration_max, h );
  input ( 'Press RETURN' )

  iteration_max = 200;
  h = 0.20;
  p01_demo ( iteration_max, h );
  input ( 'Press RETURN' )

  iteration_max = 200;
  h = 0.10;
  p01_demo ( iteration_max, h );
  input ( 'Press RETURN' )
%
%  Problem 2, unit circle with a hole.
%
  iteration_max = 200;
  h = 0.10;
  p02_demo ( iteration_max, h )

  input ( 'Press RETURN' )
%
%  Problem 3, square with a hole, uniform density.
%
  iteration_max = 200;
  h = 0.15;
  fh = @p03_fh;
  p03_demo ( iteration_max, h, fh )

  input ( 'Press RETURN' )
%
%  Problem 3, square with a hole, finer density near the hole.
%
  iteration_max = 300;
  h = 0.05;
  fh = inline ( 'min(4*sqrt(sum(p.^2,2))-1,2)', 'p' );
  p03_demo ( iteration_max, h, fh )

  input ( 'Press RETURN' )
%
%  Problem 4, hexagon with a hexagonal hole.
%
  iteration_max = 200;
  h = 0.1;
  p04_demo ( iteration_max, h );

  input ( 'Press RETURN' )
%
%  Problem 5, the horn.
%
  iteration_max = 200;
  h = 0.020;
  p05_demo ( iteration_max, h );

  input ( 'Press RETURN' )
%
%  Problem 6, the superellipse.
%  Needs MAPLE to run...
%
  if ( 0 )
    iteration_max = 200;
    h = 0.08;
    p06_demo ( iteration_max, h );
    input ( 'Press RETURN' )
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Skipping problem 6, the superellipse.\n' );
  end
%
%  Problem 7, the bicycle seat.
%  Needs MAPLE to run...
%
  if ( 0 )
    iteration_max = 200;
    h = 0.75;
    p07_demo ( iteration_max, h );
    input ( 'Press RETURN' )
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Skipping problem 7, the bicycle seat.\n' );
  end
%
%  Problem 8, the holey pie slice, uniform density.
%
  iteration_max = 200;
  h = 0.025;
  fh = @huniform;
  p08_demo ( iteration_max, h, fh );
  input ( 'Press RETURN' )
%
%  Problem 8, the holey pie slice, variable density.
%
  iteration_max = 200;
  h = 0.005;
  fh = @p08_fh;
  p08_demo ( iteration_max, h, fh );
  input ( 'Press RETURN' )
%
%  Problem 9, Jeff Borggaard's square with two hexagonal holes.
%
  iteration_max = 200;
  h = 0.08;
  p09_demo ( iteration_max, h );
  input ( 'Press RETURN' )
%
%  Problem 10, a simple square.
%
  iteration_max = 200;
  h = 0.10;
  p10_demo ( iteration_max, h )
  input ( 'Press RETURN' )
%
%  Problem 11, the L-shaped region.
%
  iteration_max = 200;
  h = 0.10;
  p11_demo ( iteration_max, h )
  input ( 'Press RETURN' )
%
%  Problem 12, John Shadid's H-shaped region.
%
  iteration_max = 200;
  h = 0.05;
  p12_demo ( iteration_max, h )
  input ( 'Press RETURN' )
%
%  Problem 13, the Sandia Fork.
%
  iteration_max = 200;
  h = 0.025;
  p13_demo ( iteration_max, h )
  input ( 'Press RETURN' )
%
%  Problem 14, Marcus Garvie's Lake Alpha with Beta island.
%
  if ( 0 )
    iteration_max = 50;
    h = 20.0;
    fh = @huniform;
    p14_demo ( iteration_max, h, fh )
    input ( 'Press RETURN' )
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Skipping problem 14.\n' );
  end
%
%  Problem 14, Marcus Garvie's Lake Alpha with Beta island.
%  You may get endless warnings from DELAUNAYN about duplicate points.
%  Shut off this annoying drivel with "warning off".
%
  if ( 0 )
    warning off
    iteration_max = 50;
    h = 10.0;
    fh = @p14_fh;
    p14_demo ( iteration_max, h, fh )
    input ( 'Press RETURN' )
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Skipping problem 14.\n' );
  end
%
%  Problem 15, Sangbum Kim's forward step problem.
%
  iteration_max = 50;
  h = 0.2;
  p15_demo ( iteration_max, h )
  input ( 'Press RETURN' )
%
%  Problem 16, Kevin Pond's elbow.
%
  iteration_max = 50;
  h = 0.05;
  p16_demo ( iteration_max, h )
  input ( 'Press RETURN' )
%
%  Problem 17, Reuleaux triangle problem.
%
  iteration_max = 50;
  h = 0.05;
  p17_demo ( iteration_max, h )
  input ( 'Press RETURN' )
%
%  Problem 18, Dumbbell problem.
%
  iteration_max = 50;
  h = 0.100;
  p18_demo ( iteration_max, h )
  input ( 'Press RETURN' )
%
%  Problem 19, Dumbbell problem.
%
  iteration_max = 200;
  h = 0.025;
  p19_demo ( iteration_max, h )
  input ( 'Press RETURN' )
%
%  Problem 20, ICAM Wright House problem.
%
  iteration_max = 2;
  h = 2.0;
  p20_demo ( iteration_max, h )
  input ( 'Press RETURN' )
%
%  Problem 21, Zhu Wang's quarter round problem.
%
  iteration_max = 200;
  h = 0.100;
  p21_demo ( iteration_max, h )
  input ( 'Press RETURN' )
%
%  Problem 22, Hans-Werner van Wyk's Big C.
%
  iteration_max = 200;
  h = 0.100;
  p22_demo ( iteration_max, h )
  input ( 'Press RETURN' )
%
%  Problem 23, Mike Schneier's nonuniform square.
%
  iteration_max = 200;
  h = 0.025;
  p23_demo ( iteration_max, h )
  input ( 'Press RETURN' )
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MESHDEMO_2D\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function res = cache(this, stuff)
% Method for retrieving/putting stuff from/into the temporary cache
% FORMAT res = cache(obj, name)
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Stefan Kiebel
% $Id: cache.m 1373 2008-04-11 14:24:03Z spm $

if ischar(stuff)
    % assume this is a retrieval
        res = getcache(this, stuff);
else
    % assume this is a put
    name = inputname(2);
        res = setcache(this, name, stuff);
end

function res = getcache(this, stuff)
if isfield(this.cache, stuff)
    eval(['res = this.cache.' stuff ';']);
else
    res = [];
end

function this = setcache(this, name, stuff)
eval(['this.cache(1).' name ' = stuff;']);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip