text:
function stats = UFstats (A, kind, nometis, skip_chol, skip_dmperm, Z)
%UFSTATS compute matrix statistics for the UF Sparse Matrix Collection
% Example:
%   stats = UFstats (A, kind, nometis, skip_chol, skip_dmperm, Z)
%
% A: a sparse matrix
% kind: a string with the Problem.kind
% nometis: if nonzero then metis(A,'col') is not used, nor is metis used in the
%       dmperm+ ordering.
% Z: empty, or a sparse matrix the same size as A.  Only used for psym and
%       nzero statistics, described below.
%
% Requires amd, cholmod, metis, RBio, and CSparse.  Computes the following
% statistics, returning them as fields in the stats struct:
%
%   nrows           number of rows
%   ncols           number of columns
%   nnz             number of entries in A
%   RBtype          Rutherford/Boeing type
%   isBinary        1 if binary, 0 otherwise
%   isReal          1 if real, 0 if complex
%   cholcand        1 if a candidate for sparse Cholesky, 0 otherwise
%   nsym            numeric symmetry (0 to 1, where 1=symmetric)
%   psym            pattern symmetry (0 to 1, where 1=symmetric)
%   nnzdiag         nnz (diag (A)) if A is square, 0 otherwise
%   nzero           nnz (Z)
%   amd_lnz         nnz(L) for chol(C(p,p)) where, C=A+A', p=amd(C)
%   amd_flops       flop count for chol(C(p,p)) where, C=A+A', p=amd(C)
%   amd_vnz         nnz in Householder vectors for qr(A(:,colamd(A)))
%   amd_rnz         nnz in R for qr(A(:,colamd(A)))
%   metis_lnz       nnz(L) for chol(C(p,p)) where, C=A+A', p=metis(C)
%   metis_flops     flop count for chol(C(p,p)) where, C=A+A', p=metis(C)
%   metis_vnz       nnz in Householder vectors for qr(A(:,metis(A,'col')))
%   metis_rnz       nnz in R for qr(A(:,metis(A,'col')))
%   nblocks         # of blocks from dmperm
%   sprank          sprank(A)
%   nzoff           # of entries not in diagonal blocks from dmperm
%   ncc             # of strongly connected components
%   dmperm_lnz      nnz(L), using dmperm plus amd or metis
%   dmperm_unz      nnz(U), using dmperm plus amd or metis
%   dmperm_flops    flop count with dperm plus
%   dmperm_vnz      nnz in Householder vectors for dmperm plus
%   dmperm_rnz      nnz in R for dmperm plus
%   posdef          1 if positive definite, 0 otherwise
%   isND	    1 if a 2D/3D problem, 0 otherwise
%
% The *_lnz, *_unz, and *_flops statistics are not computed for rectangular
% or structurally singular matrices.  nzoff and the dmperm_* stats are not
% computed for structurally singular matrices.  If a statistic is not computed,
% it is set to -2.  If an attempt to compute the statistic was made but failed,
% it is set to -1.
%
% See also UFget, UFindex, amd, metis, RBtype, cs_scc, cs_sqr, dmperm.

% Copyright 2006-2007, Timothy A. Davis

% Requires the SuiteSparse set of packages: CHOLMOD, AMD, COLAMD, RBio, CSparse;
% and METIS.

if (nargin < 3)
    nometis = 0 ;
end

%-------------------------------------------------------------------------------
% ensure the matrix is sparse
%-------------------------------------------------------------------------------

if (~issparse (A))
    A = sparse (A) ;
end

%-------------------------------------------------------------------------------
% basic stats
%-------------------------------------------------------------------------------

tic ;
[m n] = size (A) ;
stats.nrows = m ;
stats.ncols = n ;
stats.nnz = nnz (A) ;
stats.RBtype = RBtype (A) ;			% Rutherford/Boeing type
stats.isBinary = (stats.RBtype (1) == 'p') ;
stats.isReal = (stats.RBtype (1) ~= 'c') ;

fprintf ('RBtype: %s time: %g\n', stats.RBtype, toc) ;

%-------------------------------------------------------------------------------
% symmetry and Cholesky candidacy
%-------------------------------------------------------------------------------

% get the symmetry
tic ;
[s xmatched pmatched nzoffdiag nnzdiag] = spsym (A) ;

stats.cholcand = (s >= 6) ; % check if Cholesky candidate

if (m ~= n)
    stats.nsym = 0 ;
    stats.psym = 0 ;
elseif (nzoffdiag > 0)
    stats.nsym = xmatched / nzoffdiag ;
    stats.psym = pmatched / nzoffdiag ;
else
    stats.nsym = 1 ;
    stats.psym = 1 ;
end

fprintf ('cholcand: %d\n', stats.cholcand) ;
fprintf ('nsym: %g psym: %g time: %g\n', stats.nsym, stats.psym, toc) ;
tic ;

stats.nnzdiag = nnzdiag ;

if (nargin > 5)
    stats.nzero = nnz (Z) ;

    % recompute the pattern symmetry with Z included
    if (m == n)
	try
	    AZ = A+Z ;
	    if (nnz (AZ) ~= nnz (A) + nnz (Z))
		error ('A and Z overlap!')
	    end
	    [s xmatched pmatched nzoffdiag] = spsym (AZ) ;
	    clear AZ
	    if (nzoffdiag > 0)
		stats.psym = pmatched / nzoffdiag ;
	    else
		stats.psym = 1 ;
	    end
	catch
	    fprintf ('failed to compute symmetry of pattern of A+Z\n') ;
	end
    end

else
    stats.nzero = 0 ;
end

fprintf ('nsym: %g psym: %g time: %g\n', stats.nsym, stats.psym, toc) ;
tic ;

%-------------------------------------------------------------------------------
% intialize ordering statistics
%-------------------------------------------------------------------------------

% if square, Cholesky of C(p,p) where C=A+A', p = amd(C)
stats.amd_lnz = -1 ;	    % nnz (chol (C))
stats.amd_flops = -1 ;	    % flop counts for chol (C)

% if square or rectangular
stats.amd_vnz = -1 ;	    % nnz (V), upper bound on L, for A(:,colamd(A))
stats.amd_rnz = -1 ;	    % nnz (R), upper bound on U, for A(:,colamd(A))

% if square, Cholesky of C(p,p) where C=A+A', p = metis(C)
stats.metis_lnz = -1 ;	    % nnz (chol (C))
stats.metis_flops = -1 ;    % flop counts for chol (C)

% if square or rectangular
stats.metis_vnz = -1 ;	    % nnz (V), upper bound on L, for A(:,metis(A))
stats.metis_rnz = -1 ;	    % nnz (R), upper bound on U, for A(:,metis(A))

% dmperm analysis
stats.nblocks = -1 ;	    % # of blocks in block-triangular form
stats.sprank = -1 ;	    % structural rank
stats.nzoff = -1 ;	    % # of entries of A in off-diagonal blocks

% cs_scc2
stats.ncc = -1 ;	    % # of strongly connected components

% dmperm: best of amd/metis on each square block, best of colamd/metis
% on rectangular blocks
stats.dmperm_lnz = -1 ;	    % nnz (L), for square struct full rank matrices
stats.dmperm_unz = -1 ;	    % nnz (U) + nzoff, for square struct full rank mat
stats.dmperm_flops = -1 ;   % Cholesky flop count of each square block
stats.dmperm_vnz = -1 ;	    % nnz (V), upper bound on L
stats.dmperm_rnz = -1 ;	    % nnz (R), upper bound on U

stats.isND = -1 ;	    % 1 if 2D/3D problem, 0 otherwise

d = max (m,n) ;

% if the matrix has a symmetric nonzero pattern, nzoff will always be zero
if (stats.psym == 1)
    stats.nzoff = 0 ;
end

%-------------------------------------------------------------------------------
% determine if positive definite
%-------------------------------------------------------------------------------

if (~stats.cholcand)

    % not a candidate for Cholesky, so it cannot be positive definite
    stats.posdef = 0 ;

elseif (skip_chol)

    % Cholesky was skipped
    fprintf ('skip Cholesky\n') ;
    stats.posdef = -1 ;

else

    % try chol
    try
	[x, cstats] = cholmod2 (A, ones (stats.ncols,1)) ;
	rcond = cstats (1) ;
	fprintf ('rcond: %g\n', rcond) ;
	stats.posdef = (rcond > 0) ;
    catch
	% chol failed
	disp (lasterr) ;
	fprintf ('sparse Cholesky failed\n') ;
	stats.posdef = -1 ;
    end
    clear x cstats
end

fprintf ('posdef: %d time: %g\n', stats.posdef, toc) ;
tic ;

%-------------------------------------------------------------------------------
% transpose A if m < n, for ordering methods
%-------------------------------------------------------------------------------

if (m < n)
    try
	A = A' ;		    % A is now tall and thin, or square
    catch
	disp (lasterr) ;
	fprintf ('transpose failed...\n') ;
	return ;
    end
    [m n] = size (A) ;
end

if (~isreal (A))
    try
	A = spones (A) ;
    catch
	disp (lasterr) ;
	fprintf ('conversion from complex failed...\n') ;
	return ;
    end
end

fprintf ('computed A transpose if needed, time: %g\n', toc) ;
tic ;

%-------------------------------------------------------------------------------
% order entire matrix with AMD and METIS, if square
%-------------------------------------------------------------------------------

if (m == n)

    tic ;
    try
	if (stats.RBtype (2) == 'u')
	    C = A|A' ;
	else
	    C = A ;
	end
    catch
	disp (lasterr) ;
	fprintf ('A+A'' failed\n') ;
    end
    fprintf ('computed A+A'', time: %g\n', toc) ;

    % order the whole matrix with AMD
    tic ;
    try
	p = amd (C) ;
	c = symbfact (C (p,p)) ;
	stats.amd_lnz = sum (c) ;
	stats.amd_flops = sum (c.^2) ;
    catch
	disp (lasterr) ;
	fprintf ('amd failed\n') ;
    end
    clear p c
    fprintf ('AMD   lnz %d flops %g time: %g\n', ...
	stats.amd_lnz, stats.amd_flops, toc) ;

    % order the whole matrix with METIS
    tic ;
    try
	p = metis (C) ;
	c = symbfact (C (p,p)) ;
	stats.metis_lnz = sum (c) ;
	stats.metis_flops = sum (c.^2) ;
    catch
	disp (lasterr) ;
	fprintf ('metis failed\n') ;
    end
    clear p c C
    fprintf ('METIS lnz %d flops %g time: %g\n', ...
	stats.metis_lnz, stats.metis_flops, toc) ;

else

    % not computed if rectangular
    stats.amd_lnz = -2 ;
    stats.amd_flops = -2 ;
    stats.metis_lnz = -2 ;
    stats.metis_flops = -2 ;

end

%-------------------------------------------------------------------------------
% order entire matrix with COLAMD, for LU bounds
%-------------------------------------------------------------------------------

tic ;
try
    % do not ignore any rows, and do not do etree postordering
    q = colamd2mex (A, [d 10]) ;
    [vnz,rnz] = cs_sqr (A (:,q)) ;
    stats.amd_rnz = rnz ;
    stats.amd_vnz = vnz ;
catch
    disp (lasterr) ;
    fprintf ('colamd2 and cs_sqr failed\n') ;
end
clear q
fprintf ('COLAMD rnz %d vnz %d time: %g\n', stats.amd_rnz, stats.amd_vnz, toc) ;
tic ;

%-------------------------------------------------------------------------------
% order entire matrix with METIS, for LU bounds
%-------------------------------------------------------------------------------

if (~nometis)
    try
	q = metis (A, 'col') ;
	[vnz,rnz] = cs_sqr (A (:,q)) ;
	stats.metis_rnz = rnz ;
	stats.metis_vnz = vnz ;
    catch
	disp (lasterr) ;
	fprintf ('metis(A''*A) and cs_sqr failed\n') ;
    end
end
clear q
fprintf ('METIS  rnz %d vnz %d time: %g\n', ...
    stats.metis_rnz, stats.metis_vnz, toc) ;
tic ;

%-------------------------------------------------------------------------------
% strongly connected components
%-------------------------------------------------------------------------------

try
    % find the # of strongly connected components of the graph of a square A,
    % or # of connected components of the bipartite graph of a rectangular A.
    % [p,q,r,s] = cs_scc2 (A) ;
    if (m == n)
	[p r] = cs_scc (A) ;
    else
	[p r] = cs_scc (spaugment (A)) ;
    end
    stats.ncc = length (r) - 1 ;
    clear p r
catch
    disp (lasterr) ;
    fprintf ('cs_scc failed\n') ;
end

fprintf ('scc %d, time: %g\n', stats.ncc, toc) ;
tic ;

%-------------------------------------------------------------------------------
% isND
%-------------------------------------------------------------------------------

s = 0 ;
if (strfind (kind, 'structural'))
    s = 1 ;
elseif (strfind (kind, 'fluid'))
    s = 1 ;
elseif (strfind (kind, '2D'))
    s = 1 ;
elseif (strfind (kind, 'reduction'))
    s = 1 ;
elseif (strfind (kind, 'electromagnetics'))
    s = 1 ;
elseif (strfind (kind, 'semiconductor'))
    s = 1 ;
elseif (strfind (kind, 'thermal'))
    s = 1 ;
elseif (strfind (kind, 'materials'))
    s = 1 ;
elseif (strfind (kind, 'acoustics'))
    s = 1 ;
elseif (strfind (kind, 'vision'))
    s = 1 ;
elseif (strfind (kind, 'robotics'))
    s = 1 ;
end
stats.isND = s ;

fprintf ('isND %d\n', stats.isND) ;

%-------------------------------------------------------------------------------
% Dulmage-Mendelsohn permutation, and order each block
%-------------------------------------------------------------------------------

if (skip_dmperm)
    fprintf ('skip cs_dmperm, known irreducible\n') ;
else
    try
        % find the Dulmage-Mendelsohn decomposition
        [p,q,r,s,cc,rr] = cs_dmperm (A) ;
        nblocks = length (r) - 1 ;
        stats.nblocks = nblocks ;
        stats.sprank = rr(4)-1 ;
    catch
        disp (lasterr) ;
        fprintf ('cs_dmperm failed\n') ;
    end
end

fprintf ('sprank %d, time: %g\n', stats.sprank, toc) ;
fprintf ('nblocks %d\n', stats.nblocks) ;
tic

ok_square = 1 ;
ok_vnz = 1 ;

try

    mm = diff (r) ;
    nn = diff (s) ;
    square = all (mm == nn) ;

    if (~square)

	% not computed if the matrix is rectangular
	stats.dmperm_lnz = -2 ;
	stats.dmperm_unz = -2 ;
	stats.dmperm_flops = -2 ;

    end

    if (stats.sprank < min (m,n))

	% do not report DMPERM results for structurally singular matrices
	stats.nzoff = -2 ;
	stats.dmperm_lnz = -2 ;
	stats.dmperm_unz = -2 ;
	stats.dmperm_flops = -2 ;
	stats.dmperm_vnz = -2 ;
	stats.dmperm_rnz = -2 ;

    elseif (nblocks == n && m == n)

	% square triangular or diagonal
	C = A (p,q) ;
	clear p q r s

	stats.nzoff = nnz (triu (C, 1)) ;
	stats.dmperm_lnz = n ;
	stats.dmperm_unz = n + stats.nzoff ;
	stats.dmperm_flops = n ;
	stats.dmperm_vnz = n ;
	stats.dmperm_rnz = nnz (C) ;

    elseif (nblocks == 1 && m == n)

	% only one block of structural full rank, so don't redo analysis
	clear p q r s

	stats.nzoff = 0 ;
	if (stats.metis_lnz < 0 || stats.amd_lnz < stats.metis_lnz)
	    stats.dmperm_lnz = stats.amd_lnz ;
	    stats.dmperm_unz = stats.amd_lnz ;
	    stats.dmperm_flops = stats.amd_flops ;
	else
	    stats.dmperm_lnz = stats.metis_lnz ;
	    stats.dmperm_unz = stats.metis_lnz ;
	    stats.dmperm_flops = stats.metis_flops ;
	end

	if (stats.metis_vnz < 0 || stats.amd_rnz < stats.metis_rnz)
	    stats.dmperm_vnz = stats.amd_vnz ;
	    stats.dmperm_rnz = stats.amd_rnz ;
	else
	    stats.dmperm_vnz = stats.metis_vnz ;
	    stats.dmperm_rnz = stats.metis_rnz ;
	end

    else

	% analyze each block of the permuted matrix
	C = A (p,q) ;
	clear p q

	nzoff = nnz (C) ;
	lnz = 0 ;
	unz = 0 ;
	flops = 0 ;
	vnz = 0 ;
	rnz = 0 ;

	for k = 1:nblocks
	    i1 = r (k) ;
	    i2 = r (k+1) - 1 ;
	    j1 = s (k) ;
	    j2 = s (k+1) - 1 ;

	    if (i2-i1 == 1 && j2-j1 == 1)
		% singleton case
		nzoff = nzoff - 1 ;
		unz = unz + 1 ;
		lnz = lnz + 1 ;
		flops = flops + 1 ;
		rnz = rnz + 1 ;
		vnz = vnz + 1 ;
		continue ;
	    end

	    % get the kth block
	    S = C (i1:i2, j1:j2) ;
	    [ms ns] = size (S) ;
	    nzoff = nzoff - nnz (S) ;

	    if (ok_square)
		try
		    if (square)

			% all blocks are square, analyze a square block
			% best of amd and metis
			if (nometis)
			    [pblock c] = analyze (S|S', 'sym', 1) ;
			else
			    [pblock c] = analyze (S|S', 'sym', 3) ;
			end
			lnzblock = sum (c) ;
			unz = unz + lnzblock ;
			lnz = lnz + lnzblock ;
			flops = flops + sum (c.^2) ;
			clear c pblock

		    end
		catch
		    % ordering failed, but keep going to compute nzoff
		    ok_square = 0 ;
		end
	    end

	    if (ok_vnz)
		try

		    % analyze a rectangular block, or LU bounds for square block
		    if (ms < ns)
			S = S' ;
		    end
		    % best of amd and metis
		    try
			if (nometis)
			    pblock = analyze (S, 'col', 1) ;
			else
			    pblock = analyze (S, 'col', 3) ;
			end
		    catch
			pblock = colamd2mex (S, [d 10]) ;
		    end
		    [vnz2,rnz2] = cs_sqr (S (:, pblock)) ;
		    rnz = rnz + rnz2 ;
		    vnz = vnz + vnz2 ;

		catch
		    % ordering failed, but keep going to compute nzoff
		    ok_vnz = 0 ;
		end
	    end

	    clear S pblock

	end

	stats.nzoff = nzoff ;

	if (ok_square)
	    if (~square)
		stats.dmperm_unz = -2 ;
		stats.dmperm_lnz = -2 ;
		stats.dmperm_flops = -2 ;
	    else
		stats.dmperm_lnz = lnz ;
		stats.dmperm_unz = unz + nzoff ;
		stats.dmperm_flops = flops ;
	    end
	end

	if (ok_vnz)
	    stats.dmperm_vnz = vnz ;
	    stats.dmperm_rnz = rnz + nzoff ;
	end

	clear C r s

    end

    if (~ok_square)
	disp (lasterr) ;
	fprintf ('cs_dmperm (square: lnz, unz, flops) ordering failed\n') ;
    end
    if (~ok_vnz)
	disp (lasterr) ;
	fprintf ('cs_dmperm (LU bounds: vnz, rnz) ordering failed\n') ;
    end

catch
    disp (lasterr) ;
    fprintf ('cs_dmperm ordering and nzoff failed (or skipped)\n') ;
end

fprintf ('dmperm stats done, time %g\n', toc) ;
fprintf ('UFstats done\n') ;

author: fieldtrip

repo: fieldtrip


text:
function v = cvxin
superiorto('cell','struct','cvx','cvxtuple');
v = class( struct( 'active', false, 'value', [] ), 'cvxin' );

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function c = commandlines(d)
%tstoolbox/@description/commandlines
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

% C.Merkwirth,U.Parlitz,W.Lauterborn  DPI Goettingen 1998

narginchk(1,1);

c = d.commandlines;

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
classdef AlignedConstrainedCoordSysAxesEnum < matlab.mixin.SetGet
    %AlignedConstrainedCoordSysAxesEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        PosX('+X', [1;0;0], 'X');
        PosY('+Y', [0;1;0], 'Y');
        PosZ('+Z', [0;0;1], 'Z');
        NegX('-X', [-1;0;0], 'X');
        NegY('-Y', [0;-1;0], 'Y');
        NegZ('-Z', [0;0;-1], 'Z');
    end

    properties
        name(1,:) char
        vect(3,1) double
        baseAxis(1,:) char
    end

    methods
        function obj = AlignedConstrainedCoordSysAxesEnum(name, vect, baseAxis)
            obj.name = name;
            obj.vect = vect;
            obj.baseAxis = baseAxis;
        end
    end

    methods(Static)
        function listBoxStr = getListBoxStr()
            m = enumeration('AlignedConstrainedCoordSysAxesEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('AlignedConstrainedCoordSysAxesEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('AlignedConstrainedCoordSysAxesEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% DEMCLASSIFICATIONTWOIVM1 IVM for classification on a data-set sampled from a GP

% IVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'classificationTwo';
experimentNo = 1;

% load data
[X, y] = mapLoadData(dataSetName);

% Set up model
options = ivmOptions;
options.display = 2;
options.numActive = 200;
options.kern = {'rbf', 'white'};

model = ivmCreate(size(X, 1), size(y, 2), X, y, options);

if options.display > 1
  ivm3dPlot(model, 'ivmContour', i);
end
for i = 1:options.extIters;

  % Select the active set.
  model = ivmOptimiseIvm(model, options.display);
  % Plot the data.
  if options.display > 1
    ivm3dPlot(model, 'ivmContour', i);
  end
  % Optimise the kernel parameters.
  model = ivmOptimiseKernel(model, options.display, options.kernIters);
end
model = ivmOptimiseIvm(model, options.display);
if options.display > 1
  ivm3dPlot(model, 'ivmContour', i);
end
% display active points.
model = ivmOptimiseIvm(model, options.display);

% Display the final model.
ivmDisplay(model);

% Save the results.
capName = dataSetName;;
capName(1) = upper(capName(1));
[kern, noise, ivmInfo] = ivmDeconstruct(model);
save(['dem' capName num2str(experimentNo) '.mat'], ...
     'kern', ...
     'noise', ...
     'ivmInfo');

if exist('printDiagram') & printDiagram
  ivmPrintPlot(model, 'ivmContour', [], [], [], capName, experimentNo);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [detrended_data, trend] = strongDeTrend(data, thr_perc, thr_perc_global, n_sigma)
% Returns the detrend per column with a certaint percentile of values in data
% The code uses only the data filtered by the requested percentile to estimate
% an detred column by column
% With this estimation only the data within n_sigma range are used for the
% robust detrend estimation
%
% SYNTAX:
%   smean = strongDeTrend(data, thr_perc, thr_perc_global, n_sigma)
%
% INPUT:
%   data                matrix of values
%   thr_perc            percentile requested [0 1] per column
%   thr_perc_global     percentile requested [0 1] global
%   n_sigma             number of sigmas
%
% OUTPUT:
%   detrended_data   strong detrend per column

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:        Andrea Gatti
%  Contributors:      Andrea Gatti ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   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/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

    if nargin < 2 || isempty(thr_perc)
        thr_perc = 0.99;
    end

    if nargin < 3 || isempty(thr_perc_global)
        thr_perc_global = 0.97;
    end

    if nargin < 4 || isempty(n_sigma)
        n_sigma = 5;
    end

    % filter data
    detrended_data = data;
    data(abs(data) >= perc(abs(data(:)), thr_perc_global)) = nan;
    x_in = (1 : size(data,1))';
    for c = 1 : size(data, 2)
        data(:,c) = data(:,c) - movmean(data(:,c), 5, 'omitnan'); % reduce the signal with moving mean(to remove the trend)
        data((abs(data(:, c)) >= perc(abs(data(:, c)), thr_perc)), c) = nan; % keep thr_perc of data
        thr = n_sigma * std(data(:, c), 'omitnan');
        lid_ok =  abs(data(:,c)) <= thr;
        try
            warning off
            trend = Core_Utils.interp1LS(x_in(lid_ok), detrended_data(lid_ok, c), 1, x_in); % use the current filtered data to compute the trend
            warning on;
        catch
            trend = 0;
        end
        data(:,c) = detrended_data(:,c) - trend;
        thr = n_sigma * std(data(lid_ok, c), 'omitnan'); % repeat the filtering with the new detrended data
        lid_ok =  abs(data(:,c)) <= thr;
        try
            warning off
            trend = Core_Utils.interp1LS(x_in(lid_ok), detrended_data(lid_ok, c), 1, x_in); % recompute trend with the new
            warning on;
        catch
            trend = 0;
        end
        detrended_data(:,c) = detrended_data(:,c) - trend;
    end
end

author: goGPS-Project


text:
function [J, grad] = lrCostFunction(theta, X, y, lambda)
%LRCOSTFUNCTION Compute cost and gradient for logistic regression with
%regularization
%   J = LRCOSTFUNCTION(theta, X, y, lambda) computes the cost of using
%   theta as the parameter for regularized logistic regression and the
%   gradient of the cost w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%

hx = sigmoid(X * theta); %hypothesis, m * 1
% Hint: The computation of the cost function and gradients can be
%       efficiently vectorized. For example, consider the computation
%
%           sigmoid(X * theta)
%
%       Each row of the resulting matrix will contain the value of the
%       prediction for that example. You can make use of this to vectorize
%       the cost function and gradient computations.
%
J = 1 / m * sum(-y' * log(hx) - (1 - y)' * log(1 - hx)) + lambda / (2 * m) * theta(2:end)' * theta(2:end);

% Hint: When computing the gradient of the regularized cost function,
%       there're many possible vectorized solutions, but one solution
%       looks like:
%           grad = (unregularized gradient for logistic regression)
%           temp = theta;
%           temp(1) = 0;   % because we don't add anything for j = 0
%           grad = grad + YOUR_CODE_HERE (using the temp variable)
%
gradf = (1 / m) * (X(:, 1)' * (hx - y));
gradb = (1 / m) * (X(:, 2:end)' * (hx - y)) + lambda * theta(2:end) / m;

grad = [gradf;gradb];
size(grad)

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

grad = grad(:);

end

author: scruel

repo: Notes-ML-AndrewNg


text:
function varargout = panel_coordinates(varargin)
% PANEL_COORDINATES: Create a panel to add/remove/edit scouts attached to a given 3DViz figure.
%
% USAGE:  bstPanelNew = panel_coordinates('CreatePanel')

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

eval(macro_method);
end

%% ===== CREATE PANEL =====
function bstPanelNew = CreatePanel() %#ok<DEFNU>
    panelName = 'Coordinates';
    % Java initializations
    import java.awt.*;
    import javax.swing.*;
    import org.brainstorm.icon.*;
    % CONSTANTS
    TEXT_HEIGHT = java_scaled('value', 20);
    TEXT_WIDTH  = java_scaled('value', 40);
    jFontText = bst_get('Font', 11);
    % Create tools panel
    jPanelNew = gui_component('Panel');

    % ===== CREATE TOOLBAR =====
    jToolbar = gui_component('Toolbar', jPanelNew, BorderLayout.NORTH);
    jToolbar.setPreferredSize(java_scaled('dimension', 100,25));
        % Button "Select vertex"
        jButtonSelect = gui_component('ToolbarToggle', jToolbar, [], 'Select', IconLoader.ICON_SCOUT_NEW, 'Select surface point', @(h,ev)SetSelectionState(ev.getSource.isSelected()));
        % Button "View in MRI Viewer"
        gui_component('ToolbarButton', jToolbar, [], 'View/MRI', IconLoader.ICON_VIEW_SCOUT_IN_MRI, 'View point in MRI Viewer', @ViewInMriViewer);
        % Button "Remove selection"
        gui_component('ToolbarButton', jToolbar, [], 'Del', IconLoader.ICON_DELETE, 'Remove point selection', @RemoveSelection);

    % ===== Main panel =====
    jPanelMain = gui_river();
        % ===== Coordinates =====
        jPanelCoordinates = gui_river('Coordinates (millimeters)');
            % Coordinates
            jPanelCoordinates.add('br', gui_component('label', jPanelCoordinates, 'tab', ' '));
            jLabelX = gui_component('label', jPanelCoordinates, 'tab', '   X');
            jLabelY = gui_component('label', jPanelCoordinates, 'tab', '   Y');
            jLabelZ = gui_component('label', jPanelCoordinates, 'tab', '   Z');
            jLabelX.setHorizontalAlignment(javax.swing.JLabel.CENTER);
            jLabelY.setHorizontalAlignment(javax.swing.JLabel.CENTER);
            jLabelZ.setHorizontalAlignment(javax.swing.JLabel.CENTER);
            jLabelX.setPreferredSize(Dimension(TEXT_WIDTH, TEXT_HEIGHT));
            jLabelY.setPreferredSize(Dimension(TEXT_WIDTH, TEXT_HEIGHT));
            jLabelZ.setPreferredSize(Dimension(TEXT_WIDTH, TEXT_HEIGHT));
            jLabelX.setFont(jFontText);
            jLabelY.setFont(jFontText);
            jLabelZ.setFont(jFontText);
            % === MRI ===
            jPanelCoordinates.add('br', gui_component('label', jPanelCoordinates, 'tab', 'MRI: '));
            jLabelCoordMriX = JLabel('-');
            jLabelCoordMriY = JLabel('-');
            jLabelCoordMriZ = JLabel('-');
            jLabelCoordMriX.setHorizontalAlignment(javax.swing.JLabel.RIGHT);
            jLabelCoordMriY.setHorizontalAlignment(javax.swing.JLabel.RIGHT);
            jLabelCoordMriZ.setHorizontalAlignment(javax.swing.JLabel.RIGHT);
            jLabelCoordMriX.setPreferredSize(Dimension(TEXT_WIDTH, TEXT_HEIGHT));
            jLabelCoordMriY.setPreferredSize(Dimension(TEXT_WIDTH, TEXT_HEIGHT));
            jLabelCoordMriZ.setPreferredSize(Dimension(TEXT_WIDTH, TEXT_HEIGHT));
            jLabelCoordMriX.setFont(jFontText);
            jLabelCoordMriY.setFont(jFontText);
            jLabelCoordMriZ.setFont(jFontText);
            jPanelCoordinates.add('tab', jLabelCoordMriX);
            jPanelCoordinates.add('tab', jLabelCoordMriY);
            jPanelCoordinates.add('tab', jLabelCoordMriZ);
            % === SCS ===
            jPanelCoordinates.add('br', gui_component('label', jPanelCoordinates, 'tab', 'SCS: '));
            jLabelCoordScsX = JLabel('-');
            jLabelCoordScsY = JLabel('-');
            jLabelCoordScsZ = JLabel('-');
            jLabelCoordScsX.setHorizontalAlignment(javax.swing.JLabel.RIGHT);
            jLabelCoordScsY.setHorizontalAlignment(javax.swing.JLabel.RIGHT);
            jLabelCoordScsZ.setHorizontalAlignment(javax.swing.JLabel.RIGHT);
            jLabelCoordScsX.setPreferredSize(Dimension(TEXT_WIDTH, TEXT_HEIGHT));
            jLabelCoordScsY.setPreferredSize(Dimension(TEXT_WIDTH, TEXT_HEIGHT));
            jLabelCoordScsZ.setPreferredSize(Dimension(TEXT_WIDTH, TEXT_HEIGHT));
            jLabelCoordScsX.setFont(jFontText);
            jLabelCoordScsY.setFont(jFontText);
            jLabelCoordScsZ.setFont(jFontText);
            jPanelCoordinates.add('tab', jLabelCoordScsX);
            jPanelCoordinates.add('tab', jLabelCoordScsY);
            jPanelCoordinates.add('tab', jLabelCoordScsZ);
            % === WORLD ===
            jPanelCoordinates.add('br', gui_component('label', jPanelCoordinates, 'tab', 'World: '));
            jLabelCoordWrlX = JLabel('-');
            jLabelCoordWrlY = JLabel('-');
            jLabelCoordWrlZ = JLabel('-');
            jLabelCoordWrlX.setHorizontalAlignment(javax.swing.JLabel.RIGHT);
            jLabelCoordWrlY.setHorizontalAlignment(javax.swing.JLabel.RIGHT);
            jLabelCoordWrlZ.setHorizontalAlignment(javax.swing.JLabel.RIGHT);
            jLabelCoordWrlX.setPreferredSize(Dimension(TEXT_WIDTH, TEXT_HEIGHT));
            jLabelCoordWrlY.setPreferredSize(Dimension(TEXT_WIDTH, TEXT_HEIGHT));
            jLabelCoordWrlZ.setPreferredSize(Dimension(TEXT_WIDTH, TEXT_HEIGHT));
            jLabelCoordWrlX.setFont(jFontText);
            jLabelCoordWrlY.setFont(jFontText);
            jLabelCoordWrlZ.setFont(jFontText);
            jPanelCoordinates.add('tab', jLabelCoordWrlX);
            jPanelCoordinates.add('tab', jLabelCoordWrlY);
            jPanelCoordinates.add('tab', jLabelCoordWrlZ);
            % === MNI ===
            jPanelCoordinates.add('br', gui_component('label', jPanelCoordinates, 'tab', 'MNI: '));
            jLabelCoordMniX = JLabel('-');
            jLabelCoordMniY = JLabel('-');
            jLabelCoordMniZ = JLabel('-');
            jLabelCoordMniX.setHorizontalAlignment(javax.swing.JLabel.RIGHT);
            jLabelCoordMniY.setHorizontalAlignment(javax.swing.JLabel.RIGHT);
            jLabelCoordMniZ.setHorizontalAlignment(javax.swing.JLabel.RIGHT);
            jLabelCoordMniX.setPreferredSize(Dimension(TEXT_WIDTH, TEXT_HEIGHT));
            jLabelCoordMniY.setPreferredSize(Dimension(TEXT_WIDTH, TEXT_HEIGHT));
            jLabelCoordMniZ.setPreferredSize(Dimension(TEXT_WIDTH, TEXT_HEIGHT));
            jLabelCoordMniX.setFont(jFontText);
            jLabelCoordMniY.setFont(jFontText);
            jLabelCoordMniZ.setFont(jFontText);
            jPanelCoordinates.add('tab', jLabelCoordMniX);
            jPanelCoordinates.add('tab', jLabelCoordMniY);
            jPanelCoordinates.add('tab', jLabelCoordMniZ);
            % === VERTEX INDICE ===
            jLabelVertexInd = JLabel('No vertex selected');
            jLabelVertexInd.setFont(jFontText);
            jPanelCoordinates.add('br', jLabelVertexInd);

        jPanelMain.add('hfill', jPanelCoordinates);
    jPanelNew.add(jPanelMain, BorderLayout.CENTER);

    % Create the BstPanel object that is returned by the function
    % => constructor BstPanel(jHandle, panelName, sControls)
    bstPanelNew = BstPanel(panelName, ...
                           jPanelNew, ...
                           struct('jButtonSelect',     jButtonSelect, ...
                                  'jPanelCoordinates', jPanelCoordinates, ...
                                  'jLabelCoordMriX',   jLabelCoordMriX, ...
                                  'jLabelCoordMriY',   jLabelCoordMriY, ...
                                  'jLabelCoordMriZ',   jLabelCoordMriZ, ...
                                  'jLabelCoordScsX',   jLabelCoordScsX, ...
                                  'jLabelCoordScsY',   jLabelCoordScsY, ...
                                  'jLabelCoordScsZ',   jLabelCoordScsZ, ...
                                  'jLabelCoordWrlX',   jLabelCoordWrlX, ...
                                  'jLabelCoordWrlY',   jLabelCoordWrlY, ...
                                  'jLabelCoordWrlZ',   jLabelCoordWrlZ, ...
                                  'jLabelCoordMniX',   jLabelCoordMniX, ...
                                  'jLabelCoordMniY',   jLabelCoordMniY, ...
                                  'jLabelCoordMniZ',   jLabelCoordMniZ, ...
                                  'jLabelVertexInd',   jLabelVertexInd));

end

%% =================================================================================
%  === EXTERNAL PANEL CALLBACKS  ===================================================
%  =================================================================================
%% ===== UPDATE CALLBACK =====
function UpdatePanel()
    % Get panel controls
    ctrl = bst_get('PanelControls', 'Coordinates');
    if isempty(ctrl)
        return
    end
    % Get current figure
    hFig = bst_figures('GetCurrentFigure', '3D');
    % If a figure is available: get if a point select
    if ~isempty(hFig) && ishandle(hFig)
        CoordinatesSelector = getappdata(hFig, 'CoordinatesSelector');
    else
        CoordinatesSelector = [];
    end
    % Update coordinates (text fields)
    % MRI
    if ~isempty(CoordinatesSelector) && ~isempty(CoordinatesSelector.MRI)
        ctrl.jLabelCoordMriX.setText(sprintf('%3.1f', 1000 * CoordinatesSelector.MRI(1)));
        ctrl.jLabelCoordMriY.setText(sprintf('%3.1f', 1000 * CoordinatesSelector.MRI(2)));
        ctrl.jLabelCoordMriZ.setText(sprintf('%3.1f', 1000 * CoordinatesSelector.MRI(3)));
    else
        ctrl.jLabelCoordMriX.setText('-');
        ctrl.jLabelCoordMriY.setText('-');
        ctrl.jLabelCoordMriZ.setText('-');
    end
    % SCS
    if ~isempty(CoordinatesSelector) && ~isempty(CoordinatesSelector.SCS)
        ctrl.jLabelCoordScsX.setText(sprintf('%3.1f', 1000 * CoordinatesSelector.SCS(1)));
        ctrl.jLabelCoordScsY.setText(sprintf('%3.1f', 1000 * CoordinatesSelector.SCS(2)));
        ctrl.jLabelCoordScsZ.setText(sprintf('%3.1f', 1000 * CoordinatesSelector.SCS(3)));
    else
        ctrl.jLabelCoordScsX.setText('-');
        ctrl.jLabelCoordScsY.setText('-');
        ctrl.jLabelCoordScsZ.setText('-');
    end
    % World
    if ~isempty(CoordinatesSelector) && ~isempty(CoordinatesSelector.World)
        ctrl.jLabelCoordWrlX.setText(sprintf('%3.1f', 1000 * CoordinatesSelector.World(1)));
        ctrl.jLabelCoordWrlY.setText(sprintf('%3.1f', 1000 * CoordinatesSelector.World(2)));
        ctrl.jLabelCoordWrlZ.setText(sprintf('%3.1f', 1000 * CoordinatesSelector.World(3)));
    else
        ctrl.jLabelCoordWrlX.setText('-');
        ctrl.jLabelCoordWrlY.setText('-');
        ctrl.jLabelCoordWrlZ.setText('-');
    end
    % MNI
    if ~isempty(CoordinatesSelector) && ~isempty(CoordinatesSelector.MNI)
        ctrl.jLabelCoordMniX.setText(sprintf('%3.1f', 1000 * CoordinatesSelector.MNI(1)));
        ctrl.jLabelCoordMniY.setText(sprintf('%3.1f', 1000 * CoordinatesSelector.MNI(2)));
        ctrl.jLabelCoordMniZ.setText(sprintf('%3.1f', 1000 * CoordinatesSelector.MNI(3)));
    else
        ctrl.jLabelCoordMniX.setText('-');
        ctrl.jLabelCoordMniY.setText('-');
        ctrl.jLabelCoordMniZ.setText('-');
    end
    % Vertex indice
    if ~isempty(CoordinatesSelector) && ~isempty(CoordinatesSelector.iVertex)
        if ~isempty(CoordinatesSelector.Value)
            strValue = ['     Value: ' num2str(CoordinatesSelector.Value)];
        else
            strValue = '';
        end
        ctrl.jLabelVertexInd.setText(sprintf('<HTML>Vertex #%d %s', CoordinatesSelector.iVertex, strValue));
    else
        ctrl.jLabelVertexInd.setText('No vertex selected');
    end
end

%% ===== FOCUS CHANGED ======
function FocusChangedCallback(isFocused) %#ok<DEFNU>
    if ~isFocused
        RemoveSelection();
    end
end

%% ===== CURRENT FIGURE CHANGED =====
function CurrentFigureChanged_Callback() %#ok<DEFNU>
    UpdatePanel();
end

%% ===============================================================================
%  ====== POINTS SELECTION =======================================================
%  ===============================================================================
%% ===== POINT SELECTION : start/stop =====
% Manual selection of a surface point : start(1), or stop(0)
function SetSelectionState(isSelected)
    % Get panel controls
    ctrl = bst_get('PanelControls', 'Coordinates');
    if isempty(ctrl)
        return
    end
    % Get list of all figures
    hFigures = bst_figures('GetAllFigures');
    if isempty(hFigures)
        ctrl.jButtonSelect.setSelected(0);
        return
    end
    % Start selection
    if isSelected
        % Push toolbar "Select" button
        ctrl.jButtonSelect.setSelected(1);
        % Set 3DViz figures in 'SelectingCorticalSpot' mode
        for hFig = hFigures
            % Keep only figures with surfaces
            TessInfo = getappdata(hFig, 'Surface');
            if ~isempty(TessInfo)
                setappdata(hFig, 'isSelectingCoordinates', 1);
                set(hFig, 'Pointer', 'cross');
            end
        end
    % Stop selection
    else
        % Release toolbar "Select" button
        ctrl.jButtonSelect.setSelected(0);
        % Exit 3DViz figures from SelectingCorticalSpot mode
        for hFig = hFigures
            set(hFig, 'Pointer', 'arrow');
            setappdata(hFig, 'isSelectingCoordinates', 0);
        end
    end
end

%% ===== SELECT POINT =====
% Usage : SelectPoint(hFig) : Point location = user click in figure hFIg
function vi = SelectPoint(hFig, AcceptMri) %#ok<DEFNU>
    if (nargin < 2) || isempty(AcceptMri)
        AcceptMri = 1;
    end
    % Get axes handle
    hAxes = findobj(hFig, '-depth', 1, 'Tag', 'Axes3D');
    % Find the closest surface point that was selected
    [TessInfo, iTess, pout, vout, vi, hPatch] = ClickPointInSurface(hFig);
    if isempty(TessInfo)
        return
    end
    Value = [];

    % Get real location of the selected point (in SCS coord, millimeters)
    switch (TessInfo(iTess).Name)
        case 'Anatomy'
            % If MRI is not accepted
            if ~AcceptMri
                vi = [];
                return
            end
            scsLoc  = pout';
            plotLoc = scsLoc;
            iVertex = [];

            % Get MRI
            sMri = bst_memory('GetMri', TessInfo(iTess).SurfaceFile);

        case {'Scalp', 'InnerSkull', 'OuterSkull', 'Cortex', 'Other', 'FEM'}
            sSurf = bst_memory('GetSurface', TessInfo(iTess).SurfaceFile);
            scsLoc = sSurf.Vertices(vi,:);
            plotLoc = vout;
            iVertex = vi;
            % Get value
            if ~isempty(TessInfo(iTess).Data)
                if ~isempty(TessInfo(iTess).DataSource.GridAtlas)
                    Value = [];
                elseif isempty(TessInfo(iTess).DataSource.Atlas)
                    Value = TessInfo(iTess).Data(vi);
                else
                    % Look for the selected point in one of the scouts
                    iRow = panel_scout('GetScoutForVertex', TessInfo(iTess).DataSource.Atlas, vi);
                    % Return value
                    Value = TessInfo(iTess).Data(iRow(1));
                end
            end
            % Get subject
            SubjectFile = getappdata(hFig, 'SubjectFile');
            if ~isempty(SubjectFile)
                sSubject = bst_get('Subject', SubjectFile);
                % == GET MRI ==
                % If subject has a MRI defined
                if ~isempty(sSubject.iAnatomy)
                    % Load MRI
                    MriFile = sSubject.Anatomy(sSubject.iAnatomy).FileName;
                    sMri = bst_memory('LoadMri', MriFile);
                else
                    sMri = [];
                end
            else
                sMri = [];
            end
        otherwise
            scsLoc = vout;
            plotLoc = scsLoc;
            iVertex = [];
            sMri = [];
    end

    % ===== CONVERT TO ALL COORDINATES SYSTEMS =====
    % Save selected point
    CoordinatesSelector.SCS     = scsLoc;
    CoordinatesSelector.MRI     = cs_convert(sMri, 'scs', 'mri', scsLoc);
    CoordinatesSelector.MNI     = cs_convert(sMri, 'scs', 'mni', scsLoc);
    CoordinatesSelector.World   = cs_convert(sMri, 'scs', 'world', scsLoc);
    CoordinatesSelector.iVertex = iVertex;
    CoordinatesSelector.Value   = Value;
    CoordinatesSelector.hPatch  = hPatch;
    setappdata(hFig, 'CoordinatesSelector', CoordinatesSelector);

    % ===== PLOT MARKER =====
    % Remove previous mark
    delete(findobj(hAxes, '-depth', 1, 'Tag', 'ptCoordinates'));
    % Mark new point
    line(plotLoc(1)*1.005, plotLoc(2)*1.005, plotLoc(3)*1.005, ...
         'MarkerFaceColor', [1 1 0], ...
         'MarkerEdgeColor', [1 1 0], ...
         'Marker',          '+',  ...
         'MarkerSize',      12, ...
         'LineWidth',       2, ...
         'Parent',          hAxes, ...
         'Tag',             'ptCoordinates');
    % Update "Coordinates" panel
    UpdatePanel();
end

%% ===== POINT SELECTION: Surface detection =====
function [TessInfo, iTess, pout, vout, vi, hPatch] = ClickPointInSurface(hFig, SurfacesType)
    % Parse inputs
    if (nargin < 2)
        SurfacesType = [];
    end
    iTess = [];
    pout = {};
    vout = {};
    vi = {};
    hPatch = [];

    % === GET SURFACES INFORMATION ===
    % Get axes handle
    hAxes = findobj(hFig, '-depth', 1, 'Tag', 'Axes3D');
    % Get camera position
    CameraPosition = get(hAxes, 'CameraPosition');
    % Get all the surfaces in the figure
    TessInfo = getappdata(hFig, 'Surface');
    if isempty(TessInfo)
        return
    end

    % === CHECK SURFACE TYPE ===
    % Keep only surfaces that are of the required type
    if ~isempty(SurfacesType)
        iAcceptableTess = find(strcmpi({TessInfo.Name}, SurfacesType));
    else
        iAcceptableTess = 1:length(TessInfo);
    end

    % ===== GET SELECTION ON THE CLOSEST SURFACE =====
    % Get the closest point for all the surfaces and patches
    hPatch = [TessInfo(iAcceptableTess).hPatch];
    hPatch = hPatch(ishandle(hPatch));
    patchDist = zeros(1,length(hPatch));
    for i = 1:length(hPatch)
        [pout{i}, vout{i}, vi{i}] = select3d(hPatch(i));
        if ~isempty(pout{i})
            patchDist(i) = norm(pout{i}' - CameraPosition);
        else
            patchDist(i) = Inf;
        end
    end
    if all(isinf(patchDist))
        TessInfo = [];
        pout = [];
        vout = [];
        vi = [];
        return
    end
    % Find closest surface from the camera
    [minDist, iClosestPatch] = min(patchDist);
    % Keep only the point from the closest surface
    hPatch = hPatch(iClosestPatch);
    pout   = pout{iClosestPatch};
    vout   = vout{iClosestPatch};
    vi     = vi{iClosestPatch};

    % Find to which surface this tesselation belongs
    for i = 1:length(TessInfo)
        if any(TessInfo(i).hPatch == hPatch);
            iTess = i;
            break;
        end
    end
end

%% ===== REMOVE SELECTION =====
function RemoveSelection(varargin)
    % Unselect selection button
    SetSelectionState(0);
    % Find all selected points
    hCoord = findobj(0, 'Tag', 'ptCoordinates');
    % Remove coordinates from the figures
    for i = 1:length(hCoord)
        hFig = get(get(hCoord(i), 'Parent'), 'Parent');
        if ~isempty(hFig) && isappdata(hFig, 'CoordinatesSelector')
            rmappdata(hFig, 'CoordinatesSelector');
        end
    end
    % Delete selected points
    delete(hCoord);
    % Update displayed coordinates
    UpdatePanel();
end

%% ===== VIEW IN MRI VIEWER =====
function ViewInMriViewer(varargin)
    global GlobalData;
    % Get current 3D figure
    [hFig,iFig,iDS] = bst_figures('GetCurrentFigure', '3D');
    if isempty(hFig)
        return
    end
    % Get current selected point
    CoordinatesSelector = getappdata(hFig, 'CoordinatesSelector');
    if isempty(CoordinatesSelector) || isempty(CoordinatesSelector.MRI)
        return
    end
    % Get subject and subject's MRI
    sSubject = bst_get('Subject', GlobalData.DataSet(iDS).SubjectFile);
    if isempty(sSubject) || isempty(sSubject.iAnatomy)
        return
    end
    % Display subject's anatomy in MRI Viewer
    hFig = view_mri(sSubject.Anatomy(sSubject.iAnatomy).FileName);
    % Select the required point
    figure_mri('SetLocation', 'mri', hFig, [], CoordinatesSelector.MRI);
end

author: brainstorm-tools


text:
function dw2dx2 = r8poly_lagrange_2 ( npol, xpol, xval )

%*****************************************************************************80
%
%% R8POLY_LAGRANGE_2 evaluates the second derivative of the Lagrange factor.
%
%  Formula:
%
%    W(X)  = Product ( 1 <= I <= NPOL ) ( X - XPOL(I) )
%
%    W'(X) = Sum ( 1 <= J <= NPOL )
%            Product ( I /= J ) ( X - XPOL(I) )
%
%    W"(X) = Sum ( 1 <= K <= NPOL )
%            Sum ( J =/ K )
%            Product ( I /= K, J ) ( X - XPOL(I) )
%
%    For a set of points XPOL(I), 1 <= I <= NPOL, the IPOL-th Lagrange basis
%    polynomial L(IPOL)(X), has the property:
%
%      L(IPOL)( XPOL(J) ) = delta ( IPOL, J )
%
%    and may be expressed as:
%
%      L(IPOL)(X) = W(X) / ( ( X - XPOL(IPOL) ) * W'(XPOL(IPOL)) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 January 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NPOL, the number of abscissas.
%    NPOL must be at least 1.
%
%    Input, real XPOL(NPOL), the abscissas, which should be distinct.
%
%    Input, real XVAL, the point at which the Lagrange factor is to be
%    evaluated.
%
%    Output, real DW2DX2, the second derivative of W with respect to XVAL.
%
  dw2dx2 = 0.0;

  for k = 1 : npol

    for j = 1 : npol

      if ( j ~= k )
        term = 1.0;

        for i = 1: npol
          if ( i ~= j && i ~= k )
            term = term * ( xval - xpol(i) );
          end
        end

        dw2dx2 = dw2dx2 + term;

      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [R, verbose] = tapas_physio_load_other_multiple_regressors(...
    filename_other_regressors, verbose)
% reads other multiple regressors from .txt-file or as variable R from .mat
%
%   R = tapas_physio_load_other_multiple_regressors(filename_other_regressors)
%
% IN
%   filename_other_regressors   filename (.mat or .txt) of other confound
%                               regressors
%
% OUT
%   R                           [Nscans NconfoundRegressors] matrix of multiple regressors
% EXAMPLE
%   tapas_physio_load_other_multiple_regressors
%
%   See also tapas_physio_orthogonalise_physiological_regressors tapas_physio_main_create_regressors

% Author: Lars Kasper
% Created: 2013-02-21
%
% Copyright (C) 2013, Institute for Biomedical Engineering, ETH/Uni Zurich.
%
% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

[fp, fn, fs] = fileparts(filename_other_regressors);

if ~exist(filename_other_regressors, 'file')
    verbose = tapas_physio_log('No input multiple regressors found', verbose, 1);
    R = [];
else
    rp = load(filename_other_regressors);
    if strcmp('.txt', fs) % text-file
        R = rp;
    else % mat file
        R = rp.R;
    end
end
end

author: translationalneuromodeling

repo: tapas


text:
classdef IMMOEA < ALGORITHM
% <multi> <real/integer> <large/none>
% Inverse modeling based multiobjective evolutionary algorithm
% K --- 10 --- Number of reference vectors

%------------------------------- Reference --------------------------------
% R. Cheng, Y. Jin, K. Narukawa, and B. Sendhoff, A multiobjective
% evolutionary algorithm using Gaussian process-based inverse modeling,
% IEEE Transactions on Evolutionary Computation, 2015, 19(6): 838-856.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            K = Algorithm.ParameterSet(10);

            %% Generate random population
            [W,K] = UniformPoint(K,Problem.M);
            W     = fliplr(sortrows(fliplr(W)));
            Problem.N     = ceil(Problem.N/K)*K;
            Population    = Problem.Initialization();
            [~,partition] = max(1-pdist2(Population.objs,W,'cosine'),[],2);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                % Modeling and reproduction
                for k = unique(partition)'
                    Population = [Population,Operator(Problem,Population(partition==k))];
                end
                % Environmental selection
                [~,partition] = max(1-pdist2(Population.objs,W,'cosine'),[],2);
                for k = unique(partition)'
                    current = find(partition==k);
                    if length(current) > Problem.N/K
                        Del = EnvironmentalSelection(Population(current),Problem.N/K);
                        Population(current(Del)) = [];
                        partition(current(Del))  = [];
                    end
                end
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function obj = removeVirtualConstraint(obj, name)
    % Remove virtual constraints defined for the system
    %
    % Parameters:
    % name: the name of the constraint @type cellstr

    assert(ischar(name) || iscellstr(name), ...
        'The name must be a character vector or cellstr.');
    if ischar(name), name = cellstr(name); end

    for i=1:length(name)
        constr = name{i};

        if isfield(obj.VirtualConstraints, constr)
            c_obj = obj.VirtualConstraints.(constr);
            obj = removeParam(obj,c_obj.OutputParamName);
            if ~isempty(c_obj.PhaseParams)
                obj = removeParam(obj,c_obj.PhaseParamName);
            end
            obj.VirtualConstraints = rmfield(obj.VirtualConstraints,constr);
        else
            error('The virtual constraint (%s) does not exist.\n',constr);
        end
    end
end

author: ayonga

repo: frost-dev

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


text:
function a = p03_a ( m, n )

%*****************************************************************************80
%
%% P03_A returns the matrix A for problem 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Cleve Moler,
%    Numerical Computing with MATLAB,
%    SIAM, 2004,
%    ISBN13: 978-0-898716-60-3,
%    LC: QA297.M625,
%    ebook: http://www.mathworks.com/moler/chapters.html
%
%  Parameters:
%
%    Input, integer M, the number of equations.
%
%    Input, integer N, the number of variables.
%
%    Output, real A(M,N), the matrix.
%
  a = [  1.0,  2.0,  3.0; ...
         4.0,  5.0,  6.0; ...
         7.0,  8.0,  9.0; ...
        10.0, 11.0, 12.0; ...
        13.0, 14.0, 15.0 ];

  return
end


text:
% qrtimax() -  perform Quartimax rotation of rows of a data matrix.
%
% Usage: >> [Q,B] = qrtimax(data);
%        >> [Q,B] = qrtimax(data,tol,'[no]reorder');
%
% Inputs:
%        data      - input matrix
%        tol       - the termination tolerance {default: 1e-4}
%        noreorder - rotate without negation/reordering
%
% Outputs:
%        B         - B=Q*A the Quartimax rotation of A
%        Q         - the orthogonal rotation matrix
%
% Author: Sigurd Enghoff, CNL / Salk Institute, 6/18/98

% Copyright (C) Sigurd Enghoff - CNL / Salk Institute, La Jolla 6/18/98
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% Reference: Jack O. Nehaus and Charles Wrigley (1954)
% The Quartimax Method: an analytic approach to orthogonal
% simple structure, Br J Stat Psychol, 7:81-91.

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

function [Q,B] = qrtimax(A,tol,reorder)

if nargin < 1
	help qrtimax
	return
end

DEFAULT_TOL = 1e-4;
MAX_ITERATIONS = 50;

if nargin < 3
	reorder = 1;
elseif isempty(reorder) | reorder == 0
	reorder = 1; % set default
else
	reorder = strcmp('reorder',reorder);
end

if nargin < 2
	eps1 = DEFAULT_TOL;
	eps2 = DEFAULT_TOL;
else
	eps1 = tol;
	eps2 = tol;
end

% Do unto 'Q' what is done to A

Q = eye(size(A,1));

% Compute the cross-products of the rows of the squared loadings,
% i.e. the cost function.
%
%  ---  ---
%  \    \     2    2
%  /    /    f    f
%  ---  ---   ij   ik
%   i   j<k
%
%  See reference, p. 85, line 3

B = tril((A.^2)*(A.^2)');
crit = [sum(sum(B)) - trace(B) , 0];

% Initialize variables

inoim = 0;
iflip = 1;
ict = 0;

% Main iterative loop: keep looping while no two consecutive trials
% satisfy tolerance constraint AND less than MAX_ITERATIONS trials
% AND one or more rotations were performed during last trial.

while inoim < 2 & ict < MAX_ITERATIONS & iflip,
	iflip = 0;

% Run through all combinations of j and k

	for j = 1:size(A,1)-1,
		for k = j+1:size(A,1),
%
%          ---                                   ---
%          \                2     2              \       2     2  2      2   2
% fnum = 4 /    [f   f    (f   - f  )]  , fden = /    [(f   - f  )  - 4 f   f  ]
%          ---    ki  kj    ki    kj             ---     ki    kj        ki  kj
%           k                                     k
%
%  See equation (5)

			u = A(j,:) .^ 2 - A(k,:) .^2;
			v = 2 * A(j,:) .* A(k,:);
			c = sum(u .^ 2 - v .^ 2);
			d = sum(u .* v);

			fden = c;
			fnum = 2 * d;

% Skip rotation if angle is too small

			if abs(fnum) > eps1 * abs(fden)
				iflip = 1;
				angl = atan2(fnum, fden);

% Set angle of rotation according to Table I

				if fnum > 0
					if fden > 0
						angl = .25 * angl;
					else
						angl = .25 * (pi - angl);
					end
				else
					if fden > 0
						angl = .25 * (2 * pi - angl);
					else
						angl = .25 * (pi + angl);
					end
				end

% Perform rotation

				tmp    =  cos(angl) * Q(j,:) + sin(angl) * Q(k,:);
				Q(k,:) = -sin(angl) * Q(j,:) + cos(angl) * Q(k,:);
				Q(j,:) = tmp;

				tmp    =  cos(angl) * A(j,:) + sin(angl) * A(k,:);
				A(k,:) = -sin(angl) * A(j,:) + cos(angl) * A(k,:);
				A(j,:) = tmp;
			end
		end
	end

% Compute cost function.

	B = tril((A.^2)*(A.^2)');
	crit = [sum(sum(B)) - trace(B) , crit(1)];

	inoim = inoim + 1;
	ict = ict + 1;

	fprintf('#%d - crit = %g\n',ict,(crit(1)-crit(2))/crit(1));

% Check relative change of cost function (termination criterion).

	if (crit(1) - crit(2)) / crit(1) > eps2
		inoim = 0;
	end
end

% Reorder and negate if required. Determine new row order based on
% row norms and reorder accordingly. Negate those rows in which the
% accumulated sum is negative.

if reorder
	fprintf('Reordering rows...');
	[fnorm index] = sort(sum(A'.^2));
	Q = Q .* ((2 * (sum(A') > 0) - 1)' * ones(1, size(Q,2)));
	A = A .* ((2 * (sum(A') > 0) - 1)' * ones(1, size(A,2)));
	Q = Q(fliplr(index),:);
	A = A(fliplr(index),:);
	fprintf('\n');
else
	fprintf('Not reordering rows.\n');
end

B=A;

author: PatternRecognition

repo: OpenBMI


text:
%[2012]-"A new fruit fly optimization algorithm: Taking the financial
%distress model as an example"

% (9/12/2020)

function FOA = jFruitFlyOptimizationAlgorithm(feat,label,opts)
% Parameters
lb    = 0;
ub    = 1;
thres = 0.5;

if isfield(opts,'N'), N = opts.N; end
if isfield(opts,'T'), max_Iter = opts.T; end
if isfield(opts,'thres'), thres = opts.thres; end

% Objective function
fun = @jFitnessFunction;
% Number of dimensions
dim = size(feat,2);
% Initial
X = zeros(N,dim);
Y = zeros(N,dim);
for i = 1:N
  for d = 1:dim
    X(i,d) = lb + (ub - lb) * rand();
    Y(i,d) = lb + (ub - lb) * rand();
  end
end
% Compute solution
S = zeros(N,dim);
for i = 1:N
  for d = 1:dim
    % Distance between X and Y axis
    dist = sqrt(X(i,d) ^ 2 + Y(i,d) ^ 2);
    % Solution
    S(i,d) = 1 / dist;
  end
  % Boundary
  SB = S(i,:); SB(SB > ub) = ub; SB(SB < lb) = lb;
  S(i,:) = SB;
end
% Pre
fit   = zeros(1,N);
fitG  = inf;
curve = inf;
t = 1;
% Iterations
while t <= max_Iter
  % Fitness
  for i = 1:N
    % Fitness
    fit(i) = fun(feat,label,(S(i,:) > thres),opts);
    % Update better solution
    if fit(i) < fitG
      fitG = fit(i);
      Xgb  = S(i,:);
      % Update X & Y
      Xb   = X(i,:);
      Yb   = Y(i,:);
    end
  end
	for i = 1:N
    for d = 1:dim
      % Random in [-1,1]
      r1 = -1 + 2 * rand();
      r2 = -1 + 2 * rand();
      % Compute new X & Y
      X(i,d) = Xb(d) + (ub - lb) * r1;
      Y(i,d) = Yb(d) + (ub - lb) * r2;
      % Distance between X and Y axis
      dist   = sqrt((X(i,d) ^ 2) + (Y(i,d) ^ 2));
      % Solution
      S(i,d) = 1 / dist;
    end
    % Boundary
    SB = S(i,:); SB(SB > ub) = ub; SB(SB < lb) = lb;
    S(i,:) = SB;
  end
  curve(t) = fitG;
  fprintf('\nGeneration %d Best (FOA)= %f',t,curve(t))
  t = t + 1;
end
% Select features
Pos   = 1:dim;
Sf    = Pos((Xgb > thres) == 1);
sFeat = feat(:,Sf);
% Store results
FOA.sf = Sf;
FOA.ff = sFeat;
FOA.nf = length(Sf);
FOA.c  = curve;
FOA.f  = feat;
FOA.l  = label;
end

author: JingweiToo

repo: Wrapper-Feature-Selection-Toolbox


text:
function fieldstring=datastructure(structvar, substructvar, filename)
%Display nested structures in 'structvar'; 'substructvar' is used in
%      recursion, must be [] in main call. Result string is saved in
%      file 'filename'
%filename: filename to use for output
%fieldstring: output (with \n linefeeds)
%example call datastructure(saliencyData,[],'saliencyData.txt')

%store format spacing, set compact:
fs=get(0,'FormatSpacing');
format compact

%lines and double lines for display:
lines=repmat('-',1,25);
dlines=repmat('=',1,25);

%get name of variable from call
if exist('substructvar') && ~isempty(substructvar)
    fieldstring = [substructvar '\n' lines '\n'];
else
    fieldstring = [inputname(1) '\n' lines '\n'];
end

if isstruct(structvar)
    %display main structure at current level of recursion:
    fieldstring=[fieldstring ' Struct \n' evalc('disp(structvar)')];
    fieldstring=[fieldstring dlines '\n'];
    %get fields names at current level of recursion:
    fields=fieldnames(structvar);
    for k=1:length(fields)
        if isstruct(eval(['structvar.' fields{k}]))
            %recursive call to get substructures:
            fieldstring=[fieldstring datastructure(eval(['structvar.' fields{k}]),fields{k})];
        end
    end
end
%save result
if exist('filename')
     % in main call
     fid=fopen(filename,'w');
     fprintf(fid, fieldstring);
     fclose(fid);
end

%rstore format spacing
set(0,'FormatSpacing',fs);

author: opencobra

repo: cobratoolbox


text:
function [Anew Bnew] = padwithnan(A, B, dim)
% returns the two input arrays, with the smaller padded to the size of the
% larger in the particular dimension(s) with NaNs
%
% :Usage:
% ::
%
%     [Anew Bnew] = padwithnan(A, B, dim)
%

    if(nargin ~= 3)
        error('padwithnan: incorrect number of arguments');
    end

    diff = size(A,dim) - size(B,dim);

    if diff > 0
        Bnew = pwn_expand(B, dim, diff);
        Anew = A;
    elseif diff < 0
        Anew = pwn_expand(A, dim, diff);
        Bnew = B;
    else
        Anew = A;
        Bnew = B;
    end
end

function outarray = pwn_expand(inarray, dim, diff)
    sznandims = size(inarray);
    sznandims(dim) = abs(diff);
    outarray = cat(dim, inarray, NaN(sznandims));
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [ x, seed ] = triangle_sample ( a, b, c, seed )

%*****************************************************************************80
%
%% TRIANGLE_SAMPLE samples the Triangle PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, C, the parameters of the PDF.
%    A <= B <= C and A < C.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X, a sample of the PDF.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  [ cdf, seed ] = r8_uniform_01 ( seed );

  x = triangle_cdf_inv ( cdf, a, b, c );

  return
end


text:
function out=mp_log2(precision)

if nargin==0
 out=mpLog2(mp(0));
else
 out=mpLog2(mp(precision));
end


text:
function r=ellipsoidalRadius(el,elType,a,f)
%%ELLIPSOIDALRADIUS Find the distance from the center of an axis-aligned
%           oblate spheroid, as is used to approximate the shape of the
%           Earth, to a point on the surface of the spheroid as a function
%           of the (spherical) elevation angle above the equatorial plane
%           or as a function of the ellipsoidal latitude.
%
%INPUTS: el A point or a matrix of the elevation points for which one wants
%           to find radii. The points are in radians. The next parameter
%           sets the type of points.
%    elType An optional parameter specifying the type of elevation el is.
%           Possible values are:
%           0 (The default if omitted or an empty matrix is passed)
%             Geocentric elevations (as in spherical cooridinates).
%           1 Geodetic (ellipsoidal) latitudes.
%         a The semi-major axis of the reference ellipsoid. If this
%           argument is omitted, the value in Constants.WGS84SemiMajorAxis
%           is used.
%         f The flattening factor of the reference ellipsoid. If this
%           argument is omitted, the value in Constants.WGS84Flattening is
%           used.
%
%OUTPUTS: r A matrix of the radii corresponding to the values in el.
%
%The formula for the radius as a function of spherical elevation is taken
%from Section 3.2 of [1].
%
%REFERENCES:
%[1] C. Hirt and M. Rexer, "Earth2014: 1 arc-min shape, topography, bedrock
%    and ice-sheet models - available as gridded data and degree - 10,800
%    spherical harmonics," International Journal of Applied Earth
%    Observation and Geoinformation, vol. 39, pp. 103-112, Jul. 2015.
%
%June 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(elType))
   elType=0;
end

if(nargin<3||isempty(a))
   a=Constants.WGS84SemiMajorAxis;
end

if(nargin<4||isempty(f))
   f=Constants.WGS84Flattening;
end

b=a*(1-f);%Semi-minor axis
e2=(a^2-b^2)/a^2;%eccentricity squared

if(elType==0)
    %We need the square of the sine of the ellipsoidal latitude. However,
    %we only have the spherical elevation. Here, we just substituted the
    %formula from spherLat2EllipsLat into the square of the sine and
    %simplified to get the following formula in terms of the squared
    %tangent.
    spherLat=wrapRange(el,-pi/2,pi/2,true);
    sinEl2=1-1./(1+(1/(1-e2))^2*tan(spherLat).^2);
else
    sinEl2=sin(el).^2;
end

r=a*sqrt((1-e2*(2-e2)*sinEl2)./(1-e2*sinEl2));

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [gKern, gVarmeans, gVarcovars] = rbfard2VardistPsi0Gradient(rbfard2Kern, vardist, covGrad)

% RBFARD2VARDISTPSI0GRADIENT Description

% VARGPLVM
gKern = zeros(1,rbfard2Kern.nParams);
gKern(1) = covGrad*vardist.numData;

gVarmeans = zeros(1,prod(size(vardist.means)));
gVarcovars = zeros(1,prod(size(vardist.means)));

repo: vargplvm


text:
function pass = test_basic_arithmetic(pref)

% TODO: We currently only test scalar autonomous chebops.
%       We also do not check for errors boundary conditions are dealt with.

if ( nargin == 0 )
    pref = cheboppref();
end

dom = [0 3];
u = chebfun(@sin, dom);
A = chebop(@(u) u, dom);
B = chebop(@(u) u, dom);

% Basic plus:
C = A + B;
pass(1) = norm(A(u) + B(u) - C(u)) == 0;

% Basic minus:
C = A - B;
pass(2) = norm(C(u)) == 0;

% Basic uminus:
C = -B;
pass(3) = norm(B(u) + C(u)) == 0;

% Basic times (scalar):
C = 2*A;
pass(4) = norm(2*A(u) - C(u)) == 0;
C = 2.*A;
pass(5) = norm(2*A(u) - C(u)) == 0;

% Basic divide (scalar):
C = A/2;
pass(6) = norm(A(u)/2 - C(u)) == 0;
C = A./2;
pass(7) = norm(A(u)/2 - C(u)) == 0;

% eye
I = eye(A);
pass(8) = norm(u - I(u)) == 0;
pass(9) = norm((A-I)*u - (A(u)-u)) == 0;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% DEMSILHOUETTELINEAR1 Model silhouette data with independent linear models.

% FORMAT
% DESC runs a simple regression on the Agawal and Triggs data.
%
% SEEALSO : demSilhouetteGp1, demSilhouetteAverage
%
% COPYRIGHT : Neil D. Lawrence, 2008

% GP

randn('seed', 1e7)
rand('seed', 1e7)

dataSetName = 'silhouette';
experimentNo = 1;

% load data
[X, y, XTest, yTest] = mapLoadData(dataSetName);

% Set up the model
options = linearOptions;

q = size(X, 2);
d = size(y, 2);
model = linearCreate(q, d, options);

model = linearOptimise(model, X, y);
modelDisplay(model)

% Save results
capName = dataSetName;;
capName(1) = upper(capName(1));
fileBaseName = ['dem' capName 'Linear' num2str(experimentNo)];
save([fileBaseName '.mat'], 'model');

demSilhouettePlot


text:
classdef Builder_DesignVarMonitor_3D < Builder_DesignVarMonitor_Abstract

    methods (Access = protected)

        function cam = createCamera(obj)
            cam = Camera_Rotatory(obj.monitor.axes);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
matlabbatch{1}.spm.tools.physio.save_dir = {'physio_out'};
matlabbatch{1}.spm.tools.physio.log_files.vendor = 'Philips';
matlabbatch{1}.spm.tools.physio.log_files.cardiac = {'SCANPHYSLOG.log'};
matlabbatch{1}.spm.tools.physio.log_files.respiration = {'SCANPHYSLOG.log'};
matlabbatch{1}.spm.tools.physio.log_files.scan_timing = {''};
matlabbatch{1}.spm.tools.physio.log_files.sampling_interval = [];
matlabbatch{1}.spm.tools.physio.log_files.relative_start_acquisition = 0;
matlabbatch{1}.spm.tools.physio.log_files.align_scan = 'last';
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nslices = 32;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.NslicesPerBeat = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.TR = 2;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Ndummies = 0;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nscans = 180;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.onset_slice = 17;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.time_slice_to_slice = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nprep = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sync.gradient_log.grad_direction = 'y';
matlabbatch{1}.spm.tools.physio.scan_timing.sync.gradient_log.zero = 700;
matlabbatch{1}.spm.tools.physio.scan_timing.sync.gradient_log.slice = 1800;
matlabbatch{1}.spm.tools.physio.scan_timing.sync.gradient_log.vol = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sync.gradient_log.vol_spacing = [];
matlabbatch{1}.spm.tools.physio.preproc.cardiac.modality = 'PPU';
matlabbatch{1}.spm.tools.physio.preproc.cardiac.filter.no = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.min = 0.4;
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.file = 'initial_cpulse_kRpeakfile.mat';
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.max_heart_rate_bpm = 90;
matlabbatch{1}.spm.tools.physio.preproc.cardiac.posthoc_cpulse_select.off = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.respiratory.filter.passband = [0.01 2];
matlabbatch{1}.spm.tools.physio.preproc.respiratory.despike = false;
matlabbatch{1}.spm.tools.physio.model.output_multiple_regressors = 'multiple_regressors.txt';
matlabbatch{1}.spm.tools.physio.model.output_physio = 'physio.mat';
matlabbatch{1}.spm.tools.physio.model.orthogonalise = 'none';
matlabbatch{1}.spm.tools.physio.model.censor_unreliable_recording_intervals = false;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.c = 3;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.r = 4;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.cr = 1;
matlabbatch{1}.spm.tools.physio.model.rvt.yes.method = 'hilbert';
matlabbatch{1}.spm.tools.physio.model.rvt.yes.delays = 0;
matlabbatch{1}.spm.tools.physio.model.hrv.yes.delays = 0;
matlabbatch{1}.spm.tools.physio.model.noise_rois.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.movement.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.other.no = struct([]);
matlabbatch{1}.spm.tools.physio.verbose.level = 2;
matlabbatch{1}.spm.tools.physio.verbose.fig_output_file = 'PhysIO_output_level2.jpg';
matlabbatch{1}.spm.tools.physio.verbose.use_tabs = false;


text:
function add1 = cfg_example_add1
% Example script that creates an cfg_exbranch to sum two numbers. The
% inputs are entered as two single numbers, the output is just a single
% number.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_example_add1.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

%% Input Items
% Input a
input1         = cfg_entry; % This is the generic data entry item
input1.name    = 'Input a'; % The displayed name
input1.tag     = 'a';       % The name appearing in the harvested job structure. This name must be unique among all items in the val field of the superior node
input1.strtype = 'e';       % No restriction on what type of data is entered. This could be used to restrict input to real numbers, integers ...
input1.num     = [1 1];     % Number of inputs required (2D-array with exactly one row and one column)
input1.help    = {'This is input a.','This input will be added to the other input.'}; % help text displayed

% Input b
input2         = cfg_entry; % This is the generic data entry item
input2.name    = 'Input b'; % The displayed name
input2.tag     = 'b';       % The name appearing in the harvested job structure. This name must be unique among all items in the val field of the superior node
input2.strtype = 'e';       % No restriction on what type of data is entered. This could be used to restrict input to real numbers, integers ...
input2.num     = [1 1];     % Number of inputs required (2D-array with exactly one row and one column)
input2.help    = {'This is input b.','This input will be added to the other input.'}; % help text displayed

%% Executable Branch
add1      = cfg_exbranch;       % This is the branch that has information about how to run this module
add1.name = 'Add1';             % The display name
add1.tag  = 'cfg_example_add1'; % The name appearing in the harvested job structure. This name must be unique among all items in the val field of the superior node
add1.val  = {input1 input2};    % The items that belong to this branch. All items must be filled before this branch can run or produce virtual outputs
add1.prog = @cfg_example_run_add1;  % A function handle that will be called with the harvested job to run the computation
add1.vout = @cfg_example_vout_add1; % A function handle that will be called with the harvested job to determine virtual outputs
add1.help = {'Add two numbers.'};

%% Local Functions
% The cfg_example_vout_add1 function can go here, it is not useful outside
% the batch environment.
function vout = cfg_example_vout_add1(job)
% Determine what outputs will be present if this job is run. In this case,
% the structure of the inputs is fixed, and the output is always a single
% number. Note that input items may not be numbers, they can also be
% dependencies.

vout = cfg_dep;                        % The dependency object
vout.sname      = 'Add1: a + b';       % Displayed dependency name
vout.src_output = substruct('()',{1}); % The output subscript reference. This could be any reference into the output variable created during computation

author: spm