text:
% BOOTSTAT - accumulate surrogate data to assess significance by permutation of some
%              measure of two input variables.
%
%              If 'distfit','on', fits the psd with a 4th-order polynomial using the
%              data kurtosis, as in Ramberg, J.S., Tadikamalla, P.R., Dudewicz E.J.,
%              Mykkytka, E.F. "A probability distribution and its uses in fitting data."
%              Technometrics, 21:201-214, 1979.
% Usage:
%            >> [rsignif,rboot] = bootstat( { arg1 arg2 ...}, formula, varargin ...);
% Inputs:
%    arg1    - [array] 1-D, 2-D or 3-D array of values
%    arg2    - [array] 1-D, 2-D or 3-D array of values
%    formula - [string] formula to compute the given measure. Takes arguments
%                   'arg1', 'arg2' as inputs and 'res' (result, by default) as output.
%                   For data arrays of more than 1 dimension, the formula must be iterative
%                   so that shuffling can occur at each step while scanning the last
%                   array dimension.  Examples:
%                   'res = arg1 - arg2'              % difference of two 1-D data arrays
%                   'res = mean( arg1 .* arg2)'      % mean projection of two 1-D data arrays
%                   'res = res + arg1 .* conj(arg2)' % iterative, for use with 2|3-D arrays
% Optional inputs:
%   'boottype '   - ['rand'|'shuffle']
%                   'rand'  = do not shuffle data. Only flip polarity randomly (for real
%                             number) or phase (for complex numbers).
%                   'shuffle' = shuffle values of first argument (see two options below).
%                   Default.
%   'shuffledim'  - [integer] indices of dimensions to shuffle. For instance, [1 2] will
%                   shuffle the first two dimensions. Default is to shuffle along
%                   dimension 2.
%   'shufflemode' - ['swap'|'regular'] shuffle mode. Either swap dimensions (for instance
%                   swap rows then columns if dimension [1 2] are selected) or shuffle
%                   in each dimension independently (slower). If only one dimension is
%                   selected for shuffling, this option does not change the result.
%   'randmode'    - ['opposite'|'inverse'] randomize sign (or phase for complex number,
%                   or randomly set half the value to reference.
%   'alpha'       - [real] significance level (between 0 and 1) {default 0.05}.
%   'naccu'       - [integer] number of exemplars to accumulate {default 200}.
%   'bootside'    - ['both'|'upper'] side of the surrogate distribution to
%                   consider for significance. This parameter affects the size
%                   of the last dimension of the accumulation array ('accres')
%                   (size is 2 for 'both' and 1 for 'upper') {default: 'both'}.
%   'basevect'    - [integer vector] time vector indices for baseline in second dimension.
%                   {default: all time points}.
%   'rboot'       - accumulation array (from a previous call). Allows faster
%                   computation of the 'rsignif' output {default: none}.
%   'formulaout'  - [string] name of the computed variable {default: 'res'}.
%   'dimaccu'     - [integer] use dimension in result to accumulate data.
%                   For instance if the result array is size [60x50] and this value is 2,
%                   the function will consider than 50 times 60 value have been accumulated.
%
% Fitting distribution:
%   'distfit'     - ['on'|'off'] fit distribution with known function to compute more accurate
%                   limits or exact p-value (see 'vals' option). The MATLAB statistical toolbox
%                   is required. This option is currently implemented only for 1-D data.
%   'vals'        - [float array] significance values. 'alpha' is ignored and
%                   rsignif returns the p-values. Requires 'distfit' (see above).
%                   This option currently implemented only for 1-D data.
%   'correctp'    - [phat pci zerofreq] parameters for correcting for a biased probability
%                   distribution (requires 'distfit' above). See help of CORRECTFIT.
% Outputs:
%    rsignif      - significance arrays. 2 values (low high) for each point (use
%                   'alpha' to change these limits).
%    rboot        - accumulated surrogate data values.
%
% Authors: Arnaud Delorme, Bhaktivedcanta Institute, Mumbai, India, Nov 2004
%
% See also: TIMEF

% NOTE: There is an undocumented parameter, 'savecoher', [0|1]
% HELP TEXT REMOVED:  (Ex: Using option 'both', coherence during baseline would be
%                      ignored since times are shuffled during each accumulation.

% Copyright (C) 9/2002  Arnaud Delorme & Scott Makeig, SCCN/INC/UCSD
%
% 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.

% *************************************
% To fit the psd with as 4th order polynomial using the distribution kurtosis,
% Reference: Ramberg, J.S., Tadikamalla, P.R., Dudewicz E.J., Mykkytka, E.F.
% "A probability distribution and its uses in fitting data."
% Technimetrics, 1979, 21: 201-214.
% *************************************

function [accarrayout, Rbootout, Rbootout2] = bootstat(oriargs, formula, varargin)
%	nb_points, timesout, naccu, baselength, baseboot, boottype, alpha, rboot);

if nargin < 2
	help bootstat;
	return;
end

if ~ischar(formula)
	error('The second argument must be a string formula');
end

g = finputcheck(varargin, ...
                { 'dims'          'integer'  []                       []; ...
                  'naccu'         'integer'  [0 10000]                200; ...
                  'bootside'      'string'   { 'both','upper' }       'both'; ...
                  'basevect'      'integer'  []                       []; ...
				  'boottype'      'string'  { 'rand','shuffle' }      'shuffle'; ...
				  'shufflemode'   'string'  { 'swap','regular' }      'swap'; ...
				  'randmode'      'string'  { 'opposite','inverse' }  'opposite'; ...
				  'shuffledim'    'integer'  [0 Inf]                  []; ...
				  'label'         'string'   []                       formula; ...
				  'alpha'         'real'     [0 1]                    0.05; ...
				  'vals'          'real'     []                       []; ...
				  'distfit'       'string'   {'on','off' }            'off'; ...
				  'dimaccu'       'integer'  [1 Inf]                  []; ...
				  'correctp'      'real'     []                       []; ...
				  'rboot'         'real'     []                       NaN	});
if ischar(g)
	error(g);
end
if isempty(g.shuffledim) && strcmpi(g.boottype, 'rand')
    g.shuffledim = [];
elseif isempty(g.shuffledim)
    g.shuffledim = 2;
end;
unitname = '';
if 2/g.alpha > g.naccu
    if strcmpi(g.distfit, 'off') || ~((size(oriarg1,1) == 1 || size(oriarg1,2) == 1) && size(oriarg1,3) == 1)
        g.naccu = 2/g.alpha;
        fprintf('Adjusting naccu to compute alpha value');
    end
end
if isempty(g.rboot)
	g.rboot = NaN;
end

% function for bootstrap computation
% ----------------------------------
if ~iscell(oriargs) || length(oriargs) == 1,
    oriarg1 = oriargs;
    oriarg2 = [];
else
    oriarg1 = oriargs{1};
    oriarg2 = oriargs{2};
end
[nb_points times trials] = size(oriarg1);
if times == 1, disp('Warning 1 value only for shuffling dimension'); end

% only consider baseline
% ----------------------
if ~isempty(g.basevect)
    fprintf('\nPermutation statistics baseline length is %d (out of %d) points\n', length(g.basevect), times);
    arg1 = oriarg1(:,g.basevect,:);
    if ~isempty(oriarg2)
        arg2 = oriarg2(:,g.basevect,:);
    end
else
    arg1 = oriarg1;
    arg2 = oriarg2;
end

% formula for accumulation array
% ------------------------------
% if g.dimaccu is not empty, accumulate over that dimension
% of the resulting array to speed up computation
formula = [ 'res=' formula ];
g.formulapost = [ 'if index == 1, ' ...
                  '   if ~isempty(g.dimaccu), ' ...
                  '      Rbootout= zeros([ ceil(g.naccu/size(res,g.dimaccu)) size( res ) ]);' ...
                  '   else,' ...
                  '      Rbootout= zeros([ g.naccu size( res ) ]);' ...
                  '   end;' ...
                  'end,' ...
                  'Rbootout(count,:,:,:) = res;' ...
                  'count = count+1;' ...
                  'if ~isempty(g.dimaccu), ' ...
                  '   index = index + size(res,g.dimaccu);' ...
                  '   fprintf(''%d '', index-1);' ...
                  'else ' ...
                  '   index=index+1;' ...
                  '   if rem(index,10)  == 0, fprintf(''%d '', index); end;' ...
                  '   if rem(index,100) == 0, fprintf(''\n''); end;' ...
                  'end;' ];

% **************************
% case 1: precomputed values
% **************************
if ~isnan(g.rboot)
    Rbootout = g.rboot;
% ***********************************
% case 2: randomize polarity or phase
% ***********************************
elseif strcmpi(g.boottype, 'rand') && strcmpi(g.randmode, 'inverse')
    fprintf('Bootstat function: randomize inverse values\n');
    fprintf('Processing permutation statistics for %s (naccu=%d):', g.label, g.naccu);

    % compute random array
    % --------------------
    multarray = ones(size(arg1));
    totlen    = prod(size(arg1));
    if isreal(arg1),
        multarray(1:round(totlen/2)) = 0;
    end
    for shuff = 1:ndims(multarray)
         multarray = supershuffle(multarray,shuff); % initial shuffling
    end
    if isempty(g.shuffledim), g.shuffledim = 1:ndims(multarray); end
    invarg1 = 1./arg1;

    % accumulate
    % ----------
    index = 1;
    count = 1;
    while index <= g.naccu
        for shuff = g.shuffledim
            multarray = supershuffle(multarray,shuff);
        end
        tmpinds = find(reshape(multarray, 1, prod(size(multarray))));
        arg1 = oriarg1;
        arg1(tmpinds) = invarg1(tmpinds);
        eval([ formula ';' ]);
        eval( g.formulapost ); % also contains index = index+1
    end
elseif strcmpi(g.boottype, 'rand') % opposite
    fprintf('Bootstat function: randomize polarity or phase\n');
    fprintf('Processing permutation statistics for %s (naccu=%d):', g.label, g.naccu);

    % compute random array
    % --------------------
    multarray = ones(size(arg1));
    totlen    = prod(size(arg1));
    if isreal(arg1),
        multarray(1:round(totlen/2)) = -1;
    else
        tmparray            = exp(j*linspace(0,2*pi,totlen+1));
        multarray(1:totlen) = tmparray(1:end-1);
    end
    for shuff = 1:ndims(multarray)
         multarray = supershuffle(multarray,shuff); % initial shuffling
    end
    if isempty(g.shuffledim), g.shuffledim = 1:ndims(multarray); end

    % accumulate
    % ----------
    index = 1;
    count = 1;
    while index <= g.naccu
        for shuff = g.shuffledim
            multarray = supershuffle(multarray,shuff);
        end
        arg1 = arg1.*multarray;
        eval([ formula ';' ]);
        eval( g.formulapost ); % also contains index = index+1
    end
% ********************************************
% case 3: shuffle vector of only one dimension
% ********************************************
elseif length(g.shuffledim) == 1
    fprintf('Bootstat function: shuffling along dimension %d only\n', g.shuffledim);
    fprintf('Processing permutation statistics for %s (naccu=%d):', g.label, g.naccu);

    index = 1;
    count = 1;
    while index <= g.naccu
        arg1 = shuffleonedim(arg1,g.shuffledim);
        eval([ formula ';' ]);
        eval( g.formulapost );
    end
% ***********************************************
% case 5: shuffle vector along several dimensions
% ***********************************************
else
    if strcmpi(g.shufflemode, 'swap') % swap mode
        fprintf('Bootstat function: shuffling along dimension %s (swap mode)\n', int2str(g.shuffledim));
        fprintf('Processing permutation statistics for %s (naccu=%d):', g.label, g.naccu);
        index = 1;
        count = 1;
        while index <= g.naccu
            for shuff = g.shuffledim
                arg1 = supershuffle(arg1,shuff);
            end
            eval([ formula ';' ]);
            eval( g.formulapost );
        end
    else  % regular shuffling
        fprintf('Bootstat function: shuffling along dimension %s (regular mode)\n', int2str(g.shuffledim));
        fprintf('Processing permutation statistics for %s (naccu=%d):', g.label, g.naccu);
        index = 1;
        count = 1;
        while index <= g.naccu
            for shuff = g.shuffledim
                arg1 = shuffleonedim(arg1,shuff);
            end
            eval([ formula ';' ]);
            eval( g.formulapost );
        end
    end
end
Rbootout(count:end,:,:,:) = [];

% **********************
% assessing significance
% **********************

% get accumulation array
% ----------------------
accarray = Rbootout;
if ~isreal(accarray)
    accarray = sqrt(accarray .* conj(accarray)); % faster than abs()
end
% reshape the output if necessary
% -------------------------------
if ~isempty(g.dimaccu)
    if g.dimaccu+1 == 3
        accarray = permute( accarray, [1 3 2]);
    end
    accarray = reshape( accarray, size(accarray,1)*size(accarray,2), size(accarray,3) );
end
if size(accarray,1) == 1, accarray = accarray'; end; % first dim contains g.naccu

% ******************************************************
% compute thresholds on array not fitting a distribution
% ******************************************************
if strcmpi(g.distfit, 'off')

    % compute bootstrap significance level
    % ------------------------------------
    accarray  = sort(accarray,1); % always sort on naccu
    Rbootout2 = accarray;
    i         = round(size(accarray,1)*g.alpha);
    accarray1 = squeeze(mean(accarray(size(accarray,1)-i+1:end,:,:),1));
    accarray2 = squeeze(mean(accarray(1:i                     ,:,:),1));
    if abs(accarray(1,1,1) - accarray(end,1,1)) < abs(accarray(1,1,1))*1e-15
        accarray1(:) = NaN;
        accarray2(:) = NaN;
    end

else
    % *******************
    % fit to distribution
    % *******************
    sizerboot   = size (accarray);
    accarray1   = zeros(sizerboot(2:end));
    accarray2   = zeros(sizerboot(2:end));

    if ~isempty(g.vals{index})
        if ~all(size(g.vals{index}) == sizerboot(2:end) )
            error('For fitting, vals must have the same dimension as the output array (try transposing)');
        end
    end

    % fitting with Ramberg-Schmeiser distribution
    % -------------------------------------------
    if ~isempty(g.vals{index}) % compute significance for value
        for index1 = 1:size(accarrayout,1)
            for index2 = 1:size(accarrayout,2)
                accarray1(index1,index2) = 1 - rsfit(squeeze(accarray(:,index1,index2)), g.vals{index}(index1, index2));
                if length(g.correctp) == 2
                    accarray1(index1,index2) = correctfit(accarray1, 'gamparams', [g.correctp 0]); % no correction for p=0
                else
                    accarray1(index1,index2) = correctfit(accarray1, 'gamparams', g.correctp);
                end
            end
        end
    else % compute value for significance
        for index1 = 1:size(accarrayout,1)
            for index2 = 1:size(accarrayout,2)
                [p c l chi2] = rsfit(Rbootout(:),0);
                pval = g.alpha;   accarray1(index1,index2) = l(1) + (pval.^l(3) - (1-pval).^l(4))/l(2);
                pval = 1-g.alpha; accarray2(index1,index2) = l(1) + (pval.^l(3) - (1-pval).^l(4))/l(2);
            end
        end
    end

    % plot results
    % -------------------------------------
    % figure;
    % hist(abs(Rbootout)); tmpax = axis;
    % hold on;
    % valcomp = linspace(min(abs(Rbootout(:))), max(abs(Rbootout(:))), 100);
    % normy = normpdf(valcomp, mu, sigma);
    % plot(valcomp, normy/max(normy)*tmpax(4), 'r');
    % return;
end

% set output array: backward compatible
% -------------------------------------
if strcmpi(g.bootside, 'upper'); % only upper significance
    accarrayout = accarray1;
else
    if size(accarray1,1) ~= 1 && size(accarray1,2) ~= 1
        accarrayout        = accarray2;
        accarrayout(:,:,2) = accarray1;
    else
        accarrayout = [ accarray2(:) accarray1(:) ];
    end
end
accarrayout = squeeze(accarrayout);
if size(accarrayout,1) == 1 && size(accarrayout,3) == 1, accarrayout = accarrayout'; end

% better but not backward compatible
% ----------------------------------
% accarrayout = { accarray1 accarray2 };

return;

    % fitting with normal distribution (deprecated)
    % --------------------------------
    [mu sigma] = normfit(abs(Rbootout(:)));
    accarrayout = 1 - normcdf(g.vals, mu, sigma); % cumulative density distribution
                                        % formula of normal distribution
                                        % y = 1/sqrt(2) * exp( -(x-mu).^2/(sigma*sigma*2) ) / (sqrt(pi)*sigma);
% % Gamma and Beta fits:
% elseif strcmpi(g.distfit, 'gamma')
%  [phatgam pcigam] = gamfit(abs(Rbootout(:)));
%  gamy = gampdf(valcomp, phatgam(1), pcigam(2))
%  p = 1 - gamcdf(g.vals, phatgam(1), pcigam(2)); % cumulative density distribution
% elseif strcmpi(g.distfit, 'beta')
%  [phatbeta pcibeta] = betafit(abs(Rbootout(:)));
%  betay = betapdf(valcomp, phatbeta(1), pcibeta(1));
%  p = 1 - betacdf(g.vals, phatbeta(1), pcibeta(1)); % cumulative density distribution
% end

    if strcmpi(g.distfit, 'off')
        tmpsort = sort(Rbootout);
        i = round(g.alpha*g.naccu);
        sigval = [mean(tmpsort(1:i)) mean(tmpsort(g.naccu-i+1:g.naccu))];
        if strcmpi(g.bootside, 'upper'), sigval = sigval(2); end
        accarrayout = sigval;
    end

% this shuffling preserve the number of -1 and 1
% for columns and rows (assuming matrix size is multiple of 2
% -----------------------------------------------------------
function array = supershuffle(array, dim)
    if size(array, 1) == 1 || size(array,2) == 1
        array = shuffle(array);
        return;
    end
    if size(array, dim) == 1, return; end

    if dim == 1
        indrows = shuffle(1:size(array,1));
        for index = 1:2:length(indrows)-rem(length(indrows),2) % shuffle rows
            tmparray                    = array(indrows(index),:,:);
            array(indrows(index),:,:)   = array(indrows(index+1),:,:);
            array(indrows(index+1),:,:) = tmparray;
        end
    elseif dim == 2
        indcols = shuffle(1:size(array,2));
        for index = 1:2:length(indcols)-rem(length(indcols),2) % shuffle columns
            tmparray                    = array(:,indcols(index),:);
            array(:,indcols(index),:)   = array(:,indcols(index+1),:);
            array(:,indcols(index+1),:) = tmparray;
        end
    else
        ind3d = shuffle(1:size(array,3));
        for index = 1:2:length(ind3d)-rem(length(ind3d),2) % shuffle columns
            tmparray                  = array(:,:,ind3d(index));
            array(:,:,ind3d(index))   = array(:,:,ind3d(index+1));
            array(:,:,ind3d(index+1)) = tmparray;
        end
    end

% shuffle one dimension, one row/columns at a time
% -----------------------------------------------
function array = shuffleonedim(array, dim)
    if size(array, 1) == 1 || size(array,2) == 1
        array = shuffle(array, dim);
    else
        if dim == 1
            for index1 = 1:size(array,3)
                for index2 = 1:size(array,2)
                    array(:,index2,index1) = shuffle(array(:,index2,index1));
                end
            end
        elseif dim == 2
            for index1 = 1:size(array,3)
                for index2 = 1:size(array,1)
                    array(index2,:,index1) = shuffle(array(index2,:,index1));
                end
            end
        else
            for index1 = 1:size(array,1)
                for index2 = 1:size(array,2)
                    array(index1,index2,:) = shuffle(array(index1,index2,:));
                end
            end
        end
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
%{
    Tests total-variation problem

    min_x ||x||_TV
s.t.
    || A(x) - b || <= eps

The solvers solve a regularized version

see also test_sBP.m and test_sBPDN_W.m

requires the image processing toolbox for this demo
(but the TFOCS solver does not rely on this toolbox)

See also the TV examples in examples/largescale

%}

% Before running this, please add the TFOCS base directory to your path

myAwgn = @(x,snr) x + ...
        10^( (10*log10(sum(abs(x(:)).^2)/length(x(:))) - snr)/20 )*randn(size(x));

% Try to load the problem from disk
fileName = fullfile('reference_solutions','tv_problem1_smoothed_noisy');
randn('state',245);
rand('state',245);
n = 32;
n1 = n;
n2 = n-1;           % testing the code with non-square signals
N = n1*n2;
M = round(N/2);
A = randn(M,N);
if exist([fileName,'.mat'],'file')
    load(fileName);
    fprintf('Loaded problem from %s\n', fileName );
else

    % Generate a new problem

    n = max(n1,n2);
    x = phantom(n);
    x = x(1:n1,1:n2);
    x_original = x;

    mat = @(x) reshape(x,n1,n2);
    vec = @(x) x(:);

    b_original = A*vec(x_original);
    snr = 40;  % SNR in dB
    b   = myAwgn(b_original,snr);
    EPS = norm(b-b_original);

    tv = linop_TV( [n1,n2], [], 'cvx' );

    mu = .005*norm( tv(x_original) ,Inf);
    x0 = zeros(n1,n2);

    % get reference via CVX
    tic
    cvx_begin
        cvx_precision best
        variable xcvx(n1,n2)
        minimize tv(xcvx) + mu/2*sum_square(vec(xcvx)-vec(x0) )
        subject to
            norm(A*vec(xcvx) - b ) <= EPS
    cvx_end
    time_IPM = toc;
    x_ref = xcvx;
    obj_ref = tv(x_ref) + mu/2*sum_square(vec(x_ref)-vec(x0) );

    save(fileName,'x_ref','b','x_original','mu',...
        'EPS','b_original','obj_ref','x0','time_IPM','snr');
    fprintf('Saved data to file %s\n', fileName);

end

imshow( [x_original, x_ref] );

[M,N]           = size(A);
[n1,n2]         = size(x_original);
norm_x_ref      = norm(x_ref,'fro');
norm_x_orig     = norm(x_original,'fro');
er_ref          = @(x) norm(vec(x)-vec(x_ref))/norm_x_ref;
er_signal       = @(x) norm(x-x_original)/norm_x_orig;
resid           = @(x) norm(A*vec(x)-b)/norm(b);  % change if b is noisy

%% Call the TFOCS solver
er              = er_ref;  % error with reference solution (from IPM)
opts = [];
opts.restart    = 1000;
opts.errFcn     = { @(f,dual,primal) er(primal), ...
                    @(f,dual,primal) obj_ref - f  };
opts.maxIts     = 1000;

W   = linop_TV( [n1,n2] );
normW           = linop_TV( [n1,n2], [], 'norm' );
opts.normW2     = normW^2;
z0  = [];   % we don't have a good guess for the dual
tic;
[ x, out, optsOut ] = solver_sBPDN_W( A, W, b, EPS, mu, vec(x0), z0, opts );
time_TFOCS = toc;

fprintf('Solution has %d nonzeros.  Error vs. IPM solution is %.2e\n',...
    nnz(x), er(x) );

% Check that we are within allowable bounds
if out.err(end,1) < 1e-4
    disp('Everything is working');
else
    error('Failed the test');
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function [ nimg ] = num_images(  )
%NUM_IMAGES Determine number of images used in a NEB calculation.
%   N = NUM_IMAGES() determines the number of images used in an NEB
%   calculation by looking for directories 00, 01 ... N+1.

    a = dir;
    images = arrayfun(@(x) ~sum(~isstrprop(x.name,'digit')) & x.isdir , a);
    nimg = sum(images)-2;

end

author: Sable

repo: mcbench-benchmarks


text:
function c = p_polynomial_coefficients ( n )

%*****************************************************************************80
%
%% P_POLYNOMIAL_COEFFICIENTS: coefficients of Legendre polynomials P(n,x).
%
%  First terms:
%
%     1
%     0     1
%    -1/2   0      3/2
%     0    -3/2    0     5/2
%     3/8   0    -30/8   0     35/8
%     0    15/8    0   -70/8    0     63/8
%    -5/16  0    105/16  0   -315/16   0    231/16
%     0   -35/16   0   315/16   0   -693/16   0    429/16
%
%     1.00000
%     0.00000  1.00000
%    -0.50000  0.00000  1.50000
%     0.00000 -1.50000  0.00000  2.5000
%     0.37500  0.00000 -3.75000  0.00000  4.37500
%     0.00000  1.87500  0.00000 -8.75000  0.00000  7.87500
%    -0.31250  0.00000  6.56250  0.00000 -19.6875  0.00000  14.4375
%     0.00000 -2.1875   0.00000  19.6875  0.00000 -43.3215  0.00000  26.8125
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Daniel Zwillinger, editor,
%    CRC Standard Mathematical Tables and Formulae,
%    30th Edition,
%    CRC Press, 1996.
%
%  Parameters:
%
%    Input, integer N, the highest order polynomial to evaluate.
%    Note that polynomials 0 through N will be evaluated.
%
%    Output, real C(1:N+1,1:N+1), the coefficients of the Legendre polynomials
%    of degree 0 through N.  Each polynomial is stored as a row.
%
  if ( n < 0 )
    c = [];
    return
  end

  c(1:n+1,1:n+1) = 0.0;

  c(1,1) = 1.0;

  if ( n <= 0 )
    return
  end

  c(2,2) = 1.0;

  for i = 2 : n
    c(i+1,1:i-1) =                (   - i + 1 ) * c(i-1,1:i-1) / ( i );
    c(i+1,2:i+1) = c(i+1,2:i+1) + ( i + i - 1 ) * c(i  ,1:i  ) / ( i );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [success] = setupNwb(varargin)
	% Checks that NWB code is present and setup correctly.
	% Biafra Ahanonu
	% started: 2021.01.24 [14:31:24]
	% inputs
		%
	% outputs
		%

	% changelog
		% 2021.02.01 [15:19:40] - Update `_external_programs` to call ciapkg.getDirExternalPrograms() to standardize call across all functions.
		% 2021.03.26 [06:27:48] - Fix for options.defaultObjDir leading to incorrect NWB folder and cores not being generated.
	% TODO
		%

	[success] = ciapkg.nwb.setupNwb('passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [headmodel, sens] = ft_prepare_vol_sens(headmodel, sens, varargin)

% FT_PREPARE_VOL_SENS does some bookkeeping to ensure that the volume conductor model
% and the sensor array are ready for subsequent forward leadfield computations and
% takes care of some pre-computations to make the calculations more efficient.
%
% Use as
%   [headmodel, sens] = ft_prepare_vol_sens(headmodel, sens, ...)
% with input arguments
%   headmodel = structure with volume conductor definition
%   sens      = structure with gradiometer or electrode definition
%
% The headmodel structure represents a volume conductor model of the head,
% its contents depend on the type of model. It is described in more detail
% in FT_DATATYPE_HEADMODEL. The sens structure represents a electrode or
% gradiometer array. It is described in more detail in FT_DATATYPE_SENS.
%
% Additional options should be specified in key-value pairs and can be
%   'channel'  = cell-array with strings (default = 'all')
%
% The detailed behavior of this function depends on whether the input
% consists of EEG or MEG and furthermoree depends on the type of volume
% conductor model:
% - in case of EEG single and concentric sphere models, the electrodes are
%   projected onto the skin surface.
% - in case of EEG boundary element models, the electrodes are projected on
%   the surface and a blilinear interpoaltion matrix from vertices to
%   electrodes is computed.
% - in case of MEG and a localspheres model, a local sphere is determined
%   for each coil in the gradiometer definition.
%  - in case of MEG with a singleshell Nolte model, the volume conduction
%    model is initialized
% In any case channel selection and reordering will be done. The channel
% order returned by this function corresponds to the order in the 'channel'
% option, or if not specified, to the order in the input sensor array.
%
% See also FT_COMPUTE_LEADFIELD, FT_READ_HEADMODEL, FT_READ_SENS

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

if iscell(headmodel) && iscell(sens)
  % this represents combined EEG, ECoG and/or MEG
  for i=1:numel(headmodel)
    [headmodel{i}, sens{i}] = ft_prepare_vol_sens(headmodel{i}, sens{i}, varargin{:});
  end
  return
end

% get the optional input arguments
% fileformat = ft_getopt(varargin, 'fileformat');
channel = ft_getopt(varargin, 'channel', sens.label);   % cell-array with channel labels, default is all

% ensure that the sensor description is up-to-date (Aug 2011)
sens = ft_datatype_sens(sens);

% this is to support volumes saved in mat-files, particularly interpolated
if ischar(headmodel)
  vpath     = fileparts(headmodel);   % remember the path to the file
  headmodel = ft_read_headmodel(headmodel); % replace the filename with the content of the file
end

% ensure that the volume conduction description is up-to-date (Jul 2012)
headmodel = ft_datatype_headmodel(headmodel);

% determine whether the input contains EEG or MEG sensors
iseeg = ft_senstype(sens, 'eeg');
ismeg = ft_senstype(sens, 'meg');

% determine the skin compartment
if ~isfield(headmodel, 'skin_surface')
  if isfield(headmodel, 'bnd')
    headmodel.skin_surface   = find_outermost_boundary(headmodel.bnd);
  elseif isfield(headmodel, 'r') && length(headmodel.r)<=4
    [dum, headmodel.skin_surface] = max(headmodel.r);
  end
end

% determine the inner_skull_surface compartment
if ~isfield(headmodel, 'inner_skull_surface')
  if isfield(headmodel, 'bnd')
    headmodel.inner_skull_surface  = find_innermost_boundary(headmodel.bnd);
  elseif isfield(headmodel, 'r') && length(headmodel.r)<=4
    [dum, headmodel.inner_skull_surface] = min(headmodel.r);
  end
end

% FT_HEADMODELTYPE to an empty struct won't work further down
if isempty(headmodel)
  headmodel = [];
end

% this makes them easier to recognise
sens.type       = ft_senstype(sens);
headmodel.type  = ft_headmodeltype(headmodel);

if isfield(headmodel, 'unit') && isfield(sens, 'unit') && ~strcmp(headmodel.unit, sens.unit)
  ft_error('inconsistency in the units of the volume conductor and the sensor array');
end

if ismeg && iseeg
  % this is something that could be implemented relatively easily
  ft_error('simultaneous EEG and MEG not yet supported');

elseif ~ismeg && ~iseeg
  ft_error('the input does not look like EEG, nor like MEG');

elseif ismeg

  % always ensure that there is a linear transfer matrix for combining the coils into gradiometers
  if ~isfield(sens, 'tra');
    Nchans = length(sens.label);
    Ncoils = size(sens.coilpos,1);
    if Nchans~=Ncoils
      ft_error('inconsistent number of channels and coils');
    end
    sens.tra = eye(Nchans, Ncoils);
  end

  if ~ft_headmodeltype(headmodel, 'localspheres')
    % select the desired channels from the gradiometer array
    [selchan, selsens] = match_str(channel, sens.label);
    % only keep the desired channels, order them according to the users specification
    try, sens.chantype = sens.chantype(selsens,:); end
    try, sens.chanunit = sens.chanunit(selsens,:); end
    try, sens.chanpos  = sens.chanpos (selsens,:); end
    try, sens.chanori  = sens.chanori (selsens,:); end
    sens.label    = sens.label(selsens);
    sens.tra      = sens.tra(selsens,:);
  else
    % for the localspheres model it is done further down
  end

  % remove the coils that do not contribute to any channel output
  selcoil      = any(sens.tra~=0,1);
  sens.coilpos = sens.coilpos(selcoil,:);
  sens.coilori = sens.coilori(selcoil,:);
  sens.tra     = sens.tra(:,selcoil);

  switch ft_headmodeltype(headmodel)
    case {'infinite' 'infinite_monopole' 'infinite_currentdipole' 'infinite_magneticdipole'}
      % nothing to do

    case 'singlesphere'
      % nothing to do

    case 'concentricspheres'
      % nothing to do

    case 'neuromag'
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % if the forward model is computed using the external Neuromag toolbox,
      % we have to add a selection of the channels so that the channels
      % in the forward model correspond with those in the data.
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      [selchan, selsens] = match_str(channel, sens.label);
      headmodel.chansel = selsens;

    case 'localspheres'
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % If the volume conduction model consists of multiple spheres then we
      % have to match the channels in the gradiometer array and the volume
      % conduction model.
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      % the initial localspheres volume conductor has a local sphere per
      % channel, whereas it should have a local sphere for each coil
      if size(headmodel.r,1)==size(sens.coilpos,1) && ~isfield(headmodel, 'label')
        % it appears that each coil already has a sphere, which suggests
        % that the volume conductor already has been prepared to match the
        % sensor array
        return
      elseif size(headmodel.r,1)==size(sens.coilpos,1) && isfield(headmodel, 'label')
        if ~isequal(headmodel.label(:), sens.label(:))
          % if only the order is different, it would be possible to reorder them
          ft_error('the coils in the volume conduction model do not correspond to the sensor array');
        else
          % the coil-specific spheres in the volume conductor should not have a label
          % because the label is already specified for the coils in the
          % sensor array
          headmodel = rmfield(headmodel, 'label');
        end
        return
      end

      % the CTF way of representing the headmodel is one-sphere-per-channel
      % whereas the FieldTrip way of doing the forward computation is one-sphere-per-coil
      Nchans   = size(sens.tra,1);
      Ncoils   = size(sens.tra,2);
      Nspheres = size(headmodel.label);

      if isfield(headmodel, 'orig')
        % these are present in a CTF *.hdm file
        singlesphere.o(1,1) = headmodel.orig.MEG_Sphere.ORIGIN_X;
        singlesphere.o(1,2) = headmodel.orig.MEG_Sphere.ORIGIN_Y;
        singlesphere.o(1,3) = headmodel.orig.MEG_Sphere.ORIGIN_Z;
        singlesphere.r      = headmodel.orig.MEG_Sphere.RADIUS;
        % ensure consistent units
        singlesphere = ft_convert_units(singlesphere, headmodel.unit);
        % determine the channels that do not have a corresponding sphere
        % and use the globally fitted single sphere for those
        missing = setdiff(sens.label, headmodel.label);
        if ~isempty(missing)
          ft_warning('using the global fitted single sphere for %d channels that do not have a local sphere', length(missing));
        end
        for i=1:length(missing)
          headmodel.label(end+1) = missing(i);
          headmodel.r(end+1,:)   = singlesphere.r;
          headmodel.o(end+1,:)   = singlesphere.o;
        end
      end

      % make a new structure that only holds the local spheres, one per coil
      localspheres = [];
      localspheres.type = headmodel.type;
      localspheres.unit = headmodel.unit;

      % for each coil in the MEG helmet, determine the corresponding channel and from that the corresponding local sphere
      for i=1:Ncoils
        coilindex = find(sens.tra(:,i)~=0); % to which channel does this coil belong
        if length(coilindex)>1
          % this indicates that there are multiple channels to which this coil contributes,
          % which happens if the sensor array represents a synthetic higher-order gradient.
          [dum, coilindex] = max(abs(sens.tra(:,i)));
        end

        coillabel = sens.label{coilindex};               % what is the label of this channel
        chanindex = find(strcmp(coillabel, headmodel.label));  % what is the index of this channel in the list of local spheres
        localspheres.r(i,:) = headmodel.r(chanindex);
        localspheres.o(i,:) = headmodel.o(chanindex,:);
      end
      headmodel = localspheres;

      % finally do the selection of channels and coils
      % order them according to the users specification
      [selchan, selsens] = match_str(channel, sens.label);

      % first only modify the linear combination of coils into channels
      try, sens.chantype = sens.chantype(selsens,:); end
      try, sens.chanunit = sens.chanunit(selsens,:); end
      try, sens.chanpos  = sens.chanpos (selsens,:); end
      try, sens.chanori  = sens.chanori (selsens,:); end
      sens.label   = sens.label(selsens);
      sens.tra     = sens.tra(selsens,:);
      % subsequently remove the coils that do not contribute to any sensor output
      selcoil      = find(sum(sens.tra,1)~=0);
      sens.coilpos = sens.coilpos(selcoil,:);
      sens.coilori = sens.coilori(selcoil,:);
      sens.tra     = sens.tra(:,selcoil);
      % make the same selection of coils in the localspheres model
      headmodel.r = headmodel.r(selcoil);
      headmodel.o = headmodel.o(selcoil,:);

    case 'singleshell'
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % if the forward model is computed using the code from Guido Nolte, we
      % have to initialize the volume model using the gradiometer coil
      % locations
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      % compute the surface normals for each vertex point
      if ~isfield(headmodel.bnd, 'nrm')
        fprintf('computing surface normals\n');
        headmodel.bnd.nrm = surface_normals(headmodel.bnd.pos, headmodel.bnd.tri);
      end

      % estimate center and radius
      [center, radius] = fitsphere(headmodel.bnd.pos);

      % order of spherical spherical harmonics, for 'real' realistic volume conductors order=10 is o.k
      if isfield(headmodel, 'order')
        order = headmodel.order;
      else
        order = 10;
      end

      % initialize the forward calculation (only if  coils are available)
      if size(sens.coilpos,1)>0 && ~isfield(headmodel, 'forwpar')
        s = ft_scalingfactor(headmodel.unit, 'cm');
        headmodel.forwpar = meg_ini([s*headmodel.bnd.pos headmodel.bnd.nrm], s*center', order, [s*sens.coilpos sens.coilori]);
        headmodel.forwpar.scale = s;
      end

    case  'openmeeg'
      % don't do anything, h2em or h2mm generated later in ft_prepare_leadfield

    case 'duneuro'

      %compute transfer matrix
      if(~isfield(headmodel,'meg_transfer'))

        % set coils and projections
        coils = sens.coilpos;
        projections = sens.coilori;
        headmodel.driver.set_coils_and_projections(coils', projections');

        % compute transfer matrix
        cfg = [];
        cfg.solver.reduction   = headmodel.reduction;
        cfg.solver.intorderadd = headmodel.intorderadd;
        headmodel.meg_transfer = headmodel.driver.compute_meg_transfer_matrix(cfg);
      end

    case 'interpolate'
      % this is to allow moving leadfield files
      if ~exist(headmodel.filename{1}, 'file')
        for i = 1:length(headmodel.filename)
          [p, f, x] = fileparts(headmodel.filename{i});
          headmodel.filename{i} = fullfile(vpath, [f x]);
        end
      end

      matchlab = isequal(sens.label, headmodel.sens.label);
      matchpos = isequal(sens.coilpos, headmodel.sens.coilpos);
      matchtra = (~isfield(sens, 'tra') && ~isfield(headmodel.sens, 'tra')) || (isfield(sens, 'tra') && isfield(headmodel.sens, 'tra') && isequal(sens.tra, headmodel.sens.tra));

      if matchlab && matchpos && matchtra
        % the input sensor array matches precisely with the forward model
        % no further interpolation is needed
      else
        % interpolate the channels in the forward model to the desired channels
        filename = tempname;
        headmodel  = ft_headmodel_interpolate(filename, sens, headmodel);
        % update the sensor array with the one from the volume conductor
        sens = headmodel.sens;
      end % if recomputing interpolation

      % for the leadfield computations the @nifti object is used to map the image data into memory
      ft_hastoolbox('spm8up', 1);
      for i=1:length(headmodel.sens.label)
        % map each of the leadfield files into memory
        headmodel.chan{i} = nifti(headmodel.filename{i});
      end

    case 'simbio'
      ft_error('MEG not yet supported with simbio');

    otherwise
      ft_error('unsupported volume conductor model for MEG');
  end

elseif iseeg

  % the electrodes are used, the channel positions are not relevant any more
  % channel positinos need to be recomputed after projecting the electrodes on the skin
  if isfield(sens, 'chanpos'); sens = rmfield(sens, 'chanpos'); end

  % select the desired channels from the electrode array
  % order them according to the users specification
  [selchan, selsens] = match_str(channel, sens.label);
  Nchans = length(sens.label);

  sens.label     = sens.label(selsens);
  try, sens.chantype  = sens.chantype(selsens); end
  try, sens.chanunit  = sens.chanunit(selsens); end

  if isfield(sens, 'tra')
    % first only modify the linear combination of electrodes into channels
    sens.tra     = sens.tra(selsens,:);
    % subsequently remove the electrodes that do not contribute to any channel output
    selelec      = any(sens.tra~=0,1);
    sens.elecpos = sens.elecpos(selelec,:);
    sens.tra     = sens.tra(:,selelec);
  else
    % the electrodes and channels are identical
    sens.elecpos = sens.elecpos(selsens,:);
  end

  switch ft_headmodeltype(headmodel)
    case {'infinite' 'infinite_monopole' 'infinite_currentdipole'}
      % nothing to do

    case {'halfspace', 'halfspace_monopole'}
      % electrodes' all-to-all distances
      numelec = size(sens.elecpos,1);
      ref_el = sens.elecpos(1,:);
      md = dist( (sens.elecpos-repmat(ref_el,[numelec 1]))' );
      % take the min distance as reference
      md = min(md(1,2:end));
      pos = sens.elecpos;
      % scan the electrodes and reposition the ones which are in the
      % wrong halfspace (projected on the plane)... if not too far away!
      for i=1:size(pos,1)
        P = pos(i,:);
        is_in_empty = acos(dot(headmodel.ori,(P-headmodel.pos)./norm(P-headmodel.pos))) < pi/2;
        if is_in_empty
          dPplane = abs(dot(headmodel.ori, headmodel.pos-P, 2));
          if dPplane>md
            ft_error('Some electrodes are too distant from the plane: consider repositioning them')
          else
            % project point on plane
            Ppr = pointproj(P,[headmodel.pos headmodel.ori]);
            pos(i,:) = Ppr;
          end
        end
      end
      sens.elecpos = pos;

    case {'slab_monopole'}
      % electrodes' all-to-all distances
      numelc  = size(sens.elecpos,1);
      ref_elc = sens.elecpos(1,:);
      md  = dist( (sens.elecpos-repmat(ref_elc,[numelc 1]))' );
      % choose min distance between electrodes
      md  = min(md(1,2:end));
      pos = sens.elecpos;
      % looks for contacts outside the strip which are not too far away
      % and projects them on the nearest plane
      for i=1:size(pos,1)
        P = pos(i,:);
        instrip1 = acos(dot(headmodel.ori1,(P-headmodel.pos1)./norm(P-headmodel.pos1))) > pi/2;
        instrip2 = acos(dot(headmodel.ori2,(P-headmodel.pos2)./norm(P-headmodel.pos2))) > pi/2;
        is_in_empty = ~(instrip1&instrip2);
        if is_in_empty
          dPplane1 = abs(dot(headmodel.ori1, headmodel.pos1-P, 2));
          dPplane2 = abs(dot(headmodel.ori2, headmodel.pos2-P, 2));
          if dPplane1>md && dPplane2>md
            ft_error('Some electrodes are too distant from the planes: consider repositioning them')
          elseif dPplane2>dPplane1
            % project point on nearest plane
            Ppr = pointproj(P,[headmodel.pos1 headmodel.ori1]);
            pos(i,:) = Ppr;
          else
            % project point on nearest plane
            Ppr = pointproj(P,[headmodel.pos2 headmodel.ori2]);
            pos(i,:) = Ppr;
          end
        end
      end
      sens.elecpos = pos;

    case {'singlesphere', 'concentricspheres'}
      % ensure that the electrodes ly on the skin surface
      radius = max(headmodel.r);
      pos    = sens.elecpos;
      if isfield(headmodel, 'o')
        % shift the the centre of the sphere to the origin
        pos(:,1) = pos(:,1) - headmodel.o(1);
        pos(:,2) = pos(:,2) - headmodel.o(2);
        pos(:,3) = pos(:,3) - headmodel.o(3);
      end
      distance = sqrt(sum(pos.^2,2)); % to the center of the sphere
      if any((abs(distance-radius)/radius)>0.005)
        ft_warning('electrodes do not lie on skin surface -> using radial projection')
      end
      pos = pos * radius ./ [distance distance distance];
      if isfield(headmodel, 'o')
        % shift the center back to the original location
        pos(:,1) = pos(:,1) + headmodel.o(1);
        pos(:,2) = pos(:,2) + headmodel.o(2);
        pos(:,3) = pos(:,3) + headmodel.o(3);
      end
      sens.elecpos = pos;

    case {'bem', 'dipoli', 'asa', 'bemcp'}
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % do postprocessing of volume and electrodes in case of BEM model
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      % project the electrodes on the skin and determine the bilinear interpolation matrix
      if ~isfield(headmodel, 'tra') && (isfield(headmodel, 'mat') && ~isempty(headmodel.mat))
        % determine boundary corresponding with skin and inner_skull_surface
        if ~isfield(headmodel, 'skin_surface')
          headmodel.skin_surface = find_outermost_boundary(headmodel.bnd);
          fprintf('determining skin compartment (%d)\n', headmodel.skin_surface);
        end
        if ~isfield(headmodel, 'source')
          headmodel.source = find_innermost_boundary(headmodel.bnd);
          fprintf('determining source compartment (%d)\n', headmodel.source);
        end
        if size(headmodel.mat,1)~=size(headmodel.mat,2) && size(headmodel.mat,1)==length(sens.elecpos)
          fprintf('electrode transfer and system matrix were already combined\n');
        else
          fprintf('projecting electrodes on triangulated skin surface\n');
          % compute linear interpolation from triangle vertices towards electrodes
          [el, prj] = project_elec(sens.elecpos, headmodel.bnd(headmodel.skin_surface).pos, headmodel.bnd(headmodel.skin_surface).tri);
          tra       = transfer_elec(headmodel.bnd(headmodel.skin_surface).pos, headmodel.bnd(headmodel.skin_surface).tri, el);

          % replace the original electrode positions by the projected positions
          sens.elecpos = prj;

          if size(headmodel.mat,1)==size(headmodel.bnd(headmodel.skin_surface).pos,1)
            % construct the transfer from only the skin vertices towards electrodes
            interp = tra;
          else
            % construct the transfer from all vertices (also inner_skull_surface/outer_skull_surface) towards electrodes
            interp = [];
            for i=1:length(headmodel.bnd)
              if i==headmodel.skin_surface
                interp = [interp, tra];
              else
                interp = [interp, zeros(size(el,1), size(headmodel.bnd(i).pos,1))];
              end
            end
          end

          % incorporate the linear interpolation matrix and the system matrix into one matrix
          % this speeds up the subsequent repeated leadfield computations
          fprintf('combining electrode transfer and system matrix\n');

          % convert to sparse matrix to speed up the subsequent multiplication
          interp  = sparse(interp);
          headmodel.mat = interp * headmodel.mat;
          % ensure that the model potential will be average referenced
          avg = mean(headmodel.mat, 1);
          headmodel.mat = headmodel.mat - repmat(avg, size(headmodel.mat,1), 1);
        end
      end
    case  'openmeeg'
      % don't do anything, h2em or h2mm generated later in ft_prepare_leadfield

    case 'fns'
      if isfield(headmodel,'bnd')
        [el, prj] = project_elec(sens.elecpos, headmodel.bnd.pos, headmodel.bnd.tri);
        sens.tra = transfer_elec(headmodel.bnd.pos, headmodel.bnd.tri, el);
        % replace the original electrode positions by the projected positions
        sens.elecpos = prj;
      end

    case 'simbio'
      % check that the external toolbox is present
      ft_hastoolbox('simbio', 1);

      % extract the outer surface
      bnd = mesh2edge(headmodel);
      for j=1:length(sens.label)
        d = bsxfun(@minus, bnd.pos, sens.elecpos(j,:));
        [d, i] = min(sum(d.^2, 2));
        % replace the position of each electrode by the closest vertex
        sens.elecpos(j,:) = bnd.pos(i,:);
      end

      if (isfield(headmodel,'transfer') && isfield(headmodel,'elec'))
        if all(ismember(sens.label,headmodel.elec.label))
          [sensmember, senslocation] = ismember(sens.label,headmodel.elec.label);
          if (norm(sens.elecpos - headmodel.elec.elecpos(senslocation,:))<1e-8)
            headmodel.transfer = headmodel.transfer(senslocation,:);
            headmodel.elec = sens;
          else
            ft_error('Electrode positions do not fit to the given transfer matrix!');
          end
        else
          ft_error('Transfer matrix does not fit the given set of electrodes!');
        end
      else
        headmodel.transfer = sb_transfer(headmodel,sens);
        headmodel.elec = sens;
      end

    case 'duneuro'
      ft_hastoolbox('duneuro', 1);

      if(~isfield(headmodel,'eeg_transfer'))
        % set electrodes
        cfg = [];
        cfg.type = headmodel.electrodes;
        cfg.codims = headmodel.subentities;
        headmodel.driver.set_electrodes(sens.elecpos', cfg);

        % compute transfer matrix
        cfg = [];
        cfg.solver.reduction = headmodel.reduction;
        headmodel.eeg_transfer = headmodel.driver.compute_eeg_transfer_matrix(cfg);
      end

    case 'interpolate'
      % this is to allow moving leadfield files
      if ~exist(headmodel.filename{1}, 'file')
        for i = 1:length(headmodel.filename)
          [p, f, x] = fileparts(headmodel.filename{i});
          headmodel.filename{i} = fullfile(vpath, [f x]);
        end
      end

      matchlab = isequal(sens.label, headmodel.sens.label);
      matchpos = isequal(sens.elecpos, headmodel.sens.elecpos);
      matchtra = (~isfield(sens, 'tra') && ~isfield(headmodel.sens, 'tra')) || (isfield(sens, 'tra') && isfield(headmodel.sens, 'tra') && isequal(sens.tra, headmodel.sens.tra));

      if matchlab && matchpos && matchtra
        % the input sensor array matches precisely with the forward model
        % no further interpolation is needed
      else
        % interpolate the channels in the forward model to the desired channels
        filename = tempname;
        headmodel  = ft_headmodel_interpolate(filename, sens, headmodel);
        % update the sensor array with the one from the volume conductor
        sens = headmodel.sens;
      end % if recomputing interpolation

      % for the leadfield computations the @nifti object is used to map the image data into memory
      ft_hastoolbox('spm8up', 1);
      for i=1:length(headmodel.sens.label)
        % map each of the leadfield files into memory
        headmodel.chan{i} = nifti(headmodel.filename{i});
      end

    otherwise
      ft_error('unsupported volume conductor model for EEG');
  end

  % FIXME this needs careful thought to ensure that the average referencing which is
  % now done here and there, and that the linear interpolation in case of BEM are all
  % dealt with consistently always ensure that there is a linear transfer matrix for
  % rereferencing the EEG potential
  %
  % if ~isfield(sens, 'tra');
  %   sens.tra = eye(length(sens.label)) - 1/length(sens.label);
  % end

  % update the channel positions as the electrodes were projected to the skin surface
  [pos, ori, lab] = channelposition(sens);
  [selsens, selpos] = match_str(sens.label, lab);
  sens.chanpos = nan(length(sens.label),3);
  sens.chanpos(selsens,:) = pos(selpos,:);

end % if iseeg or ismeg

if isfield(sens, 'tra')
  if issparse(sens.tra) && size(sens.tra, 1)==1
    % this multiplication would result in a sparse leadfield, which is not what we want
    % the effect can be demonstrated as sparse(1)*rand(1,10), see also http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1169#c7
    sens.tra = full(sens.tra);
  elseif ~issparse(sens.tra) && size(sens.tra, 1)>1
    % the multiplication of the "sensor" leadfield (electrode or coil) with the tra matrix to get the "channel" leadfield
    % is faster for most cases if the pre-multiplying weighting matrix is made sparse
    sens.tra = sparse(sens.tra);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Ppr = pointproj(P,plane)
% projects a point on a plane
% plane(1:3) is a point on the plane
% plane(4:6) is the ori of the plane
Ppr  = [];
ori  = plane(4:6);
line = [P ori];
% get indices of line and plane which are parallel
par = abs(dot(plane(4:6), line(:,4:6), 2))<1e-14;
% difference between origins of plane and line
dp = plane(1:3) - line(:, 1:3);
% Divide only for non parallel vectors (DL)
t = dot(ori(~par,:), dp(~par,:), 2)./dot(ori(~par,:), line(~par,4:6), 2);
% compute coord of intersection point
Ppr(~par, :) = line(~par,1:3) + repmat(t,1,3).*line(~par,4:6);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function serves as a replacement for the dist function in the Neural
% Networks toolbox.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [d] = dist(x)
n = size(x,2);
d = zeros(n,n);
for i=1:n
  for j=(i+1):n
    d(i,j) = sqrt(sum((x(:,i)-x(:,j)).^2));
    d(j,i) = d(i,j);
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [exitflag, message] = lvd_executeOptimProblem(celBodyData, writeOutput, problem, recorder, callOutputFcn)
    global options_gravParamType;
    initX = []; %#ok<NASGU>

    try
        useParallel = problem.options.UseParallel;
    catch ME %#ok<NASGU>
        if(isprop(problem.options,'UseParallel') || isfield(struct(problem.options),'UseParallel'))
            useParallel = problem.options.UseParallel;
        else
            useParallel = problem.UseParallel;
        end
    end

    if(useParallel)
        pp=gcp('nocreate');

        if(not(isempty(pp)))
            fnc3 = @() setGravParamTypeForWorker(options_gravParamType);
            pp.parfevalOnAll(fnc3, 0);
        end
    end

    lvdData = problem.lvdData;
    initX = lvdData.optimizer.vars.getTotalScaledXVector();

    writeOutput('Beginning mission script optimization...','append');
    tt = tic;
    if(strcmpi(problem.solver,'fmincon'))
        [x,fval,exitflag,output,lambda,grad,hessian] = fmincon(problem);
        message = output.message;

    elseif(strcmpi(problem.solver,'patternsearch'))
        [x,fval,exitflag,output] = patternsearch(problem);
        message = output.message;

    elseif(strcmpi(problem.solver,'nomad'))
        [x,fval,exitflag,iter,nfval] = nomad(problem.objective, problem.x0, problem.lb, problem.ub, problem.options);
        message = nomadExitFlagMessageLookup(exitflag);
%         [x,fval,exitflag,iter,nfval] = nomadOpt(problem.objective, problem.x0, problem.lb, problem.ub, problem.options);

    elseif(strcmpi(problem.solver,'ipopt'))
        [x,info] = ipopt(problem.x0, problem.funcs, problem.options);
        exitflag = info.status;
        message = ipoptExitFlagMessageLookup(exitflag);

    elseif(strcmpi(problem.solver,'surrogateopt'))
        [x,fval,exitflag,output,trials] = surrogateopt(problem);
        message = output.message;

    else
        error('Unknown optimizer function: %s', problem.solver);
    end

    execTime = toc(tt);
    writeOutput(sprintf('Mission script optimization finished in %0.3f sec with exit flag "%i".', execTime, exitflag),'append');

    %%%%%%%
    % Ask if the user wants to keep the current solution or not.
    %%%%%%%
    if(callOutputFcn)
%         [~, x] = ma_OptimResultsScorecardGUI(recorder);
        out = AppDesignerGUIOutput({[],[]});
        ma_OptimResultsScorecardGUI_App(recorder, out);
        x = out.output{2};
    end

    if(~isempty(x))
        writeOutput(sprintf('Optimization results accepted: merging with mission script.'),'append');

        %%%%%%%
        % Update existing script, reprocess
        %%%%%%%
        lvdData.optimizer.vars.updateObjsWithScaledVarValues(x);
    else
        writeOutput(sprintf('Optimization results discarded: reverting to previous script.'),'append');

        lvdData.optimizer.vars.updateObjsWithScaledVarValues(initX);
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [sampleLabels, groupCodingLengths] = coding_seg(W, epsilon, affine )

% coding_seg.m
%
% Implements an agglomerative segmentation algorithm which attempts to
% minimize the number of bits needed to code the segmented vectors upto
% distortion epsilon^2. This method is described in detail in Ma, et. al.,
% Segmentation of Multivariate Mixed Data via Lossy Coding and Compression,
% PAMI 2007.
%
% Run test_coding_seg.m, included in this package, for a demonstration.
%
% Inputs:
%   W       - [ w1 | w2 | ... | wm ], the data matrix.
%               The columns are the input data vectors.
%
%   epsilon - sqrt of the allowable distortion,
%               E[ \| \hat{w_i} - w_i \|^2 ]
%
%   affine  - boolean, if true the algorithm uses the coding length for
%               nonzero-mean data, derived in Appendix II of the paper.
%
% Outputs:
%   sampleLabels - 1 x m list whose i-th entry is the group assignment of
%                   the i-th data vector w_i. Groups are indexed
%                   sequentially, starting from 1.
%
% Dependencies:
%   coding_length.m
%
% Feb '06, last revision Jan '07
%   Questions? John Wright -- jnwright@uiuc.edu

% Copyright 2007, University of Illinois. All rights reserved.

    VERBOSE = 0;

    epsilon2 = epsilon^2;

    % n - dimension
    % m - number of samples
    [n,m] = size(W);

    % initially treat each sample as its own group
    sampleLabels = 1:m;
    curGroupCount = m;

    % compute the coding length of each of the initial groups
    squaredNormW = sum(W.^2,1);
    if affine
        groupCodingLengths = n * log2(1 + (squaredNormW/epsilon2))/ 2 + log2(m);
    else
        groupCodingLengths = (n + 1) * log2(1 + (n / epsilon2)*squaredNormW)/ 2 + log2(m);
    end

    if VERBOSE
        disp('   Computing initial table.');
    end

    % entropyChange is a matrix whose ij-th element is the entropy lost by
    % merging groups i and j into one group.
    delta_L_matrix = ones(m);

    for groupIndex1 = 1 : m-1
        for groupIndex2 = groupIndex1+1 : m
            L1 = groupCodingLengths(groupIndex1);
            L2 = groupCodingLengths(groupIndex2);
            L_union = coding_length(W(:, [groupIndex1 groupIndex2]),epsilon2,m,affine);
            delta_L_matrix(groupIndex1, groupIndex2) = L_union - L1 - L2;
        end;
    end;

    if VERBOSE
        disp('   Starting merging process');
    end

    while curGroupCount > 1

        if mod( curGroupCount, 50 ) == 0 && VERBOSE
            disp(['   Group count: ' num2str(curGroupCount)]);
        end

        % Find the best two groups to merge together
        [temp minIndex]=min(delta_L_matrix);
        [min_delta_L minIndex2]=min(temp);
        minIndex1=minIndex(minIndex2);

        % If this does not decrease the coding length, we are done
        if ( min_delta_L > 0 )
            break;
        end;

        % Merging groups with labels minIndex1 and minIndex2
        if (minIndex1 > minIndex2)
            temp = minIndex1;
            minIndex1 = minIndex2;
            minIndex2 = temp;
        end

        if (minIndex2 ~= curGroupCount)
            % First make the group to merge have the last label

            % Change the sample labels of the group minIndex2 to be the last
            % label.
            sampleLabels(sampleLabels == curGroupCount) = -1;
            sampleLabels(sampleLabels == minIndex2) = curGroupCount;
            sampleLabels(sampleLabels == -1) = minIndex2;

            % Swap the entropy of groups minIndex2 and the last group.
            [groupCodingLengths(minIndex2), groupCodingLengths(curGroupCount)] = ...
                swap(groupCodingLengths(minIndex2), groupCodingLengths(curGroupCount));

            delta_L_matrix = swap_matrix(delta_L_matrix, minIndex2, curGroupCount);

        end

        % Now we can assume the last group is the group to be merged
        sampleLabels(sampleLabels == curGroupCount) = minIndex1;
        groupCodingLengths(curGroupCount) = [];
        groupCodingLengths(minIndex1) = coding_length(W(:, sampleLabels == minIndex1),epsilon2,m,affine);
        delta_L_matrix = delta_L_matrix(1:end-1, 1:end-1);

        curGroupCount = curGroupCount - 1;

        for groupIndex2 = [1:minIndex1-1 minIndex1+1:curGroupCount],
            i = min(minIndex1, groupIndex2); j = max(minIndex1, groupIndex2);
                L1 = groupCodingLengths(i);
                L2 = groupCodingLengths(j);
                L_union = coding_length(W(:, (sampleLabels == i) | (sampleLabels == j)) ,epsilon2,m,affine);
                delta_L_matrix(i, j) = L_union - L1 - L2;
        end;
    end;

    %codingLength = sum(groupCodingLengths);

    if VERBOSE
        disp(['   Final group count: ' num2str(curGroupCount)]);
    end

end
% If A is an upper triangular matrix where A(m,n) = f(x_m,x_n) for m < n
% relabel so that the roles of x_i and x_j are swapped.
function A = swap_matrix(A, i, j)
[A(1:i-1, i), A(1:i-1, j)] = swap(A(1:i-1, i), A(1:i-1, j));
[A(i, i+1:j-1), A(i+1:j-1, j)] = swap(A(i, i+1:j-1), A(i+1:j-1, j));
[A(i, j+1:end), A(j, j+1:end)] = swap(A(i, j+1:end), A(j, j+1:end));
end

% swap the values of x and y
function [y, x] = swap(x,y)
end

author: SuTanTank

repo: VideoStitchingViaShakinessRemoving


text:
mex abs.c
mex max.c
mex min.c
mex minus.c
mex plus.c
mex rdivide.c
mex times.c


text:
function stack_samples = LDRStackSubSampling(stack, stack_exposure, nSamples, sampling_strategy, outliers_percentage )
%
%       stack_samples = LDRStackSubSampling(stack, stack_exposure, nSamples, sampling_strategy, outliers_percentage )
%
%       This function subsamples a stack
%
%        Input:
%           -stack: a stack of LDR images. If the stack is a single or
%           double values are assumed to be in [0,1]
%           -nSamples: number of samples for computing the CRF
%           -sampling_strategy: how to select samples:
%               -'Grossberg': picking samples according to Grossberg and
%               Nayar algorithm (CDF based)
%               -'RandomSpatial': picking random samples in the image
%               -'RegularSpatial': picking regular samples in the image
%               -outliers_percentage:
%
%        Output:
%           -stack_samples: sub-sampled stack
%
%     Copyright (C) 2016  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

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

if(isempty(stack))
    error('LDRStackSubSampling: a stack cannot be empty!');
end

if(~exist('nSamples', 'var'))
    nSamples = 256;
end

if(~exist('sampling_strategy', 'var'))
    sampling_strategy = 'Grossberg';
end

[~, sort_index] = sort(stack_exposure, 'descend');
sort_index = sort_index';

%stack sub-sampling
switch sampling_strategy
    case 'Grossberg'
        stack_hist = ComputeLDRStackHistogram(stack);
        stack_samples = GrossbergSampling(stack_hist, nSamples);

    case 'RandomSpatial'
        stack_samples = SpatialSampling(stack, sort_index, nSamples, sampling_strategy);
        stack_samples = round(stack_samples * 255);

    case 'RegularSpatial'
        stack_samples = SpatialSampling(stack, sort_index, nSamples, sampling_strategy);
        stack_samples = round(stack_samples * 255);

    otherwise
        error('LDRStackSubSampling: A sampling strategy neeeds to be specified.');
end

if(outliers_percentage > 0.0)
    t_min = outliers_percentage;
    t_max = 1.0 - t_min;
    stack_samples(stack_samples < (t_min * 255.0) | stack_samples > (t_max * 255)) = -1.0;
end

end

author: banterle


text:
function varargout=minspane(varargin)
%
% For calling details please see v_minspane.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_minspane.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: minspane.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_minspane(varargin{:});
else
    v_minspane(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
% COMPUTE_ERSP_TIMES - computes the widest possible ERSP/ITC time window,
%        which depends on requested ERSP/ITC parameters such as epoch limits,
%        frequency range, wavelet parameters, sampling rate and frequency
%        resolution that are used by TIMEF.
%        This helper function is called by POP_PRECLUST & STD_ERSP.
% Example:
%    [time_range, winsize] = compute_ersp_times(cycles,  ALLEEG(seti).srate, ...
%                              [ALLEEG(seti).xmin ALLEEG(seti).xmax]*1000, freq(1),padratio);
%
% Authors: Hilit Serby & Arnaud Delorme, SCCN, INC, UCSD, Feb 03, 2005

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

function [time_range, winsize] = compute_ERSP_times(cycles, srate, epoch_lim, lowfreq, padratio)

if cycles == 0 %FFT option
    if ~exist('padratio')
        error('You must enter padratio value for FFT ERSP');
    end
    lowfreq = lowfreq*padratio;
    t = 1/lowfreq;%time window in sec
    winsize = t*srate;%time window in points
    %time window in points (must be power of 2) for FFT
    winsize =pow2(nextpow2(winsize));
    %winsize =2^round(log2(winsize));
else %wavelet
    t = cycles(1)/lowfreq; %time window in sec
    winsize  = round(t*srate); %time window in points
end

time_range(1) = epoch_lim(1) + .5*t*1000;
time_range(2) = epoch_lim(2) - .5*t*1000;


text:
function geometry_test205 ( )

%*****************************************************************************80
%
%% TEST205 tests TMAT_MXP2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2005
%
%  Author:
%
%    John Burkardt
%
  n = 4;
  dim_num = 3;

  point = [ ...
    1.0, 0.0, 0.0; ...
    0.0, 1.0, 0.0; ...
    0.0, 0.0, 1.0; ...
    1.0, 1.0, 1.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST205\n' );
  fprintf ( 1, '  TMAT_MXP2 applies a geometric transformation\n' );
  fprintf ( 1, '  matrix to a set of points.\n' );

  r8mat_transpose_print ( 3, n, point, '  Points:' );
%
%  Initialization of transformation matrix.
%
  a = tmat_init ( );

  r8mat_print ( 4, 4, a, '  Initial transformation matrix:' );
%
%  Rotation about an axis.
%
  angle = 30.0;
  axis1 = 'x';
  b = tmat_rot_axis ( a, angle, axis1 );

  point2 = tmat_mxp2 ( b, n, point );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Rotation about %s\n', axis1 );
  fprintf ( 1, '  by %f\n' , angle );

  r8mat_transpose_print ( 3, n, point2, ' ' );
%
%  Rotation about a vector.
%
  angle = 30.0;
  axis(1:3) = [ 1.0, 2.0, 3.0 ];

  b = tmat_rot_vector ( a, angle, axis );

  point2 = tmat_mxp2 ( b, n, point );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Rotation about %f  %f  %f\n', axis(1:3) );
  fprintf ( 1, '  of %f\n', angle );

  r8mat_transpose_print ( 3, n, point2, ' ' );
%
%  Scaling.
%
  v(1:3) = [ 2.0, 0.5, 10.0 ];
  b = tmat_scale ( a, v );

  point2 = tmat_mxp2 ( b, n, point );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Scaling by %f  %f  %f\n', v(1:3) );

  r8mat_transpose_print ( 3, n, point2, ' ' );
%
%  Shear.
%
  axis2 = 'xy';
  s = 0.5;
  b = tmat_shear ( a, axis2, s );

  point2 = tmat_mxp2 ( b, n, point );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  %s\n', axis2 );
  fprintf ( 1, ' shear coefficient of %f\n', s );

  r8mat_transpose_print ( 3, n, point2, ' ' );
%
%  Translation.
%
  v(1:3) = [ 1.0, 2.0, 3.0 ];
  b = tmat_trans ( a, v );

  point2 = tmat_mxp2 ( b, n, point );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Translation by %f  %f  %f\n', v(1:3) );

  r8mat_transpose_print ( 3, n, point2, ' ' );

  return
end


text:
function value = r4_int ( x )

%*****************************************************************************80
%
%% R4_INT returns the integer part of an R4 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 October 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the integer part of X.
%
  persistent npart
  persistent scale
  persistent xbig
  persistent xmax

  if ( isempty ( npart ) )
    ibase = i4_mach ( 10 );
    xmax = 1.0 / r4_mach ( 4 );
    xbig = min ( i4_mach ( 9 ), xmax );
    expo = floor ( log ( xbig ) / log ( ibase ) - 0.5 );
    scale = ibase ^ expo;
    npart = floor ( log ( xmax ) / log ( scale ) + 1.0 );
  end

  if ( x < - xmax )

    value = x;

  elseif ( x < - xbig )

    xscl = - x;

    for i = 1 : npart
      xscl = xscl / scale;
    end

    value = 0.0;
    for i = 1 : npart
      xscl = xscl * scale;
      ipart = ceil ( xscl );
      part = ipart;
      xscl = xscl - part;
      value = value * scale + part;
    end

    value = - value;

  else if ( x < 0 )

    value = ceil ( x );

  elseif ( x < + xbig )

    value = floor ( x );

  elseif ( x < + xmax )

    xscl = x;

    for i = 1 : npart
      xscl = xscl / scale;
    end

    value = 0.0;
    for i = 1 : npart
      xscl = xscl * scale;
      ipart = floor ( xscl );
      part = ipart;
      xscl = xscl - part;
      value = value * scale + part;
    end

  else

    value = x;

  end

  return
end


text:
% clc; clear all; close all;

%% hydro data
hydro = struct();
hydro = readAQWA(hydro, 'WEC3.AH1', 'WEC3.LIS');
hydro = radiationIRF(hydro,100,[],[],[],[]);
hydro = radiationIRFSS(hydro,[],[]);
hydro = excitationIRF(hydro,100,[],[],[],[]);
writeBEMIOH5(hydro)

%% Plot hydro data
% plotBEMIO(hydro)

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
classdef Tx < adi.AD9144.Tx
    % adi.DAQ2.Tx Transmit data from the DAQ2 evaluation platform
    %   The adi.DAQ2.Tx System object is a signal source that can
    %   send complex data to the DAQ2.
    %
    %   tx = adi.DAQ2.Tx;
    %   tx = adi.DAQ2.Tx('uri','192.168.2.1');
    %
    %   User Guide
    %
    %   See also adi.AD9144.Tx, adi.DAQ2.Rx
    methods
        %% Constructor
        function obj = Tx(varargin)
            % Returns the matlabshared.libiio.base object
            coder.allowpcode('plain');
            obj = obj@adi.AD9144.Tx(varargin{:});
        end
    end

end

author: analogdevicesinc


text:
function [rglm,yclean] = spm_rglm (y,X,m,priors,verbose)
% Fit a Robust GLM
% FORMAT [rglm,yclean] = spm_rglm (y,X,m,priors,verbose)
%
% The noise is modelled with a Mixture of Zero-Mean Gaussians
%
% y          [N x 1] data vector
% X          [N x p] design matrix
% m          Number of mixture components
% priors     .alpha      [1 x 1] weight precision (default=0.001)
%            .mean_err   [m x 1] vector of mean error SD
%            .std_err    [m x 1] vector of dev of error SD
% verbose    0/1 to printout inner workings (default=0)
%
% rglm       Returned model
% yclean     'Clean' data
%
% -------------------------------------------------------
% The fields in rglm are:
%
% m                The number of error components
% fm               The negative free energy
% loops            Number of iterations used
%
%                  In the field priors:
%
% lambda_0         Dirichlet parameters for mixing coeffs
% b_0,c_0          Gamma parameters for precisions
%
%                  In the field posts:
%
% lambda           Dirichlet parameters  for mixing coeffs
% b,c              Gamma parameters for precisions
% w_mean           Mean estimated regression coefficients
% w_cov            Covariance of regression coefficients
% pi               mixing coefficients (lambda/sum(lambda))
% variances        variances (1./(b.*c))
% gamma            the responsilities of each noise component
%___________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_rglm.m 1276 2008-03-28 18:29:19Z guillaume $

if nargin < 4 | isempty(priors)
    mean_alpha=0.001;
    % Set variance priors
    b_0=1000*ones(m,1);
    c_0=0.001*ones(m,1);
else
    mean_alpha=priors.alpha;
    for mm=1:m,
        b_0(mm)=(priors.std_err(mm)^2)/priors.mean_err(mm);
        c_0(mm)=priors.mean_err(mm)/b_0(mm);
    end
end

if (m==1)
    rglm=spm_glm(y,X,mean_alpha);
    return
end

if nargin < 5 | isempty(verbose)
    verbose=0;
end

y=y(:);
N=length(y);
p=size(X,2);
spX=issparse(X);

% Initialise regression coefficients to maximum likelihood solution
if spX
    iX=inv(X'*X);
    w_mean=iX*X'*y;
else
    iX=inv(X'*X);
    w_mean = pinv(X)*y;
end
y_pred = full(X*w_mean);
err=y-y_pred;
v=mean((y-y_pred).^2);
w_cov = v*iX;

% Set mixing priors
lambda_0=5;

% Cluster on absolute difference from mean
zmix=spm_kmeans1(abs(err-mean(err)),m);
% Posterior for mixers
lambda=100*zmix.pi+lambda_0;
zmean=[zmix.m].^2;
% Posterior for precisions
var_precision=1/std(zmean)^2;
for s=1:m,
      % Set so that b*c=precision and  b^2*c=var_precision
      precision=1/zmean(s);
      b(s)=var_precision/precision;
      c(s)=(precision^2)/var_precision;
end

if verbose
    disp('Init');
    for s=1:m,
        disp(sprintf('State %d mix=%1.2f var=%1.2f',s,lambda(s)/sum(lambda),1/(b(s)*c(s))));
    end
end

lik=[];
tol=0.0001;
max_loops=100;
WLOOPS=5;
for loops=1:max_loops,

    % E-step
    lambda_tot=sum(lambda);
    ypred=full(X*w_mean);
    ypred2=ypred.^2;
    y2=y.^2;
    xt=X';
    y_err=full(sum(xt.*(w_cov*xt)));

    tv=y2-2*ypred.*y+y_err'+ypred2;
    tv=tv';
    for s=1:m,
        log_tilde_pi(s)=psi(lambda(s))-psi(lambda_tot);
        log_tilde_beta(s)=psi(c(s))+log(b(s));
        tilde_pi(s)=exp(log_tilde_pi(s));
        tilde_beta(s)=exp(log_tilde_beta(s));
        mean_beta(s)=c(s)*b(s);
        tilde_var(s,:)=tv;
        gamma(s,:)=tilde_pi(s)*(tilde_beta(s)^0.5)*exp(-0.5*mean_beta(s)*tv);
    end
    gamma_n=sum(gamma);
    for s=1:m,
        if mean(gamma_n) > eps
            % If component still exists
            gamma(s,:)=gamma(s,:)./gamma_n;
        end
    end

    % M-step
    % Part I
    for s=1:m,
        pi_bar(s)=mean(gamma(s,:));
        N_bar(s)=N*pi_bar(s);
        mean_bar(s)=mean(gamma(s,:)'.*y);
        dg=sparse(1:N,1:N,gamma(s,:));
        x_bar(s,:)=mean(dg*X);
        var_bar_mu(s)=mean(gamma(s,:).*tilde_var(s,:));
    end

    % CALCULATE THE FREE ENERGY
    avg_likelihood=sum(N_bar.*(log_tilde_pi+0.5*log_tilde_beta));
    fit=-0.5*N*sum(mean_beta.*var_bar_mu);
    % Ensure that 0 log 0 = 0
    ent_s=sum(sum(-gamma.*log(gamma+eps)));
    avg_likelihood=avg_likelihood+fit+ent_s;
    lambda_p=lambda_0*ones(1,m);
    kl_dir=spm_kl_dirichlet(lambda,lambda_p,log_tilde_pi);
    kl_gamm=0;
    for s=1:m,
        kl_gamm=kl_gamm+spm_kl_gamma(b(s),c(s),b_0(s),c_0(s));
    end
    kl_weights=spm_kl_normald(w_mean,w_cov,zeros(1,p),(1/mean_alpha)*eye(p));
    fm= avg_likelihood - kl_dir - kl_gamm - kl_weights;

    % Convergence criterion
    oldlik=lik;
    lik=fm;
    if (mod(loops-1,WLOOPS)==0)
        if (loops>1)
            if abs((lik-oldlik)/lik) < tol
                break;
            end
        end
    end

    % M-Step: Part II
    for s=1:m,
        % Mixers
        lambda(s)=N_bar(s)+lambda_0;
        % Precisions
        b(s)=1/(0.5*N*var_bar_mu(s)+1/b_0(s));
        c(s)=0.5*N_bar(s)+c_0(s);
        mean_beta(s)=c(s)*b(s);
    end

    if mod(loops,WLOOPS)==0
        % Regression coefficients
        cc=sparse(p,p);
        cw=zeros(p,1);
        for s=1:m,
            dg=sparse(1:N,1:N,gamma(s,:));
            cc=cc+mean_beta(s)*X'*dg*X;
            cw=cw+mean_beta(s)*X'*dg*y;
        end
        cc=cc+mean_alpha*speye(p);
        w_cov=inv(cc);
        w_mean=w_cov*cw;
    end

    if verbose
        disp(sprintf('It=%d, L_AV =%1.2f, KL Mix=%1.2f, KL Prec=%1.2f, KL-W=%1.2f, Fm=%1.2f',loops,avg_likelihood,kl_dir,kl_gamm,kl_weights,fm));
        disp(' ');
        disp(sprintf('Iteration number=%d',loops));
        for s=1:m,
            var=1/(b(s)*c(s));
            disp(sprintf('State %d pi_bar=%1.2f var=%1.2f',s,pi_bar(s),var));
        end
    end

end

% Put variables into data structure
rglm.m=m;
rglm.fm=fm;
rglm.kl_dir=kl_dir;
rglm.kl_gamm=kl_gamm;
rglm.kl_w=kl_weights;
rglm.mean_alpha=mean_alpha;
rglm.priors.lambda_0=lambda_0;
rglm.priors.c_0=c_0;
rglm.priors.b_0=b_0;
for k=1:m,
  rglm.posts.lambda(k)=lambda(k);
end
for k=1:m,
  rglm.posts.c(k)=c(k);
  rglm.posts.b(k)=b(k);
end
rglm.posts.w_mean=w_mean;
rglm.posts.w_cov=w_cov;

for k=1:m,
  rglm.posts.pi(k)=lambda(k)/sum(lambda);
end
rglm.posts.variances=1./(b.*c);
rglm.posts.gamma=gamma;
rglm.loops=loops;

% Get 'clean' data
if m > 1
    [tmp,outlier_class]=min(rglm.posts.pi);
    e=y-ypred;
    outlier_error=rglm.posts.gamma(outlier_class,:)'.*e;
    yclean=ypred+e-outlier_error;
else
    yclean=y;
end

author: spm


text:
%DEMOS_PHRIWARE shows some demos of the pHRIWARE software package
%
% Opens a menu to choose some demonstrations, which give step-by-step
% explanations in the Command Window of using certain functionalities
% of pHRIWARE. Based on the demo format of RTB by Peter Corke.
%
% Copyright (C) Bryan Moutrie, 2013-014
% Licensed under the GNU General Public License, see file for statement
%
% This file modifies file(s) from The Robotics Toolbox for MATLAB (RTB)
% by Peter Corke (www.petercorke.com), see file for statement
%
% See also pHRIWARE

% LICENSE STATEMENT:
%
% This file is part of pHRIWARE.
%
% pHRIWARE 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.
%
% pHRIWARE 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 pHRIWARE.  If not, see <http://www.gnu.org/licenses/>.
%
% RTB LIBRARY:
%
% Copyright (C) 1993-2014, by Peter I. Corke
% http://www.petercorke.com
% Released under the GNU Lesser General Public license,
% Modified 24/7/2014 (code adapted from rtbdemo.m, uses runscript.m)

demos = {
    '*HAL basics', 'hal';
    '*Using HAL', 'hal2';
    'Load HAT collision model', 'hat';
    'Create a CollisionModel', 'collisionmodel';
    '*collisions basics', 'collisions';
    '*collisions with dynamic objects', 'collisions2';
    'Exit', '';
    };

while true
    selection = menu('*Requires Robotics Toolbox', demos{:,1});

    if strcmp(demos{selection,1}, 'Exit')
        % quit now
        delete( get(0, 'Children') );
        break;
    else
        % run the appropriate script
        script = demos{selection,2};
        runscript(['demo_', script])
    end
end

author: Allopart

repo: rbpf-gmapping


text:
function [sizeROI] = getSize(ROIonly,pixelW,sliceS)
% -------------------------------------------------------------------------
% function [sizeROI] = getSize(ROIonly,pixelW,sliceS)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes the size (longest diameter) of the region of
% interest (ROI) of an input volume.
% -------------------------------------------------------------------------
% INPUTS:
% - ROIonly: 3D array, with voxels outside the ROI set to NaNs.
% - pixelW: Pixel width, or in-plane resolution, in mm.
% - sliceS: Slice spacing, in mm.
% -------------------------------------------------------------------------
% OUTPUTS:
% - sizeROI: Longest diameter of the ROI, in mm.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: May 2015
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package 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 package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

% IMPORTANT: SMALL ERROR DUE TO Z EXTREMA DEFINED BY CENTER OF VOXELS, AND
% NOT BY EXTREME FACE DISTANCES OF VOXELS (as in regionprops.m). ANOTHER
% WAY WOULD BE TO NOT CONSIDER FACE EXTREMA AND ONLY CENTER OF VOXEL. ONE
% OR THE OTHER SOLUTION SHOULD BE USED, BUT NOT A HYBRID SOLUTION LIKE
% RIGHT NOW.
% --> TO BE CORRECTED ON NEXT RELEASE (message written on Sep 12, 2016)

mask = ~isnan(ROIonly); % Find mask covering the ROI
vectX = zeros(1,8*size(mask,3));
vectY = zeros(1,8*size(mask,3));
vectZ = zeros(1,8*size(mask,3));

for i = 1:size(mask,3)
    temp = regionprops(mask(:,:,i),'Extrema');
    try
        temp = temp.Extrema;
        temp = temp';
        vectX((i-1)*8 + 1:i*8) = temp(1,:) * pixelW;
        vectY((i-1)*8 + 1:i*8) = temp(2,:) * pixelW;
        vectZ((i-1)*8 + 1:i*8) = (i-1) * sliceS;
    catch % Will always work, except when the first slice contains all zeros
        vectX((i-1)*8 + 1:i*8) = vectX(((i-1)-1)*8 + 1:(i-1)*8);
        vectY((i-1)*8 + 1:i*8) = vectY(((i-1)-1)*8 + 1:(i-1)*8);
        vectZ((i-1)*8 + 1:i*8) = vectZ(((i-1)-1)*8 + 1:(i-1)*8);
    end
end

max = 0;
for i = 1:8*size(mask,3) - 1
    for j = i + 1:8*size(mask,3)
       dist = (vectX(i)-vectX(j))^2 + (vectY(i)-vectY(j))^2 + (vectZ(i)-vectZ(j))^2;
       if dist > max
           max = dist;
       end
    end
end

sizeROI = sqrt(max);

end

author: mvallieres

repo: radiomics


text:
function score = visualize(board,goal,solutionVector)
solutionVector = validateSolutionVector(solutionVector, board);
boardToScore  = runSolutionVector(solutionVector, board, goal, true);
score = sum(boardToScore(:)) + sum(solutionVector);


text:
function p = predict(Theta1, Theta2, X)
%PREDICT Predict the label of an input given a trained neural network
%   p = PREDICT(Theta1, Theta2, X) outputs the predicted label of X given the
%   trained weights of a neural network (Theta1, Theta2)

% Useful values
m = size(X, 1);
num_labels = size(Theta2, 1);

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

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned neural network. You should set p to a
%               vector containing labels between 1 to num_labels.
%
% Hint: The max function might come in useful. In particular, the max
%       function can also return the index of the max element, for more
%       information see 'help max'. If your examples are in rows, then, you
%       can use max(A, [], 2) to obtain the max for each row.
%

Z = zeros(size(X, 1), 1);

X = [ones(m, 1) X];

A2 = sigmoid(X * Theta1');

A2 = [ones(m, 1), A2];

A3 = sigmoid(A2 * Theta2');

[Z, p] = max(A3, [], 2);

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

end

author: UtkarshPathrabe


text:
function RDot = rotationDerivativeFromAngVel(omega,R)

    % ROTATIONDERIVATIVEFROMANGVEL computes the derivative of a rotation matrix
    %                              given the angular velocity. It makes use of
    %                              the following convention: if the rotation of
    %                              a body b is expressed in the WORLD FRAME w
    %                              (i.e., R = w_R_b), then this function is
    %                              expecting the angular velocity of the body
    %                              to be expressed in the WORLD FRAME, i.e.
    %                              omega = w_omega.
    %
    % USAGE: please note that this function has been designed for being inserted
    %        in a Simulink model.
    %
    % FORMAT: RDot = rotationDerivativeFromAngVel(omega,R)
    %
    % INPUT:  - omega = [3 * 1] angular velocity
    %         - R = [3 * 3] rotation matrix
    %
    % OUTPUT: - RDot = [3 * 3] rotation matrix derivative
    %
    % Authors: Daniele Pucci, Marie Charbonneau, Gabriele Nava
    %
    %          all authors are with the Italian Istitute of Technology (IIT)
    %          email: name.surname@iit.it
    %
    % Genoa, Dec 2017
    %

    %% --- Initialization ---

    % Rotation matrix derivative (with correction to keep the integration
    % inside the space of rotation matrices). Omega is assumed to be w.r.t.
    % the inertial frame, i.e. w_omega.
    kCorr   = 1;
    RDot    = wbc.skew(omega)*R +kCorr*(eye(3)-R*transpose(R))*R;
end

author: robotology

repo: whole-body-controllers


text:
% test for computation of curvature on meshes

name = 'mushroom';
name = 'fandisk';
name = 'horse';
name = 'nefertiti';
name = 'elephant-50kv';
name = 'bunny';
name = 'skull';
name = 'lion-head';
name = 'vase-lion';
name = 'david50kf';
name = 'david_head';
name = 'armadillo';
name = 'aphro';
name = 'gargoyle';
name = 'tyra';
name = 'screwdriver';
name = 'hand';
options.name = name;

path(path, '../toolbox_graph_data/off/');

rep = ['results/curvature/' name '/'];
if not(exist(rep))
    mkdir(rep);
end

[vertex,face] = read_mesh(name);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% local covariance analysis
options.covariance_smoothing = 15;
[C,U,D] = compute_mesh_local_covariance(vertex,face,vertex,options);

% options for display
tau = 1.2;
options.normal_scaling = 1.5;

options.normal = squeeze(U(:,2,:));
clf;
options.face_vertex_color = perform_saturation( -D(2,:)' - D(3,:)',tau);
plot_mesh(vertex,face, options);
shading interp; camlight; colormap jet(256);
saveas(gcf, [rep name '-covariance.png'], 'png');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% curvature

options.curvature_smoothing = 10;
[Umin,Umax,Cmin,Cmax,Cmean,Cgauss,Normal] = compute_curvature(vertex,face,options);

options.normal = [];
clf;
options.face_vertex_color = perform_saturation(Cmax,tau);
plot_mesh(vertex,face, options);
shading interp; camlight; colormap jet(256);
saveas(gcf, [rep name '-cmax.png'], 'png');
clf;
options.face_vertex_color = perform_saturation(Cmin,tau);
plot_mesh(vertex,face, options);
shading interp; camlight; colormap jet(256);
saveas(gcf, [rep name '-cmin.png'], 'png');
clf;
options.face_vertex_color = perform_saturation(Cmean,tau);
plot_mesh(vertex,face, options);
shading interp; camlight; colormap jet(256);
saveas(gcf, [rep name '-cmean.png'], 'png');
clf;
options.face_vertex_color = perform_saturation(Cgauss,tau);
plot_mesh(vertex,face, options);
shading interp; camlight; colormap jet(256);
saveas(gcf, [rep name '-cgauss.png'], 'png');
clf;
options.face_vertex_color = perform_saturation(abs(Cmin)+abs(Cmax),tau);
plot_mesh(vertex,face, options);
shading interp; camlight; colormap jet(256);
saveas(gcf, [rep name '-cabs.png'], 'png');

author: gpeyre

repo: matlab-toolboxes


text:
function [tests,pass,perf]=test_snip(varargin)

inputs={'verbose'};
verbose=0;
for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};']);
    end
end

%Test the examples
test_string={'Fs=360;err=snip(''mitdb/100'',''100cut'',Fs*60,Fs*2*60);'};
clean_up={['delete([pwd filesep ''100cut*'']);']};
[tests,pass,perf]=test_wrapper(test_string,clean_up,verbose);

author: ikarosilva

repo: wfdb-app-toolbox


text:
function z=rpsola(x,y,t)
  % The solution of a Riemann problem (continuous)
  z=zeros(size(x));
  lne=0.25*t;
  lm=-0.5*t;
  lM=0.5*t;

  ind=find((x>=lM)&(y>=lM));
  z(ind)=0.5;

  ind=find((x<lM)&(x>=lne)&(y>=lne)&(y>=x));
  z(ind)=x(ind)/t;

  ind=find((x>=lne)&(y<lM)&(y>=lne)&(y<x));
  z(ind)=y(ind)/t;

  ind=find(((x>=lne)&(y<lne))|((x<lne)&(y>=lne)));
  z(ind)=0.25;

  ind=find((x<lne)&(x>=lm)&(x>=y));
  z(ind)=x(ind)/t;

  ind=find((y<lne)&(y>=lm)&(x<y));
  z(ind)=y(ind)/t;

  ind=find((x<=lm)&(y<=lm));
  z(ind)=-0.5;

repo: Aero-matlab


text:
function test_issue2075

% WALLTIME 00:10:00
% MEM 4gb
% DEPENDENCY nanmean

% A user reported an issue related to nanmean, using MATLAB 2022a.
% Operating system unknown. This script tries to reproduce it. On our end,
% running on Linux we cannot reproduce an error (the exact error is as of
% yet still underspecified). The purpose of this function is to document
% it, being aware that the test functions are not necessarily executed in
% their automatic overnight batch in the reported matlab version

[ftver, ftpath] = ft_version;

global ft_default
ft_default.toolbox.signal = 'compat';
ft_default.toolbox.stats= 'compat';
restoredefaultpath;
addpath(ftpath);
ft_defaults;
ft_preproc_bandpassfilter(randn(1,1000), 1000, [230 0.5], [], 'fir');
which nanmean

ft_default.toolbox.signal = 'matlab'; % rather than compat
ft_default.toolbox.stats= 'matlab';
restoredefaultpath;
addpath(ftpath);
ft_defaults;
ft_preproc_bandpassfilter(randn(1,1000), 1000, [230 0.5], [], 'fir');
which nanmean


text:
% Test file for trigtechh/mtimes.m

function pass = test_mtimes(pref)

% Get preferences.
if ( nargin < 1 )
    pref = trigtech.techPref();
end

testclass = trigtech();

% Generate a few random points to use as test values.
seedRNG(6178);
x = 2 * rand(100, 1) - 1;

% A random number to use as an arbitrary scalar multiplier.
alpha = randn() + 1i*randn();

%%
% Check operation in the face of empty arguments.

f = testclass.make(@(x) sin(10*pi*x), [], pref);
g = testclass.make();
pass(1) = isempty(f*[]) && isempty([]*f) && isempty(2*g) && isempty(g*2);

%%
% Check operation for scalar TRIGTECH objects.

f = testclass.make(@(x) sin(10*pi*x), [], pref);
g1 = alpha*f;
g2 = f*alpha;
pass(2) = isequal(g1, g2);
g_exact = @(x) alpha*sin(10*pi*x);
pass(3) = norm(feval(g1, x) - g_exact(x), inf) < ...
    100*vscale(g1)*eps;

g = 0*f;
pass(4) = all(g.coeffs == 0);

%%
% Check operation for array-valued TRIGTECH objects.

f = testclass.make(@(x) [sin(10*pi*x) cos(20*pi*x) cos(sin(pi*x))], [], pref);
g1 = alpha*f;
g2 = f*alpha;
pass(5) = isequal(g1, g2);
g_exact = @(x) alpha*[sin(10*pi*x) cos(20*pi*x) cos(sin(pi*x))];
err = abs(feval(g1, x) - g_exact(x));
pass(6) = max(err(:)) < 100*max(vscale(g1)*eps);

g = 0*f;
pass(7) =  all(g.coeffs == 0);

A = randn(3, 3);
g = f*A;
g_exact = @(x) [sin(10*pi*x) cos(20*pi*x) cos(sin(pi*x))]*A;
err = abs(feval(g, x) - g_exact(x));
pass(8) = max(err(:)) < 100*max(vscale(g)*eps);

f = testclass.make(@(x) [exp(1i*11*pi*x) cos(20*pi*x) cos(sin(pi*x))], [], pref);
g = f*A;
g_exact = @(x) [exp(1i*11*pi*x) cos(20*pi*x) cos(sin(pi*x))]*A;
err = abs(feval(g, x) - g_exact(x));
pass(9) = max(err(:)) < 100*max(vscale(g)*eps);

f = testclass.make(@(x) [exp(1i*11*pi*x) cos(20*pi*x) cos(sin(pi*x))], [], pref);
A = randn(3, 3) + 1i*randn(3, 3);
g = f*A;
g_exact = @(x) [exp(1i*11*pi*x) cos(20*pi*x) cos(sin(pi*x))]*A;
err = abs(feval(g, x) - g_exact(x));
pass(10) = max(err(:)) < 100*max(vscale(g)*eps);

%%
% Verify error handling and corner cases.

% Multiply non-scalar double and TRIGTECH.
try
    f = testclass.make(@(x) exp(cos(pi*x)));
    disp([1 2 3]*f);
    pass(11) = false;
catch ME
    pass(11) = strcmp(ME.identifier, 'CHEBFUN:TRIGTECH:mtimes:size') ...
        && strcmp(ME.message, 'Inner matrix dimensions must agree.');
end

% Multiply TRIGTECH and non-scalar double with mismatching dimensions.
try
    f = testclass.make(@(x) [sin(10*pi*x) cos(20*pi*x)]);
    disp(f*[1 ; 2 ; 3]);
    pass(12) = false;
catch ME
    pass(12) = strcmp(ME.identifier, 'CHEBFUN:TRIGTECH:mtimes:size2') ...
        && strcmp(ME.message, 'Inner matrix dimensions must agree.');
end

% Using * for multiplication of two TRIGTECH objects.
try
    g = testclass.make(@(x) cos(20*pi*x));
    disp(f*g);
    pass(13) = false;
catch ME
    pass(13) = strcmp(ME.message, 'Use .* to multiply TRIGTECH objects.');
end

% Using * to multiply a TRIGTECH and something else.
try
    disp(f*uint8(128));
    pass(14) = false;
catch ME
    pass(14) = strcmp(ME.message, ...
        'mtimes does not know how to multiply a TRIGTECH and a uint8.');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun