text:
function test
load trimesh2d;
zfe   = zeros(numel(xfe),1);
nodes = [xfe,yfe,zfe];
triangles = trife;
filename  = 'testSTL.stl';
writeSTL(nodes,triangles,filename);
%% 100k nodes and triangles test(~17MB stl file)
nodesOneM          = rand(100000,3);
trianglesOneM      = ones(100000,3);
trianglesOneM(:,2) = 2;
trianglesOneM(:,3) = 3;
filenameOneM       = 'test100k.stl';
a                  = clock;
writeSTL(nodesOneM,trianglesOneM,filenameOneM);
b                  = clock;
time               = b-a;
    %machine performance
    %http://www.mathworks.com/help/matlab/ref/bench.html(LU 0.54, FFT 0.651,
    %output for my machine:ODE 0.2263, Sparse 0.7735, 2D 1.4288, 3D 1.24)

end


text:
% ATANH  Inverse hyperbolic tangent.
%    ATANH(X) is the inverse hyperbolic tangent of the elements of X.
%
%    See also TANH.
%
%    Reference page in Doc Center
%       doc atanh
%
%    Other functions named atanh
%
%       codistributed/atanh    gpuArray/atanh    sym/atanh    ts/atanh
%

author: jmaih


% =========================================================================
% This code is part of the Matlab-based toolbox
% LagLDDDM - A Lagrangian Gauss--Newton--Krylov Solver for Mass- and
%                        Intensity-Preserving Diffeomorphic Image Registration
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM
%
% Setup NIREP problem (data needs to be obtained
% separately)
%
% =========================================================================

dataset='na05';
example = ['3D-nirep-',dataset];
checkSetupDataFile; if OK, return; end;

% set view options and interpolation options
[viewer,viewOptn] = viewImage('reset','viewImage','imgmontage','colormap','gray(256)','direction','-zyx');

% setup interpolation scheme
imgOptn  = {'imgModel','splineInterMex','regularizer','moments','theta',.01};

% setup transformation used in the parametric part
traOptn  = {'trafo','affine3D'};

% setup distance measure
disOptn  = {'distance','SSD'};

% initialize the regularizer for the non-parametric part
regOptn = {'regularizer','mbHyperElastic','alpha',1,'alphaLength',1,'alphaArea',.1,'alphaVolume',2};

FAIRmessage(mfilename)
load('na01-128x150x128');
load('na01-128x150x128-labels');
load([dataset,'-128x150x128']);
load([dataset,'-128x150x128-labels']);
dataT = double(dataT);
dataR = double(dataR);
dataTl = double(dataTl);
dataRl = double(dataRl);
m = size(dataT);

dataT = 256.*dataT;
dataR = 256.*dataR;
max(dataT(:))

%omega = [0,m(1),0,m(2),0,m(3)];
omega = [0,20,0,23.4375,0,20];
ML = getMultilevel({dataT,dataR},omega,m,'fig',2);
save(outfile,'dataT','dataR','dataTl','dataRl','omega','m','ML');
save(outfile,'-append','viewOptn','imgOptn','traOptn','disOptn','regOptn');
checkSetupDataFile;

% xc       = getCellCenteredGrid(omega,m);
% viewData = @(I) viewImage(imgModel(I,omega,xc),omega,m);

% FAIRfigure(1,'figname',mfilename); clf;
% subplot(1,2,1); viewData(dataT); title('template');
% subplot(1,2,2); viewData(dataR); title('reference');


text:
function [obj, objGrad] = pathObjective(u)
% [obj, objGrad] = pathObjective(u)
%
% Computes the objective function (and gradients) for the simple pendulum
%

obj = u.^2;

if nargout == 2  % Analytic gradients
    nTime = length(u);

    objGrad = zeros(4,nTime); %4 = [time + angle + rate + torque];

    objGrad(4,:) = 2*u;  %gradient obj wrt u

end

end

author: MatthewPeterKelly

repo: OptimTraj


text:
function Y = vanherk(X,N,TYPE,varargin)
%  VANHERK    Fast max/min 1D filter
%
%    Y = VANHERK(X,N,TYPE) performs the 1D max/min filtering of the row
%    vector X using a N-length filter.
%    The filtering type is defined by TYPE = 'max' or 'min'. This function
%    uses the van Herk algorithm for min/max filters that demands only 3
%    min/max calculations per element, independently of the filter size.
%
%    If X is a 2D matrix, each row will be filtered separately.
%
%    Y = VANHERK(...,'col') performs the filtering on the columns of X.
%
%    Y = VANHERK(...,'shape') returns the subset of the filtering specified
%    by 'shape' :
%        'full'  - Returns the full filtering result,
%        'same'  - (default) Returns the central filter area that is the
%                   same size as X,
%        'valid' - Returns only the area where no filter elements are outside
%                  the image.
%
%    X can be uint8 or double. If X is uint8 the processing is quite faster, so
%    dont't use X as double, unless it is really necessary.
%

% Initialization
[direc, shape] = parse_inputs(varargin{:});
if strcmp(direc,'col')
   X = X';
end
if strcmp(TYPE,'max')
   maxfilt = 1;
elseif strcmp(TYPE,'min')
   maxfilt = 0;
else
   error([ 'TYPE must be ' char(39) 'max' char(39) ' or ' char(39) 'min' char(39) '.'])
end

% Correcting X size
fixsize = 0;
addel = 0;
if mod(size(X,2),N) ~= 0
   fixsize = 1;
   addel = N-mod(size(X,2),N);
   if maxfilt
      f = [ X zeros(size(X,1), addel) ];
   else
      f = [X repmat(X(:,end),1,addel)];
   end
else
   f = X;
end
lf = size(f,2);
lx = size(X,2);
clear X

% Declaring aux. mat.
g = f;
h = g;

% Filling g & h (aux. mat.)
ig = 1:N:size(f,2);
ih = ig + N - 1;

g(:,ig) = f(:,ig);
h(:,ih) = f(:,ih);

if maxfilt
   for i = 2 : N
      igold = ig;
      ihold = ih;

      ig = ig + 1;
      ih = ih - 1;

      g(:,ig) = max(f(:,ig),g(:,igold));
      h(:,ih) = max(f(:,ih),h(:,ihold));
   end
else
   for i = 2 : N
      igold = ig;
      ihold = ih;

      ig = ig + 1;
      ih = ih - 1;

      g(:,ig) = min(f(:,ig),g(:,igold));
      h(:,ih) = min(f(:,ih),h(:,ihold));
   end
end
clear f

% Comparing g & h
if strcmp(shape,'full')
   ig = [ N : 1 : lf ];
   ih = [ 1 : 1 : lf-N+1 ];
   if fixsize
      if maxfilt
         Y = [ g(:,1:N-1)  max(g(:,ig), h(:,ih))  h(:,end-N+2:end-addel) ];
      else
         Y = [ g(:,1:N-1)  min(g(:,ig), h(:,ih))  h(:,end-N+2:end-addel) ];
      end
   else
      if maxfilt
         Y = [ g(:,1:N-1)  max(g(:,ig), h(:,ih))  h(:,end-N+2:end) ];
      else
         Y = [ g(:,1:N-1)  min(g(:,ig), h(:,ih))  h(:,end-N+2:end) ];
      end
   end

elseif strcmp(shape,'same')
   if fixsize
      if addel > (N-1)/2
         %disp('hoi')
         ig = [ N : 1 : lf - addel + floor((N-1)/2) ];
         ih = [ 1 : 1 : lf-N+1 - addel + floor((N-1)/2)];
         if maxfilt
            Y = [ g(:,1+ceil((N-1)/2):N-1)  max(g(:,ig), h(:,ih)) ];
         else
            Y = [ g(:,1+ceil((N-1)/2):N-1)  min(g(:,ig), h(:,ih)) ];
         end
      else
         ig = [ N : 1 : lf ];
         ih = [ 1 : 1 : lf-N+1 ];
         if maxfilt
            Y = [ g(:,1+ceil((N-1)/2):N-1)  max(g(:,ig), h(:,ih))  h(:,lf-N+2:lf-N+1+floor((N-1)/2)-addel) ];
         else
            Y = [ g(:,1+ceil((N-1)/2):N-1)  min(g(:,ig), h(:,ih))  h(:,lf-N+2:lf-N+1+floor((N-1)/2)-addel) ];
         end
      end
   else % not fixsize (addel=0, lf=lx)
      ig = [ N : 1 : lx ];
      ih = [ 1 : 1 : lx-N+1 ];
      if maxfilt
         Y = [  g(:,N-ceil((N-1)/2):N-1) max( g(:,ig), h(:,ih) )  h(:,lx-N+2:lx-N+1+floor((N-1)/2)) ];
      else
         Y = [  g(:,N-ceil((N-1)/2):N-1) min( g(:,ig), h(:,ih) )  h(:,lx-N+2:lx-N+1+floor((N-1)/2)) ];
      end
   end

elseif strcmp(shape,'valid')
   ig = [ N : 1 : lx];
   ih = [ 1 : 1: lx-N+1];
   if maxfilt
      Y = [ max( g(:,ig), h(:,ih) ) ];
   else
      Y = [ min( g(:,ig), h(:,ih) ) ];
   end
end

if strcmp(direc,'col')
   Y = Y';
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [direc, shape] = parse_inputs(varargin)
direc = 'lin';
shape = 'same';
flag = [0 0]; % [dir shape]

for i = 1 : nargin
   t = varargin{i};
   if strcmp(t,'col') & flag(1) == 0
      direc = 'col';
      flag(1) = 1;
   elseif strcmp(t,'full') & flag(2) == 0
      shape = 'full';
      flag(2) = 1;
   elseif strcmp(t,'same') & flag(2) == 0
      shape = 'same';
      flag(2) = 1;
   elseif strcmp(t,'valid') & flag(2) == 0
      shape = 'valid';
      flag(2) = 1;
   else
      error(['Too many / Unkown parameter : ' t ])
   end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function definput = arg_wfbtcommon(definput)

% definput.keyvals.J = 1;
% definput.flags.treetype = {'full','dwt'};
% Frequency vs. natral ordering of the output subbands
definput.flags.forder = {'freq','nat'};

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function output = callopticsdp(interfacedata)

% Retrieve needed data
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
x0      = interfacedata.x0;
ub      = interfacedata.ub;
lb      = interfacedata.lb;

model = yalmip2opticsdp(interfacedata);
if options.showprogress;showprogress(['Calling ' interfacedata.solver.tag],options.showprogress);end

if options.savedebug
    save csdpdebug model
end

solvertime = tic;
model.ops = csdpset;
[y,fvals,exitflag,stats,X] = csdp(model.f,model.A,model.b,model.lb,model.ub,model.sdcone,model.y0,model.ops);
solvertime = toc(solvertime);

% Create dual variable in internal format
if options.saveduals
    if any(K.l)
        top = 1;
        D_struc = X{1};
    else
        top = 0;
        D_struc = [];
    end
    if any(K.s)
        for j = 1:length(K.s)
            D_struc = [D_struc;X{j+top}(:)];
        end
    end
else
    D_struc = [];
end

x = y;  % Our notation do not coincide ...
switch exitflag
    case 0
       problem = 0;
    case 1
        problem = 2;
    case 2
        problem = 1;
    case {3,5,6,7,8,9,10}
        problem = 4;
    case {4,-27}
        problem = 3;
    case -50
        problem = 16;
    case 11
        problem = 7;
    otherwise
        problem = -1;
end

if options.savesolveroutput
	solveroutput.y = y;
    solveroutput.fvals = fvals;
    solveroutput.exitflag = exitflag;
    solveroutput.stats = stats;
    solveroutput.X = X;
else
	solveroutput = [];
end

if options.savesolverinput
	solverinput.model = model;
else
	solverinput = [];
end

% Standard interface
output = createOutputStructure(x,D_struc,[],problem,interfacedata.solver.tag,solverinput,solveroutput,solvertime);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function foff=comp_warpedfoff(fc,bw,fs,L,freqtoscale,scaletofreq,do_symmetric)
%COMP_WARPEDFOFF  foff for warped filters

fcwasnegative = fc < 0;

if fcwasnegative && do_symmetric
   fc = -fc;
   fcscale = freqtoscale(fc);
   foff = -floor(scaletofreq(fcscale+.5*bw)/fs*L)+1;
else
   fcscale = freqtoscale(fc);
   foff = floor(scaletofreq(fcscale-.5*bw)/fs*L)+1;
end


text:
function p3 = max(p1)
%PREAL/MAX Overloaded MAX function for class preal.
% Note: Unlike datafun/max preal/max returns the largest single element of
% the array, resulting in one scalar.

global useUnitsFlag

if ~(useUnitsFlag) % If physunits is disabled...
    p3=max(double(p1)); % ... treat as double.
    return
end

p1=preal(p1);
p3=p1(1);
for k=2:numel(p1)
    if p1(k)>p3
        p3=p1(k);
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function res = spm_eeg_specest_mtmspec(S, data, time)
% Plugin for spm_eeg_tf using SPM implementation of multitaper method
% FORMAT res = spm_eeg_specest_mtmspec(S, data, time)
%
% S                     - input structure
% fields of S:
%    S.bandwidth   - time bandwidth parameter determining the degree of
%                    spectral smoothing (typically 3 or 4).
%    S.frequencies - vector of frequencies
%    S.timeres     - time resolution in ms (length of the sliding time-window)
%    S.timestep    - time step (in ms) to slide the time-window by.
%
% Output:
%  res -
%   If no input is provided the plugin returns a cfg branch for itself
%
%   If input is provided:
%      res.fourier - the complex output of wavelet transform (in the case
%                    of single taper)
%      res.pow     - power (in case of multiple tapers, phase is not computed)
%      res.time    - time axis
%      res.freq    - frequency axis
%______________________________________________________________________________________
% Copyright (C) 2010 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak based on the code contributed by Krish Singh
% $Id: spm_eeg_specest_mtmspec.m 4021 2010-07-28 12:43:16Z vladimir $

%-This part if for creating a config branch that plugs into spm_cfg_eeg_tf
% Any parameters can be specified and they are then passed to the plugin
% when it's called.
%--------------------------------------------------------------------------
if nargin == 0
    timeres = cfg_entry;
    timeres.tag = 'timeres';
    timeres.name = 'Time resolution';
    timeres.strtype = 'r';
    timeres.num = [1 1];
    timeres.val = {400};
    timeres.help = {'Length of the sliding time window (in ms)'};

    timestep = cfg_entry;
    timestep.tag = 'timestep';
    timestep.name = 'Time step';
    timestep.strtype = 'r';
    timestep.num = [1 1];
    timestep.val = {50};
    timestep.help = {'Step to slide the time window by (in ms)'};

    bandwidth = cfg_entry;
    bandwidth.tag = 'bandwidth';
    bandwidth.name = 'Time bandwidth';
    bandwidth.strtype = 'n';
    bandwidth.num = [1 1];
    bandwidth.val = {3};
    bandwidth.help = {'Time bandwidth parameter (e.g. 3 or 4)'};

    mtmspec = cfg_branch;
    mtmspec.tag = 'mtmspec';
    mtmspec.name = 'SPM multitaper';
    mtmspec.val = {timeres, timestep, bandwidth};

    res = mtmspec;

    return
elseif nargin < 3
    error('Three input arguments are required');
end

%-Defaults
%--------------------------------------------------------------------------
if ~isfield(S, 'timeres')
    S.timeres = 400;
end

if ~isfield(S, 'timestep')
    S.timestep = 50;
end

if ~isfield(S, 'bandwidth')
    S.bandwidth = 3;
end

%-Data dimensions
%--------------------------------------------------------------------------
fsample = 1./diff(time(1:2));

timeres  = 1e-3*S.timeres;
timestep = 1e-3*S.timestep;

if timestep>timeres
    error('Time resolution should exceed time step');
end

%-Do the spectral analysis
%--------------------------------------------------------------------------
res = [];
[p, f, t] = spm_mmtspec(data', fsample, S.frequencies, timeres, timestep, S.bandwidth);

if size(data, 1) == 1
    res.pow = shiftdim(p, -1);
else
    res.pow = permute(p, [3 1 2]);
end

res.freq = f;

dt = (time(end)-time(1))./length(t);

res.time = (time(1)+dt/2):dt:time(end);

author: spm


text:
function SCKS = spm_SCK(SCKS)
% FORMAT SCKS = spm_SCK(SCKS)
%__________________________________________________________________________
% Square-root Cubature Kalman Filters [2] & Square-root Rauch-Tang-Striebel
% Smoother (SCKF-SCKS [1]).
%==========================================================================
% This function performs joint estimation of the states, input and parameters
% of a model that is described as a stochastic continuous-discrete
% state-space in terms of nonlinear blind deconvolution. The state equations
% must have the form of ordinary differential equations, where the
% discretization is performed through local-linearization scheme [3].
% Additionally, the parameter noise covariance is estimated online via
% stochastic Robbins-Monro approximation method [4], and the measurement noise
% covariance is estimated using a combined variational Bayesian (VB)
% approach with a nonlinear filter/smoother [5].
%__________________________________________________________________________
%
% SCKS.M  - model structure (based on DEM [6] in SPM8 toolbox)
% SCKS.Y  - response variable, output or data
%__________________________________________________________________________
%
% generative model:
%--------------------------------------------------------------------------
%   M(1).f  = dx/dt = f(x,v,P)    {inline function, string or m-file}
%   M(1).g  = y(t)  = g(x,v,P)    {inline function, string or m-file}
%
%   M(1).xP = state error covariance matrix
%   M(1).uP = input error variance
%   M(1).wP = parameter error covariance matrix
%
%   M(1).pE = prior expectation of p model-parameters
%   M(1).pC = prior covariances of p model-parameters
%   M(1).ip = parameter indices
%   M(1).cb = constrain on parameters [lower, upper];
%
%   M(1).Q  = precision components on observation noise
%   M(1).V  = fixed precision (input noise)
%   M(1).W  = precision on state noise (approximated by annealing)
%
%   M(i).m  = number of inputs v(i + 1);
%   M(1).n  = number of states x(i);
%   M(1).l  = number of output v(i);
%
%   M(1).Qf      = form of measurement noise cov estimate:
%                  'auto'[default],'min','mean'
%   M(1).E.nN    = number of SCKF-SCKS algorithm iterations
%   M(1).E.Itol  = tolerance value for SCKF-SCKS convergence
%   M(1).E.nD    = number of integration step between observations
%   M(1).VB.N    = number of VB algorithm iterations
%   M(1).VB.Itol = tolerance value for VB convergence
%   M(1).VB.l    = VB scaling factor;
%
%   conditional moments of model-states - q(u)
%--------------------------------------------------------------------------
%   qU.x  = Conditional expectation of hidden states (backward estimate)
%   qU.v  = Conditional expectation of input (backward estimate)
%   qU.z  = Conditional prediction error
%   qU.S  = Conditional covariance: cov(x) (states - backward estimate)
%   qU.C  = Conditional covariance: cov(u) (input - backward estimate)
%
% conditional moments of model-parameters - q(p)
%--------------------------------------------------------------------------
%   qP.P    = Conditional expectation
%   qP.C    = Conditional covariance
%
%      F    = negative log-likelihood
%__________________________________________________________________________
% Copyright (c) Brno University of Technology (2010)...
% Martin Havlicek 05-12-2010
%
% References:
% [1] Havlicek M et al (2011)
% [2] Arasaratnam, I., Haykin, S. (2009) Cubature Kalman Filters. IEEE
%     Transactions on Automatic Control 54, 1254-1269.
% [3] Jimenez, J.C. (2002) A simple algebraic expression to evaluate the
%     local linearization schemes for stochastic differential equations*
%     1. Applied Mathematics Letters 15, 775-780.
% [4] Van der Merwe, R., 2004. Sigma-point Kalman filters for probabilistic
%     inference in dynamic state-space models. Ph.D.thesis, Oregon Graduate
%     Institute of Science and Technology.
% [5] Sarkka, S., Hartikainen, J. (2011?) Extension of VB-AKF to Estimation
%     of Full Covariance and Non-Linear Systems. In Press.
% [6] Friston, K.J., et al. (2008) DEM: a variational treatment of dynamic
%     systems. Neuroimage 41, 849-885.
%__________________________________________________________________________
% Copyright (C) - Martin Havlicek

% Martin Havlicek
% $Id: spm_SCK.m 4628 2012-01-27 20:51:41Z karl $
% check model specification
%--------------------------------------------------------------------------
M  = SCKS.M;
M  = spm_DEM_M_set(M);

% get integration step dt:
%--------------------------------------------------------------------------
nD = M(1).E.nD;
dt = 1/nD;

% INITIALISATION:
% =========================================================================

% interpolate observation according to integration step
%--------------------------------------------------------------------------
y     = SCKS.Y;            % observations
if size(y,1) > size(y,2)   % check the dimensions
    y = y';
end

% interpolate if dt < 1:
%--------------------------------------------------------------------------
y    = spm_interp(y',1/dt)';
if size(y,1) > size(y,2)   % check dimensions again
    y = y';
end
T    = size(y,2);          % number of time points

% initial condition:
%--------------------------------------------------------------------------
x     = M(1).x;            % states
u     = M(2).v;            % input
pE    = spm_vec(M(1).pE);  % all model parameter
ip    = M(1).ip;           % parameter indices to be estimated
theta = pE(ip);            % selected parameters

try cb  = M(1).cb;               catch, cb = []; end; % parameter constraints
try tE  = spm_vec(SCKS.pP.P{1}); catch, tE = []; end; % true parameters for display (if available)

% covariances (square-roots)
%--------------------------------------------------------------------------
sR      = cell(1,T);
[sR{:}] = deal(sparse(real(chol(inv(M(1).V)))*dt));   % observation noise variance
sQ      = sparse(real(chol(inv(M(1).W)))*dt);         % hidden state noise variance
if ~isempty(M(2).v)
    sV  = sparse(real(chol(inv(M(2).V)))*dt);         % input noise variance
else
    sV  = [];
end

% process error covariances (square-roots)
%--------------------------------------------------------------------------
Sx = sparse(real(chol(M(1).xP))*dt);
if ~isempty(M(2).v)
    Su = sparse(real(chol(M(1).uP))*dt);
else
    Su = [];
end
if ~isempty(ip)
    Sw = sparse(real(chol(M(1).wP(ip,ip)))*dt);
    sW = sparse(real(chol(M(1).pC(ip,ip)))*dt);  % parameter noise variance
    dv = diag(sW);
else
    Sw = [];
    sW = [];
end

% number of states, inputs and parameters:
%--------------------------------------------------------------------------
nx      = size(Sx,1);         % number of states
nu      = size(sV,1);         % number of states
nw      = size(sW,1);         % number of parameters
no      = size(sR{1},1);      % number of observations
noises  = nx + nu + nw + no;  % number of noise components

% concatenate state vector and square-root error covariance:
%--------------------------------------------------------------------------
xc      = [x(:); u(:); theta(:)];
xx      = zeros(nx+nu+nw,T);
xx(:,1) = xc;
Sc      = cell(1,T);
[Sc{:}] = deal(sparse(nx+nu+nw,nx+nu+nw));
Sc{1}   = blkdiag(Sx,Su,Sw);

% get vector indices for components of concatenated state vector
%--------------------------------------------------------------------------
xmask   = [ones(1,nx),ones(1,nu)*2,ones(1,nw)*3,ones(1,no)*4];
xind    = find(xmask==1);
uind    = find(xmask==2);
wind    = find(xmask==3);
clear xmask;

% setting for VB: observation noise estimation:
%--------------------------------------------------------------------------
if ~isempty(M(1).Q)
    try, iter0  = M(1).VB.N;   catch, iter0 = 3;          end
    try, lambda = M(1).VB.l;   catch, lambda = 1-exp(-2); end
    NU      = 6;
    V       = diag(repmat(1e-4,1,no));
    [sR{:}] = deal(sqrt(1./(NU-no-1).*V));
    B       = sqrt(lambda)*eye(no);
    k       = size(sR{1},1);
    iter    = iter0;
    MSE0    = zeros(no,1);
    RR0     = zeros(no,T);
    VBrun   = [];
else
    iter0   = 1;
    iter    = iter0;
    RR      = [];
    VBrun   = [];
end

% Pre-calculate cubature points:
%--------------------------------------------------------------------------
n          = nx + nu + nw + noises;       % total state vector dimension
nPts       = 2*n;                         % number of cubature points
CubPtArray = sqrt(n)*[eye(n) -eye(n)];    % cubature points array

% augment paramter matrix by number of cubature points:
%--------------------------------------------------------------------------
pE   = pE(:,ones(1,nPts));

% prepare matrix template for integration by Local linearization scheme:
%--------------------------------------------------------------------------
EXPm = repmat({[ones(nx),2*ones(nx,1);zeros(1,nx+1)]},1,nPts);
EXPm = sparse(blkdiag(EXPm{:}));
xt   = repmat([zeros(1,nx) 1],1,nPts)';

% =========================================================================
% Iteration scheme:
% =========================================================================
% get maximum number of iterations and tolerance:
%--------------------------------------------------------------------------
try, ItolVB = M(1).VB.Itol;  catch,  ItolVB = 1e-4;      end
try, Itol   = M(1).E.Itol;   catch,  Itol   = 1e-3;      end
try, RUN    = M(1).E.nN;     catch,  RUN    = 32;        end
try, ap     = M(1).E.RM;     catch,  ap     = [1e3 1e6]; end   % Robins-Monro approximation parameters

MLdiff0  = 1e-1;
mloglik0 = 0;
ML       = [];
VBrun    = RUN;
t0       = tic;
% =========================================================================
% Iteration loop (until convergence)
% =========================================================================

for run = 1:RUN
    t1  = tic;
    mloglik = -log(2.*pi).*(T/dt);
    % =====================================================================
    %   Forward pass:
    % =====================================================================
     for t = 2:T

        sQ = diag(diag((1/sqrt(0.9995)-1)*Sc{t-1}(xind,xind)));
        Sa = blkdiag(Sc{t-1},sQ,sV,sW,sR{t-1});
        xa = [xc;zeros(noises,1)];
        Xi =  xa(:,ones(1,nPts)) + Sa*CubPtArray;

        %------------------------------------------------------------------
        % PREDICTION STEP:
        %------------------------------------------------------------------
        xPred(uind,:) = Xi(uind,:) + Xi(uind+nx+nu+nw,:);  % add input noise
        xPred(wind,:) = Xi(wind,:) + Xi(wind+nx+nu+nw,:);  % add parameter noise

        % parameter constraint:
        %------------------------------------------------------------------
        if ~isempty(cb) && ~isempty(ip)
            xPred(wind,:) = min(cb(:,2*ones(1,nPts)),xPred(wind,:)); % upper constrain
            xPred(wind,:) = max(cb(:,1*ones(1,nPts)),xPred(wind,:)); % lower constrain
        end
        pE(ip,:) = xPred(wind,:);

        % propagation of cubature points through nonlinear function:
        %------------------------------------------------------------------
        f             = M(1).f(Xi(xind,:),xPred(uind,:),pE);

        % integration by local-linearization scheme:
        %------------------------------------------------------------------
        dfdx          = spm_diff_all(M(1).f,Xi(xind,:),xPred(uind,:),pE,1);
        dx            = expmall(dfdx,f,dt,EXPm)*xt;
        xPred(xind,:) = Xi(xind,:) + reshape(dx(~xt),nx,nPts) + Xi(xind+nx+nu+nw,:);

        % mean prediction:
        %------------------------------------------------------------------
        x1      = sum(xPred,2)/nPts;
        X       = (xPred-x1(:,ones(1,nPts)))/sqrt(nPts);
        [foo,S] = qr([X]',0);
        S       = S';

        % in the case of VB observation noise estimation:
        %------------------------------------------------------------------
        if ~isempty(M(1).Q) && iter~=1
            NU   = lambda.*(NU-k-1)+k+1;
            V    = B*V*B';
            NU   = NU + 1;
            V0   = V;
        end
        x10    = x1;
        S0     = S;
        yPred0 = M(1).g(xPred(xind,:),xPred(uind,:),pE);

        %------------------------------------------------------------------
        % UPDATE STEP:
        %------------------------------------------------------------------

        % VB estimation of sR (iteratively)
        %------------------------------------------------------------------
        for it = 1:iter

           % VB part - update of square-root measurement noise cov
           %---------------------------------------------------------------
            if ~isempty(M(1).Q) && iter ~= 1
                sR{t} = real(chol(1/(NU - no - 1)*V));
            end

            % propagate cubature points through observation function:
            %--------------------------------------------------------------
            yPred    = yPred0 + sR{t}*CubPtArray(end-no+1:end,:);
            y1       = sum(yPred,2)/nPts;

            Y        = (yPred-y1(:,ones(1,nPts)))/sqrt(nPts);
            [foo,Sy] = qr([Y]',0);
            Sy       = Sy';

            Pxy      = X*Y';                % cross covariance
            K        = (Pxy/Sy')/Sy;        % Kalman gain
            resid    =  y(:,t) - y1;        % innovations

            % state (input,parameter) estimates:
            %--------------------------------------------------------------
            xc = x10 + K*(resid);

            % check parameter constraints:
            %--------------------------------------------------------------
            if ~isempty(cb) && ~isempty(ip)
                xc(wind) = min(cb(:,2),xc(wind)); % upper constrain
                xc(wind) = max(cb(:,1),xc(wind)); % lower constrain
            end
            % estimate of process error covarinace:
            %--------------------------------------------------------------
            [foo,S] = qr([(X - K*Y)]',0);
            S       = S';

            % VB part
            %--------------------------------------------------------------
            if ~isempty(M(1).Q) && iter~=1
                Sa     = blkdiag(S,sQ,sV,sW,sR{t});
                Xi     = repmat([xc;zeros(noises,1)],1,nPts) + Sa*CubPtArray;
                yPreds = M(1).g(Xi(xind,:),Xi(uind,:),pE);  % no additive noise!
                D      = repmat(y(:,t),1,nPts)-yPreds;
                D      = D*D'/nPts;
                V      = V0 + D;
            end

        end
        Sc{t}   = S;
        xx(:,t) = xc;

        % Maximum log-Likelihood
        %------------------------------------------------------------------
         mloglik = mloglik - log(det(Sy*Sy')) - resid'/(Sy*Sy')*resid;

        % Robins-Monro stochastic approximation for of parameters noise cov
        %------------------------------------------------------------------
         if ~isempty(ip)
             subKG = K(wind,:);
             dv    = sqrt((1-1/ap(1))*(dv.^2) + 1/ap(1)*diag(subKG*(subKG*resid*resid')'));
             sW    = diag(dv);
             ap(1) = min(ap(1)+1,ap(2));
         end
         if ~isempty(M(1).Q) && iter~=1
             RR(:,t) = diag(sR{t});
         end
    end
    xxf = xx;
    Sf  = Sc;
    %----------------------------------------------------------------------
    % END of forward pass
    % ---------------------------------------------------------------------

    % =====================================================================
    %   Backward pass:
    % =====================================================================
    for t = T-1:-1:1
        % VB part:
        %------------------------------------------------------------------
        if ~isempty(M(1).Q) && iter~=1
            NU    = lambda.*(NU-k-1)+k+1;
            V     = B*V*B';
            NU    = NU + 1;
            V0    = V;
            sR{t} = real(chol(1/(NU-no-1)*V));
        end

        % Square-root Cubature Rauch-Tung-Striebel smoother
        %------------------------------------------------------------------

        % evaluate cubature points:
        %------------------------------------------------------------------
        Sa = blkdiag(Sc{t},sQ,sV,sW,sR{t});
        xa = [xx(:,t);zeros(noises,1)];
        Xi = xa(:,ones(1,nPts)) + Sa*CubPtArray;

        xPred(uind,:) = Xi(uind,:) + Xi(uind+nx+nu+nw,:);
        xPred(wind,:) = Xi(wind,:) + Xi(wind+nx+nu+nw,:);

        % check parameter constraints:
        %------------------------------------------------------------------
        if ~isempty(cb) && ~isempty(ip)
            xPred(wind,:) = min(cb(:,2*ones(1,nPts)),xPred(wind,:)); % upper constrain
            xPred(wind,:) = max(cb(:,1*ones(1,nPts)),xPred(wind,:)); % lower constrain
        end

        pE(ip,:)      = xPred(wind,:);
        % propagate cubature points through nonlinear function:
        %------------------------------------------------------------------
        f             = M(1).f(Xi(xind,:),xPred(uind,:),pE);
        dfdx          = spm_diff_all(M(1).f,Xi(xind,:),xPred(uind,:),pE,1);
        dx            = expmall(dfdx,f,dt,EXPm)*xt;
        xPred(xind,:) = Xi(xind,:) + reshape(dx(~xt),nx,nPts) + Xi(xind+nx+nu+nw,:) ;

        x1       = sum(xPred,2)/nPts;
        X        = (xPred-x1(:,ones(1,nPts)))/sqrt(nPts) + eps;
        x01      = xx(:,t);
        X01      = (Xi([xind,uind,wind],:) - repmat(x01,1,nPts))/sqrt(nPts);
        [foo,S]  = qr([X]',0);
        S        = S';

        Pxy      = X01*X';      % cross covariance
        K        = (Pxy/S')/S;  % Kalman gain

        % smoothed estimate of the states (input, parameters)
        % and process error covariance:
        %------------------------------------------------------------------
        xx(:,t)  = xx(:,t) + K*(xx(:,t+1) - x1);
        [foo,S]  = qr([X01 - K*X, K*Sc{t+1}]',0);
        S        = S';
        Sc{t}    = S;

        % check parameter constraint:
        %------------------------------------------------------------------
        if ~isempty(cb) && ~isempty(ip)
            xx(wind,t) = min(cb(:,2),xx(wind,t)); % upper constrain
            xx(wind,t) = max(cb(:,1),xx(wind,t));
        end

        % VB part (smoothing):
        %------------------------------------------------------------------
        if ~isempty(M(1).Q) && iter~=1
            Sa     = blkdiag(S,sQ,sV,sW,sR{t});
            Xi     = repmat([xx(:,t);zeros(noises,1)],1,nPts) + Sa*CubPtArray;
            yPreds = M(1).g(Xi(xind,:),Xi(uind,:),pE);  % no additive noise!
            D      = repmat(y(:,t),1,nPts)-yPreds;
            D      = D*D'/nPts;
            V      = V0 + D;
        end

        if ~isempty(M(1).Q) && iter~=1
            RR(:,t) = diag(sR{t});
        end

    end
    xxb = xx;
    Sb  = Sc;
    %----------------------------------------------------------------------
    % END of backward pass
    %----------------------------------------------------------------------

    str{1} = sprintf('SCKS: %i (1:%i)',run,iter);

    % iteration condition for measurement noise estimate:
    % iterate until stabilization of sR estimate
    %----------------------------------------------------------------------
    if ~isempty(M(1).Q) && iter0 ~= 1

        MSE     = mean((RR -(RR0)).^2,2);
        RR0     = RR;
        MSEdiff = abs(MSE - MSE0);
        MSE0    = MSE;

        if MSEdiff < ItolVB  % (till it gets stable)
            switch(lower(M(1).Qf))
                case('all')
                    % take all
                case('mean')
                    sR      = cell(1,T);
                    [sR{:}] = deal(diag(mean(RR,2)));
                case('min')
                    RRs     = sort(RR,2,'descend');
                    sR      = cell(1,T);
                    [sR{:}] = deal(diag(mean(RRs(:,round(T*0.90):end),2)));
                case('auto')
                    dlim    = min(RR,[],2);
                    ulim    = max(RR,[],2);
                    if all(ulim./dlim<4)
                        % take all
                    else
                        RRs     = sort(RR,2,'descend');
                        sR      = cell(1,T);
                        [sR{:}] = deal(diag(mean(RRs(:,round(T*0.90):end),2)));
                    end
            end
            iter0    = 1;
            iter     = iter0;
            mloglik0 = 0;
            VBrun    = run;
        end
    end

    % log-likelihood difference:
    %----------------------------------------------------------------------
    MLdiff(run) = abs(mloglik-mloglik0);
    ML(run)     = mloglik;
    if mloglik > 0
        mloglik     = mloglik - 5000;
        MLdiff(run) = abs(mloglik-mloglik0);
        ML(run)     = mloglik;
    end

    timed  = toc(t1);
    str{2} = sprintf('F:%.4e',ML(end));
    str{3} = sprintf('dF:%.4e',MLdiff(end));
    str{4} = sprintf('(%.2e sec)',timed);
    fprintf('%-16s%-16s%-16s%-16s\n',str{:})

    % plot estimates:
    %----------------------------------------------------------------------
    doplotting(M,xxf,xxb,Sf,Sb,ML,T,wind,ip,run,RR,VBrun);

    % stopping condition:
    %----------------------------------------------------------------------
    if RUN > 1 && (~isempty(ip) || ~isempty(M(1).Q))
        if run == 2
            MLdiff0 = MLdiff(run);
        elseif run > 2
            if MLdiff0 < MLdiff(run),
                MLdiff0 = MLdiff(run);
            end
        end
        if (((MLdiff(run)/MLdiff0)<Itol || run==RUN) || (isempty(ip)&&MLdiff(run)<Itol)) && iter0 == 1,

            yy       = M(1).g(xx(xind,:),xx(uind,:),pE(:,1));
            res      = y - yy;

            try SCKS = rmfield(SCKS,'qU'); end
            try SCKS = rmfield(SCKS,'qP'); end
            try SCKS = rmfield(SCKS,'qH'); end

            % save results in structure:
            %--------------------------------------------------------------
            SCKS.qU.v{1} = yy(:,1:nD:end);
            SCKS.qU.x{1} = xxb(xind,1:nD:end);
            SCKS.qU.v{2} = xxb(uind,1:nD:end);
            SCKS.qU.z{1} = res(:,1:nD:end);
            if ~isempty(ip)
                qP           = SCKS.M(1).pE;
                qP(ip)       = mean(xx(wind,:),2);
                SCKS.qP.P{1} = qP;
                SCKS.qP.P{2} = [];
            end
            SCKS.F       = ML;

            for i = 1:nD:T
                j = 1 + (i - 1)/nD;
                SCKS.qU.S{j} = Sb{i}(xind,xind);
                SCKS.qU.C{j} = Sb{i}(uind,uind);
                if ~isempty(ip)
                    SCKS.qP.p{j} = xxb(wind,:);
                    SCKS.qP.c{j} = Sb{i}(wind,wind);
                    SCKS.qP.ip   = ip;
                end
            end
            return
        end

        mloglik0 = mloglik;
        xc       = [xx([xind,uind],1); mean(xx(wind,:),2)];
        xx(:,1)  = xc;

    else

        pE(ip,1) = mean(xx(wind,:),2);
        yy       = M(1).g(xx(xind,:),xx(uind,:),pE(:,1));
        res      = y - yy;

        try SCKS = rmfield(SCKS,'qU'); end
        try SCKS = rmfield(SCKS,'qP'); end
        try SCKS = rmfield(SCKS,'qH'); end

        % save results in structure:
        %------------------------------------------------------------------
        SCKS.qU.v{1} = yy(:,1:nD:end);
        SCKS.qU.x{1} = xxb(xind,1:nD:end);
        SCKS.qU.v{2} = xxb(uind,1:nD:end);
        SCKS.qU.z{1} = res(:,1:nD:end);
        if ~isempty(ip)
            qP           = SCKS.M(1).pE;
            qP(ip)       = mean(xx(wind,:),2);
            SCKS.qP.P{1} = qP;
            SCKS.qP.P{2} = [];
        end
        SCKS.F       = ML;

        for i = 1:nD:T
            j = 1 + (i - 1)/nD;
            SCKS.qU.S{j} = Sb{i}(xind,xind);
            SCKS.qU.C{j} = Sb{i}(uind,uind);
            if ~isempty(ip)
                SCKS.qP.p{j} = xxb(wind,:);
                SCKS.qP.c{j} = Sb{i}(wind,wind);
                SCKS.qP.ip   = ip;
            end
        end
        return
    end
end

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

%--------------------------------------------------------------------------
% Plot estimates at each iteration:
%--------------------------------------------------------------------------
function doplotting(M,xxf,xxb,Sf,Sb,ML,T,wind,ip,run,RR,VBrun)

% Initialize display:
%--------------------------------------------------------------------------
spm_figure('GetWin','SCKF-SCKS estimates');
set(gcf,'Renderer','painter'); clf

% Hidden states
%--------------------------------------------------------------------------
for p=1:2
    subplot(3,3,[1:3]+3*(p-1)),
    hax = gca;
    si    = spm_invNcdf(1 - 0.05);
    s     = [];
    if p == 1,
        xxfig = xxf;
        Sfig  = Sf;
        tit   = 'SCKF - forward pass';
    else
        xxfig = xxb;
        Sfig  = Sb;
        tit   = 'SCKS - backward pass';
    end
    for i = 1:T
        s = [s abs(diag(Sfig{i}))];
    end

    % conditional covariances
    %----------------------------------------------------------------------
    j           = [1:size(xxfig(:,:),1)];
    ss          = si*s(j,:);
    [ill indss] = sort(full(mean(ss,2)),'descend');

    pf = plot(1:T,xxfig,'linewidth',1.5);
    set(hax,'xlim',[1,T],'nextplot','add')
    box(hax,'on')
    for ic = 1:size(xxfig,1)
        col0 = get(pf(indss(ic)),'color');
        col  = (ones(1,3)-col0)*0.65 + col0;
        fill([(1:T) fliplr(1:T)],[(xxfig(indss(ic),:) + ss(indss(ic),:)) fliplr((xxfig(indss(ic),:) - ss(indss(ic),:)))],...
            'r',...
            'FaceColor',col,...
            'EdgeColor',col);
        hold on;
        COL{ic} = col0;
    end
    for ic = 1:size(xxfig,1)
        plot(xxfig(indss(ic),:),'color',COL{ic},'linewidth',0.75);
    end
    title(tit);
    grid(hax,'on');
    axis(hax,'tight');
end

% Parameter estimates
%--------------------------------------------------------------------------
subplot(3,3,7)
h = plot([1:length(ML)],ML);
if ~isempty(M(1).Q)
    AYlim = get(gca,'Ylim');
    if VBrun>=run
        bkg = ones(1,run)*max(AYlim);
    else
        bkg = [ones(1,VBrun)*max(AYlim),ones(1,abs(VBrun-run))*min(AYlim)];
    end
    a = area(bkg,min(AYlim));
    axis([1 run+1 AYlim(1) AYlim(2)]); hold on;
    set(a(1),'FaceColor',ones(1,3)*0.6,'EdgeColor',ones(1,3)*0.6)
    h = plot([1:length(ML)],ML);
    axis([1 run+1 AYlim(1) AYlim(2)]); hold on;
    if VBrun>=run
        text(run/2+1,mean(AYlim),'VB-SCKS','HorizontalAlignment','center','VerticalAlignment','top');
    else
        text(VBrun/2+1,mean(AYlim),'VB-SCKS','HorizontalAlignment','center','VerticalAlignment','top');
        text(VBrun+(abs(VBrun-run)/2)+1,mean(AYlim),'SCKS','HorizontalAlignment','center','VerticalAlignment','top');
    end
    set(gca,'Layer','top')
    hold off
end
set(h,'color','k','Marker','o','Markersize',4,'MarkerFaceColor','k','linewidth',1);
title('Log-Likelihood');

if ~isempty(ip)
    subplot(3,3,8)
    b1 = bar(ip,mean(xxb(wind,:),2)','FaceColor','k');
    set(b1,'BarWidth',0.5); hold on;
    title('Parameters');
end

% Standard deviation of residuals
%--------------------------------------------------------------------------
if ~isempty(M(1).Q)
    subplot(3,3,9)
    if VBrun >= run
        plot(RR'); hold on;
    else
        plot(RR'); hold on;
        switch(lower(M(1).Qf))
            case('all')
                sR  = RR';
            case('mean')
                sR  = repmat(mean(RR,2),1,T)';
            case('min')
                RRs = sort(RR,2,'descend');
                sR  = repmat(mean(RRs(:,round(T*0.90):end),2),1,T)';
            case('auto')
                dlim    = min(RR,[],2);
                ulim    = max(RR,[],2);
                if all(ulim./dlim<4)
                    sR  = RR';
                else
                    RRs  = sort(RR,2,'descend');
                    sR   = repmat(mean(RRs(:,round(T*0.90):end),2),1,T)';
                end
        end
        plot([1:T],sR,'r','linewidth',2);
    end
    axis(gca,'tight');
    title('Standard deviation of residuals');
    hold off;
end
drawnow


text:
function [fStdDevB, fStdDevMc, fBValue, fMc, vBValues] = calc_BootstrapB(mCatalog, nNumberRuns, nMinNum, nCalculateMc, fBinning)
    % Computes standard deviation of b-value and Mc by bootstrapping the dataset.
    %
    % [fStdDevB, fStdDevMc, fBValue, fMc, vBValues] = calc_BootstrapB(mCatalog, nNumberRuns, nMinNum, nCalculateMc, fBinning)
    %
    %
    % Input parameters:
    %   mCatalog          Earthquake catalog to be used
    %   nNumberRuns       Number of simulation runs (Bootstrap)
    %   nMinNum           Minimum number of events > Mc for computing a b-value (after bottstrapping sample)
    %   nCalculateMC      Method to determine the magnitude of completeness (see also: help calc_Mc)
    %   fBinning          Magnitude binning of the catalog
    %
    % Output parameters:
    %   fStdDevB          Standard deviation of the computed b-values as the second moment of the b-value distribution
    %   fStdDevMc         Standard deviation of the computed Mc as the second moment of the Mc distribution
    %   fBValue           Mean b-value of bootstrapped samples
    %   fMc               Mean Mc of bootstrapped samples
    %   vBValues          Vector of all bootstrapped b-values
    %
    % Danijel Schorlemmer
    % June 18, 2003

    report_this_filefun();

    % Get number of events in catalog
    nLength = mCatalog.Count;
    % Init container
    mResult = [];
    % Bootstrap loop
    for nRuns = 1:nNumberRuns
        % Iniy the bootstrapped catalog
        mLoopCatalog = [];
        % Get the random selection of events (multiples allowed)
        vRnd = ceil(rand(nLength,1) * nLength);
        % Create the bootstrapped catalog
        mLoopCatalog = mCatalog.subset(vRnd);
        % Reduce bootstrapped catalog to all events with M >= Mc
        fMc = calc_Mc(mLoopCatalog, McMethods.McBestCombo, fBinning);
        vSel = mLoopCatalog(:,6) >= fMc;
        mLoopCatalog = mLoopCatalog(vSel,:);
        % If enough events remain, compute b-value
        if length(mLoopCatalog(:,1)) >= nMinNum
            % Calculate b-value from bootstrapped catalog
            [fBValue] =  calc_bmemag(mLoopCatalog, fBinning);
        else
            % Not enough events available
            fBValue = nan;
        end
        % Store the results
        mResult = [mResult; fBValue fMc];
    end
    % Return values
    fBValue = mean(mResult(:,1), 'omitnan');
    fMc = mean(mResult(:,2), 'omitnan');
    % Compute the standard deviation of Mc as the second moment of the Mc distribution
    vSel = ~isnan(mResult(:,2));
    vDist = mResult(vSel,2);
    fStdDevMc = std(vDist,1,'omitnan');
    % Compute the standard deviation of b as the second moment of the b distribution
    vSel = ~isnan(mResult(:,1));
    vBValues = mResult(vSel,1);
    fStdDevB = std(vBValues,1,'omitnan');
end

author: CelsoReyes


text:
function [ y, next ] = p17_equil ( neqn, next )

%*****************************************************************************80
%
%% P17_EQUIL returns equilibrium solutions of problem p17.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, integer NEXT, the index of the previous
%    equilibrium, which should be 0 on first call.
%
%    Output, real Y(NEQN), the "next" equilibrium solution, if any.
%
%    Output, integer NEXT, the index of the current equilibrium,
%    or 0 if there are no more.
%
  if ( next == 0 )
    next = 1;
    y(1:neqn,1) = 0.0;
  else
    next = 0;
    y = [];
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function X = tenrand(varargin)
%TENRAND Uniformly distributed pseudo-random tensor.
%
%   X = TENRAND(SZ) forms a tensor of size SZ with pseudo-random
%   values drawn from a uniform distribution on the unit interval.
%
%   TENRAND(SZ) is equivalent to TENSOR(RAND(SZ(1),SZ(2),...),SZ).
%
%   See also TENSOR, SPTENRAND, RAND.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

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

if nargin == 1
    sz = varargin{1};
else
    sz = cell2mat(varargin);
end

data = rand([sz 1 1]);
X = tensor(data,sz);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function [ fea, out ] = ex_laplace2( varargin )
%EX_LAPLACE2 2D Laplace equation on a circle with nonzero boundary conditions.
%
%   [ FEA, OUT ] = EX_LAPLACE2( VARARGIN ) Laplace equation on a circle with
%   source term 0, nonzero boundary conditions and exact solution u=r^3*sin(3*th).
%   Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       radi        scalar {1}             Radius of circle
%       hmax        scalar {0.1}           Max grid cell size
%       sfun        string {sflag1}        Shape function
%       iphys       scalar 0/{1}           Use physics mode to define problem    (=1)
%                                          or directly define fea.eqn/bdr fields (=0)
%       iplot       scalar 0/{1}           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { ...
  'radi',     1; ...
  'hmax',     0.1; ...
  'refsol',   'hypot(x,y)^3*sin(3*atan2(y,x))'; ...
  'sfun',     'sflag1'; ...
  'icub',     2; ...
  'iphys',    1; ...
  'iplot',    1; ...
  'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

% Geometry definition.
gobj = gobj_circle( [0 0], opt.radi );
fea.geom.objects = { gobj };

% Grid generation.
fea.grid = gridgen(fea,'hmax',opt.hmax,'fid',fid,'dprim',false);
n_bdr    = max(fea.grid.b(3,:));        % Number of boundaries.

% Problem definition.
fea.sdim  = { 'x' 'y' };                % Coordinate names.
if ( opt.iphys==1 )

  fea = addphys(fea,@poisson);          % Add Poisson equation physics mode.
  fea.phys.poi.sfun = { opt.sfun };     % Set shape function.
  fea.phys.poi.eqn.coef{3,4} = { 0 };   % Set source term to zero.
  fea.phys.poi.bdr.coef{1,end} = repmat({opt.refsol},1,n_bdr);
  fea = parsephys(fea);                 % Check and parse physics modes.

else

  fea.dvar  = { 'u' };                  % Dependent variable name.
  fea.sfun  = { opt.sfun  };            % Shape function.

  % Define equation system.
  fea.eqn.a.form = { [2 3;2 3] };       % First row indicates test function space   (2=x-derivative + 3=y-derivative),
                                        % second row indicates trial function space (2=x-derivative + 3=y-derivative).
  fea.eqn.a.coef = { 1 };               % Coefficient used in assembling stiffness matrix.

  fea.eqn.f.form = { 1 };               % Test function space to evaluate in right hand side (1=function values).
  fea.eqn.f.coef = { 0 };               % Coefficient used in right hand side.

  % Define boundary conditions.
  fea.bdr.d     = cell(1,n_bdr);
 [fea.bdr.d{:}] = deal(opt.refsol);     % Assign reference solution to all boundaries (Dirichlet).

  fea.bdr.n     = cell(1,n_bdr);        % No Neumann boundaries ('fea.bdr.n' empty).

end

% Parse and solve problem.
fea       = parseprob(fea);             % Check and parse problem struct.
fea.sol.u = solvestat(fea,'fid',fid,'icub',opt.icub);   % Call to stationary solver.

% Postprocessing.
s_err = ['abs(',opt.refsol,'-u)'];
if ( opt.iplot>0 )
  figure
  subplot(3,1,1)
  postplot(fea,'surfexpr','u','axequal','on')
  title('Solution u')
  subplot(3,1,2)
  postplot(fea,'surfexpr',opt.refsol,'axequal','on')
  title('Exact solution')
  subplot(3,1,3)
  postplot(fea,'surfexpr',s_err,'axequal','on')
  title('Error')
end

% Error checking.
if ( size(fea.grid.c,1)==4 )
  xi = [0;0];
else
  xi = [1/3;1/3;1/3];
end
err = evalexpr0(s_err,xi,1,1:size(fea.grid.c,2),[],fea);
ref = evalexpr0('u',xi,1,1:size(fea.grid.c,2),[],fea);
err = sqrt(sum(err.^2)/sum(ref.^2));

if( ~isempty(fid) )
  fprintf(fid,'\nL2 Error: %f\n',err)
  fprintf(fid,'\n\n')
end

out.err  = err;
out.pass = out.err<5e-3;
if ( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
% This is the main extended Kalman filter SLAM loop. This script calls all the required
% functions in the correct order.
%
% You can disable the plotting or change the number of steps the filter
% runs for to ease the debugging. You should however not change the order
% or calls of any of the other lines, as it might break the framework.
%
% If you are unsure about the input and return values of functions you
% should read their documentation which tells you the expected dimensions.

% Turn off pagination:
more off;

% clear all variables and close all windows
clear all;
close all;

% Make tools available
addpath('tools');

% Read world data, i.e. landmarks. The true landmark positions are not given to the robot
landmarks = read_world('../data/world.dat');
% load landmarks;
% Read sensor readings, i.e. odometry and range-bearing sensor
data = read_data('../data/sensor_data.dat');
%load data;

INF = 1000;
% Get the number of landmarks in the map
N = size(landmarks,2);

% observedLandmarks is a vector that keeps track of which landmarks have been observed so far.
% observedLandmarks(i) will be true if the landmark with id = i has been observed at some point by the robot
observedLandmarks = repmat(false,1,N);

% Initialize belief:
% mu: 2N+3x1 vector representing the mean of the normal distribution
% The first 3 components of mu correspond to the pose of the robot,
% and the landmark poses (xi, yi) are stacked in ascending id order.
% sigma: (2N+3)x(2N+3) covariance matrix of the normal distribution
mu = repmat([0.0], (2*N+3), 1);
robSigma = zeros(3);
robMapSigma = zeros(3,2*N);
mapSigma = INF*eye(2*N);
sigma = [[robSigma robMapSigma];[robMapSigma' mapSigma]];

% toogle the visualization type
%showGui = true;  % show a window while the algorithm runs
showGui = false; % plot to files instead

% Perform filter update for each odometry-observation pair read from the
% data file.
for t = 1:size(data.timestep, 2)
%for t = 1:80

    % Perform the prediction step of the EKF
    [mu, sigma] = prediction_step(mu, sigma, data.timestep(t).odometry);

    % Perform the correction step of the EKF
    [mu, sigma, observedLandmarks] = correction_step(mu, sigma, data.timestep(t).sensor, observedLandmarks);

    %Generate visualization plots of the current state of the filter
    plot_state(mu, sigma, landmarks, t, observedLandmarks, data.timestep(t).sensor, showGui);
    disp("Current state vector:")
    disp("mu = "), disp(mu)
endfor

disp("Final system covariance matrix:"), disp(sigma)
% Display the final state estimate
disp("Final robot pose:")
disp("mu_robot = "), disp(mu(1:3)), disp("sigma_robot = "), disp(sigma(1:3,1:3))

author: kiran-mohan

repo: SLAM-Algorithms-Octave


text:
function [out] = melt_2(p1,S,dt)
%melt_2

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Snowmelt at a constant rate
% Constraints:  f <= S/dt
% @(Inputs):    p1   - melt rate [mm/d]
%               S    - current storage [mm]
%               dt   - time step size [d]

out = min(p1,S/dt);

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function [K] = ku0u0(x, xp, hyp, ubar, ubarp, dt, i)

logsigma = hyp(1);
logtheta = hyp(2);

a1 = hyp(3);
a2 = hyp(4);
a3 = hyp(5);
a4 = hyp(6);

n_x = size(x,1);
n_xp = size(xp,1);

x = repmat(x,1,n_xp);
xp = repmat(xp',n_x,1);

ubar = repmat(ubar,1,n_xp);
ubarp = repmat(ubarp',n_x,1);

switch i

case 0

K=exp(1).^(logsigma+(-8).*logtheta+(-1/2).*exp(1).^((-1).*logtheta).*(x+( ...
  -1).*xp).^2).*(exp(1).^(2.*logtheta).*(exp(1).^(2.*logtheta).*(exp(1).^( ...
  2.*logtheta).*(exp(1).^(2.*logtheta)+a1.*dt.*exp(1).^logtheta.*(a1.*dt.* ...
  ubar.*ubarp+(-1).*(ubar+(-1).*ubarp).*(x+(-1).*xp))+(-1).*a1.^2.*dt.^2.* ...
  ubar.*ubarp.*(x+(-1).*xp).^2)+a2.*dt.*exp(1).^logtheta.*((-2).*exp(1).^( ...
  2.*logtheta)+exp(1).^logtheta.*(3.*a1.*dt.*(ubar+(-1).*ubarp)+2.*(x+(-1) ...
  .*xp)).*(x+(-1).*xp)+(-1).*a1.*dt.*(ubar+(-1).*ubarp).*(x+(-1).*xp).^3)+ ...
  a2.^2.*dt.^2.*(3.*exp(1).^(2.*logtheta)+(-6).*exp(1).^logtheta.*(x+(-1) ...
  .*xp).^2+(x+(-1).*xp).^4))+(-1).*a1.*a3.*dt.^2.*exp(1).^(2.*logtheta).*( ...
  ubar+ubarp).*(3.*exp(1).^(2.*logtheta)+(-6).*exp(1).^logtheta.*(x+(-1).* ...
  xp).^2+(x+(-1).*xp).^4)+a3.^2.*dt.^2.*(15.*exp(1).^(3.*logtheta)+(-45).* ...
  exp(1).^(2.*logtheta).*(x+(-1).*xp).^2+15.*exp(1).^logtheta.*(x+(-1).* ...
  xp).^4+(-1).*(x+(-1).*xp).^6))+a4.*dt.*exp(1).^(2.*logtheta).*((-2).* ...
  a2.*dt.*(15.*exp(1).^(3.*logtheta)+(-45).*exp(1).^(2.*logtheta).*(x+(-1) ...
  .*xp).^2+15.*exp(1).^logtheta.*(x+(-1).*xp).^4+(-1).*(x+(-1).*xp).^6)+ ...
  exp(1).^logtheta.*(6.*exp(1).^(3.*logtheta)+(-3).*exp(1).^(2.*logtheta) ...
  .*(5.*a1.*dt.*(ubar+(-1).*ubarp)+4.*(x+(-1).*xp)).*(x+(-1).*xp)+(-1).* ...
  a1.*dt.*(ubar+(-1).*ubarp).*(x+(-1).*xp).^5+2.*exp(1).^logtheta.*(x+(-1) ...
  .*xp).^3.*(5.*a1.*dt.*(ubar+(-1).*ubarp)+x+(-1).*xp)))+a4.^2.*dt.^2.*( ...
  105.*exp(1).^(4.*logtheta)+(-420).*exp(1).^(3.*logtheta).*(x+(-1).*xp) ...
  .^2+210.*exp(1).^(2.*logtheta).*(x+(-1).*xp).^4+(-28).*exp(1) ...
  .^logtheta.*(x+(-1).*xp).^6+(x+(-1).*xp).^8));

case 1 % logsigma

K=exp(1).^(logsigma+(-8).*logtheta+(-1/2).*exp(1).^((-1).*logtheta).*(x+( ...
  -1).*xp).^2).*(exp(1).^(2.*logtheta).*(exp(1).^(2.*logtheta).*(exp(1).^( ...
  2.*logtheta).*(exp(1).^(2.*logtheta)+a1.*dt.*exp(1).^logtheta.*(a1.*dt.* ...
  ubar.*ubarp+(-1).*(ubar+(-1).*ubarp).*(x+(-1).*xp))+(-1).*a1.^2.*dt.^2.* ...
  ubar.*ubarp.*(x+(-1).*xp).^2)+a2.*dt.*exp(1).^logtheta.*((-2).*exp(1).^( ...
  2.*logtheta)+exp(1).^logtheta.*(3.*a1.*dt.*(ubar+(-1).*ubarp)+2.*(x+(-1) ...
  .*xp)).*(x+(-1).*xp)+(-1).*a1.*dt.*(ubar+(-1).*ubarp).*(x+(-1).*xp).^3)+ ...
  a2.^2.*dt.^2.*(3.*exp(1).^(2.*logtheta)+(-6).*exp(1).^logtheta.*(x+(-1) ...
  .*xp).^2+(x+(-1).*xp).^4))+(-1).*a1.*a3.*dt.^2.*exp(1).^(2.*logtheta).*( ...
  ubar+ubarp).*(3.*exp(1).^(2.*logtheta)+(-6).*exp(1).^logtheta.*(x+(-1).* ...
  xp).^2+(x+(-1).*xp).^4)+a3.^2.*dt.^2.*(15.*exp(1).^(3.*logtheta)+(-45).* ...
  exp(1).^(2.*logtheta).*(x+(-1).*xp).^2+15.*exp(1).^logtheta.*(x+(-1).* ...
  xp).^4+(-1).*(x+(-1).*xp).^6))+a4.*dt.*exp(1).^(2.*logtheta).*((-2).* ...
  a2.*dt.*(15.*exp(1).^(3.*logtheta)+(-45).*exp(1).^(2.*logtheta).*(x+(-1) ...
  .*xp).^2+15.*exp(1).^logtheta.*(x+(-1).*xp).^4+(-1).*(x+(-1).*xp).^6)+ ...
  exp(1).^logtheta.*(6.*exp(1).^(3.*logtheta)+(-3).*exp(1).^(2.*logtheta) ...
  .*(5.*a1.*dt.*(ubar+(-1).*ubarp)+4.*(x+(-1).*xp)).*(x+(-1).*xp)+(-1).* ...
  a1.*dt.*(ubar+(-1).*ubarp).*(x+(-1).*xp).^5+2.*exp(1).^logtheta.*(x+(-1) ...
  .*xp).^3.*(5.*a1.*dt.*(ubar+(-1).*ubarp)+x+(-1).*xp)))+a4.^2.*dt.^2.*( ...
  105.*exp(1).^(4.*logtheta)+(-420).*exp(1).^(3.*logtheta).*(x+(-1).*xp) ...
  .^2+210.*exp(1).^(2.*logtheta).*(x+(-1).*xp).^4+(-28).*exp(1) ...
  .^logtheta.*(x+(-1).*xp).^6+(x+(-1).*xp).^8));

case 2 % logtheta

K=exp(1).^(logsigma+(-8).*logtheta+(-1/2).*exp(1).^((-1).*logtheta).*(x+( ...
  -1).*xp).^2).*(exp(1).^(3.*logtheta).*(exp(1).^(2.*logtheta).*(6.* ...
  a2.^2.*dt.^2.*(exp(1).^logtheta+(-1).*(x+(-1).*xp).^2)+exp(1) ...
  .^logtheta.*(4.*exp(1).^(2.*logtheta)+3.*a1.*dt.*exp(1).^logtheta.*(a1.* ...
  dt.*ubar.*ubarp+(-1).*(ubar+(-1).*ubarp).*(x+(-1).*xp))+(-2).*a1.^2.* ...
  dt.^2.*ubar.*ubarp.*(x+(-1).*xp).^2)+a2.*dt.*((-6).*exp(1).^(2.* ...
  logtheta)+2.*exp(1).^logtheta.*(3.*a1.*dt.*(ubar+(-1).*ubarp)+2.*(x+(-1) ...
  .*xp)).*(x+(-1).*xp)+(-1).*a1.*dt.*(ubar+(-1).*ubarp).*(x+(-1).*xp).^3)) ...
  +2.*exp(1).^logtheta.*(exp(1).^(2.*logtheta).*(exp(1).^(2.*logtheta)+ ...
  a1.*dt.*exp(1).^logtheta.*(a1.*dt.*ubar.*ubarp+(-1).*(ubar+(-1).*ubarp) ...
  .*(x+(-1).*xp))+(-1).*a1.^2.*dt.^2.*ubar.*ubarp.*(x+(-1).*xp).^2)+a2.* ...
  dt.*exp(1).^logtheta.*((-2).*exp(1).^(2.*logtheta)+exp(1).^logtheta.*( ...
  3.*a1.*dt.*(ubar+(-1).*ubarp)+2.*(x+(-1).*xp)).*(x+(-1).*xp)+(-1).*a1.* ...
  dt.*(ubar+(-1).*ubarp).*(x+(-1).*xp).^3)+a2.^2.*dt.^2.*(3.*exp(1).^(2.* ...
  logtheta)+(-6).*exp(1).^logtheta.*(x+(-1).*xp).^2+(x+(-1).*xp).^4))+(-6) ...
  .*a1.*a3.*dt.^2.*exp(1).^(2.*logtheta).*(ubar+ubarp).*(exp(1).^logtheta+ ...
  (-1).*(x+(-1).*xp).^2)+15.*a3.^2.*dt.^2.*(3.*exp(1).^(2.*logtheta)+(-6) ...
  .*exp(1).^logtheta.*(x+(-1).*xp).^2+(x+(-1).*xp).^4)+(-2).*a1.*a3.* ...
  dt.^2.*exp(1).^logtheta.*(ubar+ubarp).*(3.*exp(1).^(2.*logtheta)+(-6).* ...
  exp(1).^logtheta.*(x+(-1).*xp).^2+(x+(-1).*xp).^4))+2.*exp(1).^(2.* ...
  logtheta).*(exp(1).^(2.*logtheta).*(exp(1).^(2.*logtheta).*(exp(1).^(2.* ...
  logtheta)+a1.*dt.*exp(1).^logtheta.*(a1.*dt.*ubar.*ubarp+(-1).*(ubar+( ...
  -1).*ubarp).*(x+(-1).*xp))+(-1).*a1.^2.*dt.^2.*ubar.*ubarp.*(x+(-1).*xp) ...
  .^2)+a2.*dt.*exp(1).^logtheta.*((-2).*exp(1).^(2.*logtheta)+exp(1) ...
  .^logtheta.*(3.*a1.*dt.*(ubar+(-1).*ubarp)+2.*(x+(-1).*xp)).*(x+(-1).* ...
  xp)+(-1).*a1.*dt.*(ubar+(-1).*ubarp).*(x+(-1).*xp).^3)+a2.^2.*dt.^2.*( ...
  3.*exp(1).^(2.*logtheta)+(-6).*exp(1).^logtheta.*(x+(-1).*xp).^2+(x+(-1) ...
  .*xp).^4))+(-1).*a1.*a3.*dt.^2.*exp(1).^(2.*logtheta).*(ubar+ubarp).*( ...
  3.*exp(1).^(2.*logtheta)+(-6).*exp(1).^logtheta.*(x+(-1).*xp).^2+(x+(-1) ...
  .*xp).^4)+a3.^2.*dt.^2.*(15.*exp(1).^(3.*logtheta)+(-45).*exp(1).^(2.* ...
  logtheta).*(x+(-1).*xp).^2+15.*exp(1).^logtheta.*(x+(-1).*xp).^4+(-1).*( ...
  x+(-1).*xp).^6))+2.*a4.*dt.*exp(1).^(2.*logtheta).*((-2).*a2.*dt.*(15.* ...
  exp(1).^(3.*logtheta)+(-45).*exp(1).^(2.*logtheta).*(x+(-1).*xp).^2+15.* ...
  exp(1).^logtheta.*(x+(-1).*xp).^4+(-1).*(x+(-1).*xp).^6)+exp(1) ...
  .^logtheta.*(6.*exp(1).^(3.*logtheta)+(-3).*exp(1).^(2.*logtheta).*(5.* ...
  a1.*dt.*(ubar+(-1).*ubarp)+4.*(x+(-1).*xp)).*(x+(-1).*xp)+(-1).*a1.*dt.* ...
  (ubar+(-1).*ubarp).*(x+(-1).*xp).^5+2.*exp(1).^logtheta.*(x+(-1).*xp) ...
  .^3.*(5.*a1.*dt.*(ubar+(-1).*ubarp)+x+(-1).*xp)))+(exp(1).^(2.*logtheta) ...
  .*(exp(1).^(2.*logtheta).*(exp(1).^(2.*logtheta).*(exp(1).^(2.*logtheta) ...
  +a1.*dt.*exp(1).^logtheta.*(a1.*dt.*ubar.*ubarp+(-1).*(ubar+(-1).*ubarp) ...
  .*(x+(-1).*xp))+(-1).*a1.^2.*dt.^2.*ubar.*ubarp.*(x+(-1).*xp).^2)+a2.* ...
  dt.*exp(1).^logtheta.*((-2).*exp(1).^(2.*logtheta)+exp(1).^logtheta.*( ...
  3.*a1.*dt.*(ubar+(-1).*ubarp)+2.*(x+(-1).*xp)).*(x+(-1).*xp)+(-1).*a1.* ...
  dt.*(ubar+(-1).*ubarp).*(x+(-1).*xp).^3)+a2.^2.*dt.^2.*(3.*exp(1).^(2.* ...
  logtheta)+(-6).*exp(1).^logtheta.*(x+(-1).*xp).^2+(x+(-1).*xp).^4))+(-1) ...
  .*a1.*a3.*dt.^2.*exp(1).^(2.*logtheta).*(ubar+ubarp).*(3.*exp(1).^(2.* ...
  logtheta)+(-6).*exp(1).^logtheta.*(x+(-1).*xp).^2+(x+(-1).*xp).^4)+ ...
  a3.^2.*dt.^2.*(15.*exp(1).^(3.*logtheta)+(-45).*exp(1).^(2.*logtheta).*( ...
  x+(-1).*xp).^2+15.*exp(1).^logtheta.*(x+(-1).*xp).^4+(-1).*(x+(-1).*xp) ...
  .^6))+a4.*dt.*exp(1).^(2.*logtheta).*((-2).*a2.*dt.*(15.*exp(1).^(3.* ...
  logtheta)+(-45).*exp(1).^(2.*logtheta).*(x+(-1).*xp).^2+15.*exp(1) ...
  .^logtheta.*(x+(-1).*xp).^4+(-1).*(x+(-1).*xp).^6)+exp(1).^logtheta.*( ...
  6.*exp(1).^(3.*logtheta)+(-3).*exp(1).^(2.*logtheta).*(5.*a1.*dt.*(ubar+ ...
  (-1).*ubarp)+4.*(x+(-1).*xp)).*(x+(-1).*xp)+(-1).*a1.*dt.*(ubar+(-1).* ...
  ubarp).*(x+(-1).*xp).^5+2.*exp(1).^logtheta.*(x+(-1).*xp).^3.*(5.*a1.* ...
  dt.*(ubar+(-1).*ubarp)+x+(-1).*xp)))+a4.^2.*dt.^2.*(105.*exp(1).^(4.* ...
  logtheta)+(-420).*exp(1).^(3.*logtheta).*(x+(-1).*xp).^2+210.*exp(1).^( ...
  2.*logtheta).*(x+(-1).*xp).^4+(-28).*exp(1).^logtheta.*(x+(-1).*xp).^6+( ...
  x+(-1).*xp).^8)).*((-8)+(1/2).*exp(1).^((-1).*logtheta).*(x+(-1).*xp) ...
  .^2)+28.*a4.^2.*dt.^2.*exp(1).^logtheta.*(15.*exp(1).^(3.*logtheta)+( ...
  -45).*exp(1).^(2.*logtheta).*(x+(-1).*xp).^2+15.*exp(1).^logtheta.*(x+( ...
  -1).*xp).^4+(-1).*(x+(-1).*xp).^6)+a4.*dt.*exp(1).^(3.*logtheta).*(6.* ...
  exp(1).^(3.*logtheta)+(-30).*a2.*dt.*(3.*exp(1).^(2.*logtheta)+(-6).* ...
  exp(1).^logtheta.*(x+(-1).*xp).^2+(x+(-1).*xp).^4)+2.*exp(1).^logtheta.* ...
  (9.*exp(1).^(2.*logtheta)+(-3).*exp(1).^logtheta.*(5.*a1.*dt.*(ubar+(-1) ...
  .*ubarp)+4.*(x+(-1).*xp)).*(x+(-1).*xp)+(x+(-1).*xp).^3.*(5.*a1.*dt.*( ...
  ubar+(-1).*ubarp)+x+(-1).*xp))+(-3).*exp(1).^(2.*logtheta).*(5.*a1.*dt.* ...
  (ubar+(-1).*ubarp)+4.*(x+(-1).*xp)).*(x+(-1).*xp)+(-1).*a1.*dt.*(ubar+( ...
  -1).*ubarp).*(x+(-1).*xp).^5+2.*exp(1).^logtheta.*(x+(-1).*xp).^3.*(5.* ...
  a1.*dt.*(ubar+(-1).*ubarp)+x+(-1).*xp)));

case 3 % a1

K=dt.*exp(1).^(logsigma+(-5).*logtheta+(-1/2).*exp(1).^((-1).*logtheta).*( ...
  x+(-1).*xp).^2).*(exp(1).^(2.*logtheta).*(exp(1).^logtheta.*(2.*a1.*dt.* ...
  ubar.*ubarp.*(exp(1).^logtheta+(-1).*(x+(-1).*xp).^2)+(-1).*exp(1) ...
  .^logtheta.*(ubar+(-1).*ubarp).*(x+(-1).*xp))+(-1).*a2.*dt.*(ubar+(-1).* ...
  ubarp).*((-3).*exp(1).^logtheta+(x+(-1).*xp).^2).*(x+(-1).*xp))+(-1).* ...
  a3.*dt.*exp(1).^logtheta.*(ubar+ubarp).*(3.*exp(1).^(2.*logtheta)+(-6).* ...
  exp(1).^logtheta.*(x+(-1).*xp).^2+(x+(-1).*xp).^4)+a4.*dt.*(ubar+(-1).* ...
  ubarp).*((-15).*exp(1).^(2.*logtheta)+10.*exp(1).^logtheta.*(x+(-1).*xp) ...
  .^2+(-1).*(x+(-1).*xp).^4).*(x+(-1).*xp));

case 4 % a2

K=dt.*exp(1).^(logsigma+(-6).*logtheta+(-1/2).*exp(1).^((-1).*logtheta).*( ...
  x+(-1).*xp).^2).*(exp(1).^(2.*logtheta).*(exp(1).^logtheta.*((-2).*exp( ...
  1).^(2.*logtheta)+exp(1).^logtheta.*(3.*a1.*dt.*(ubar+(-1).*ubarp)+2.*( ...
  x+(-1).*xp)).*(x+(-1).*xp)+(-1).*a1.*dt.*(ubar+(-1).*ubarp).*(x+(-1).* ...
  xp).^3)+2.*a2.*dt.*(3.*exp(1).^(2.*logtheta)+(-6).*exp(1).^logtheta.*(x+ ...
  (-1).*xp).^2+(x+(-1).*xp).^4))+(-2).*a4.*dt.*(15.*exp(1).^(3.*logtheta)+ ...
  (-45).*exp(1).^(2.*logtheta).*(x+(-1).*xp).^2+15.*exp(1).^logtheta.*(x+( ...
  -1).*xp).^4+(-1).*(x+(-1).*xp).^6));

case 5 % a3

K=dt.^2.*exp(1).^(logsigma+(-6).*logtheta+(-1/2).*exp(1).^((-1).*logtheta) ...
  .*(x+(-1).*xp).^2).*((-1).*a1.*exp(1).^(2.*logtheta).*(ubar+ubarp).*(3.* ...
  exp(1).^(2.*logtheta)+(-6).*exp(1).^logtheta.*(x+(-1).*xp).^2+(x+(-1).* ...
  xp).^4)+2.*a3.*(15.*exp(1).^(3.*logtheta)+(-45).*exp(1).^(2.*logtheta).* ...
  (x+(-1).*xp).^2+15.*exp(1).^logtheta.*(x+(-1).*xp).^4+(-1).*(x+(-1).*xp) ...
  .^6));

case 6 % a4

K=dt.*exp(1).^(logsigma+(-8).*logtheta+(-1/2).*exp(1).^((-1).*logtheta).*( ...
  x+(-1).*xp).^2).*(exp(1).^(2.*logtheta).*((-2).*a2.*dt.*(15.*exp(1).^( ...
  3.*logtheta)+(-45).*exp(1).^(2.*logtheta).*(x+(-1).*xp).^2+15.*exp(1) ...
  .^logtheta.*(x+(-1).*xp).^4+(-1).*(x+(-1).*xp).^6)+exp(1).^logtheta.*( ...
  6.*exp(1).^(3.*logtheta)+(-3).*exp(1).^(2.*logtheta).*(5.*a1.*dt.*(ubar+ ...
  (-1).*ubarp)+4.*(x+(-1).*xp)).*(x+(-1).*xp)+(-1).*a1.*dt.*(ubar+(-1).* ...
  ubarp).*(x+(-1).*xp).^5+2.*exp(1).^logtheta.*(x+(-1).*xp).^3.*(5.*a1.* ...
  dt.*(ubar+(-1).*ubarp)+x+(-1).*xp)))+2.*a4.*dt.*(105.*exp(1).^(4.* ...
  logtheta)+(-420).*exp(1).^(3.*logtheta).*(x+(-1).*xp).^2+210.*exp(1).^( ...
  2.*logtheta).*(x+(-1).*xp).^4+(-28).*exp(1).^logtheta.*(x+(-1).*xp).^6+( ...
  x+(-1).*xp).^8));

otherwise

        K = zeros(n_x, n_xp);
end

if K == 0

    K = zeros(n_x, n_xp);

end

end

author: maziarraissi

repo: HPM

save_path: github-repos/MATLAB/maziarraissi-HPM


text:
function r = subsref(a,s)
%SUBSREF      Implements subscripted references for Taylor
%

% written  05/21/09     S.M. Rump
% modified 07/07/10     S.M. Rump  treatment of "end"-index
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  INTLAB_TAYLOR_ORDER = getappdata(0,'INTLAB_TAYLOR_ORDER');

  while 1
    if ~isa(a,'taylor')                  % index reference a.x(i) etc.
      r = subsref(a,s(1));
    elseif strcmp(s(1).type,'()')        % index reference a(i)
      setappdata(0,'INTLAB_TAYLOR_END',0);  % reset INTLAB variable
      index = reshape(1:prod(a.size),a.size);
      index = index(s(1).subs{:});
      r.size = size(index);
      r.t = a.t(:,index(:));
      r = class(r,'taylor');
    elseif strcmp(s(1).type,'{}')        % Taylor coefficient reference a{i}
      INTLAB_TAYLOR_END = getappdata(0,'INTLAB_TAYLOR_END');
      if INTLAB_TAYLOR_END
        setappdata(0,'INTLAB_TAYLOR_END',0); % reset INTLAB variable
        error('"end" cannot be used in an index expression when accessing derivatives by {}')
      end
      setappdata(0,'INTLAB_TAYLOR_END',0);   % reset INTLAB variable
      index = ( 1:INTLAB_TAYLOR_ORDER+1 );
      index = index(s(1).subs{:}+1);
      r = a.t(index,:).';
    elseif strcmp(s(1).type,'.')         % index reference a.t
      if strcmp(s(1).subs,'t')
        r = a.t.';
      elseif strcmp(s(1).subs,'mid')
        r = mid(a);
      else
        error('invalid subscript reference for taylor')
      end
    else
      error('invalid index reference for taylor')
    end
    if length(s)==1
      if rndold
        setround(rndold)
      end
      return
    end
    s = s(2:end);
    a = r;
  end

repo: OpenMAS


text:
function [lambda_vec, error_train, error_val] = ...
    validationCurve(X, y, Xval, yval)
%VALIDATIONCURVE Generate the train and validation errors needed to
%plot a validation curve that we can use to select lambda
%   [lambda_vec, error_train, error_val] = ...
%       VALIDATIONCURVE(X, y, Xval, yval) returns the train
%       and validation errors (in error_train, error_val)
%       for different values of lambda. You are given the training set (X,
%       y) and validation set (Xval, yval).
%

% Selected values of lambda (you should not change this)
lambda_vec = [0 0.001 0.003 0.01 0.03 0.1 0.3 1 3 10]';

% You need to return these variables correctly.
error_train = zeros(length(lambda_vec), 1);
error_val = zeros(length(lambda_vec), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return training errors in
%               error_train and the validation errors in error_val. The
%               vector lambda_vec contains the different lambda parameters
%               to use for each calculation of the errors, i.e,
%               error_train(i), and error_val(i) should give
%               you the errors obtained after training with
%               lambda = lambda_vec(i)
%
% Note: You can loop over lambda_vec with the following:
%
%       for i = 1:length(lambda_vec)
%           lambda = lambda_vec(i);
%           % Compute train / val errors when training linear
%           % regression with regularization parameter lambda
%           % You should store the result in error_train(i)
%           % and error_val(i)
%           ....
%
%       end
%
%

for i = 1:length(lambda_vec)
    lambda = lambda_vec(i);
    [theta] = trainLinearReg(X, y, lambda);
    error_train(i) = linearRegCostFunction(X, y, theta, 0);
    error_val(i) = linearRegCostFunction(Xval, yval, theta, 0);
end

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

end

author: zlotus


text:
filename='Gripping_quad_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'nonadjoint_compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
constraint_case = 'INEQUALITY';
optimizer = 'PROJECTED GRADIENT'; incrementFactor = 1;
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.5;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [mtype, mkind, skind] = RBtype (A)                                 %#ok
%RBTYPE determine the Rutherford/Boeing type of a sparse matrix
% Usage:
%   [mtype mkind skind] = RBtype (A)
%
% A must be a sparse matrix.  RBtype determines the Rutherford/Boeing type of A.
% Very little memory is used (just size(A,2) integer workspace), so this can
% succeed where a test such as nnz(A-A')==0 will fail.
%
%       mkind:  R: (0), A is real, and not binary
%               P: (1), A is binary (all values or 0 or 1)
%               C: (2), A is complex
%               I: (3), A is integer
%
%       skind:  R: (-1), A is rectangular
%               U: (0), A is unsymmetric (not S, H, or Z)
%               S: (1), A is symmetric (nnz(A-A.') is 0)
%               H: (2), A is Hermitian (nnz(A-A') is 0)
%               Z: (3), A is skew symmetric (nnz(A+A.') is 0)
%
% mtype is a 3-character string, where mtype(1) is the mkind
% ('R', 'P', or 'C').  mtype(2) is the skind ('R', 'U', 'S', 'H', or 'Z'),
% and mtype(3) is 'A'.
%
% Example:
%   load west0479
%   A = west0479 ;
%   RBtype (A)
%   RBtype (spones (A))
%   RBtype (2*spones (A))
%   C = A+A' ;
%   RBtype (C)
%
% See also RBread, RBwrite.

% Copyright 2007, Timothy A. Davis, University of Florida

error ('RBtype mexFunction not found') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [K, P] = whiteVardistPsi1Compute(whitekern, vardist, Z)

% WHITEVARDISTPSI1COMPUTE one line description
% FORMAT
% DESC description
% RETURN K : description
% RETURN P : description
% ARG whiteKern : the kernel structure associated with the white kernel.
% ARG vardist : description
% ARG Z : description
%
%
% SEEALSO : others
%
%
% COPYRIGHT : Michalis K. Titsias, 2009
%

% VARGPLVM

  K = zeros(size(vardist.means,1), size(Z,1));

  P = [];

author: SheffieldML

repo: vargplvm


text:
classdef NSLS < ALGORITHM
% <multi> <real/integer>
% Multiobjective optimization framework based on nondominated sorting and
% local search

%------------------------------- Reference --------------------------------
% B. Chen, W. Zeng, Y. Lin, and D. Zhang, A new local search-based
% multiobjective optimization algorithm, IEEE Transactions on Evolutionary
% Computation, 2015, 19(1): 50-73.
%------------------------------- 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)
            %% Generate random population
            Population = Problem.Initialization();

            %% Optimization
            while Algorithm.NotTerminated(Population)
                Offspring  = Operator(Problem,Population);
                Population = EnvironmentalSelection([Population,Offspring],Problem.N);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%PRM Probabilistic RoadMap navigation class
%
% A concrete subclass of the abstract Navigation class that implements the
% probabilistic roadmap navigation algorithm over an occupancy grid.  This
% performs goal independent planning of roadmaps, and at the query stage
% finds paths between specific start and goal points.
%
% Methods::
%
% plan         Compute the roadmap
% path         Compute a path to the goal
% visualize    Display the obstacle map (deprecated)
% plot         Display the obstacle map
% display      Display the parameters in human readable form
% char         Convert to string
%
% Example::
%
%        load map1              % load map
%        goal = [50,30];        % goal point
%        start = [20, 10];      % start point
%        prm = PRM(map);        % create navigation object
%        prm.plan()             % create roadmaps
%        prm.path(start, goal)  % animate path from this start location
%
% References::
%
% - Probabilistic roadmaps for path planning in high dimensional configuration spaces,
%   L. Kavraki, P. Svestka, J. Latombe, and M. Overmars,
%   IEEE Transactions on Robotics and Automation, vol. 12, pp. 566-580, Aug 1996.
% - Robotics, Vision & Control, Section 5.2.4,
%   P. Corke, Springer 2011.
%
% See also Navigation, DXform, Dstar, PGraph.

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

% Peter Corke 8/2009.

classdef PRM < Navigation

    properties
        npoints         % number of sample points
        distthresh      % distance threshold, links between vertices
                        % must be less than this.

        graph           % graph Object representing random nodes

        vgoal           % index of vertex closest to goal
        vstart          % index of vertex closest to start
        localGoal       % next vertex on the roadmap
        localPath       % set of points along path to next vertex
        gpath           % list of vertices between start and goal
    end

    methods

        % constructor
        function prm = PRM(varargin)
            %PRM.PRM Create a PRM navigation object
            %
            % P = PRM(MAP, options) is a probabilistic roadmap navigation
            % object, and MAP is an occupancy grid, a representation of a
            % planar world as a matrix whose elements are 0 (free space) or 1
            % (occupied).
            %
            % Options::
            %  'npoints',N      Number of sample points (default 100)
            %  'distthresh',D   Distance threshold, edges only connect vertices closer
            %                   than D (default 0.3 max(size(occgrid)))
            %
            % Other options are supported by the Navigation superclass.
            %
            % See also Navigation.Navigation.

            % invoke the superclass constructor, it handles some options
            prm = prm@Navigation(varargin{:});

            % create an empty 2D graph
            prm.graph = PGraph(2);

            % parse out PRM specific options and save in the navigation object
            opt.npoints = 100;
            opt.distthresh = 0.3*max(size(prm.occgrid));
            [opt,args] = tb_optparse(opt, varargin);
            prm.npoints = opt.npoints;
            prm.distthresh = opt.distthresh;
        end

        function plan(prm)
        %PRM.plan Create a probabilistic roadmap
        %
        % P.plan() creates the probabilistic roadmap by randomly
        % sampling the free space in the map and building a graph with
        % edges connecting close points.  The resulting graph is kept
        % within the object.

        % build a graph over the free space
            prm.message('create the graph');

            prm.graph.clear();  % empty the graph
            create_roadmap(prm);  % build the graph
        end

        function p = path(prm, start, goal)
        %PRM.path Find a path between two points
        %
        % P.path(START, GOAL) finds and displays a path from START to GOAL
        % which is overlaid on the occupancy grid.
        %
        % X = P.path(START) returns the path (2xM) from START to GOAL.

            if nargin < 3
                error('must specify start and goal');
            end

            % set the goal coordinate
            prm.goal = goal;

            % invoke the superclass path function, which iterates on our
            % next method
            if nargout == 0
                path@Navigation(prm, start);
            else
                p = path@Navigation(prm, start);
            end
        end

        % Handler invoked by Navigation.path() to start the navigation process
        %
        %   - find a path through the graph
        %   - determine vertices closest to start and goal
        %   - find path to first vertex
        function navigate_init(prm, start)

            % find the vertex closest to the goal
            prm.vgoal = prm.graph.closest(prm.goal);

            % find the vertex closest to the start
            prm.vstart = prm.graph.closest(start);

            % are the vertices connected?
            if prm.graph.component(prm.vstart) ~= prm.graph.component(prm.vgoal)
                error('PRM:plan:nopath', 'PRM: start and goal not connected: rerun the planner');
            end

            % find a path through the graph
            prm.message('planning path through graph');
            prm.graph.goal(prm.vgoal);   % set the goal
            prm.gpath = prm.graph.path(prm.vstart);

            % the path is a list of nodes from vstart to vgoal
            % discard the first vertex, since we plan a local path to it
            prm.gpath = prm.gpath(2:end);

            % start the navigation engine with a path to the nearest vertex
            prm.graph.highlight_node(prm.vstart);

            prm.localPath = bresenham(start, prm.graph.coord(prm.vstart));
            prm.localPath = prm.localPath(2:end,:);
        end

        % Invoked for each step on the path by path() method.
        function n = next(prm, p)

            if all(p(:) == prm.goal)
                n = [];     % signal that we've arrived
                return;
            end

            % we take the next point from the localPath
            if numrows(prm.localPath) == 0
                % local path is consumed, move to next vertex
                if isempty(prm.gpath)
                    % we have arrived at the goal vertex
                    % make the path from this vertex to the goal coordinate
                    prm.localPath = bresenham(p, prm.goal);
                    prm.localPath = prm.localPath(2:end,:);
                    prm.localGoal = [];
                else
                    % set local goal to next vertex in gpath and remove it from the list
                    prm.localGoal = prm.gpath(1);
                    prm.gpath = prm.gpath(2:end);

                    % compute local path to the next vertex
                    prm.localPath = bresenham(p, prm.graph.coord(prm.localGoal));
                    prm.localPath = prm.localPath(2:end,:);
                    prm.graph.highlight_node(prm.localGoal);
                end
            end

            n = prm.localPath(1,:)';     % take the first point
            prm.localPath = prm.localPath(2:end,:); % and remove from the path
        end

        function s = char(prm)
        %PRM.char  Convert to string
        %
        % P.char() is a string representing the state of the PRM
        % object in human-readable form.
        %
        % See also PRM.display.

            % invoke the superclass char() method
            s = char@Navigation(prm);

            % add PRM specific stuff information
            s = char(s, sprintf('  graph size: %d', prm.npoints));
            s = char(s, sprintf('  dist thresh: %f', prm.distthresh));
            s = char(s, char(prm.graph) );
        end

        function plot(prm, varargin)
        %PRM.plot Visualize navigation environment
        %
        % P.plot() displays the occupancy grid with an optional distance field.
        %
        % Options::
        %  'goal'            Superimpose the goal position if set
        %  'nooverlay'       Don't overlay the PRM graph

            opt.nooverlay = false;
            [opt,args] = tb_optparse(opt, varargin);

            % display the occgrid
            plot@Navigation(prm, args{:});

            if ~opt.nooverlay
                hold on
                prm.graph.plot()%varargin{:});
                hold off
            end
        end

    end % method

    methods (Access='protected')
    % private methods
        % create the roadmap
        function create_roadmap(prm)

            for j=1:prm.npoints
                % pick a point not in obstacle
                while true
                    x = prm.randi(numcols(prm.occgrid));
                    y = prm.randi(numrows(prm.occgrid));
                    if prm.occgrid(y,x) == 0
                        break;
                    end
                end
                new = [x; y];

                vnew = prm.graph.add_node(new);

                [d,v] = prm.graph.distances(new);
                % test neighbours in order of increasing distance
                for i=1:length(d)
                    if v(i) == vnew
                        continue;
                    end
                    if d(i) > prm.distthresh
                        continue;
                    end
                    if ~prm.testpath(new, prm.graph.coord(v(i)))
                        continue;
                    end
                    prm.graph.add_edge(vnew, v(i));
                end
            end
        end

        % test the path from p1 to p2 is entirely in free space
        function c = testpath(prm, p1, p2)
            p = bresenham(p1, p2);

            for pp=p'
                if prm.occgrid(pp(2), pp(1)) > 0
                    c = false;
                    return;
                end
            end
            c = true;
        end

    end % private methods

end % classdef

author: Allopart

repo: rbpf-gmapping


text:
function [ px, pxp ] = least_val2 ( nterms, b, c, d, x )

%*****************************************************************************80
%
%% LEAST_VAL2 evaluates a least squares polynomial defined by LEAST_SET.
%
%  Discussion:
%
%    This routine also computes the derivative of the polynomial.
%
%    The least squares polynomial is assumed to be defined as a sum
%
%      P(X) = SUM ( I = 1 to NTERMS ) D(I) * P(I-1,X)
%
%    where the orthogonal basis polynomials P(I,X) satisfy the following
%    three term recurrence:
%
%      P(-1,X) = 0
%      P(0,X) = 1
%      P(I,X) = ( X - B(I-1) ) * P(I-1,X) - C(I-1) * P(I-2,X)
%
%    Therefore, the least squares polynomial can be evaluated as follows:
%
%    If NTERMS is 1, then the value of P(X) is D(1) * P(0,X) = D(1).
%
%    Otherwise, P(X) is defined as the sum of NTERMS > 1 terms.  We can
%    reduce the number of terms by 1, because the polynomial P(NTERMS,X)
%    can be rewritten as a sum of polynomials;  Therefore, P(NTERMS,X)
%    can be eliminated from the sum, and its coefficient merged in with
%    those of other polynomials.  Repeat this process for P(NTERMS-1,X)
%    and so on until a single term remains.
%    P(NTERMS,X) of P(NTERMS-1,X)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NTERMS, the number of terms in the least squares
%    polynomial.  NTERMS must be at least 1.  The value of NTERMS
%    may be reduced from the value given to LEAST_SET.
%    This will cause LEAST_VAL to evaluate the least squares polynomial
%    of the lower degree specified.
%
%    Input, real B(NTERMS), C(NTERMS), D(NTERMS), the information
%    computed by LEAST_SET.
%
%    Input, real X, the point at which the least squares polynomial
%    is to be evaluated.
%
%    Output, real PX, PXP, the value and derivative of the least
%    squares polynomial at X.
%
  px = d(nterms);
  pxp = 0.0;
  pxm1 = 0.0;
  pxpm1 = 0.0;

  for i = nterms-1 : -1 : 1

    pxm2 = pxm1;
    pxpm2 = pxpm1;
    pxm1 = px;
    pxpm1 = pxp;

    if ( i == nterms-1 )
      px = d(i) + ( x - b(i) ) * pxm1;
      pxp = pxm1 + ( x - b(i) ) * pxpm1;
    else
      px = d(i) + ( x - b(i) ) * pxm1 - c(i+1) * pxm2;
      pxp = pxm1 + ( x - b(i) ) * pxpm1 - c(i+1) * pxpm2;
    end

  end

  return
end


text:
function integral = ball01_monomial_integral ( e )

%*****************************************************************************80
%
%% BALL01_MONOMIAL_INTEGRAL returns monomial integrals in the unit ball.
%
%  Discussion:
%
%    The integration region is
%
%      X^2 + Y^2 + Z^2 <= 1.
%
%    The monomial is F(X,Y,Z) = X^E(1) * Y^E(2) * Z^E(3).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 January 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gerald Folland,
%    How to Integrate a Polynomial Over a Sphere,
%    American Mathematical Monthly,
%    Volume 108, Number 5, May 2001, pages 446-448.
%
%  Parameters:
%
%    Input, integer E(3), the exponents of X, Y and Z in the
%    monomial.  Each exponent must be nonnegative.
%
%    Output, real INTEGRAL, the integral.
%
  if ( any ( e(1:3) < 0 ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BALL01_MONOMIAL_INTEGRAL - Fatal error!\n' );
    fprintf ( 1, '  All exponents must be nonnegative.\n' );
    error ( 'BALL01_MONOMIAL_INTEGRAL - Fatal error!' );
  end
%
%  Integrate over the surface.
%
  if ( all ( e(1:3) == 0 ) )

    integral = 2.0 * sqrt ( pi ^ 3 ) / gamma ( 1.5 );

  elseif ( any ( mod ( e(1:3), 2 ) == 1 ) )

    integral = 0.0;

  else

    integral = 2.0;

    for i = 1 : 3
      integral = integral * gamma ( 0.5 * ( e(i) + 1 ) );
    end

    integral = integral / gamma ( 0.5 * ( sum ( e(1:3) + 1 ) ) );

  end
%
%  The surface integral is now adjusted to give the volume integral.
%
  r = 1.0;
  s = sum ( e(1:3) ) + 3;

  integral = integral * r ^ s / s;

  return
end


text:
function test_spm12

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

% currently (Jan, 2017) SPM12 support for:
% - ft_volumerealign
% - ft_volumedownsample (incl smoothing)
% - private/mni2tal
% - private/tal2mni
% - private/volumesmooth
% - ft_read_mri (nifti_spm, default remains spm8) - based on function deps, not real data
% - private/volumewrite_spm (default remains spm8) - based on function deps, not real data
% - utilities/private/sn2individual (called by ft_warp_apply) - based on function deps, not real data
% - utilities/private/individual2sn (called by ft_warp_apply) - forced use of spm8 since spm_invdef is missing in spm12 (in toolbox/OldSeg)
% to do (function-specific intelligence required):
% - ft_volumenormalise
% - ft_volumesegment

mrifile1 = dccnpath('/home/common/matlab/fieldtrip/data/test/latest/mri/freesurfer/T1.mgz');
mri1     = ft_read_mri(mrifile1);
mri1.coordsys = 'tal'; % I don't think this is technically correct, it should be fsaverage aka mni305

mrifile2 = dccnpath('/home/common/matlab/fieldtrip/data/test/latest/mri/nifti/single_subj_T1.nii');
mri2     = ft_read_mri(mrifile2);
mri2.coordsys = 'mni';

mrifile3 = dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/Subject01.mri');
mri3     = ft_read_mri(mrifile3);
mri3.coordsys = 'ctf';

%----------------------------- SPM12 -----------------------------------

%ft_convert_coordsys: currently uses OldNorm sub-toolbox (i.e. SPM8)
ft_hastoolbox('spm12',1)
c2a = ft_convert_coordsys(mri3, 'tal', 1);
c2b = ft_convert_coordsys(mri3, 'tal', 2);

try
  rmpath(spm('dir'));
end

%ft_volumerealign: coregistration (used in human ecog tutorial)
cfg               = [];
cfg.method        = 'spm';
cfg.spmversion    = 'spm12';
cfg.spm.cost_fun  = 'nmi';
r12 = ft_volumerealign(cfg, mri1, mri2);

try
  rmpath(spm('dir'));
end

%ft_volumedownsample
tmp            = randn(256,256,256);
mri1.pow        = tmp;
mri1.pow(1)     = 0;
cfg            = [];
cfg.downsample = 2;
cfg.spmversion = 'spm12';
cfg.smooth     = 'no';
d12            = ft_volumedownsample(cfg, mri1);
cfg.smooth     = 5;
d12s           = ft_volumedownsample(cfg, mri1);

try
  rmpath(spm('dir'));
end

%mni2tal and tal2mni
[v, p] = ft_version;
cd(fullfile(p, 'private'));

inpoints  = randn(100,3);
outpoints = mni2tal(inpoints);

try
  rmpath(spm('dir'));
end

inpoints = tal2mni(outpoints);

%----------------------------- SPM8 COMPAT -----------------------------

try
  rmpath(spm('dir'));
end

%ft_volumenormalise
cfg            = [];
cfg.nonlinear  = 'no';
cfg.spmversion = 'spm8';
n8 = ft_volumenormalise(cfg, mri1);

%ft_warp_apply with spm12
elecpos = ft_warp_apply(n8.params, [4 4 4; 1 1 1], 'individual2sn');

try
  rmpath(spm('dir'));
end

%ft_volumesegment
cfg             = [];
cfg.spmversion  = 'spm8';
s8 = ft_volumesegment(cfg, mri1);

%-----------------------------------------------------------------------
% Notes from test_old_spm8:
%
% the following m-files contain a line with 'spm_' and have to be dealt with
%
%./external/bemcp/Makefile
%./spmcalls.txt
%./private/tal2mni.m
%./private/prepare_mesh_segmentation.m
%./private/prepare_dipole_grid.m
%./private/mni2tal.m
%./multivariate/toolboxes/gerven/bayesbrain/examples/html/fmri_bn_demo.html
%./multivariate/toolboxes/gerven/bayesbrain/examples/fmri_bn_demo.m
%./testing/test_bug62.m
%
%%done
%./private/avw_img_read.m
%./fileio/private/avw_img_read.m
%./fileio/read_mri.m
%./fileio/ft_read_mri.m
%
%%done and changed
%./ft_volumesegment.m
%./ft_volumenormalise.m
%./ft_volumedownsample.m
%./private/volumewrite_spm.m
%------------------------------------------------------------------------


text:
function [obSlice,reflections] = mrReflectObl(obSlice,obSize,reflections,which,update)
%
%NAME:	 [obSlice,reflections] = mrReflectObl(obSlice,obSize,reflections,which,update)
%AUTHOR:  Poirson
%DATE:    08.09.96
%HISTORY  11.08.96 SPG added bool update
%BUGS:

global obwin volslimin2 volslimax2;

row = obSize(1);
col = obSize(2);

if which == 1
	% First entry tell us left/right flipping
	reflections(1) = reflections(1) * (-1.0);
	tmp = fliplr(reshape(obSlice,row,col));

elseif which == 2
	% Second entry tell us up/down flipping
	reflections(2) = reflections(2) * (-1.0);
	tmp = flipud(reshape(obSlice,row,col));
else
	disp('mrReflectObl: invalid case');
end
 obSlice = tmp(:)';

if (update == 1)
 figure(obwin);
 myShowImageVol(obSlice',obSize,max(obSlice)*get(volslimin2,'value'),max(obSlice)*get(volslimax2,'value'));
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function plotClass(X, label)
% Plot 2d/3d samples of different classes with different colors.
% Written by Mo Chen (sth4nth@gmail.com).
[d,n] = size(X);
if nargin == 1
    label = ones(n,1);
end
assert(n == length(label));

color = 'brgmcyk';
m = length(color);
c = max(label);

figure(gcf);
clf;
hold on;
switch d
    case 2
        view(2);
        for i = 1:c
            idc = label==i;
%             plot(X(1,label==i),X(2,label==i),['.' color(i)],'MarkerSize',15);
            scatter(X(1,idc),X(2,idc),36,color(mod(i-1,m)+1));
        end
    case 3
        view(3);
        for i = 1:c
            idc = label==i;
%             plot3(X(1,idc),X(2,idci),X(3,idc),['.' idc],'MarkerSize',15);
            scatter3(X(1,idc),X(2,idc),X(3,idc),36,color(mod(i-1,m)+1));
        end
    otherwise
        error('ERROR: only support data of 2D or 3D.');
end
axis equal
grid on
hold off

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function normals = faceNormal(nodes, faces)
%FACENORMAL Compute normal vector of faces in a 3D mesh.
%
%   NORMALS = faceNormal(VERTICES, FACES)
%   VERTICES is a set of 3D points (as a N-by-3 array), and FACES is either
%   a N-by-3 index array or a cell array of indices. The function computes
%   the normal vector of each face.
%   The orientation of the normal is defined by the sign of cross product
%   between vectors joining vertices 1 to 2 and 1 to 3.
%
%
%   Example
%     [v e f] = createIcosahedron;
%     normals1 = faceNormal(v, f);
%     centros1 = faceCentroids(v, f);
%     figure; drawMesh(v, f);
%     hold on; axis equal; view(3);
%     drawVector3d(centros1, normals1);
%
%     pts = rand(50, 3);
%     hull = minConvexHull(pts);
%     normals2 = faceNormal(pts, hull);
%
%   See also
%   meshes3d, drawMesh, convhull, convhulln, drawVector3d

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2006-07-05
% Copyright 2006 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas)

% deprecation warning
warning('geom3d:deprecated', ...
    [mfilename ' is deprecated, use ''meshFaceNormals'' instead']);

if isnumeric(faces)
    % compute vector of first edges
	v1 = nodes(faces(:,2),1:3) - nodes(faces(:,1),1:3);
    v2 = nodes(faces(:,3),1:3) - nodes(faces(:,1),1:3);

%     % normalize vectors
%     v1 = normalizeVector3d(v1);
%     v2 = normalizeVector3d(v2);

    % compute normals using cross product (nodes have same size)
	normals = normalizeVector3d(cross(v1, v2, 2));

else
    % initialize empty array
    normals = zeros(length(faces), 3);

    for i = 1:length(faces)
        face = faces{i};
        % compute vector of first edges
        v1 = nodes(face(2),1:3) - nodes(face(1),1:3);
        v2 = nodes(face(3),1:3) - nodes(face(1),1:3);

%         % normalize vectors
%         v1 = normalizeVector3d(v1);
%         v2 = normalizeVector3d(v2);

        % compute normals using cross product
        normals(i, :) = normalizeVector3d(cross(v1, v2, 2));
    end
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function hx = happly (v, beta, x)
%HAPPLY apply Householder reflection to a vector
% Example:
%   hx = happly (v,beta,x) ;        % computes hx = x - v * (beta * (v' *x)) ;
% See also: testall

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

hx = x - v * (beta * (v' *x)) ;