text:
function model = findSExRxnInd(model, nRealMet, printLevel)
% Returns a model with boolean vectors indicating internal vs external (exchange/demand/sink) reactions.
% Finds the reactions in the model which export/import from the model boundary
%
% e.g. Exchange reactions,
% Demand reactions,
% Sink reactions
%
% USAGE:
%
%    model = findSExRxnInd(model, nRealMet, printLevel)
% model.SIntRxnBool         Boolean of reactions heuristically though to be mass balanced.
%
% INPUT:
%    model:         structure with:
%
%                     * S - m x n stoichiometric matrix
%
% OPTIONAL INPUT:
%    model:         structure with:
%                     * model.biomassRxnAbbr - abbreviation of biomass reaction
%    nRealMet:      specified in case extra rows in `S` which dont correspond to metabolties
%    printLevel:    verbose level
%
% OUTPUT:
%    model:         structure with:
%
%                     * .SIntRxnBool - Boolean of reactions heuristically though to be mass balanced.
%                     * .SIntMetBool - Boolean of metabolites heuristically though to be involved in mass balanced reactions.
%                     * .SOnlyIntMetBool - Boolean of metabolites heuristically though only to be involved in mass balanced reactions.
%                     * .SExMetBool - Boolean of metabolites heuristically though to be involved in mass imbalanced reactions.
%                     * .SOnlyExMetBool - Boolean of metabolites heuristically though only to be involved in mass imbalanced reactions.
%                     * .biomassBool - Boolean of biomass reaction
%                     * .DMRxnBool - Boolean of demand reactions. Prefix `DM_` (optional field)
%                     * .SinkRxnBool - Boolean of sink reactions. Prefix `sink_` (optional field)
%                     * .ExchRxnBool - Boolean of exchange reactions. Prefix `EX_` or `Exch_` or `Ex_` (optional field)
%
% .. Author: -  Ronan Fleming

if ~exist('printLevel','var')
    printLevel=0;
end

[nMet,nRxn]=size(model.S);

if ~exist('nRealMet','var')
    nRealMet=length(model.mets);
    if nMet~=nRealMet
        if printLevel>0
            fprintf('%s\n','Detected extra rows of S without corresponding metabolite abbreviations.')
        end
    end
else
    if isempty(nRealMet)
        nRealMet=length(model.mets);
    end
end

%locate biomass reaction if there is one
biomassBool=false(nRxn,1);
if ~isfield(model,'c')
    model.c=zeros(nMet,1);
end
if ~isfield(model,'biomassRxnAbbr')
    bool=model.c~=0;
    if nnz(bool)==1
        model.biomassRxnAbbr=model.rxns{model.c~=0};
        if printLevel>0
            fprintf('%s%s\n','Assuming biomass reaction is: ', model.biomassRxnAbbr);
        end
        biomassBool(bool)=1;
    else
        if nnz(bool)==0
            if printLevel>0
                fprintf('%s\n','No model.biomassRxnAbbr? Give abbreviation of biomass reaction if there is one.');
            end
        else
            warning('More than one biomass reaction?');
        end
    end
else
    biomassBool=contains(model.rxns,model.biomassRxnAbbr);%finds a subset of the abbreviation
    if any(biomassBool)
        if nnz(biomassBool)>1
            ind = find(biomassBool);
            for p=1:length(ind)
                if printLevel>0
                    fprintf('%s%s\n','Found multiple possible biomass reactions: ', model.rxns{ind(p)});
                end
            end
        else
            fprintf('%s%s\n','Found biomass reaction: ', model.rxns{biomassBool});
        end
    else
        if printLevel>0
            fprintf('%s\n','Assuming no biomass reaction.');
        end
    end
end

%grab any possible biomass reactions
biomassBool = biomassBool | contains(model.rxns,'biomass');
biomassBool = biomassBool | contains(model.rxns,'Biomass');

model.biomassBool=biomassBool;

SExRxnBoolOneCoefficient=false(nRxn,1);
for n=1:nRxn
    %find reactions with only one coefficient
    %or no coefficient at all - Ronan May 29th 2011
    if nnz(model.S(1:nRealMet,n))<=1
        SExRxnBoolOneCoefficient(n,1)=1;
        if printLevel>2
            if nonzeros(model.S(1:nRealMet,n))>0
                fprintf('%s\t%s\n','Positive coefficient:',model.rxns{n});
            else
                fprintf('%s\t%s\n','Negative coefficient:',model.rxns{n});
                %                 fprintf('%s%s%s%s%s\n','''',model.rxns{n},''',0 ,0 ''',model.mets{find(model.S(1:nRealMet,n)~=0)},''',0 ,0 ;');
            end
        end
    end
end

% models with typical HMR subsystems - heuristic
if isfield(model,'subSystems')
    model.ExchRxnBool=cellfun(@(x) any(ismember({'Exchange reactions','Artificial reactions','Pool reactions'},x)),model.subSystems);
    if isfield(model,'rxnComps')
        model.ExchRxnBool=model.ExchRxnBool | strcmp('x',model.rxnComps);
    end
    % models with typical COBRA abbreviations - heuristic
    model.ExchRxnBool=strncmp('EX_', model.rxns, 3)==1 | strncmp('Exch_', model.rxns, 5)==1 | strncmp('Ex_', model.rxns, 5)==1 | biomassBool | model.ExchRxnBool;
else
    % models with typical COBRA abbreviations - heuristic
    model.ExchRxnBool=strncmp('EX_', model.rxns, 3)==1 | strncmp('Exch_', model.rxns, 5)==1 | strncmp('Ex_', model.rxns, 5)==1 | biomassBool;
end
%demand reactions going out of model
model.DMRxnBool=strncmp('DM_', model.rxns, 3)==1;
%sink reactions going into or out of model
model.SinkRxnBool=strncmp('sink_', model.rxns, 5) | strncmp('Sink_', model.rxns, 5)==1;

%input/output
SExRxnBoolHeuristic = model.ExchRxnBool | model.DMRxnBool | model.SinkRxnBool;

%remove ATP demand as it is usually mass balanced
bool=strcmp('ATPM',model.rxns);
if any(bool)
    if printLevel>0
        fprintf('%s\n','ATP maintenance reaction is not considered an exchange reaction by default. It should be mass balanced:')
        formulas = printRxnFormula(model,{'ATPM'});
    end
    model.DMRxnBool(bool)=0;
end
bool=strcmp('DM_atp(c)',model.rxns);
if any(bool)
    if printLevel>0
        fprintf('%s\n','ATP demand reaction is not considered an exchange reaction by default. It should be mass balanced')
        formulas = printRxnFormula(model,{'DM_atp(c)'});
    end
    model.DMRxnBool(bool)=0;
end
bool=strcmp('DM_atp_c_',model.rxns);
if any(bool)
    if printLevel>0
        fprintf('%s\n','ATP demand reaction is not considered an exchange reaction by default. It should be mass balanced:')
        formulas = printRxnFormula(model,{'DM_atp_c_'});
    end
    model.DMRxnBool(bool)=0;
end
%remove atp demand
SExRxnBoolHeuristic(bool)=0;

diffBool= ~SExRxnBoolHeuristic & SExRxnBoolOneCoefficient;
if any(diffBool)
    if printLevel>0
        fprintf('%s\n','Exchanges that would otherwise have been missed without abbreviation prefix search:')
        fprintf('%s\t%s\t%s\t\t%s\t\t%s\n','Coefficient','Metabolite','#','Reaction','#')
    end
    for n=1:nRxn
        if diffBool(n)
            objMetInd=find(model.S(:,n));
            for m=1:length(objMetInd)
                Sij=full(model.S(objMetInd(m),n));
                if length(model.mets{objMetInd(m)})<4
                    if printLevel>0
                        fprintf('%g\t\t\t%s\t\t\t%i\t%s\t\t%i\n',Sij,model.mets{objMetInd(m)},objMetInd(m),model.rxns{n},n)
                    end
                else
                    if length(model.mets{objMetInd(m)})<8
                        if printLevel>0
                            fprintf('%g\t\t\t%s\t\t%i\t%s\t\t%i\n',Sij,model.mets{objMetInd(m)},objMetInd(m),model.rxns{n},n)
                        end
                    else
                        if length(model.mets{objMetInd(m)})<12
                            if printLevel>0
                                fprintf('%g\t\t\t%s\t%i\t%s\t\t%i\n',Sij,model.mets{objMetInd(m)},objMetInd(m),model.rxns{n},n)
                            end
                        end
                    end
                end
            end
        end
    end
end

% %dont check if there are coupling constraints
% %(E. coli E matrix specific)
% if ~isfield(model,'A')
%     diffBool= SExRxnBoolHeuristic & ~SExRxnBoolOneCoefficient;
%     if any(diffBool)
%         if printLevel>0
%             fprintf('%s\n','Exchanges missed by prefix search:')
%             fprintf('%s\t%s\n','#', 'Exchange')
%         end
%         for n=1:length(diffBool)
%             if diffBool(n)
%                 equation=printRxnFormula(model,model.rxns(n),0);
%                 if printLevel>0
%                     fprintf('%i\t%s\t%s\n',n,model.rxns{n},equation{1});
%                 end
%             end
%         end
%     end
% end

%amalagamate all exchanges
SExRxnBool= SExRxnBoolHeuristic | SExRxnBoolOneCoefficient;
model.SIntRxnBool=~SExRxnBool;
%rows corresponding to internal reactions
boolMet=true(nMet,1);
%first pair
model.SIntMetBool = getCorrespondingRows(model.S,boolMet,model.SIntRxnBool,'inclusive');
model.SOnlyExMetBool = getCorrespondingRows(model.S,boolMet,~model.SIntRxnBool,'exclusive');
%second pair
model.SOnlyIntMetBool = getCorrespondingRows(model.S,boolMet,model.SIntRxnBool,'exclusive');
model.SExMetBool = getCorrespondingRows(model.S,boolMet,~model.SIntRxnBool,'inclusive');
%sanity check
if nnz(model.SIntMetBool)+nnz(model.SOnlyExMetBool) ~= nnz(model.SIntMetBool)+nnz(model.SOnlyExMetBool)
    error('Inconsistency in metabolite counts')
end

author: opencobra

repo: cobratoolbox


text:
function y=expSlope1(x)
%%EXPSLOPE1  Evaluate the function (exp(x)-1)/x avoiding problems at the
%            point x=0, where the limit equals 1. This function arises in
%            range enclosure algorithms for interval arithmetic.
%
%INPUTS: x A vector of matrix of real values at which the function should
%          be evaluated.
%
%OUTPUTS: y The value(s) of the function (exp(x)-1)/x at the point(s) in x.
%
%For values of x below 0.1, a Taylor series expansion about x=0 is used.
%For values of x above 0.1, the function is implemented as in Chapter
%1.14.1 of [2]. Though the method in [2] can go to lower values of x, it
%still has a singularity at x=0. The Taylor series does not.
%
%This is one of the slope function in Table 10.6 of Section 10.6.2 of [1].
%
%REFERENCES:
%[1] IEEE Standard for Interval Arithmetic," in IEEE Std 1788-2015,
%    pp.1-97, June 30 2015.
%[2] N. J. Higham, Accuracy and Stability of Numerical Algorithms.
%    Philadelphia: SIAM, 1996.
%
%November 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

y=zeros(size(x));

selSmall=x<0.1;
xSmall=x(selSmall);
xLarge=x(~selSmall);

if(~isempty(xSmall))
    %Use a Taylor-series approximation about x=0. THis avoud the
    %singularity.
    expanOrder=20;
    ySmall=ones(size(xSmall));

    xPow=xSmall;
    kFactorial=1;
    for k=1:expanOrder
       kFactorial=kFactorial*k;

       ySmall=ySmall+xPow/((k+1)*kFactorial);

       xPow=xPow.*xSmall;
    end
else
    ySmall=[];
end

if(~isempty(xLarge))
    yLarge=exp(xLarge);
    yLarge=(yLarge-1)/log(yLarge);
else
    yLarge=[];
end

y(selSmall)=ySmall;
y(~selSmall)=yLarge;
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function y = isfinite(x)
    %ISFINITE True for finite elements of symbolic arrays.
    %   ISFINITE(X) returns an array that contains 1's where
    %   the elements of X are finite and 0's where they are not.
    %   For example, isfinite(sym([pi NaN Inf -Inf])) is [1 0 0 0].
    %
    %   For any X, exactly one of ISFINITE(X), ISINF(X), or ISNAN(X)
    %   is 1 for each element.
    %
    %   See also SYM/ISINF, SYM/ISNAN.

    y = ~isinf(x) & ~isnan(x);
end

author: ayonga

repo: frost-dev

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


text:
% CORRIMAGE - compute correlation image between an event and amplitude
%               and phase in the time-frequency domain.
%
% Usage:
%     corrimage(data, sortvar, times);
%     [times, freqs, alpha, sigout, limits ] = corrimage(data, sortvar, ...
%                                    times, 'key1', val1, 'key2', val2, ...)
%
% Inputs:
%   data    - [float array] data of size (times,trials)
%   sortvar - [float array] sorting variable of size (trials)
%   times   - [float array] time indices for every data point. Size of
%             (times). Note: same input as the times vector for erpimage.
%
% Optional input:
%   'mode'     - ['phase'|'amp'] compute correlation of event values with
%                phase ('phase') or amplitude ('amp') of signal at the
%                given time-frequency points. Default is 'amp'.
%   'freqs'    - [min nfreqs max] build a frequency vector of size nfreqs
%                with frequency spaced in a log scale. Then compute
%                correlation at these frequency value.
%                Default is 50 points in a log scale between 2.5Hz to 50Hz.
%   'times'    - [float vector] compute correlation at these time value (ms).
%                (uses closest times points in 'times' input and return
%                them in the times output).
%                Default is 100 steps between min time+5% of time range and
%                max time-5%. Enter only 2 values [N X] to generate N
%                time points and trim by X %. If N is negative, uses it as
%                a subsampling factor [-3 5] trims times by 5% and subsample
%                by 3 (by subsampling one obtains a regularly spaced times).
%   'trim'     - [low high] low and high percentile of sorted sortvar values
%                to retain. i.e. [5 95] remove the 5 lowest and highest
%                percentile of sortvar values (and associated data) before
%                computing statistics. Default is [0 100].
%   'align'    - [float] same as 'align' parameter of ERPIMAGE. This
%                parameter is used to constrain the 'times' parameter so
%                correlation with data trials containing 0-values (as a
%                result of data alignment) are avoided: computing these
%                correlations would produce spurious significant results.
%                Default is no alignment.
%   'method'   - ['erpimage'|timefreq'] use either the ERPIMAGE function
%                of the TIMEFREQ function to compute spectral decomposition.
%                Default is 'timefreq' for speed reasons (note that both
%                methods should return the same results).
%   'erpout'   - [min max] regress out ERP using the selected time-window [min
%                max] in ms (the ERP is subtracted from the whole time period
%                but only regressed out in the selected time window).
%   'triallimit' - [integer array] specify trial boundaries for subjects. For
%                instance [1 200 400] indicates 2 subjects, trials 1 to 199
%                for subject 1 and trials 200 to 399 for subject 2. This is
%                currently only used for regressing erp out.
%
% Processing options:
%   'erpimopt' - [cell array] erpimage additional options (number of cycle ...).
%   'tfopt'    - [cell array] timefreq additional options (number of cycle ...).
%
% Plotting options:
%   'plotvals' - [cell array of output values] enter output values here
%                { times freqs alpha sigout} to replot them.
%   'nofig'    - ['on'|'off'] do not create figure.
%   'cbar'     - ['on'|'off'] plot color bar. Default: 'on'.
%   'smooth'   - ['on'|'off'] smooth significance array. Default: 'on'.
%   'plot'     - ['no'|'alpha'|'sigout'|'sigoutm'|'sigoutm2'] plot -10*log(alpha)
%                values ('alpha'); output signal (slope or ITC) ('sigout'),
%                output signal masked by significance ('sigoutm') or the last
%                2 option combined ('sigoutm2'). 'no' prevent the function
%                from plotting. In addition, see pmask. Default is 'sigoutmasked'.
%   'pmask'    - [real] maximum p value to show in plot. Default is 0.00001
%                (0.001 taking into account multiple comparisons (100)). Enter
%                0.9XXX to get the higher tail or a negative value (e.e., -0.001
%                to get both tails).
%   'vert'     - [real array] time vector for vertivcal lines.
%   'limits'   - [min max] plotting limits.
%
% Outputs:
%   times    - [float vector] vector of times (ms)
%   freqs    - [float vector] vector of frequencies (Hz)
%   alpha    - [float array] array (freqs,times) of p-values.
%   sigout   - [float array] array (freqs,times) of signal out (coherence
%              for phase and slope for amplitude).
%
% Important note: the 'timefreq' method may truncate the time window to
%                 compute the spectral decomposition at the lowest freq.
%
% Author: Arnaud Delorme & Scott Makeig, SCCN UCSD,
%         and CNL Salk Institute, 18 April 2003
%
% See also: ERPIMAGE

%123456789012345678901234567890123456789012345678901234567890123456789012

% Copyright (C) 2002 Arnaud Delorme
%
% 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.

function [time, freq, alpha, sigout, limits, tf, sortvar] = corrimage(data, sortvar, timevect, varargin)

if nargin < 3
    help corrimage;
    return;
end

if isempty(timevect), timevect = [1 2]; end

% check inputs
% ------------
g = finputcheck(varargin, { 'freqs'    'real'   [0 Inf]    [2.5 50 50];
                            'times'    'real'   []         [100 5]; % see function at the end
                            'mode'     'string' { 'phase','amp' }  'amp';
                            'vert'     'real'   []         [];
                            'align'    { 'real','cell' }   []         [];
                            'plotvals' 'cell'   []         {};
                            'pmask'    'real'   []         0.00001;
                            'triallimit' 'integer'   []       [];
                            'trim'     'real'   [0 100]    [0 100];
                            'limits'   'real'   []         [];
                            'method'   'string' { 'erpimage','timefreq' }        'timefreq';
                            'plot'     'string' { 'no','alpha','sigout','sigoutm','sigoutp','sigoutm2' }  'sigoutm';
                            'nofig'    'string' { 'on','off' } 'off';
                            'cbar'     'string' { 'on','off' } 'on';
                            'smooth'   'string' { 'on','off' } 'off';
                            'erpout'   'real'   []             [];
                            'tfopt'    'cell'   []             {};
                            'erpimopt' 'cell'   []             {} });
if ischar(g), error(g); end

fprintf('Generating %d frequencies in log scale (ignore message on linear scale)\n', g.freqs(2));
g.freqs = logscale(g.freqs(1), g.freqs(3), g.freqs(2));

frames = length(timevect);
if size(data,1) == 1
    data = reshape(data, frames, size(data,2)*size(data,3)/frames);
end

% trim sortvar values
% -------------------
[sortvar sortorder] = sort(sortvar);
data = data(:, sortorder);
len = length(sortvar);
lowindex  = round(len*g.trim(1)/100)+1;
highindex = round(len*g.trim(2)/100);
sortvar = sortvar(lowindex:highindex);
data    = data(:, lowindex:highindex);
if lowindex ~=1 || highindex ~= length(sortorder)
    fprintf('Actual percentiles %1.2f-%1.2f (indices 1-%d -> %d-%d): event vals min %3.2f; max %3.2f\n', ...
             100*(lowindex-1)/len, 100*highindex/len, len, lowindex, highindex, min(sortvar), max(sortvar));
end

% assign subject number for each trial
% ------------------------------------
if ~isempty(g.triallimit)
    alltrials = zeros(1,len);
    for index = 1:length(g.triallimit)-1
        alltrials([g.triallimit(index):g.triallimit(index+1)-1]) = index;
    end
    alltrials = alltrials(sortorder);
    alltrials = alltrials(lowindex:highindex);
end

%figure; hist(sortvar)

% constraining time values depending on data alignment
% ----------------------------------------------------
srate = 1000*(length(timevect)-1)/(timevect(end)-timevect(1));
if ~isempty(g.align)

    if iscell(g.align)
        fprintf('Realigned sortvar plotted at %g ms.\n',g.align{1});
        shifts = round((g.align{1}-g.align{2})*srate/1000); % shifts can be positive or negative
        g.align = g.align{1};
    else
        if isinf(g.align)
            g.align = median(sortvar);
        end
        fprintf('Realigned sortvar plotted at %g ms.\n',g.align);

        % compute shifts
        % --------------
        shifts = round((g.align-sortvar)*srate/1000); % shifts can be positive or negative
    end

    %figure; hist(shifts)
    minshift = min(shifts); % negative
    maxshift = max(shifts); % positive
    if minshift > 0, error('minshift has to be negative'); end
    if maxshift < 0, error('maxshift has to be positive'); end

    % realign data for all trials
    % ---------------------------
    aligndata=zeros(size(data))+NaN; % begin with matrix of zeros()
    for t=1:size(data,2), %%%%%%%%% foreach trial %%%%%%%%%
        shft = shifts(t);
        if shft>0, % shift right
            aligndata(shft+1:frames,t)=data(1:frames-shft,t);
        elseif shft < 0 % shift left
            aligndata(1:frames+shft,t)=data(1-shft:frames,t);
        else % shft == 0
            aligndata(:,t) = data(:,t);
        end
    end % end trial
    data = aligndata(maxshift+1:frames+minshift,:);
    if any(any(isnan(data))), error('NaNs remaining after data alignment'); end
    timevect = timevect(maxshift+1:frames+minshift);

    % take the time vector subset
    % ---------------------------
    if isempty(timevect), error('Shift too big, empty time vector');
    else fprintf('Time vector truncated for data alignment between %1.0f and %1.0f\n', ...
                 min(timevect), max(timevect));
    end;
end

% regress out the ERP from the data (4 times to have residual ERP close to 0)
% ---------------------------------------------------------------------------
if ~isempty(g.erpout)
    %data = erpregout(data);
    %data = erpregout(data);
    %data = erpregout(data);
    disp('Regressing out ERP');
    erpbeg = mean(data,2);
    if ~isempty(g.triallimit)
        for index = 1:length(g.triallimit)-1
            trials = find(alltrials == index);
            %[data(:,trials) erp factors]= erpregout(data(:,trials), [timevect(1) timevect(end)], [300 400]);
            erpstart = mean(data(:,trials),2);
            data(:,trials) = erpregout(data(:,trials));
            data(:,trials) = erpregout(data(:,trials));
            data(:,trials) = erpregout(data(:,trials));
            data(:,trials) = erpregout(data(:,trials));
            erpend = mean(data(:,trials),2);
            fprintf([ '***********************************************\n' ...
                      'Ratio in ERP after regression (compare to before) is %3.2f\n' ...
                      '***********************************************\n' ], mean(erpend./erpstart));
        end
    end
    erpend = mean(data,2);
    fprintf([ '***********************************************\n' ...
              'Ratio in grand ERP after regression (compare to before) is %3.2f\n' ...
              '***********************************************\n' ], mean(erpend./erpbeg));

    if 0
    %trials = find(alltrials == 1);
    data2 = data;
    dasf

    [data(:,trials) erp factors]= erpregout(data(:,trials), [timevect(1) timevect(end)], [300 400]);

    figure; subplot(1,2,1); erpimage(data, sortvar, timevect, '', 300, 500, 'erp');

    figure;
    subplot(1,2,1); erpimage(data2, sortvar, timevect, '', 300, 500, 'erp');
    subplot(1,2,2); erpimage(data , sortvar, timevect, '', 300, 500, 'erp');

    figure;
    subplot(1,2,1); erpimage(data2(:,trials), sortvar(trials), timevect, '', 0, 0, 'erp');
    subplot(1,2,2); erpimage(data (:,trials), sortvar(trials), timevect, '', 0, 0, 'erp');

    figure;
    for index = 1:length(g.triallimit)-1
        subplot(3,5,index);
        trials = find(alltrials == index);
        erpimage(data(:,trials), sortvar(trials), timevect, '', 50, 100, 'erp');
    end

    disp('Regressing out ERP');
    if ~isempty(g.triallimit)
        for index = 1:length(g.triallimit)-1
            trials = find(alltrials == index);
            [data(:,trials) erp factors]= erpregout(data(:,trials), [timevect(1) timevect(end)], [300 400]);
            data(:,trials) = erpregout(data(:,trials));
            data(:,trials) = erpregout(data(:,trials));
            data(:,trials) = erpregout(data(:,trials));
            data(:,trials) = erpregout(data(:,trials));
            data(:,trials) = erpregout(data(:,trials));
            data(:,trials) = erpregout(data(:,trials));
            data(:,trials) = erpregout(data(:,trials));
            %data(:,trials) = erpregout(data(:,trials), [timevect(1) timevect(end)], g.erpout);
            %data(:,trials) = erpregout(data(:,trials), [timevect(1) timevect(end)], g.erpout);
            %data(:,trials) = erpregout(data(:,trials), [timevect(1) timevect(end)], g.erpout);
            %data(:,trials) = erpregout(data(:,trials), [timevect(1) timevect(end)], g.erpout);
            %data(:,trials) = erpregout(data(:,trials), [timevect(1) timevect(end)], g.erpout);
            %data(:,trials) = erpregout(data(:,trials), [timevect(1) timevect(end)], g.erpout);
            %data(:,trials) = erpregout(data(:,trials), [timevect(1) timevect(end)], g.erpout);
            %data(:,trials) = erpregout(data(:,trials), [timevect(1) timevect(end)], g.erpout);
        end
    else
        %data = erpregout(data, [timevect(1) timevect(end)], g.erpout);
    end
    end
end

% generate time points
% --------------------
g.times = gettimes(timevect, g.times);
data = reshape(data, 1, size(data,2)*size(data,1));

% time frequency decomposition
% ----------------------------
if strcmpi(g.method, 'timefreq') && isempty(g.plotvals)
    data = reshape(data, length(data)/length(sortvar), length(sortvar));
    [tf, g.freqs, g.times] = timefreq(data, srate, 'freqs', g.freqs, 'timesout', g.times, ...
                                      'tlimits', [min(timevect) max(timevect)], 'wavelet', 3);
    outvar   = sortvar;
end

% compute correlation
% -------------------
if strcmpi(g.mode, 'phase') && isempty(g.plotvals)
    for freq = 1:length(g.freqs)
        fprintf('Computing correlation with phase %3.2f Hz ----------------------\n', g.freqs(freq));
        for time = 1:length(g.times)
            if strcmpi(g.method, 'erpimage')
                [outdata,outvar,outtrials,limits,axhndls,erp, ...
                 amps,cohers,cohsig,ampsig,outamps,phsangls,phsamp,sortidx] =erpimage(data,sortvar,timevect, ...
                        '', 0,0,g.erpimopt{:}, 'phasesort', [g.times(time) 0 g.freqs(freq)], 'noshow', 'yes');
            else
                phsangls = angle( squeeze(tf(freq, time, :))' );
            end

            % computing ITCs
            [ ITC(freq, time) alpha(freq, time) ] = ...
                    bootcircle(outvar/mean(outvar), exp(j*phsangls), 'naccu', 250);

            % accumulate 200 values, fitted with a normal distribution
            % --------------------------------------------------------
            %cmplx = outvar .* exp(j*phsangls)/mean(outvar);
            %ITC(freq, time) = mean(cmplx);
            %alpha(freq,time) = bootstat(outvar/mean(outvar), exp(j*phsangls), 'res = mean(arg1 .* arg2)', ...
            %                         'naccu', 100, 'vals', abs(ITC(freq, time)), 'distfit', 'norm');
        end
    end
    sigout = ITC;
elseif isempty(g.plotvals)
    for freq = 1:length(g.freqs)
        fprintf('Computing correlation with amplitude %3.2f Hz ----------------------\n', g.freqs(freq));
        for time = 1:length(g.times)
            if strcmpi(g.method, 'erpimage')
                [outdata,outvar,outtrials,limits,axhndls,erp, ...
                 amps,cohers,cohsig,ampsig,outamps,phsangls,phsamp,sortidx] =erpimage(data,sortvar,timevect, ...
                             '', 0,0,g.erpimopt{:}, 'ampsort', [g.times(time) 0 g.freqs(freq)], 'noshow', 'yes');
            else
                phsamp = abs(squeeze(tf(freq, time, :)))';
            end

            % computing ITCs
            [ypred alpha(freq, time) Rsq slope(freq, time)] = myregress(outvar, 20*log10(phsamp));
        end
    end;
    sigout = slope;
else
    g.times = g.plotvals{1};
    g.freqs = g.plotvals{2};
    alpha   = g.plotvals{3};
    sigout  = g.plotvals{4};
end

% plot correlation
% ----------------
if ~strcmp('plot', 'no')
    if ~isreal(sigout),
        if strcmpi(g.plot, 'sigoutp')
            sigoutplot = angle(sigout);
        else
            sigoutplot = abs(sigout);
        end
    else sigoutplot = sigout;
    end
    sigouttmp = sigoutplot;
    if strcmpi(g.smooth, 'on')
        tmpfilt = gauss2d(3,3,.3,.3);
        tmpfilt = tmpfilt/sum(sum(tmpfilt));
        alpha = conv2(alpha, tmpfilt, 'same');
    end

    % mask signal out
    if g.pmask > 0.5
        indices = find( alpha(:) < g.pmask);
        sigouttmp = sigoutplot;
        sigouttmp(indices) = 0;
    elseif g.pmask < 0 % both sides, i.e. 0.01
        sigouttmp = sigoutplot;
        indices = intersect_bc(find( alpha(:) > -g.pmask), find( alpha(:) < 1+g.pmask));
        sigouttmp(indices) = 0;
    else
        sigouttmp = sigoutplot;
        indices = find( alpha(:) > g.pmask);
        sigouttmp(indices) = 0;
    end

    if strcmpi(g.nofig, 'off'), figure; end
    switch g.plot
     case 'alpha',   limits = plotfig(g.times, g.freqs, -log10(alpha), g);
     case 'sigout',  limits = plotfigsim(g.times, g.freqs, sigoutplot, g);
     case { 'sigoutm' 'sigoutp' }, limits = plotfigsim(g.times, g.freqs, sigouttmp, g);
     case 'sigoutm2', limits = subplot(1,2,1); plotfigsim(g.times, g.freqs, sigoutplot, g);
                      limits = subplot(1,2,2); plotfigsim(g.times, g.freqs, sigouttmp, g);
    end
end

time = g.times;
freq = g.freqs;

return;

% formula for the normal distribution
% -----------------------------------
function y = norm(mu, sigma, x);

y = 1/sqrt(2) * exp( -(x-mu).^2/(sigma*sigma*2) ) / (sqrt(pi)*sigma);

% get time points
% ---------------
function val = gettimes(timevect, timevar);
    if length(timevar) == 2

        if timevar(1) > 0
            % generate linearly space vector
            % ------------------------------
            npoints = timevar(1);
            trim    = timevar(2);
            if length(timevect)-2*round(trim/100*length(timevect)) < npoints
                npoints = length(timevect)-round(trim/100*length(timevect));
            end
            fprintf('Generating %d times points trimmed by %1.1f percent\n', npoints, trim);
            timer = max(timevect) - min(timevect);
            maxt = max(timevect)-timer*trim/100;
            mint = min(timevect)+timer*trim/100;
            val = linspace(mint,maxt, npoints);
        else
            % subsample data
            % --------------
            nsub    = -timevar(1);
            trim    = timevar(2);
            len     = length(timevect);
            trimtimevect = timevect(round(trim/100*len)+1:len-round(trim/100*len));
            fprintf('Subsampling by %d and trimming data by %1.1f percent (%d points)\n', nsub, trim, round(trim/100*len));
            val = trimtimevect(1:nsub:end);
        end
    else
        val = timevar;
    end;

    % find closet points in data
    oldval = val;
    for index = 1:length(val)
        [dum ind] = min(abs(timevect-val(index)));
        val(index) = timevect(ind);
    end
    if length(val) < length(unique(val))
        disp('Warning: duplicate times, reduce the number of output times');
    end
    if all(oldval == val)
        disp('Debug msg: Time value unchanged by finding closest in data');
    end;

% get log scale (for frequency)
% -------------
function val = logscale(a,b,n);
    %val = [5 7 9]; return;
    val = linspace(log(a), log(b), n);
    val = exp(val);

% plot figure
% -----------
function limits = plotfig(times, freqs, vals, g)
    icadefs;
    imagesc(times, [1:size(vals,1)], vals);

    colormap(DEFAULT_COLORMAP);
    ticks = linspace(1, size(vals,1), length(freqs));
    ticks = ticks(1:4:end);
    set(gca, 'ytick', ticks);
    set(gca, 'yticklabel', num2str(freqs(1:4:end)', 3))

    xlabel('Time (ms)'); ylabel('Freq (Hz)');
    if ~isempty(g.limits), caxis(g.limits); end
    limits = caxis;
    if ~isempty(g.vert)
        for vert = g.vert(:)'
            hold on; plot([vert vert], [0.001 500], 'k', 'linewidth', 2);
        end
    end
    if strcmpi(g.cbar,'on')
        cbar;
    end

% plot figure with symmetrical phase
% ---------------------------------
function limits = plotfigsim(times, freqs, vals, g)
    icadefs;
    imagesc(times, [1:size(vals,1)], vals);

    colormap(DEFAULT_COLORMAP);
    ticks = linspace(1, size(vals,1), length(freqs));
    ticks = ticks(1:4:end);
    set(gca, 'ytick', ticks);
    set(gca, 'yticklabel', num2str(freqs(1:4:end)', 3))

    if ~isempty(g.limits)
        caxis(g.limits);
        limits = g.limits;
    else
        clim = max(abs(caxis));
        limits = [-clim clim];
        caxis(limits);
    end
    xlabel('Time (ms)'); ylabel('Freq (Hz)');
    if ~isempty(g.vert)
        for vert = g.vert(:)'
            hold on; plot([vert vert], [0.01 500], 'k', 'linewidth', 2);
        end
    end
    if strcmpi(g.cbar,'on')
        cbar;
    end
    return;

% -----------------------------------------
% plot time-freq point (when clicking on it
% -----------------------------------------
function plotpoint(data, sortvar, timevect, freq, timepnts);

figure;
subplot(2,2,1);
erpimage(act_all(:,:),sortvar_all,timepnts, ...
   '', 300,10,'erp', 'erpstd', 'caxis',[-1.0 1.0], 'srate',256,'align',352, 'yerplabel', '', erpimopt{:}, ...
                                                   'phasesort', [500 0 5]); % aligned to median rt

% plot in polar coordinates phase versus RT
% -----------------------------------------
phaseang2 = [phsangls phsangls-2*pi]; phaseang2 = movav(phaseang2,[], 100);
outvar2   = [outvar   outvar];   outvar2   = movav(outvar2,[], 100);
phaseang2 = phaseang2(length(phsangls)/2-50:length(phsangls)+length(phsangls)/2-50);
outvar2   = outvar2  (length(phsangls)/2-50:length(phsangls)+length(phsangls)/2-50);
subplot(2,2,3);
polar(phsangls, outvar, '.');
hold on; polar(phaseang2, outvar2, 'r');

% computing ITC
% -------------
cmplx = outvar .* exp(j*phsangls);
ITCval = mean(cmplx);
angle(ITCval)/pi*180;
abs(ITCval)/mean(outvar);
text(-1300,-1400,sprintf('ITC value: amplitude %3.4f, angle %3.1f\n', abs(ITCval)/mean(outvar), angle(ITCval)/pi*180));

% accumulate 200 values
% ---------------------
alpha = 0.05;
if exist('accarray') ~= 1, accarray = NaN; end
[sigval accarray] = bootstat(outvar/mean(outvar), phsangls, 'res = mean(arg1 .* exp(arg2*complex(0,1)))', ...
                      'accarray', accarray, 'bootside', 'upper', 'alpha', alpha);
text(-1300,-1600,sprintf('Threshold for significance at %d percent(naccu=200): %3.4f\n', alpha*100, sigval));
title(sprintf('Clust %d corr. theta phase at 500 ms and RT', clust));

% amplitude sorting
% -----------------
figure;
[outdata,outvar,outtrials,limits,axhndls,erp, ...
          amps,cohers,cohsig,ampsig,outamps,phsangls,phsamp,sortidx]=erpimage(act_all(:,:),sortvar_all,timepnts, ...
   '', 0,0,'erp', 'erpstd', 'caxis', 0.5, 'cbar', ...
   'srate',256,'align',352, 'yerplabel', '', erpimopt{:}, 'ampsort', [500 0 5]); % aligned to median rt
close;
subplot(2,2,2);
erpimage(act_all(:,:),sortvar_all,timevect, ...
   '', 300,10,'erp', 'erpstd', 'caxis', 0.5, 'cbar', ...
   'srate',256,'align',352, 'yerplabel', '', erpimopt{:}, 'ampsort', [500 0 5]); % aligned to median rt

% compute correlation
% -------------------
[ypred alpha Rsq] = myregress(outvar, phsamp);
subplot(2,2,4);
plot(outvar, phsamp, '.');
hold on;
plot(outvar, ypred, 'r');
xlabel('Reaction time');
ylabel('Amplitude');
title(sprintf('Theta amp. at 500 ms vs RT (p=%1.8f, R^2=%3.4f)', alpha, Rsq));
set(gcf, 'position', [336 485 730 540], 'paperpositionmode', 'auto');

setfont(gcf, 'fontsize', 10);
eval(['print -djpeg clust' int2str(clust) 'corrthetart.jpg']);

% set up for command line call
% copy text from plotcorrthetaart
%timevect = times;
sortvar = sortvar_all;
data = act_all;
time = 1
freq = 1
g.times = 0;
g.freqs = 5;
g.erpimopt = {};

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function center = triangle_incenter_2d ( t )

%*****************************************************************************80
%
%% TRIANGLE_INCENTER_2D computes the incenter of a triangle in 2D.
%
%  Discussion:
%
%    The incenter of a triangle is the center of the inscribed circle.
%
%    The inscribed circle of a triangle is the largest circle that can
%    be drawn inside the triangle.
%
%    The inscribed circle is tangent to all three sides of the triangle.
%
%    The angle bisectors of the triangle intersect at the center of the
%    inscribed circle.
%
%    In geometry, the incenter is often represented by "I".
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Adrian Bowyer and John Woodwark,
%    A Programmer's Geometry,
%    Butterworths, 1983.
%
%  Parameters:
%
%    Input, real T(2,3), the triangle vertices.
%
%    Output, real CENTER(2,1), the incenter.
%
  dim_num = 2;
%
%  Compute the length of each side.
%
  a = sqrt ( sum ( ( t(1:dim_num,1) - t(1:dim_num,2) ).^2 ) );
  b = sqrt ( sum ( ( t(1:dim_num,2) - t(1:dim_num,3) ).^2 ) );
  c = sqrt ( sum ( ( t(1:dim_num,3) - t(1:dim_num,1) ).^2 ) );

  perimeter = a + b + c;

  if ( perimeter == 0.0 )
    center(1:dim_num,1) = t(1:dim_num,1);
  else
    center(1:dim_num,1) = ( b * t(1:dim_num,1) ...
                          + c * t(1:dim_num,2) ...
                          + a * t(1:dim_num,3) ) / perimeter;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Check whether a graph is weighted, i.e edges have weights.
%
% INPUTS: edge list, m x 3, m: number of edges, [node 1, node 2, edge weight]
% OUTPUTS: Boolean variable, 0 or 1
%
% GB: last updated, Sep 23, 2012

function S=isWeighted(el)

S=true;

if numel( find(el(:,3)==1) ) == size(el,1); S=false; end

author: aeolianine

repo: octave-networks-toolbox


text:
function avgNEES = ANEES(trans_err, rot_err, err_sigma)

    stateErr = [rot_err;trans_err];
    stateVar = err_sigma.^2;
    avgNEES = 0;
    stepNum = size(stateErr, 2);
    for i = 1:stepNum
        avgNEES = avgNEES + (1/stepNum)*stateErr(:,i)'*inv(diag(stateVar(:,i)))*stateErr(:,i);
    end
end

repo: MSCKF


text:
function t_pf_radial(quiet)
%T_PF_RADIAL  Tests for distribution power flow solvers.

if nargin < 1
    quiet = 0;
end

t_begin(329, quiet);

if quiet
    verbose = 0;
else
    verbose = 1;
end
if have_feature('octave')
    if have_feature('octave', 'vnum') >= 4
        file_in_path_warn_id = 'Octave:data-file-in-path';
    else
        file_in_path_warn_id = 'Octave:load-file-in-path';
    end
    s1 = warning('query', file_in_path_warn_id);
    warning('off', file_in_path_warn_id);
end
mpopt = mpoption('out.all', 0, 'verbose', verbose);

%% define named indices into bus, gen, branch matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;

%% Test Distribution Power Flow
mpopt = mpoption(mpopt, 'pf.alg', 'NR', 'pf.radial.max_it', 100);
casefile = {
    'case4_dist'
    'case18'
    'case22'
    'case33bw'
    'case69'
    'case85'
    'case141'
    };
method = {
    'PQSUM' 'Power Summation'
    'ISUM'  'Current Summation'
    'YSUM'  'Admittance Summation'
    };

% Original test cases (no PV buses)
iter_expected = {
           'none'    'NR' 'PQSUM' 'ISUM' 'YSUM'
     'case4_dist'      3      12     12     11
         'case18'      4       9     12     10
         'case22'      3       3      6      6
       'case33bw'      3       5      8      8
         'case69'      4       5      8      8
         'case85'      4       5      9      9
        'case141'      3       4      7      7
        };
iter_nopv = zeros(size(casefile,1),size(method,1)+1);
for i = 1:length(casefile)
    % which row in iter_expected is for current casefile
    row = find(strcmp(casefile{i}, iter_expected(:,1)) == 1);
    mpc = loadcase(casefile{i});
    % solve it with NR
    r = runpf(mpc, mpopt);
    iter_nopv(i,1) = r.iterations;
    % which columng in iter_expected is for current method
    col = find(strcmp('NR', iter_expected(1,:)) == 1);
    t = ['Newton''s method, ' casefile{i} ' : '];
    t_is(r.iterations, iter_expected{row,col}, 6, [t 'iterations']);
    for j = 1:size(method,1)
        mpopt1 = mpoption(mpopt,'pf.alg',method{j,1});
        r1 = runpf(mpc,mpopt1);
        iter_nopv(i,j+1) = r1.iterations;
        % which columng in iter_expected is for current method
        col = find(strcmp(method{j,1}, iter_expected(1,:)) == 1);
        t = [method{j,2} ', ' casefile{i} ' : '];
        t_ok(r1.success, [t 'success']);
        t_is(r1.iterations, iter_expected{row,col}, 6, [t 'iterations']);
        t_is(r1.bus, r.bus, 6, [t 'bus']);
        t_is(r1.gen, r.gen, 6, [t 'gen']);
        t_is(r1.branch, r.branch, 6, [t 'branch']);
    end
end

% Test cases with added PV buses
clear iter_expected
iter_expected{1} = {
% pf.radial.vcorr = 0
           'none'    'NR' 'PQSUM' 'ISUM' 'YSUM'
     'case4_dist'      3      12     12     11
         'case18'      4      20     16     13
         'case22'      4      11     25     25
       'case33bw'      4       8     16     16
         'case69'      4      12     27     27
         'case85'      5       9     17     17
        'case141'      5       9     18     18
};
iter_expected{2} = {
% pf.radial.vcorr = 1
           'none'    'NR' 'PQSUM' 'ISUM' 'YSUM'
     'case4_dist'      3       7      6      7
         'case18'      4      10     12     12
         'case22'      4      12     14     14
       'case33bw'      4       8     11     11
         'case69'      4      13     17     17
         'case85'      5       9     13     13
        'case141'      5       9     12     12
};
iter_pv = zeros(size(casefile,1),size(method,1)+1,2);
for i = 1:length(casefile)
    % which row in iter_expected is for current casefile
    row = find(strcmp(casefile{i}, iter_expected{1}(:,1)) == 1);
    mpc = loadcase(casefile{i});
    pv = mpc.bus(:,BUS_TYPE) == PV;
    if all(~pv)
        %%% ADD PV BUSES %%%
        N = 5; % How many PV buses should I add?
        % Solve the case without PV buses with Newton method
        r = runpf(mpc, mpopt);
        % Find the last N buses with lowest voltage
        [Vm, B] = sort(r.bus(:,VM));
        B = B(1:N);
        Vm = Vm(1:N);
        % add PV generators at buses in B
        % set VG 0.05 pu bigger then voltage at buses in B
        mpc.gen = repmat(mpc.gen,N+1,1);
        mpc.gen(2:end,GEN_BUS) = mpc.bus(B, BUS_I);
        mpc.gen(2:end,VG) = Vm + 0.05;
        mpc.bus(B,BUS_TYPE) = 2;
        if isfield(mpc,'gencost')
            mpc.gencost = repmat(mpc.gencost,N+1,1);
        end
        %%% END OF ADD PV BUSES %%%
    end
    r = runpf(mpc, mpopt);
    iter_pv(i,1,1) = r.iterations;
    iter_pv(i,1,2) = r.iterations;
    % which columng in iter_expected is for current method
    col = find(strcmp('NR', iter_expected{1}(1,:)) == 1);
    t = ['Newton''s method, ' casefile{i} ' : '];
    t_is(r.iterations, iter_expected{1}{row,col}, 6, [t 'iterations']);
    for j = 1:size(method,1)
        for vcorr = 0:1
            mpopt1 = mpoption(mpopt,'pf.alg',method{j,1},'pf.radial.vcorr',vcorr);
            r1 = runpf(mpc,mpopt1);
            iter_pv(i,j+1,vcorr+1) = r1.iterations;
            % which columng in iter_expected is for current method
            col = find(strcmp(method{j,1}, iter_expected{vcorr+1}(1,:)) == 1);
            t = [method{j,2} ', vcorr = ' num2str(vcorr) ', ' casefile{i} ' : '];
            t_ok(r1.success, [t 'success']);
            t_is(r1.iterations, iter_expected{vcorr+1}{row,col}, 6, [t 'iterations']);
            t_is(r1.bus, r.bus, 6, [t 'bus']);
            t_is(r1.gen, r.gen, 6, [t 'gen']);
            t_is(r1.branch, r.branch, 6, [t 'branch']);
        end
    end
end

t_end;

if verbose > 0
    fprintf('\nITERATIONS: Original test cases (no PV buses), pf.radial.vcorr = 0\n');
    print_iterations(casefile,method,iter_nopv)
    fprintf('\nITERATIONS: Test cases with added PV buses, pf.radial.vcorr = 0\n');
    print_iterations(casefile,method,iter_pv(:,:,1))
    fprintf('\nITERATIONS: Test cases with added PV buses, pf.radial.vcorr = 1\n');
    print_iterations(casefile,method,iter_pv(:,:,2))
end

if have_feature('octave')
    warning(s1.state, file_in_path_warn_id);
end

function print_iterations(casefile,method,iterations)
fprintf('%22s','NR');
for j = 1:size(method,1)
    fprintf('%7s',method{j,1});
end
fprintf('\n');
for i = 1:length(casefile)
    fprintf('%15s',casefile{i});
    for j = 1:size(method,1)+1
        fprintf('%7i',iterations(i,j));
    end
    fprintf('\n');
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%  doinvers
% This file calculates orintation of the stress tensor
% based on Gephard's algorithm.
% stress tensor orientation. The actual calculation is done
% using a call to a fortran program.
%
% Stefan Wiemer 03/96

global mi mif1 mif2  hndl3 a newcat2 mi2
global tmp cumu2
report_this_filefun(mfilename('fullpath'));
think

% select the gridpoints

% do the loop

i2 = 0;
clear howmany
for i= 1:length(newgri(:,1))
    x = newgri(i,1);y = newgri(i,2);
    allcount = allcount + 1.;
    i2 = i2+1;

    % calculate distance from center point and sort wrt distance
    l = sqrt(((xsecx' - x)).^2 + ((xsecy + y)).^2) ;
    [s,is] = sort(l);
    b = newa(is(:,1),:) ;       % re-orders matrix to agree row-wise
    b = b(1:ni,:);

    tmp = [b(:,10:14)];
    comm = ['!mkdir /home/david/tmp/in' num2str(i2)]
    eval(comm)
    comm = ['save /home/david/tmp/in' num2str(i2) '/data.inp tmp -ascii']
    eval(comm)
    comm = ['save /home/david/tmp/in' num2str(i2) '/loc.inp  x y  -ascii']
    eval(comm)
end
save tmpinv.mat

com =['gps        '
    'ugle       '
    'dutton     '
    'kiska      '
    'moment     '
    %'nordic     '
    'kanaga     '
    'spurr      '
    'chaos      '
    'megathrust '
    %'pele       '
    %'geoid      '
    'model      '
    'marvin     '];

cd /home/stefan/ZMAP/invers
i = 1;

while i <  length(newgri(:,1))

    for k = 1:length(com(:,1))

        comm = ['!/bin/rm /home/stefan/ZMAP/invers/howmany.dat '];
        eval(comm)
        comm = ['!rsh ' com(k,:) ' ps -axuw | grep fmsiWindow | cut -c10-15 ',...
            ' > /home/stefan/ZMAP/invers/howmany.dat '];
        eval(comm)
        pause(3)
        comm =['load /home/stefan/ZMAP/invers/howmany.dat '];
        eval(comm,'disp('' no'')')

        if ~exist('howmany','var'); howmany = 0; end
        if com(k,:) == 'moment     ' ; howmany(1:2) = [] ; end
        if length(howmany(:,1)) < 2
            comm = ['! rsh ' com(k,:) ' /home/stefan/ZMAP/invers/invshell ',...
                num2str(length(tmp(:,1))) ' ' num2str(i)  ' &']
            eval(comm)
            i = i+1
            clear howmany
            pause(3)
        end   % if  length
        clear howmany
    end   % for k

    % wait for 1 minutes
    pause(10)

end  % while i

author: CelsoReyes


text:
function [ns_RESULT, nsNeuralInfo] = ns_GetNeuralInfo(hFile, EntityID);

%ns_GetNeuralInfo   Retrieves information for neural event entities
%
%   Usage:
%      [ns_RESULT, nsNeuralInfo] = ns_GetNeuralInfo(hFile, EntityID)
%
%   Description:
%       Retrieves information on Neural Event entity EntityID from the
%       file referenced by hFile.  The information is returned in the
%       nsNeuralInfo.
%
%   Parameters:
%       hFile	    Handle/Indentification number to an open file.
%       EntityID	Identification number of the entity in the data file.
%
%   Return Values:
%       nsNeuralInfo    ns_NEURALINFO structure to receive the Neural
%                       Event information.
%       ns_RESULT   This function returns ns_OK if the file is successfully
%                   opened. Otherwise one of the following error codes is
%                   generated:
%
%                       ns_BADFILE	    Invalid file handle passed to
%                                       function
%                       ns_BADENTITY	Invalid or inappropriate entity
%                                       identifier specified
%                       ns_FILEERROR	File access or read error
%
%   Copyright (C) 2003 Neuroshare Project
%   Author: Almut Branner
%   Last modification: 10/24/2003

[ns_RESULT, nsNeuralInfo] = mexprog(12, hFile, EntityID - 1);

SourceEntityID = [nsNeuralInfo.SourceEntityID] + 1;
SourceUnitID = [nsNeuralInfo.SourceUnitID];

ind = find(SourceUnitID == 1);
SourceUnitID(ind) = 255;

ind = find((SourceUnitID > 1) & (SourceUnitID < 255));
SourceUnitID(ind) = log2(SourceUnitID(ind));

for i = 1:length(nsNeuralInfo)
    nsNeuralInfo(i).SourceEntityID = SourceEntityID(i);
    nsNeuralInfo(i).SourceUnitID = SourceUnitID(i);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [km] = au2km(au)
% Convert length from astronomical units to kilometers.
% Chad A. Greene 2012
km = au*1.49597870691e+8;

author: Sable

repo: mcbench-benchmarks


text:
function midi = readmidi(filename, rawbytes)
% midi = readmidi(filename, rawbytes)
% midi = readmidi(filename)
%
% Read MIDI file and store in a Matlab structure
% (use midiInfo.m to see structure detail)
%
% Inputs:
%  filename - input MIDI file
%  rawbytes - 0 or 1: Include raw bytes in structure
%             This info is redundant, but can be
%             useful for debugging. default=0
%

% Copyright (c) 2009 Ken Schutte
% more info at: http://www.kenschutte.com/midi

if (nargin<2)
  rawbytes=0;
end

fid = fopen(filename);
[A count] = fread(fid,'uint8');
fclose(fid);

if (rawbytes) midi.rawbytes_all = A; end

% realtime events: status: [F8, FF].  no data bytes
%clock, undefined, start, continue, stop, undefined, active
%sensing, systerm reset

% file consists of "header chunk" and "track chunks"
%   4B  'MThd' (header) or 'MTrk' (track)
%   4B  32-bit unsigned int = number of bytes in chunk, not
%       counting these first 8

% HEADER CHUNK --------------------------------------------------------
% 4B 'Mthd'
% 4B length
% 2B file format
%    0=single track, 1=multitrack synchronous, 2=multitrack asynchronous
%    Synchronous formats start all tracks at the same time, while asynchronous formats can start and end any track at any time during the score.
% 2B track cout (must be 1 for format 0)
% 2B num delta-time ticks per quarter note
%

if ~isequal(A(1:4)',[77 84 104 100])  % double('MThd')
    error('File does not begin with header ID (MThd)');
end

header_len = decode_int(A(5:8));
if (header_len == 6)
else
    error('Header length != 6 bytes.');
end

format = decode_int(A(9:10));
if (format==0 || format==1 || format==2)
     midi.format = format;
else
    error('Format does not equal 0,1,or 2');
end

num_tracks = decode_int(A(11:12));
if (format==0 && num_tracks~=1)
    error('File is format 0, but num_tracks != 1');
end

time_unit = decode_int(A(13:14));
if (bitand(time_unit,2^15)==0)
  midi.ticks_per_quarter_note = time_unit;
else
  error('Header: SMPTE time format found - not currently supported');
end

if (rawbytes)
  midi.rawbytes_header = A(1:14);
end

% end header parse ----------------------------------------------------

% BREAK INTO SEPARATE TRACKS ------------------------------------------
% midi.track(1).data = [byte byte byte ...];
% midi.track(2).date = ...
% ...
%
% Track Chunks---------
% 4B 'MTrk'
% 4B length (after first 8B)
%
ctr = 15;
for i=1:num_tracks

  if ~isequal(A(ctr:ctr+3)',[77 84 114 107])  % double('MTrk')
    error(['Track ' num2str(i) ' does not begin with track ID=MTrk']);
  end
  ctr = ctr+4;

  track_len = decode_int(A(ctr:ctr+3));
  ctr = ctr+4;

  % have track.rawbytes hold initial 8B also...
  track_rawbytes{i} = A((ctr-8):(ctr+track_len-1));

  if (rawbytes)
    midi.track(i).rawbytes_header = A(ctr-8:ctr-1);
  end

  ctr = ctr+track_len;
end
% ----------------------------------------------------------------------

% Events:
%  - meta events: start with 'FF'
%  - MIDI events: all others

% MIDI events:
%  optional command byte + 0,1,or 2 bytes of parameters
%  "running mode": command byte omitted.
%
% all midi command bytes have MSB=1
% all data for inside midi command have value <= 127 (ie MSB=0)
% -> so can determine running mode
%
% meta events' data may have any values (meta events have to set
% len)
%

% 'Fn' MIDI commands:
%  no chan. control the entire system
%F8 Timing Clock
%FA start a sequence
%FB continue a sequence
%FC stop a sequence

% Meta events:
%  1B 0xFF
%  1B event type
%  1B length of additional data
%  ?? additional data
%

% "channel mode messages"
% have same code as "control change": 0xBn
%  but uses reserved controller numbers 120-127
%

%Midi events consist of an optional command byte
% followed by zero, one or two bytes of parameters.
% In running mode, the command can be omitted, in
% which case the last MIDI command specified is
% assumed.  The first bit of a command byte is 1,
% while the first bit of a parameter is always 0.
%   In addition, the last 4 bits of a command
%   indicate the channel to which the event should
%   be sent; therefore, there are 6 possible
%   commands (really 7, but we will discuss the x'Fn'
%   commands later) that can be specified.  They are:

% parse tracks -----------------------------------------
for i=1:num_tracks

  track = track_rawbytes{i};

  if (rawbytes); midi.track(i).rawbytes = track; end

  msgCtr = 1;
  ctr=9;  % first 8B were MTrk and length
  while (ctr < length(track_rawbytes{i}))

    clear currMsg;
    currMsg.used_running_mode = 0;
    % note:
    %  .used_running_mode is necessary only to
    %  be able to reconstruct a file _exactly_ from
    %  the 'midi' structure.  this is helpful for
    %  debugging since write(read(filename)) can be
    %  tested for exact replication...
    %

    ctr_start_msg = ctr;

    [deltatime,ctr] = decode_var_length(track, ctr);

    % ?
    %if (rawbytes)
    %  currMsg.rawbytes_deltatime = track(ctr_start_msg:ctr-1);
    %end

    % deltaime must be 1-4 bytes long.
    % could check here...

    % CHECK FOR META EVENTS ------------------------
    % 'FF'
    if track(ctr)==255

      type = track(ctr+1);

      ctr = ctr+2;

      % get variable length 'length' field
      [len,ctr] = decode_var_length(track, ctr);

      % note: some meta events have pre-determined lengths...
      %  we could try verifiying they are correct here.

      thedata = track(ctr:ctr+len-1);
      chan = [];

      ctr = ctr + len;

      midimeta = 0;

    else
      midimeta = 1;
      % MIDI EVENT ---------------------------

      % check for running mode:
      if (track(ctr)<128)

	% make it re-do last command:
	%ctr = ctr - 1;
	%track(ctr) = last_byte;
	currMsg.used_running_mode = 1;

	B = last_byte;
	nB = track(ctr); % ?

      else

	B  = track(ctr);
	nB = track(ctr+1);

	ctr = ctr + 1;

      end

      % nibbles:
      %B  = track(ctr);
      %nB = track(ctr+1);

      Hn = bitshift(B,-4);
      Ln = bitand(B,15);

      chan = [];

      msg_type = midi_msg_type(B,nB);

      % DEBUG:
      if (i==2)
	  if (msgCtr==1)
	    disp(msg_type);
	  end
      end

      switch msg_type

       case 'channel_mode'

	% UNSURE: if all channel mode messages have 2 data byes (?)
	type = bitshift(Hn,4) + (nB-120+1);
	thedata = track(ctr:ctr+1);
	chan = Ln;

	ctr = ctr + 2;

	% ---- channel voice messages:
       case 'channel_voice'

	type = bitshift(Hn,4);
	len = channel_voice_msg_len(type); % var length data:
	thedata = track(ctr:ctr+len-1);
	chan = Ln;

	% DEBUG:
	if (i==2)
	  if (msgCtr==1)
	    disp([999  Hn type])
	  end
	end

	ctr = ctr + len;

       case 'sysex'

	% UNSURE: do sysex events (F0-F7) have
	%  variable length 'length' field?

	[len,ctr] = decode_var_length(track, ctr);

	type = B;
	thedata = track(ctr:ctr+len-1);
	chan = [];

	ctr = ctr + len;

       case 'sys_realtime'

	% UNSURE: I think these are all just one byte
	type = B;
	thedata = [];
	chan = [];

      end

      last_byte = Ln + bitshift(Hn,4);

    end % end midi event 'if'

    currMsg.deltatime = deltatime;
    currMsg.midimeta = midimeta;
    currMsg.type = type;
    currMsg.data = thedata;
    currMsg.chan = chan;

    if (rawbytes)
      currMsg.rawbytes = track(ctr_start_msg:ctr-1);
    end

    midi.track(i).messages(msgCtr) = currMsg;
    msgCtr = msgCtr + 1;

  end % end loop over rawbytes
end % end loop over tracks

function val=decode_int(A)

val = 0;
for i=1:length(A)
  val = val + bitshift(A(length(A)-i+1), 8*(i-1));
end

function len=channel_voice_msg_len(type)

if     (type==128); len=2;
elseif (type==144); len=2;
elseif (type==160); len=2;
elseif (type==176); len=2;
elseif (type==192); len=1;
elseif (type==208); len=1;
elseif (type==224); len=2;
else
  disp(type); error('bad channel voice message type');
end

%
% decode variable length field (often deltatime)
%
%  return value and new position of pointer into 'bytes'
%
function [val,ptr] = decode_var_length(bytes, ptr)

keepgoing=1;
val = 0;
while (keepgoing)
  % check MSB:
  %  if MSB=1, then delta-time continues into next byte...
  if(~bitand(bytes(ptr),128))
    keepgoing=0;
  end
  % keep appending last 7 bits from each byte in the deltatime:
  val = val*128 + rem(bytes(ptr), 128);
  ptr=ptr+1;
end

%
% Read first 2 bytes of msg and
%  determine the type
%  (most require only 1st byte)
%
% str is one of:
%  'channel_mode'
%  'channel_voice'
%  'sysex'
%  'sys_realtime'
%
function str=midi_msg_type(B,nB)

Hn = bitshift(B,-4);
Ln = bitand(B,7);

% ---- channel mode messages:
%if (Hn==11 && nB>=120 && nB<=127)
if (Hn==11 && nB>=122 && nB<=127)
  str = 'channel_mode';

  % ---- channel voice messages:
elseif (Hn>=8 && Hn<=14)
  str = 'channel_voice';

  %  ---- sysex events:
elseif (Hn==15 && Ln>=0 && Ln<=7)
  str = 'sysex';

  % system real-time messages
elseif (Hn==15 && Ln>=8 && Ln<=15)
  % UNSURE: how can you tell between 0xFF system real-time
  %   message and 0xFF meta event?
  %   (now, it will always be processed by meta)
  str = 'sys_realtime';

else
  % don't think it can get here...
  error('bad midi message');
end

author: kts

repo: matlab-midi

save_path: github-repos/MATLAB/kts-matlab-midi


text:
% This file is part of the project NILM-Eval (https://github.com/beckel/nilm-eval).
% Licence: GPL 2.0 (http://www.gnu.org/licenses/gpl-2.0.html)
% Copyright: ETH Zurich, 2014
% Author: Romano Cicchetti

function [best_clustersOfFsm] = find_best_variation(events_FSMs_matrix, clusterOfEvents, valid_cluster_permutations, fsm_idx, events, maxSequenceLength)

    best_deviation_from_median = intmax;
    best_clustersOfFsm = valid_cluster_permutations(1,:);

    for variation = 1:size(valid_cluster_permutations)
        clustersOfFsm = valid_cluster_permutations(variation, :);
        for i=1:length(clustersOfFsm)
            eventsOfFsm = clusterOfEvents == clustersOfFsm(1,i);
            events_FSMs_matrix(eventsOfFsm, fsm_idx) = i;
        end

        % build event sequences of a finte state machine
        sequencesOfFsm = buildSequences(events_FSMs_matrix, events, fsm_idx, maxSequenceLength);
        events_FSMs_matrix(eventsOfFsm, fsm_idx) = 0;
        if isempty(sequencesOfFsm)
            continue;
        end

        % compute median duration of all event sequences that have a high
        % probability of belonging to the shortest path
        duration_of_sequences = events(sequencesOfFsm(:,2),4) - events(sequencesOfFsm(:,1),4);
        if nnz(sequencesOfFsm(:,3) == 1) > 1
            median_duration_of_sequences = median(duration_of_sequences(sequencesOfFsm(:,3) == 1,1));
        else
            median_duration_of_sequences = median(duration_of_sequences(:,1));
        end

        % build directed graph with event sequences as nodes and
        % compute the shortest path
        [directed_graph] = buildGraph(sequencesOfFsm, duration_of_sequences, median_duration_of_sequences);
        num_of_sequences = size(sequencesOfFsm,1);
        [~, path, ~] = graphshortestpath(directed_graph, 1, num_of_sequences, 'Method', 'Bellman-Ford');
        median_duration_of_sequences_on_path = median(duration_of_sequences(path, 1));
        deviation_from_median = abs((duration_of_sequences(path, 1) - median_duration_of_sequences_on_path)/median_duration_of_sequences_on_path);
        quality = deviation_from_median .* log(deviation_from_median);
        quality(isnan(quality)) = -10;

        if sum(quality) < best_deviation_from_median
            best_clustersOfFsm = valid_cluster_permutations(variation, :);
            best_deviation_from_median = deviation_from_median;
        end
    end

end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function [e1,e2] = ExplNoise(t)
%% --------------------------------------------------------------
% ExplNoise:
% Function to generate Exploration Noise
% ------------------------------------------------------------------------
e1 = 10*sum(sin([38.1558   76.5517   79.5200   18.6873   48.9764   44.5586   64.6313   70.9365   75.4687   27.6025]*t));
e2 = 20*sum(sin([17.9703   15.5098  -33.7388  -38.1002   100   45.9744  -15.9614    8.5268  -27.6188   25.1267]*t));
end

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
% jCell:  Tools for operating on cell arrays of column vectors
%
% Basic mathematical operations
%   cellabs    - Absolute value of each element in a cell array.
%   cellmax    - Maximum of each element in a cell array.
%   cellmin    - Minimum of each element in a cell array.
%   cellsum    - Sum of each element a cell array, possibly weighted.
%   cellmean   - Mean value of each element a cell array, possibly weighted.
%   cellstd    - Standard deviation of each element a cell array, possibly weighted.
%   cellmed    - Median value of each element a cell array.
%   cellreal   - Real part of each element in a cell array.
%   cellimag   - Imaginary part of each element in a cell array.
%   cellpair   - Complex pairing for elements in a cell array.
%   celllog10  - Base ten logarithm of each element in a cell array.
%   celladd    - Addition acting on each element in a cell array.
%   cellmult   - Multiplication acting on each element in a cell array.
%   celldiv    - Division acting on each element in a cell array.
%   celldot    - Dot product for arrays of column vectors.
%
% Reshaping, indexing, and sizes
%   cell2col   - Converts cell arrays of numeric arrays into 'column-appended' form.
%   col2cell   - Converts 'column-appended' data into cell arrays of numeric arrays.
%   cellindex  - Applies a cell array of indices to a cell array of column vectors.
%   cellchunk  - Converts cell array data into uniform length 'chunks'.
%   cellength  - Length of each element in a cell array.
%   cellsize   - Size of each element in a cell array along specified dimension.
%   cellget    - Indexes a cell array of numerical arrays by ID number.
%   cellimit   - Limits the ranges of times in a cell array of numerical arrays.
%
% Data processing
%   cellstrip  - Strips INF values from the beginnings or ends of cell arrays.
%   cellsplit  - Splits cell arrays of numeric arrays at data gaps.
%   cellprune  - Removes all empty cells, or cells less than a specified length.
%   cellpack   - Removes all INF values from cell arrays of numeric arrays.
%   cellfill   - Fills missing data marked by INFs in a cell array.
%   cellgrid   - Interpolate a cell array of numeric arrays onto a regular grid.
%   cellfirst  - Returns the first element of each entry in a cell array.
%
% Plotting
%   cellplot   - Rapidly plot all element of a cell array of numeric arrays.
%
% See also jVarfun.
help jCell

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function truefalse = ch_is_format_code ( c )

%*****************************************************************************80
%
%% CH_IS_FORMAT_CODE returns TRUE if a character is a FORTRAN format code.
%
%  Discussion:
%
%    The format codes accepted here are not the only legal format
%    codes in FORTRAN90.  However, they are more than sufficient
%    for my needs!
%
%  Table:
%
%    A  Character
%    B  Binary digits
%    D  Real number, exponential representation
%    E  Real number, exponential representation
%    F  Real number, fixed point
%    G  General format
%    I  Integer
%    L  Logical variable
%    O  Octal digits
%    Z  Hexadecimal digits
%    *  Free format
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, the character to be analyzed.
%
%    Output, logical TRUEFALSE, is TRUE (1) if C is a FORTRAN format code
%     and FALSE (0) otherwise.
%
  FALSE = 0;
  TRUE = 1;

      if ( ch_eqi ( c, 'A' ) )
    truefalse = TRUE;
  elseif ( ch_eqi ( c, 'B' ) )
    truefalse = TRUE;
  elseif ( ch_eqi ( c, 'D' ) )
    truefalse = TRUE;
  elseif ( ch_eqi ( c, 'E' ) )
    truefalse = TRUE;
  elseif ( ch_eqi ( c, 'F' ) )
    truefalse = TRUE;
  elseif ( ch_eqi ( c, 'G' ) )
    truefalse = TRUE;
  elseif ( ch_eqi ( c, 'I' ) )
    truefalse = TRUE;
  elseif ( ch_eqi ( c, 'L' ) )
    truefalse = TRUE;
  elseif ( ch_eqi ( c, 'O' ) )
    truefalse = TRUE;
  elseif ( ch_eqi ( c, 'Z' ) )
    truefalse = TRUE;
  elseif ( c == '*' )
    truefalse = TRUE;
  else
    truefalse = FALSE;
  end

  return
end


text:
function pcomb = stat_stouffer_pcomb(p)
% Stouffer et al's (1949) unweighted method for combination of
% independent p-values via z's
% From: http://imaging.mrc-cbu.cam.ac.uk/statswiki/FAQ/CombiningPvalues

if isempty(p)
    error('pfast was passed an empty array of p-values')
    pcomb=1;
else
    pcomb = (1-erf(sum(sqrt(2) * erfinv(1-2*p))/sqrt(2*length(p))))/2;
end

repo: REST


text:
% Fig. 6.10   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%

clear all;
%close all;
clf

num=10;
den=conv([1 0],[1 0.4 4]);
w=logspace(-1,2,100);
[m,p]=bode(num,den,w);

figure(1)
loglog(w,m);
xlabel('\omega (rad/sec)');
ylabel('Magnitude');
title('Fig. 6.10 Bode plot for a TF with a complex pole :(a) magnitude');
bodegrid;
%pause;
figure(2)
semilogx(w,p);
grid;
xlabel('\omega (rad/sec)');
ylabel('Phase');
title('Fig. 6.10 (b) phase');
bodegrid;


text:
function taulist = InverseDynamics(thetalist, dthetalist, ddthetalist, ...
                                   g, Ftip, Mlist, Glist, Slist)
% *** CHAPTER 8: DYNAMICS OF OPEN CHAINS ***
% Takes thetalist: n-vector of joint variables,
%       dthetalist: n-vector of joint rates,
%       ddthetalist: n-vector of joint accelerations,
%       g: Gravity vector g,
%       Ftip: Spatial force applied by the end-effector expressed in frame
%             {n+1},
%       Mlist: List of link frames {i} relative to {i-1} at the home
%              position,
%       Glist: Spatial inertia matrices Gi of the links,
%       Slist: Screw axes Si of the joints in a space frame, in the format
%              of a matrix with the screw axes as the columns.
% Returns taulist: The n-vector of required joint forces/torques.
% This function uses forward-backward Newton-Euler iterations to solve the
% equation:
% taulist = Mlist(thetalist) * ddthetalist + c(thetalist, dthetalist) ...
%           + g(thetalist) + Jtr(thetalist) * Ftip
% Example Input (3 Link Robot):
%
% clear; clc;
% thetalist = [0.1; 0.1; 0.1];
% dthetalist = [0.1; 0.2; 0.3];
% ddthetalist = [2; 1.5; 1];
% g = [0; 0; -9.8];
% Ftip = [1; 1; 1; 1; 1; 1];
% M01 = [[1, 0, 0, 0]; [0, 1, 0, 0]; [0, 0, 1, 0.089159]; [0, 0, 0, 1]];
% M12 = [[0, 0, 1, 0.28]; [0, 1, 0, 0.13585]; [-1, 0 ,0, 0]; [0, 0, 0, 1]];
% M23 = [[1, 0, 0, 0]; [0, 1, 0, -0.1197]; [0, 0, 1, 0.395]; [0, 0, 0, 1]];
% M34 = [[1, 0, 0, 0]; [0, 1, 0, 0]; [0, 0, 1, 0.14225]; [0, 0, 0, 1]];
% G1 = diag([0.010267, 0.010267, 0.00666, 3.7, 3.7, 3.7]);
% G2 = diag([0.22689, 0.22689, 0.0151074, 8.393, 8.393, 8.393]);
% G3 = diag([0.0494433, 0.0494433, 0.004095, 2.275, 2.275, 2.275]);
% Glist = cat(3, G1, G2, G3);
% Mlist = cat(3, M01, M12, M23, M34);
% Slist = [[1; 0; 1;      0; 1;     0], ...
%        [0; 1; 0; -0.089; 0;     0], ...
%        [0; 1; 0; -0.089; 0; 0.425]];
% taulist = InverseDynamics(thetalist, dthetalist, ddthetalist, g, ...
%                         Ftip, Mlist, Glist, Slist)
%
% Output:
% taulist =
%   74.6962
%  -33.0677
%   -3.2306

n = size(thetalist, 1);
Mi = eye(4);
Ai = zeros(6, n);
AdTi = zeros(6, 6, n + 1);
Vi = zeros(6, n + 1);
Vdi = zeros(6, n + 1);
Vdi(4: 6, 1) = -g;
AdTi(:, :, n + 1) = Adjoint(TransInv(Mlist(:, :, n + 1)));
Fi = Ftip;
taulist = zeros(n, 1);
for i=1: n
    Mi = Mi * Mlist(:, :, i);
    Ai(:, i) = Adjoint(TransInv(Mi)) * Slist(:, i);
    AdTi(:, :, i) = Adjoint(MatrixExp6(VecTose3(Ai(:, i) ...
                    * -thetalist(i))) * TransInv(Mlist(:, :, i)));
    Vi(:, i + 1) = AdTi(:, :, i) * Vi(:, i) + Ai(:, i) * dthetalist(i);
    Vdi(:, i + 1) = AdTi(:, :, i) * Vdi(:, i) ...
                    + Ai(:, i) * ddthetalist(i) ...
                    + ad(Vi(:, i + 1)) * Ai(:, i) * dthetalist(i);
end
for i = n: -1: 1
    Fi = AdTi(:, :, i + 1)' * Fi + Glist(:, :, i) * Vdi(:, i + 1) ...
         - ad(Vi(:, i + 1))' * (Glist(:, :, i) * Vi(:, i + 1));
    taulist(i) = Fi' * Ai(:, i);
end
end

author: ShuoYangRobotics

repo: QuadrupedSim


text:
% Lifted matrices for MPC on horizon N
function [Ab Bb] = createMPCmatrices(A,B,N,includex0)
    if(~exist('includex0','var'))
        includex0 = 0; %first block rows of Ab and Bb corresponding to x0 removed
    end

    n = size(A,1);
    Ab = zeros((N+1)*n, n);
    Ab(1:n,:) = eye(n,n);
    for i = 2:N+1
        if( size(A,1) > 1000 )
            i
        end
        %Ab((i-1)*n+1:i*n,:) = A^(i-1);
        Ab((i-1)*n+1:i*n,:) = Ab((i-2)*n+1:(i-1)*n,:)*A;
    end
    if(includex0 == 0)
        Ab = Ab(n+1:end,:); % Seems to take a lot of time if A is huge, possible to speedup
    end

    for q = 1:length(B)
        m = size(B{q},2);
        Bb{q} = zeros((N+1)*n, N*m);
        for i = 2:N+1
            Bb{q}((i-1)*n+1:i*n,:) = A * Bb{q}((i-2)*n+1:(i-1)*n,:);
            Bb{q}((i-1)*n+1:n*i,(i-2)*m+1:m*(i-1)) = B{q};
        end
        if(includex0 == 0)
            Bb{q} = Bb{q}(n+1:end,:);
        end
    end

end

author: arbabiha


text:
% eegplugin_mffmatlabio() - plugin for importing and exporting MFF files
%
% Usage:
%   >> eegplugin_mffmatlabio(fig, trystrs, catchstrs);
%
% Inputs:
%   fig        - [integer] eeglab figure.
%   trystrs    - [struct] "try" strings for menu callbacks.
%   catchstrs  - [struct] "catch" strings for menu callbacks.
%
% Main files
% ----------
% eegplugin_mffmatlabio.m - create menu in the EEGLAB interface
% pop_mffimport.m     - import MFF file into EEGLAB (graphic interface)
% pop_mffexport.m     - export EEGLAB structure to MFF file/folder (graphic interface)
%
% Support files
% -------------
% mff_import.m          - import MFF file/folder to EEGLAB structure
% mff_importcategories.m   - import 'categories.xml' file
% mff_importcoordinates.m  - import 'coordinates.xml' file
% mff_importepochs.m       - import 'epochs.xml' file
% mff_importevents.m       - import 'eventsxxxx.xml' file(s)
% mff_importinfo.m         - import 'info.xml' file
% mff_importinfon.m        - import 'info1.xml' file
% mff_importpnsset.m       - import PNS file
% mff_importsensorlayout.m - import 'sensorlayout.xml' file
% mff_importsignal.m       - import 'signal1.bin' file
% mff_importsubject.m      - import subject file
% mff_export.m             - export MFF file/folder from EEGLAB structure
% mff_createmff.m          - create empty MFF file/folder
% mff_exportcategories.m   - export 'categories.xml' file
% mff_exportcoordinates.m  - export 'coordinates.xml' file
% mff_exportepochs.m       - export 'epochs.xml' file
% mff_exportevents.m       - export 'events.xml' file
% mff_exportinfo.m         - export 'info.xml' file
% mff_exportinfon.m        - export 'info1.xml' file
% mff_exportpnsset.m       - export PNS file
% mff_exportsensorlayout.m - export 'sensorlayout.xml' file
% mff_exportsignal.m       - export 'signal1.bin' file
% mff_exportsubject.m      - export subject file
% mff_decodetime.m         - convert MFF time to Matlab time
% mff_encodetime.m         - convert Matlab time to MFF time
% mff_getobj.m             - convert Java object to Matlab
% mff_setobj.m             - convert Matlab object to Java
% pop_mffimport.m          - EEGLAB import GUI
% pop_mffexport.m          - EEGLAB export GUI
% eegplugin_mffmatlabio.m  - EEGLAB startup function
% eeg_compare.m            - Function to compare EEGLAB sturtures
% mff_fileio_read_data.m   - File-IO function to read data
% mff_fileio_read_header.m - File-IO function to read header
% mff_fileio_read_event.m  - File-IO function to read event
% mff_fileio_write.m       - File-IO function to write data and events

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

function versionstr = eegplugin_mffmatlabio(fig, trystrs, catchstrs)

    %global EEG
    versionstr = '3.4';
    if nargin < 3
        disp('eegplugin_mffmatlabio requires 3 arguments');
        return
    end

    % add amica folder to path
    % -----------------------
    if ~exist('mff_import')
        p = which('eegplugin_mffmatlabio');
        p = p(1:findstr(p,'eegplugin_mffmatlabio.m')-1);
        addpath(p);
    end

    % find tools menu
    % ---------------
    menui = findobj(fig, 'tag', 'import data');
    menue = findobj(fig, 'tag', 'export');

    % menu callback commands
    % ----------------------
    comload    = [  trystrs.no_check '[EEG, LASTCOM] = pop_mffimport;' catchstrs.store_and_hist ];
    comwrite   = [  trystrs.no_check 'LASTCOM = pop_mffexport(EEG);' catchstrs.store_and_hist ];

    % create menus (CHANGING THESE MENUS AFFECTS THE MAIN eeglab.m FUNCTION)
    % ------------
    submenu = uimenu( menui, 'Label', 'Import EGI/Philips .mff file', 'separator', 'on', 'CallBack', comload);
    submenu = uimenu( menue, 'Label', 'Export EGI/Philips .mff file', 'separator', 'on', 'CallBack', comwrite);


text:
function p = predict(theta, X)
%PREDICT Predict whether the label is 0 or 1 using learned logistic
%regression parameters theta
%   p = PREDICT(theta, X) computes the predictions for X using a
%   threshold at 0.5 (i.e., if sigmoid(theta'*x) >= 0.5, predict 1)

m = size(X, 1); % Number of training examples

% You need to return the following variables correctly
p = zeros(m, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters.
%               You should set p to a vector of 0's and 1's
%

f=@(n) round(n);

s = sigmoid(X * theta);

p = f(s);

% =========================================================================

end

author: vkosuri

repo: CourseraMachineLearning


text:
function [C,k] = str2cell_startup(str, delimiters, options)

% Option tells weather to keep leading whitespaces.
% (Trailing whitespaces are always removed)
if ~exist('options','var')
    options = '';
end

if ~strcmpi(options, 'keepblanks')
    str = strtrim(str);
end
str = deblank(str);

if ~exist('delimiters','var') || isempty(delimiters)
    delimiters{1} = sprintf('\n');
elseif ~iscell(delimiters)
    foo{1} = delimiters;
    delimiters = foo;
end

% Get indices of all the delimiters
k=[];
for kk=1:length(delimiters)
    k = [k, find(str==delimiters{kk})];
end
j = find(~ismember(1:length(str),k));

% The following line seems to hurt performance a little bit. It was
% meant to preallocate to speed things up but it does not seem to do that.
% C = repmat({blanks(max(diff([k,length(str)])))}, length(k)+1, 1);
C = {};
ii=1; kk=1;
while ii<=length(j)
    C{kk} = str(j(ii));
    ii=ii+1;
    jj=2;
    while (ii<=length(j)) && ((j(ii)-j(ii-1))==1)
        C{kk}(jj) = str(j(ii));
        jj=jj+1;
        ii=ii+1;
    end
    C{kk}(jj:end)='';
    kk=kk+1;
end
C(kk:end) = [];

author: BUNPC


text:
classdef MimMouseCursorStatus < handle
    % MimMouseCursorStatus. Data representing the voxel under the mouse cursor
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        GlobalCoordX
        GlobalCoordY
        GlobalCoordZ

        ImageExists

        ImageValue
        RescaledValue
        RescaleUnits
        OverlayValue
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function idx = scimat_world2index(x, scimat, CHOP)
% SCIMAT_WORLD2INDEX  Convert real world coordinates to image indices for
% the SCIMAT image struct that we use in Gerardus.
%
%   Function SCIMAT_WORLD2INDEX() converts the coordinates of a voxel given
%   as real world coordinates [x, y, z, t] into index coordinates
%   [row, column, slice, frame].
%
%      [x, y, z, t] -> [r, c, s, f]
%
%   This agrees with Matlab's convention that images are expected to be
%   (r, c, s) <-> (y, x, z), but point coordinates are given in the
%   (x, y, z)-order.
%
%   This function can also be applied to images that are not 4D, and in
%   that case, index and real world coordinates will have the same number
%   of elements as dimensions the image has.
%
%   The relation between indices IDX and real world coordinates X is
%
%     X = s.*(IDX-1)*R + t
%
%   where s is the voxel size, R the rotation matrix, and t the
%   image offset.
%
%   For points that are not within the data volume, the returned
%   indices are "NaN".
%
%   Note also that the indices are not rounded, to allow for sub-pixel
%   accuracy. If integer indices are required, then just use round(idx).
%
% IDX = SCIMAT_WORLD2INDEX(X, SCIMAT)
%
%   X is a 3-column matrix where each row contains the real world
%   (x,y,z)-coordinates of a point.
%
%   IDX has the same size as X, and the voxel indices in
%   (row, column, slice)-order, that corresponds to (y, x, z)-order.
%
%   SCIMAT is a struct with the image space metadata, i.e. spacing, offset
%   and orientation (see "help scimat" for details). SCIMAT.data (the fild
%   that contains the image itself) is not used by the function, and thus
%   can be present or absent. Note that Matlab will pass SCIMAT.data by
%   reference, so passing the whole image does not require more memory or
%   slow the function down.
%
% IDX = SCIMAT_WORLD2INDEX(..., CHOP)
%
%   CHOP is a flag to convert points outside the image volume to NaNs. By
%   default, CHOP=true.
%
%
% Example:
%
% >> idx = scimat_world2index([.01, .011, .02], scimat)
%
% idx =
%
%     55   189   780
%
% See also: scimat, scimat_index2world, scimat_load, scimat_im2scimat.

% Authors: Ramon Casero <rcasero@gmail.com>,
% Benjamin Villard <b.016434@gmail.com>,
% Christopher Kelly  <christopher.kelly28@googlemail.com>
% Copyright © 2009-2015 University of Oxford
% Version: 0.5.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
narginchk(2, 3);
nargoutchk(0, 1);

% defaults
if (nargin < 3 || isempty(CHOP))
    CHOP = true;
end
if (~isfield(scimat, 'rotmat'))
    scimat.rotmat = [];
end

% extract parameters
xmin = [scimat.axis.min];
dx = [scimat.axis.spacing];
n = [scimat.axis.size];
orig = xmin + dx/2;
R = scimat.rotmat;

% number of dimensions
D = length(scimat.axis);

%% convert real world coordinates to indices

% remove offset
idx = x - repmat(orig([2 1 3:end]), size(x, 1), 1);

% apply inverse rotation only to the spatial coordinates
if (~isempty(R))
    idx(:, 1:size(R, 1)) = idx(:, 1:size(R, 1)) * R';
end

% (x, y) => (y, x)
idx = idx(:, [2 1 3:end]);

% i = x / dx + 1
idx = idx ./ repmat(dx, size(idx, 1), 1) + 1;

% find which coordinates are outside the volume
if CHOP
    for I = 1:D
        idx(idx(:, I) < 0.5 | idx(:, I) > n(I)+0.5, I) = NaN;
    end
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function vw = UpdateMapWindow(vw)
%
% vw = UpdateMapWindow
%
% Sets the MapWindow sliders to correspond to current parameter map (if any
% is present).
%
% Ress, 6/03
% ras, 11/04 -- doesn't always throw out the prev. set points
% jw,  2/10  -- avoid error if map range is not defined (e.g., due to all map
%                               values being NaNs)
if isempty(vw.map), return, end

parMap = viewGet(vw,'map');

if isfield(vw, 'ui') && isfield(vw.ui,'windowHandle'); % test for non-hidden vw
    % Find min/max and reset sliders
	minVal = inf;
	maxVal = -inf;
    for iScan = 1:viewGet(vw, 'numScans')
        if length(parMap) >= iScan && ~isempty(parMap{iScan})
            % the 'mrvMinmax' function ignores NaNs and Infs, this helps
            % prevent slider-setting errors
            rng = mrvMinmax( parMap{iScan}(:) );

            % deal with NaNs to avoid error when setting map window limits
            if isempty(rng), rng = [0 1]; end

            minVal = min(minVal, rng(1));
            maxVal = max(maxVal, rng(2));
        end
    end

    % deal w/ min/max being switched
    if (minVal>maxVal), tmp=maxVal; maxVal=minVal; minVal=tmp; end

    % if new min/max are outside current set pts, keep same set pts;
    % otherwise, reset points to minVal & maxVal
    setVals = viewGet(vw, 'mapWin');
    setMin = setVals(1); setMax = setVals(2);
%    if setMin < minVal
        setMin = minVal;
%    end
%    if setMax > maxVal
        setMax = maxVal;
%    end

	vw = resetSlider(vw, vw.ui.mapWinMin, minVal, maxVal, setMin);
	vw = resetSlider(vw, vw.ui.mapWinMax, minVal, maxVal, setMax);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% Test file for chebtech1/coeffs2vals.m

function pass = test_coeffs2vals(varargin)

% Set a tolerance (pref.chebfuneps doesn't matter)
tol = 100*eps;

%%
% Test that a single coefficient is converted correctly
c = sqrt(2);
v = chebtech1.coeffs2vals(c);
pass(1) = (c == v);

%%
% Simple data (even case)
c = (6:-1:1).';
% Exact values
vTrue = [ -3*sqrt(6)/2-5/sqrt(2)+2*sqrt(3)+7 ; 4 - sqrt(2)/2 ; -3*sqrt(6)/2+5/sqrt(2)-2*sqrt(3)+7 ; 3*sqrt(6)/2-5/sqrt(2)-2*sqrt(3)+7 ; 4 + sqrt(2)/2 ; 3*sqrt(6)/2+5/sqrt(2)+2*sqrt(3)+7];

%%
% Test real branch
v = chebtech1.coeffs2vals(c);
pass(2) = norm(v - vTrue, inf) < tol;
pass(3) = ~any(imag(v));

%%
% Test imaginary branch
v = chebtech1.coeffs2vals(1i*c);
pass(4) = norm(v - 1i*vTrue, inf) < tol;
pass(5) = ~any(real(v));

%%
% Test general branch
v = chebtech1.coeffs2vals((1+1i)*c);
pass(6) = norm(v - (1+1i)*vTrue, inf) < tol;

%%
% Test for array input
v = chebtech1.coeffs2vals([c, -c]);
pass(7) = norm(v(:,1) - vTrue, inf) < tol && ...
          norm(v(:,2) + vTrue, inf) < tol;

%%
% Simple data (odd case)
c = (5:-1:1).';
% Exact values
vTrue = [ 11/2+sqrt(5)-2*sqrt((5+sqrt(5))/2)-sqrt((5-sqrt(5))/2) ; 11/2-sqrt(5)-2*sqrt((5-sqrt(5))/2)+sqrt((5+sqrt(5))/2) ; 3 ; 11/2-sqrt(5)+2*sqrt((5-sqrt(5))/2)-sqrt((5+sqrt(5))/2) ; 11/2+sqrt(5)+2*sqrt((5+sqrt(5))/2)+sqrt((5-sqrt(5))/2) ];

%%
% Test real branch
v = chebtech1.coeffs2vals(c);
pass(8) = norm(v - vTrue, inf) < tol;
pass(9) = ~any(imag(v));

%%
% Test imaginary branch
v = chebtech1.coeffs2vals(1i*c);
pass(10) = norm(v - 1i*vTrue, inf) < tol;
pass(11) = ~any(real(v));

%%
% Test general branch
v = chebtech1.coeffs2vals((1+1i)*c);
pass(12) = norm(v - (1+1i)*vTrue, inf) < tol;

%%
% Test for array input
v = chebtech1.coeffs2vals([c, -c]);
pass(13) = norm(v(:,1) - vTrue, inf) < tol && ...
          norm(v(:,2) + vTrue, inf) < tol;

%%
% Test for symmetry preservation
c = kron(ones(10,1),eye(2));
v = chebtech1.coeffs2vals(c);
pass(14) = norm(v(:,1) - flipud(v(:,1)), inf) == 0 && ...
          norm(v(:,2) + flipud(v(:,2)), inf) == 0;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function gui_edit_channelflag( DataFile )
% GUI_EDIT_CHANNELFLAG: Edit the ChannelFlag array of a MEG/EEG recordings file, with the associated Channel file.
%
% USAGE: gui_edit_channelflag( DataFile )

% @=============================================================================
% 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, 2008-2018

global GlobalData;
ChannelFile = [];

%% ===== LOAD CHANNEL/DATA FILE =====
% Progress bar
bst_progress('start', 'Channel editor', 'Initialization...');
% Get the DataFile in the Brainstorm database
[sStudy, iStudy] = bst_get('DataFile', DataFile);
% Try to find file in stat files
if isempty(sStudy)
    [sStudy, iStudy] = bst_get('StatFile', DataFile);
end
% If DataFile in database
if ~isempty(sStudy)
    % Get channel file
    Channel = bst_get('ChannelForStudy', iStudy);
    if ~isempty(Channel)
        ChannelFile = Channel.FileName;
    end
% DataFile not in database: try absolute reference
elseif file_exist(DataFile)
    % Try to find a channel file in the DataFile directory
    dirFiles = dir(bst_fullfile(bst_fileparts(DataFile), 'channel*.mat'));
    if ~isempty(dirFiles)
        ChannelFile = bst_fullfile(bst_fileparts(DataFile), dirFiles(1).name);
    end
end
% No ChannelFile
if isempty(ChannelFile)
    bst_error('Channel file not found.', 'Edit ChannelFlag');
    return
end

% ===== START CHANNEL EDITOR =====
% Check if a "Channel Editor" panel already exists
panel = bst_get('Panel', 'ChannelEditor');
% If panel exists
if ~isempty(panel)
    % Close it
    gui_hide(panel);
end
% Create new panel
newpanel = panel_channel_editor('CreatePanel', ChannelFile, DataFile);
if isempty(newpanel)
    bst_progress('stop');
    return
end
% Show panel
bstContainer = gui_show(newpanel, 'JavaWindow', ['Edit ChannelFlag: ', DataFile], 0, 0, 1);
if isempty(bstContainer)
    return;
end

% Get current layout
[jBstArea, FigArea, nScreens, jFigArea, jInsets] = gui_layout('GetScreenBrainstormAreas');
% Place this figure on the screen to take all the possible figure space
jFrame = bstContainer.handle{1};
jFrame.setLocation(jFigArea.getX(), jFigArea.getY() - jInsets.bottom);
jFrame.setSize(jFigArea.getWidth(), jFigArea.getHeight());

% Set selected sensors
if ~isempty(GlobalData.DataViewer.SelectedRows)
    panel_channel_editor('SetChannelSelection', GlobalData.DataViewer.SelectedRows, 0);
end
bst_progress('stop');

author: brainstorm-tools


text:
function err = getHcurlerror3NE2(node,elem,curlE,Eh,markedElem)
%% GETHCURLERROR3NE2 Hcurl norm of approximation error for the quadratic (1st type) Nedelect element in 3-D.
%
% err = getHcurlerror3NE1(node,elem,curlE,Eh,markedElem);
%
% NOTE: it is identical to getHcurlerror3NE since the added basis has no
% contribution to the curl part.
%
% Example
%
%     node = [-1,-1,-1; 1,-1,-1; 1,1,-1; -1,1,-1; -1,-1,1; 1,-1,1; 1,1,1; -1,1,1];
%     elem = [1,2,3,7; 1,6,2,7; 1,5,6,7; 1,8,5,7; 1,4,8,7; 1,3,4,7];
%     maxIt = 4;
%     HcurlErr = zeros(maxIt,1);
%     N = zeros(maxIt,1);
%     for k = 1:maxIt
%         [node,elem] = uniformbisect3(node,elem);
%         [elem2dof,dofSign,edge] = dof3edge(elem);
%         pde = Maxwelldata2;
%         uI = edgeinterpolate1(pde.exactu,node,edge);
%         HcurlErr(k) = getHcurlerror3NE1(node,elem,pde.curlu,uI);
%         N(k) = length(uI);
%     end
%     r = showrate(N,HcurlErr,1,'b-+');
%     legend('||u-u_I||_{curl}',['N^{' num2str(r) '}'],'LOCATION','Best');
%
% See also getHcurlerror3NE1, getHcurlerror3NE2, getL2error3NE
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Sort elem to ascend ordering
elem = sort(elem,2);

%% Construct Data Structure
[Dlambda,volume] = gradbasis3(node,elem);
% elem2dof
[elem2edge,edge] = dof3edge(elem);
% elem2face = dof3RT0(elem);
[elem2face,face] = dof3face(elem);
NT = size(elem,1);  NE = size(edge,1); NF = size(face,1);
elem2dof = [elem2edge elem2edge+NE elem2face+2*NE elem2face+2*NE+NF];
% local indices
locBasesIdx = [1 2 0; 1 3 0; 1 4 0; 2 3 0; 2 4 0; 3 4 0; ... % phi
               1 2 0; 1 3 0; 1 4 0; 2 3 0; 2 4 0; 3 4 0; ... % psi
               3 2 4; 3 1 4; 2 1 4; 2 1 3; ...
               4 2 3; 4 1 3; 4 1 2; 3 1 2]; % chi

%% compute H1 error element-wise using quadrature rule with order quadOrder
err = zeros(NT,1);
[lambda,w] = quadpts3(4);
nQuad = size(lambda,1);
for p = 1:nQuad
    % quadrature points in the x-y-z coordinate
    pxy = lambda(p,1)*node(elem(:,1),:) ...
        + lambda(p,2)*node(elem(:,2),:) ...
        + lambda(p,3)*node(elem(:,3),:) ...
        + lambda(p,4)*node(elem(:,4),:);
    if isnumeric(curlE) % a constant vector
        curlEp = repmat(curlE,NT,1);
    else % function handel
        curlEp = curlE(pxy);
    end
    curlEhp = zeros(NT,3);
    % compute Ehp at quadrature points
    for k = 1:20
        k1 = locBasesIdx(k,1);
        k2 = locBasesIdx(k,2);
        k3 = locBasesIdx(k,3);
        if k<=6
            % curl phi = 2*Dlambda_i cross Dlambda_j;
            curlBasis_k = 2*cross(Dlambda(:,:,k1),Dlambda(:,:,k2),2);
        elseif k<=12
            % curl psi = 0;
            curlBasis_k = 0;
        else % chi = lambda_{i1}phi_{i2,i3}
            % curl chi =  Dlambda_{i1}cross phi_{i2,i3} + lambda_{i1}curl phi_{i2,i3}
            curlBasis_k = cross(Dlambda(:,:,k1),lambda(p,k2)*Dlambda(:,:,k3) ...
                                   -lambda(p,k3)*Dlambda(:,:,k2),2) ...
                        + lambda(p,k1)*2*cross(Dlambda(:,:,k2),Dlambda(:,:,k3),2);
        end
        curlEhp = curlEhp + repmat(Eh(elem2dof(:,k)),1,3).*curlBasis_k;
    end
    err = err + w(p)*volume.*sum((curlEp - curlEhp).^2,2);
end
% modify the error
err(isnan(err)) = 0; % remove the singular part
if (nargin == 5) && ~isempty(markedElem)
    err = err(markedElem); % L2 error on some marked region
end
err = sqrt(sum(err));
%% TODO write more M-lint

repo: Aero-matlab


text:
function [ n_data, mu, sigma, b, x, fx ] = ...
  truncated_normal_b_pdf_values ( n_data )

%*****************************************************************************80
%
%% TRUNCATED_NORMAL_B_PDF_VALUES: values of the Truncated Normal B PDF.
%
%  Discussion:
%
%    The Normal distribution, with mean Mu and standard deviation Sigma,
%    is truncated to the interval (-oo,B].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real MU, the mean of the distribution.
%
%    Output, real SIGMA, the standard deviation of the distribution.
%
%    Output, real B, the upper truncation limit.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 11;

  b_vec = [ ...
     150.0, ...
     150.0, ...
     150.0, ...
     150.0, ...
     150.0, ...
     150.0, ...
     150.0, ...
     150.0, ...
     150.0, ...
     150.0, ...
     150.0 ];

  fx_vec = [ ...
    0.01507373507401876, ...
    0.01551417047139894, ...
    0.01586560931024694, ...
    0.01612150073158793, ...
    0.01627701240029317, ...
    0.01632918226724295, ...
    0.01627701240029317, ...
    0.01612150073158793, ...
    0.01586560931024694, ...
    0.01551417047139894, ...
    0.01507373507401876 ];

  mu_vec = [ ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0 ];

  sigma_vec = [ ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0  ];

  x_vec = [ ...
     90.0, ...
     92.0, ...
     94.0, ...
     96.0, ...
     98.0, ...
    100.0, ...
    102.0, ...
    104.0, ...
    106.0, ...
    108.0, ...
    110.0 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    b = 0.0;
    mu = 0.0;
    sigma = 0.0;
    x = 0.0;
    fx = 0.0;
  else
    b = b_vec(n_data);
    mu = mu_vec(n_data);
    sigma = sigma_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function cnk = combin ( n, k )

%*****************************************************************************80
%
%% COMBIN computes the combinatorial coefficient C(N,K).
%
%  Discussion:
%
%    Real arithmetic is used, and C(N,K) is computed directly, via
%    Gamma functions, rather than recursively.
%
%    C(N,K) is the number of distinct combinations of K objects
%    chosen from a set of N distinct objects.  A combination is
%    like a set, in that order does not matter.
%
%    C(N,K) = N% / ( (N-K)% * K% )
%
%  Example:
%
%    The number of combinations of 2 things chosen from 5 is 10.
%
%    C(5,2) = ( 5 * 4 * 3 * 2 * 1 ) / ( ( 3 * 2 * 1 ) * ( 2 * 1 ) ) = 10.
%
%    The actual combinations may be represented as:
%
%      (1,2), (1,3), (1,4), (1,5), (2,3),
%      (2,4), (2,5), (3,4), (3,5), (4,5).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the value of N.
%
%    Input, integer K, the value of K.
%
%    Output, real CNK, the value of C(N,K)
%
  if ( n < 0 )

    cnk = 0.0;

  elseif ( k == 0 )

    cnk = 1.0;

  elseif ( k == 1 )

    cnk = n;

  elseif ( 1 < k && k < n - 1 )

    facn = gammaln ( n + 1 );
    fack = gammaln ( k + 1 );
    facnmk = gammaln ( n - k + 1 );

    cnk = round ( exp ( facn - fack - facnmk ) );

  elseif ( k == n - 1 )

    cnk = n;

  elseif ( k == n )

    cnk = 1.0;

  else

    cnk = 0.0;

  end

  return
end