text:
function K=kronSym(A,B)
%%KRONSYM Take the symmetric Kronecker product of the matrices A and B. The
%         standard Kronecker product of two real, square matrices A and B,
%         one can write the following relation with respect to any real
%         square matrix S:
%         kron(A,B)*vec(S)==vec(B*S*A')
%         However, the symmetric Kronecker product is such that for any two
%         real square matrices A and B and a square SYMMETRIC matrix S, one
%         can write
%         kronSym(A,B)*vech(S,sqrt(2))==vech((1/2)*(B*S*A'+A*S*B'),sqrt(2))
%         Whereas is A and B are nXn, the standard kronecker product is
%         n^2Xn^2, the symmetric Kronecker product is
%         (n*(n+1)/2)X(n*(n+1)/2).
%
%INPUTS: A, B Two real nXn matrices. They need not be square.
%
%OUTPUTS: K The (n*(n+1)/2)X(n*(n+1)/2) symmetric Kronecker product of A
%           and B.
%
%Symmetric Kronecker products are discussed in the appendix of [1], where
%they play a role in the implementation of a semidefinite programming
%algorithm.
%
%The symmetric Kronecker product is implemented by noting that in the
%identity
%kronSym(A,B)*vech(S,sqrt(2))==vech((1/2)*(B*S*A'+A*S*B'),sqrt(2)), if S is
%a symmetric matrix where on or below the main diagonal there is only a
%single nonzero element, then one effectively selects a column of
%kronSym(A,B). To get an unscaled column of kronSym(A,B), the nonzero
%element should be 1 if it is on the main diagonal and it should be
%1/sqrt(2) for elements below (and above) the main diagonal due to the
%sqrt(2) scaling in the vech command. Basically, a simple implementation
%would be
% n=size(A,1);
% prodDim=n*(n+1)/2;
%
% K=zeros(prodDim,prodDim);
% SelMat=zeros(n,n);
% for curEl=1:prodDim
%     [row,col]=vechInd2Sub(n,curEl);
%     if(row==col)
%         SelMat(row,col)=1;
%     else%row>col
%         SelMat(row,col)=1/sqrt(2);
%         SelMat(col,row)=1/sqrt(2);
%     end
%     K(:,curEl)=vech((1/2)*(B*SelMat*A'+A*SelMat*B'),sqrt(2));
%
%     SelMat(row,col)=0;
%     SelMat(col,row)=0;
% end
%However, that is rather slow. Thus, this implementation tries to directly
%write the elements of the output matrix K.
%
%REFERENCES:
%[1] F. Alizadeh, J.-P. A. Haeberly, and M. L. Overton, "Primal-dual
%    interior-point methods for semidefinite programming: Convergence
%    rates, stability and numerical results," SIAM Journal on Optimization,
%    vol. 8, no. 3, pp. 746-768, 1998.
%
%February 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%The columns of K are indexed using (i1,i2) and the rows are indexed using
%(j1,j2).
n=size(A,1);
prodDim=n*(n+1)/2;

K=zeros(prodDim,prodDim);
sqrt2=sqrt(2);
dblSqrt2=2*sqrt2;
for i1=1:n
    col=i1+(i1-1)*(2*n-i1)/2;
    for j1=1:n
        row=j1+(j1-1)*(2*n-j1)/2;

        K(row,col)=(A(j1,i1)*B(j1,i1)+B(j1,i1)*A(j1,i1))/2;
        for j2=(j1+1):n
            row=row+1;

            K(row,col)=(A(j1,i1)*B(j2,i1)+B(j1,i1)*A(j2,i1))/sqrt2;
        end
    end

    for i2=(i1+1):n
      col=col+1;
      for j1=1:n
         row=j1+(j1-1)*(2*n-j1)/2;

         K(row,col)=(A(j1,i1)*B(j1,i2)+A(j1,i2)*B(j1,i1)+B(j1,i1)*A(j1,i2)+B(j1,i2)*A(j1,i1))/dblSqrt2;
         for j2=(j1+1):n
            row=row+1;

            K(row,col)=(A(j1,i1)*B(j2,i2)+A(j1,i2)*B(j2,i1)+B(j1,i1)*A(j2,i2)+B(j1,i2)*A(j2,i1))/2;
         end
      end
    end
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [params, names] = expvarMeanExtractParam(model)

% EXPVARMEANEXTRACTPARAM Extract weights and biases from an EXPVARMEAN mapping.
% FORMAT
% DESC returns a vector of all the parameters from the mapping
% associated with the mean function of a exponential variational approximation.
% ARG model : the mappings from which we wish to extract the parameters.
% and biases.
% RETURN params : vector of all the parameters.
% RETURN names : optional additional returned cell array of the
% names of the parameters.
%
% SEEALSO : expvarMeanCreate, expvarMeanExpandParam, modelExtractParam
%
% COPYRIGHT : Neil D. Lawrence, 2006

% SHEFFIELDML

if nargout > 1
  [params, names] = kernExtractParam(model.kern);
else
  params = kernExtractParam(model.kern);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function v = sr_clamp(v, vLowerB, vUpperB)

% SC_CLAMP: clamp value v

v = max(v, vLowerB);
v = min(v, vUpperB);

end

repo: SelfExSR


text:
function [y,mask] = vl_nndropout(x,varargin)
%VL_NNDROPOUT CNN dropout.
%   [Y,MASK] = VL_NNDROPOUT(X) applies dropout to the data X. MASK
%   is the randomly sampled dropout mask. Both Y and MASK have the
%   same size as X.
%
%   VL_NNDROPOUT(X, 'rate', R) sets the dropout rate to R.
%
%   [DZDX] = VL_NNDROPOUT(X, DZDY, 'mask', MASK) computes the
%   derivatives of the blocks projected onto DZDY. Note that MASK must
%   be specified in order to compute the derivative consistently with
%   the MASK randomly sampled in the forward pass. DZDX and DZDY have
%   the same dimesnions as X and Y respectivey.
%
%   Note that in the original paper on dropout, at test time the
%   network weights for the dropout layers are scaled down to
%   compensate for having all the neurons active. In this
%   implementation the dropout function itself already does this
%   compensation during training. So at test time no alterations are
%   required.

% Copyright (C) 2014-15 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

opts.rate = 0.5 ;
opts.mask = [] ;

backMode = numel(varargin) > 0 && ~isstr(varargin{1}) ;
if backMode
  dzdy = varargin{1} ;
  opts = vl_argparse(opts, varargin(2:end)) ;
else
  opts = vl_argparse(opts, varargin) ;
end

% determine mask
mask = opts.mask ;
scale = single(1 / (1 - opts.rate)) ;
if backMode && isempty(mask)
  warning('vl_nndropout: when using in backward mode, the mask should be specified') ;
end
if isempty(mask)
  if isa(x,'gpuArray')
    mask = scale * single(gpuArray.rand(size(x)) >= opts.rate) ;
  else
    mask = scale * single(rand(size(x)) >= opts.rate) ;
  end
end

% do job
if ~backMode
  y = mask .* x ;
else
  y = mask .* dzdy ;
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function p=t2p(t,n,df,units)
% t2p - get p-value from t-statistical values
% p=t2p(t,n,df,units);
%  t     = t-statistical values
%  n     = number of comparisons for Bonferroni correction [1]
%  df    = degrees of freedom [Inf]
%  units = 'p', 'log10p' [p]

% 08/2005 SOD: created it.

% input checks and defaults
if ~exist('t','var') || isempty(t),   error('Need t');   end
if ~exist('n','var') || isempty(n),   n = 1;             end
if ~exist('df','var') || isempty(df), df = Inf;          end
if ~exist('units','var') || isempty(units), units = 'p'; end

% preserve the sign of t in p values
tsign = sign(t);
tsign(tsign==0)=1;
t = abs(t);

% Find the upper tail probs of the two-tailed t-distribution:
if df==Inf,
  p=0.5*erfc(t/sqrt(2));
else
  p=0.5*betainc(df./(df+t.^2),df/2,0.5);
end;

% Bonferroni
p=(n+1).*p;
p(p>1)=1; % otherwise log10p is funny

% convert
switch units
 case 'log10p',
  p(p==0) = 10^-50; % remove 0
  p       = -log10(abs(p)) .* tsign;
 case 'p',
  % preserve the sign of t in p values
  p = p.*tsign;
 otherwise,
  disp(sprintf('Unknown unit: %s',units));
end;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% module = noise_module_laplace(noise_module, varargin)

function module = noise_module_laplace(varargin)

options = struct('init', struct(), ...
                 'update_u', true);
[options, errmsg] = argparse(options, varargin{:});
error(errmsg);

U = [];
LogU = [];

module.initialize = @initialize;
module.update = @update;
module.get_struct = @get_struct;

  function [Tau] = initialize(M, N)

  U = nan(M,N);

  % Parse custom initialization
  init = struct('U', 1);
  [init, errmsg] = argparse(init, options.init);
  error(errmsg);

  % Initialize
  U(:,:) = init.U;
  Tau = U;
  end

  function S = get_struct()
  S.module = 'noise_module_independent_t';
  S.U = U;
  end

  function [Tau, LogTau, KL_Tau] = update(iter, E2, Obs)

  % OBS is the number of observations for each element of E2.

  [M,N] = size(E2);

  %
  % Update U
  %

  if index_selected(iter, options.update_u)

    % Update the distribution

    % a_U = bsxfun(@plus, nu, Obs) / 2;
    % b_U = bsxfun(@plus, nu, E2) / 2;
    % U = a_U ./ b_U;
    U(Obs) = E2(Obs) .^ (-0.5);
    LogU = nan(M,N); %psi(a_U) - log(b_U);
  end

% $$$   if iter == 4
% $$$     figure
% $$$     imagesc(U);
% $$$     error('debuggin')
% $$$   end
% $$$   if iter == 1
% $$$     figure
% $$$     imagesc(U);
% $$$     error('debuggin')
% $$$   end

  %minU = min(U(:));
  %maxU = max(U(:));

  % KL-term: <log q(U)> - <log p(U)>
  KL_U = NaN;
% $$$     0;
% $$$     KL_U = KL_U - sum(gamma_entropy(a_U(Obs),b_U(Obs)));
% $$$     switch options.nu
% $$$
% $$$      case 'pooled'
% $$$       KL_U = KL_U - sum(gamma_logpdf(U(Obs), nu/2, nu/2, LogU(Obs)));
% $$$
% $$$      case 'separate_rows'
% $$$       for m=1:M
% $$$         Nmv = Obs(m,:);
% $$$         KL_U = KL_U - sum(gamma_logpdf(U(m,Nmv),LogU(m,Nmv),nu(m)/2,nu(m)/2));
% $$$       end
% $$$
% $$$      case 'separate_columns'
% $$$       error('Not implemented yet')
% $$$
% $$$      otherwise
% $$$       error('Unknown model for nu');
% $$$
% $$$     end
% $$$
% $$$   else
% $$$
% $$$     KL_U = NaN;
% $$$
% $$$   end

  %
  % Compute the joint noise stuff
  %

  Tau = U;
  LogTau = LogU;
  KL_Tau = KL_U;

  end

end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
% Copyright (C) 1993-2014, by Joern Malzahn
%
% 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

%%begin

%% Getting started with symbolics and code generation
% This is a brief example about how we can derive symbolic robot model
% expressions and how we can generate robot specific functions as well as
% real-time capable Simulink blocks using the |CodeGenerator| class. The
% example uses a reduced version of the Puma 560 arm with the first 3
% links.
%
% A requirement for this demo is that we have the Mathworks Symbolic Toolbox
% installed besides the Robotics Toolbox.
%

%% Instantiate a |CodeGenerator| class object
% We start off with the instantiation of a |CodeGenerator| class object.
% First, we load the |SerialLink| object for which we intend to generate
% code.
mdl_planar3

%%
% After that, we find a |SerialLink| object named p3 in the workspace. This
% object is used to instantiate the CodeGenerator.
cGen = CodeGenerator(p3)

%% Code generation
% By default |CodeGenerator| class objects are configured to generate:
%
% * symbolic expressions
% * m-code
% * Simulink blocks
%
% and they document the CodeGeneration progress on the Matlab console.
% We may modify this behaviour by passing extra arguments to the
% |CodeGenerator| constructor. (Type |help CodeGenerator| for details)
%
% Now let's generate code for the forward kinematics of our reduced Puma
% 560.
symExp = cGen.genfkine

%%
% The text output to the console may be disabled
cGen.verbose = false;
%
% or logged to disk by specifying a log file name
cGen.logfile = 'roblog.txt'

%%
% The output variable |symExp| now contains the symbolic expression for the
% forward kinematics. This expression is the same as would be obtained by
% the following code.
symp3 = p3.sym;
q = symp3.gencoords;
symExpDir = symp3.fkine(q)

% So we have basically two ways for deriving symbolic expressions using the
% Robotics Toolbox.

%%
% The difference is that in addition the functional output the symbolic
% expression has now been saved to disk along with the generated m-code and
% Simulink blocks. The storage directory is given in |cGen.basepath|, which
% we now add to our search path.
addpath(cGen.basepath)
ls(cGen.basepath)
%%
% The m-code is contained in a specialized robot class.
ls(cGen.robjpath)

%% Using the generated m-code
% The |mdl_puma560_3| robot definition script defines some special joint
% configurations:
%
% * qz         zero joint angle configuration
% * qr         vertical 'READY' configuration
% * qstretch   arm is stretched out in the X direction
% * qn         arm is at a nominal non-singular configuration
%
% The use of the symbolic expressions and generated code will be
% exemplified in the following based on the zero joint angle configuration.
%
%
qz
%%
% With the generic version of the fkine function from the |SerialLink|
% class we would compute the forward kinematics as follows:
tic; Tz1 = p3.fkine(qz); t1 = toc
%%
% In order to use the generated robot specific m-functions we add them to
% the search path and instantiate a new robot object.
addpath(cGen.basepath)
specRob = eval(cGen.getrobfname)
tic; Tz2 = specRob.fkine(qz); t2 = toc

%% Speedup
% The specialized robot version of fkine runs a little faster
% because it only performs the computations necessary for the specific robot.
% The speedup of the generated robot specific m-code becomes even more appearent if we
% repeat the comparison of the execution times for dynamics
% functions such as:
%
% * gravload -> cGen.gengravload
% * inertia  -> cGen.geninertia
% * coriolis -> cGen.gencoriolis
% * invdyn   -> cGen.geninvdyn
%
% This way the specialized m-code can be used to decrease simulation times.
%
%%
% We obtain the exact solution without floating point notation if we use
% the symbolic expression as follows:
tic; Tz1 = subs(symExp, {'q1', 'q2', 'q3'},qz); toc
%%
% This is however more time consuming. Most probably we might use the
% symbolic expressions for algorithm development, controller design,
% stability proofs as well as analysis, system identification or teaching.
%
% It is also possible to get the symbolic expressions for the homogenous
% transformations of up to each individual joint. This has been found to be
% useful for example for during derivation of analytical inverse kinematics
% solutions. See the documentation of genfkine for details.
%

%% C-Code generation
% Since Release 9.9 the RTB is able to also generate ready to use C-code.
% You can enable C-code generation by activating the CodeGenerator property
% flag |cGen.genccode|:
cGen.genccode = true;
%%
% Now all higher level generator methods (|cGen.genfkine|, |cGen.geninvdyn|
% etc. ...) also produce .c and .h files. They are written to the directory
% specified by the |cGen.ccodepath| property. You can use the C-files in
% your projects outside the MATLAB world. The header files are documented
% and compatible with Doxygen.
%
% Instead of using the higher level generator methods, you can also
% directly call the C-code generation routine for the model code of your
% choice. In the following we complement the previously generated m-functions
% for the forward kinematics by their C-equivalent:
cGen.genccodefkine;
disp('Generated C-headers:')
ls(fullfile(cGen.ccodepath,'include'))
disp('Generated C-definitions:')
ls(fullfile(cGen.ccodepath,'src'))

%% Generating C-MEX functions
% We can use the generated C-code outside the MATLAB world and use it in
% arbitrary C-applications. In addition we can also benefit from it inside
% the MATLAB world by means of C-MEX functions. The automated generation of
% C-MEX functions is controlled by the CodeGenerator flag properties
% |cGen.genmex| and |cGen.compilemex|:
cGen.genmex = true;
%%
% Now all higher level generator methods (|cGen.genfkine|, |cGen.geninvdyn|
% etc. ...) also produce C-MEX files. The MEX files are stored in the class
% directory |cGen.robjpath| of the specialized robot object also incorporating the
% m-functions we generated before.
%%
% By default the flag compilemex is active. This means that the
% CodeGenerator always compiles the generated MEX function after generation.
% We require an installed C-compiler and our MATLAB MEX environment being
% configured properly. See the MATLAB documentation ond MEX files for
% details. In order to proceed with this demo in the case where we do not
% have this prerequisites, we now deactivate the automatic generation:
cGen.compilemex = false;
%%
% Nevhertheless, we can create the C-MEX code for the forward kinematics
% and inspect the |cGen.robjpath| directory:
cGen.genmexfkine
disp('Robot object directory with new MEX source file fkine.c:')
ls(cGen.robjpath)
%%
% The readily compiled MEX functions will shadow the previously generated
% m-functions. The function calls as such remain identical. Using the
% specialized robot object with MEX files we experience an additional and
% substantial computation speed up compared to the robot specific m-code
% as well as the generic rne functions (both, m and MEX version).
%

%% Inheritance
% Even though we have not yet generated robot specific code for |SerialLink|
% metods other than |fkine|, we can still use all functionality of
% |SerialLink| objects with our new specialized robot object which inherits
% from |SerialLink|.
J01 = p3.jacob0(qz)
J02 = specRob.jacob0(qz)

%% A look at the generated Simulink blocks
% The Simulink blocks are stored in a Simulink library file. By opening the
% generated Simulink library we can investigate the already optimized robot
% specific code within the blocks.
% The usage of these blocks is also accompanied with a noticable speedup
% compared to the blocks based on generic |SerialLink| objects.
eval(cGen.slib);
snapnow;
%%
% Beyond the speedup for simulations all blocks in the generated library
% may be directly compiled for real-time operating systems such as xPC-Target or
% dSpace systems for model based control of real hardware setups.
% This way we avoid tedious and error prone reimplementation of the model
% on the target hardware.
%

bdclose(cGen.slib);

%% Further information
% For further information on symbolics and code generation see the
% documentation of the |SerialLink| and |CodeGenerator| class.
%
% All generated functions come with their own description headers so that information
% about their usage can be found by typing |help funname|.
%
% A list of all available methods provide by the |CodeGenerator| class can be
% displayed.
methods CodeGenerator

%%
% The same applies to the configurable properties:
properties CodeGenerator

%% Cleanup
% If we whish to clean our disk from all the generated code, we can simply
% remove it from the search path
rmpath(cGen.basepath)
%%
% and purge everything.
cGen.purge(1)
snapnow

author: petercorke

repo: robotics-toolbox-matlab


text:
function tbxStruct=demos
% DEMOS Demo list for The N-way Toolbox for MATLAB

% Copyright (C) 1995-2006  Rasmus Bro & Claus Andersson
% Copenhagen University, DK-1958 Frederiksberg, Denmark, rb@life.ku.dk
%
% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 2 of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
% You should have received a copy of the GNU General Public License along with
% this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
% Street, Fifth Floor, Boston, MA  02110-1301, USA.

if nargout==0,
    demo toolbox 'The N-way Toolbox for MATLAB';
    return;
end

tbxStruct.Name='The N-way Toolbox for MATLAB';
tbxStruct.Type='Toolbox';
tbxStruct.Help= {
' The N-way Toolbox for MATLAB contains state-of-the-art'
' tools for doing multi-way analysis in MATLAB.'};
tbxStruct.DemoList={
'PARAFAC Demo', 'parademo', '',
'Tucker Demo','tuckdemo', '',};

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function rot = rotation(T)
% convertes orthogonal rank 2 tensors into rotations
%
% Syntax
%   rot = rotation(T)
%
% Input
%  T - orthogonal rank 2 @tensor
%
% Output
%  rot - @rotation
%

rot = rotation.nan(length(T));

switch T.rank

  case 2
    for i = 1:length(T)
      rot(i) = mat2quat(T.M(:,:,i));
    end

    rot = rot .* sign(det(T));

  otherwise
    error('Tensor needs to have rank 2!')
end

author: mtex-toolbox

repo: mtex

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


text:
function measure_types = nst_measure_types()
% NST_MEASURE_TYPES return an enumeration of measure types (eg wavelength, Hb)
%
% MEASURE_TYPES = NST_MEASURE_TYPES()
%    MEASURE_TYPES: struct with numerical fields listing all available
%                   channel types:
%                   - MEASURE_TYPES.WAVELENGTH
%                   - MEASURE_TYPES.Hb

measure_types.WAVELENGTH = 1;
measure_types.HB = 2;
end

author: brainstorm-tools


text:
function Z = iminpaint(im,M)
% Z = iminpaint(im,M)
%
% Inputs:
%   im  w by h by c image
%   M  w by h mask
% Outputs:
%   Z  w by h by c image with mask in-painted with Laplacian
%
  if islogical(M)
    assert(numel(M) == size(im,1)*size(im,2));
    if any(size(M)==1)
      M = reshape(M,size(im,1),size(im,2));
    end
  end
  L = fd_laplacian(size(im(:,:,1)));
  b = find(M);
  X = reshape(im,[],size(im,3));
  Z = min_quad_with_fixed(-L,[],b,X(b,:),[],[],struct('force_Aeq_li',true));
  Z = reshape(Z,size(im));
end

author: alecjacobson

repo: gptoolbox


text:
%% BRENTMIN: Brent's minimization method in one dimension
function [xmin,fmin,funccount,varargout] = ...
                                  brentmin(xlow,xupp,Nitmax,tol,f,nout,varargin)
% code taken from
%    § 10.2 Parabolic Interpolation and Brent's Method in One Dimension
%    Press, Teukolsky, Vetterling & Flannery
%    Numerical Recipes in C, Cambridge University Press, 2002
%
% [xmin,fmin,funccout,varargout] = BRENTMIN(xlow,xupp,Nit,tol,f,nout,varargin)
%    Given a function f, and given a search interval this routine isolates
%    the minimum of fractional precision of about tol using Brent's method.
%
% INPUT
% -----
% xlow,xupp:  search interval such that xlow<=xmin<=xupp
% Nitmax:     maximum number of function evaluations made by the routine
% tol:        fractional precision
% f:          [y,varargout{:}] = f(x,varargin{:}) is the function
% nout:       no. of outputs of f (in varargout) in addition to the y value
%
% OUTPUT
% ------
% fmin:      minimal function value
% xmin:      corresponding abscissa-value
% funccount: number of function evaluations made
% varargout: additional outputs of f at optimum
%
% Copyright (c) by Hannes Nickisch 2010-01-10.

if nargin<6, nout = 0; end
varargout = cell(nout,1);

% tolerance is no smaller than machine's floating point precision
tol = max(tol,eps);

% Evaluate endpoints
fa = f(xlow,varargin{:});
fb = f(xupp,varargin{:});
funccount = 2; % number of function evaluations
% Compute the start point
seps = sqrt(eps);
c = 0.5*(3.0 - sqrt(5.0));% golden ratio
a = xlow; b = xupp;
v = a + c*(b-a);
w = v; xf = v;
d = 0.0; e = 0.0;
x = xf; [fx,varargout{:}] = f(x,varargin{:});
funccount = funccount + 1;

fv = fx; fw = fx;
xm = 0.5*(a+b);
tol1 = seps*abs(xf) + tol/3.0;
tol2 = 2.0*tol1;

% Main loop
while ( abs(xf-xm) > (tol2 - 0.5*(b-a)) )
    gs = 1;
    % Is a parabolic fit possible
    if abs(e) > tol1
        % Yes, so fit parabola
        gs = 0;
        r = (xf-w)*(fx-fv);
        q = (xf-v)*(fx-fw);
        p = (xf-v)*q-(xf-w)*r;
        q = 2.0*(q-r);
        if q > 0.0,  p = -p; end
        q = abs(q);
        r = e;  e = d;

        % Is the parabola acceptable
        if ( (abs(p)<abs(0.5*q*r)) && (p>q*(a-xf)) && (p<q*(b-xf)) )

            % Yes, parabolic interpolation step
            d = p/q;
            x = xf+d;

            % f must not be evaluated too close to ax or bx
            if ((x-a) < tol2) || ((b-x) < tol2)
                si = sign(xm-xf) + ((xm-xf) == 0);
                d = tol1*si;
            end
        else
            % Not acceptable, must do a golden section step
            gs=1;
        end
    end
    if gs
        % A golden-section step is required
        if xf >= xm, e = a-xf;    else e = b-xf;  end
        d = c*e;
    end

    % The function must not be evaluated too close to xf
    si = sign(d) + (d == 0);
    x = xf + si * max( abs(d), tol1 );
    [fu,varargout{:}] = f(x,varargin{:});
    funccount = funccount + 1;

    % Update a, b, v, w, x, xm, tol1, tol2
    if fu <= fx
        if x >= xf, a = xf; else b = xf; end
        v = w; fv = fw;
        w = xf; fw = fx;
        xf = x; fx = fu;
    else % fu > fx
        if x < xf, a = x; else b = x; end
        if ( (fu <= fw) || (w == xf) )
            v = w; fv = fw;
            w = x; fw = fu;
        elseif ( (fu <= fv) || (v == xf) || (v == w) )
            v = x; fv = fu;
        end
    end
    xm = 0.5*(a+b);
    tol1 = seps*abs(xf) + tol/3.0; tol2 = 2.0*tol1;

    if funccount >= Nitmax
        % typically we should not get here
        % warning(sprintf(['Maximum number of iterations (%d) exceeded:', ...
        %                  'precision is not guaranteed'],Nitmax))
        % fprintf('[%1.3f,%1.3f,%1.3f]\n',xlow,xf,xupp)
        break
    end
end % while

% check that endpoints are less than the minimum found
if ( (fa < fx) && (fa <= fb) )
    xf = xlow; fx = fa;
elseif fb < fx
    xf = xupp; fx = fb;
end
fmin = fx;
xmin = xf;

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
 function [xfbp, sino] = em_fbp(sg, ig, yi, ci, ri, varargin)
%function [xfbp, sino] = em_fbp(sg, ig, yi, ci, ri, [options])
% Emission FBP reconstruction from Poisson measurements
% model: Y_i ~ Poisson(c_i [G x]_i + r_i)
% in
%	sg		sino_geom()
%	ig		image_geom()
%	yi		transmission sinogram
%	ci		calibration factors
%	ri		background (randoms, scatter, crosstalk, etc)
%	ci,ri:		optional (can use empty matrices)
%	yi,ci,ri	must have identical dimensions
% option
%	'kernel'	apodization filter kernel (default: [1/3 1/3 1/3])
% out
%	x [np]		image estimate
%	sino [nb,na]	filtered sinogram
%
% Copyright Apr 2000, Jeff Fessler, The University of Michigan

if nargin < 3, ir_usage, end

if ~isvar('ci') || isempty(ci)
	ci = sg.ones;
end
if ~isvar('ri') || isempty(ri)
	ri = sg.zeros;
end

arg.kernel = ones(3,1)/3;
arg = vararg_pair(arg, varargin);

eml_check(yi, ci, ri, 'fbp');
[nb na nz] = size(yi);

tmp = fbp2(sg, ig);
xfbp = ig.zeros('nz', nz);
for iz=1:nz
	proj = (yi(:,:,iz) - ri(:,:,iz)) ./ ci(:,:,iz);
	if any(size(arg.kernel) ~= 1)
		proj = conv2(proj, arg.kernel, 'same'); % filter
	end
	xfbp(:,:,iz) = fbp2(proj, tmp);
	xfbp(:,:,iz) = xfbp(:,:,iz) .* ig.mask;
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%% elephant
% Below is a demonstration of the features of the |elephant| function

%%
clear; close all; clc;

%% Syntax
% |[F,V]=elephant;|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function suppData_ClinicalRF(pathSUPP,nBoot,nSplit,seeds,matlabPATH)

startpath = pwd;

cd(pathSUPP), load('training'), load('testing')
nameOutcomes = {'Locoregional','Distant','Death'}; nOutcomes = numel(nameOutcomes);
testSplit = 1/3; % Proportion of test cases in the splits
comb =  {{[1,2,3],...      % Loco        : Age, Subtype, T_Stage
        [1,2,5],...        % Loco        : Age, Subtype, TNM_Stage
        [1,2,3,4]},...     % Loco        : Age, Subtype, T_Stage, N_Stage
        {[1,2,4],...       % Distant     : Age, Subtype, N_Stage
        [1,2,5],...        % Distant     : Age, Subtype, TNM_Stage
        [1,2,3,4]},...     % Distant     : Age, Subtype, T_Stage, N_Stage
        {[1,2,3],...       % Death       : Age, Subtype, T_Stage
        [1,2,4],...        % Death       : Age, Subtype, N_Stage
        [1,2,5],...        % Death       : Age, Subtype, TNM_Stage
        [1,2,3,4]}};       % Death       : Age, Subtype, T_Stage, N_Stage
testCost = 0.4:0.1:2.5; % Emphasis factor on positive instances during random forest training
time = 60;
mkdir('batchLog_RFclinical'), cd('batchLog_RFclinical'), pathBatch = pwd;

% PRODUCE BATCH COMPUTATIONS
seed = seeds(1); batch = 0;
save('workspace','pathSUPP','nameOutcomes','comb','testCost','nSplit','testSplit','nBoot','seed'), pause(2);
for o = 1:nOutcomes
    batch = batch + 1;
    nameScript = ['batch',num2str(batch),'_script.m'];
    fid = fopen(nameScript,'w');
    fprintf(fid,'load(''workspace'')\n');
    fprintf(fid,['[indClinic,bestCost] = estimateClinicalRF(pathSUPP,nameOutcomes{',num2str(o),'},comb{',num2str(o),'},testCost,nSplit,testSplit,nBoot,seed);\n']);
    fprintf(fid,['save(''result_batch',num2str(batch),''',''indClinic'',''bestCost'')\n']);
    fprintf(fid,['system(''touch batch',num2str(batch),'_end'');\n']);
    fprintf(fid,'clear all');
    fclose(fid);
    system([matlabPATH,' -nojvm -nodisplay -nodesktop -nosplash < ',nameScript,' >& ',nameScript(1:end-1),'log &']);
end

% WAITING LOOP
waitBatch(pathBatch,time,nOutcomes)
delete('workspace.mat')

% GROUP RESULTS
cd(pathBatch)
batch = 0; parameters = struct;
for o = 1:nOutcomes
    batch = batch + 1;
    load(['result_batch',num2str(batch)]) % Variables 'indClinic' and 'bestCost' gets out of there
    parameters.clinical.(nameOutcomes{o}) = indClinic;
    parameters.cost.(nameOutcomes{o}) = bestCost;
    delete(['result_batch',num2str(batch),'.mat'])
    clear indClinic bestCost
end
cd(pathSUPP), save('clinicalPARAM','parameters')

% COMPUTE THE RANDOM FORESTS
for o = 1:nOutcomes
    if strcmp(nameOutcomes{o},'DeathSign')
        outcome = training.outcomes.Death;
    else
        outcome = training.outcomes.(nameOutcomes{o});
    end
    indClinic = parameters.clinical.(nameOutcomes{o});
    cost = parameters.cost.(nameOutcomes{o});
    tableTrain = training.clinical.table(:,indClinic);
    cat = logical(training.clinical.categories(indClinic));
    rng(seeds(2)), [RF] = trainRF_table(tableTrain,outcome,cat,nBoot,cost);
    RF = compact(RF); % Compact version
    save(['RFclinical_',nameOutcomes{o}],'RF')
end

% TEST THE RANDOM FORESTS
cd(pathSUPP), load('clinicalPARAM')
for o = 1:nOutcomes
    if strcmp(nameOutcomes{o},'DeathSign')
        outcome = testing.outcomes.Death;
        time = testing.timeToEvents.Death;
    else
        outcome = testing.outcomes.(nameOutcomes{o});
        time = testing.timeToEvents.(nameOutcomes{o});
    end
    censoring = 1 - outcome;
    results = struct;
    RF = load(['RFclinical_',nameOutcomes{o}]); RF = struct2cell(RF); RF = RF{1};
    indClinic = parameters.clinical.(nameOutcomes{o});
    tableTest = testing.clinical.table(:,indClinic);
    [prob] = predictRF(tableTest,RF);
    results.probResponse = prob;
    [AUC,sensitivity,specificity,accuracy] = calcPerformMetrics(prob,outcome,0.5);
    CI = calcCI(prob,time,censoring);
    results.AUC = AUC; results.Sensitivity = sensitivity; results.Specificity = specificity; results.Accuracy = accuracy;
    results.CI = CI;
    save(['testResultsRFclinical_',nameOutcomes{o}],'results'), clear results
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function gauss_seidel_test02 ( )

%*****************************************************************************80
%
%% GAUSS_SEIDEL_TEST02 tests GAUSS_SEIDEL2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GAUSS_SEIDEL_TEST02:\n' );

  it_num = 400;
  n = 20;

  x_exact = ( 1 : n )';

  l = sparse ( 2:n,   1:n-1, -1.0, n, n );
  d = sparse ( 1:n,   1:n,    2.0, n, n );
  u = sparse ( 1:n-1, 2:n,   -1.0, n, n );

  a = l + d + u;

  b = a * x_exact;

  x = zeros ( n, 1 );
  x_plot(1:n,1) = x;

  step = 1 : it_num + 1;
  e = nan ( it_num+1, 1 );
  xm = nan ( it_num+1, 1 );

  e(1,1) = ( norm ( a * x - b ) ).^2;

  for it = 1 : it_num

    x_new = gauss_seidel2 ( n, l, d, u, b, x );

    e(it+1,1) = ( norm ( a * x_new - b ) ).^2;
    x_plot(1:n,it+1) = x_new;
%
%  Display the error.
%
    figure ( 1 )
    plot ( step, log ( e ), 'm-*' )
    title ( 'Log (Error^2)' )
    xlabel ( 'Step' )
    ylabel ( 'Error' )
    grid
%
%  Display the motion.
%
    xm(it,1) = sum ( ( x_new(:) - x(:) ).^2 ) / n;

    figure ( 2 )
    plot ( step, log ( xm ), 'm-*' )
    title ( 'Log (Average generator motion)' )
    xlabel ( 'Step' )
    ylabel ( 'Energy' )
    grid
%
%  Update the solution
%
    x = x_new;

  end
%
%  Plot the evolution of the locations of the generators.
%
  figure ( 3 )

  y = ( 0 : it_num );
  for k = 1 : n
    plot ( x_plot(k,1:it_num+1), y )
    hold on;
  end
  grid on
  hold off;

  title ( 'Generator evolution.' );
  xlabel ( 'Generator positions' );
  ylabel ( 'Iterations' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function result = out_of_FOV_mask(img, geo, dia)
% ROI_mask=Cylindrical Mask to Input image matrix
% Date: 2015-10-22
% Author: Yi Du (yi.du@hotmail.com)
% Modified: Ander Biguri

[Length, Width, ~] = size(img);

% diameter -> radius
Rad = dia*0.5;

[XX, YY] = meshgrid((1:Width)-(Width+1)/2-geo.offOrigin(2)./geo.dVoxel(2), (1:Length)-(Length+1)/2-geo.offOrigin(1)./geo.dVoxel(1));
XX= XX.^2;
YY = YY.^2;
R2 = Rad^2;

mask = (XX+YY)-R2;
mask(mask>0)=0;
mask(mask<0)=1;

result = bsxfun(@times,img,mask);

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function SavePTS(pts, ptsPath)
fPTS = fopen(ptsPath, 'w');
numPoints = size(pts, 1);
fprintf(fPTS, 'version: 1.0\r\nn_points: %d\r\n{\r\n', numPoints);
for j = 1 : numPoints
    fprintf(fPTS, '%f %f\r\n', pts(j, 1), pts(j, 2));
end
fprintf(fPTS, '}\r\n');
fclose(fPTS);
end

author: jiankangdeng

repo: MenpoBenchmark


text:
function c = gamvbcost(x,logx,apost,bpost,aprior,bprior)
% function c = gamvbcost(x,logx,apost,bpost,aprior,bprior)
%
% Calculates c = -KL(q||p) = <log p(X)> - <log q(X)>
% where the expectation <.> is taken over q(X). This is used for
% calculating the lower bound of the marginal loglikelihood in VB
% models.
%
% Let X be a Gamma distributed variable:
% p(X) = G(aprior,bprior)
% q(X) = G(apost,bpost)
% <X> = x
% <log X> = logx

% Cost from prior
c = aprior*log(bprior) - gammaln(aprior) + (aprior-1)*logx - bprior*x;
% Cost from q-posterior
c = c + gamentropy(apost,bpost);


text:
function [data] = ft_datatype_raw(data, varargin)

% FT_DATATYPE_RAW describes the FieldTrip MATLAB structure for raw data
%
% The raw datatype represents sensor-level time-domain data typically
% obtained after calling FT_DEFINETRIAL and FT_PREPROCESSING. It contains
% one or multiple segments of data, each represented as Nchan X Ntime
% arrays.
%
% An example of a raw data structure with 151 MEG channels is
%
%          label: {151x1 cell}      the channel labels (e.g. 'MRC13')
%           time: {1x266 cell}      the timeaxis [1*Ntime double] per trial
%          trial: {1x266 cell}      the numeric data [151*Ntime double] per trial
%     sampleinfo: [266x2 double]    the begin and endsample of each trial relative to the recording on disk
%      trialinfo: [266x1 double]    optional trigger or condition codes for each trial
%            hdr: [1x1 struct]      the full header information of the original dataset on disk
%           grad: [1x1 struct]      information about the sensor array (for EEG it is called elec)
%           cfg: [1x1 struct]       the configuration used by the function that generated this data structure
%
% Required fields:
%   - time, trial, label
%
% Optional fields:
%   - sampleinfo, trialinfo, grad, elec, opto, hdr, cfg
%
% Deprecated fields:
%   - fsample
%
% Obsoleted fields:
%   - offset
%
% Revision history:
%
% (2011/latest) The description of the sensors has changed, see FT_DATATYPE_SENS
% for further information.
%
% (2010v2) The trialdef field has been replaced by the sampleinfo and
% trialinfo fields. The sampleinfo corresponds to trl(:,1:2), the trialinfo
% to trl(4:end).
%
% (2010v1) In 2010/Q3 it shortly contained the trialdef field which was a copy
% of the trial definition (trl) is generated by FT_DEFINETRIAL.
%
% (2007) It used to contain the offset field, which correcponds to trl(:,3).
% Since the offset field is redundant with the time axis, the offset field is
% from now on not present any more. It can be recreated if needed.
%
% (2003) The initial version was defined
%
% See also FT_DATATYPE, FT_DATATYPE_COMP, FT_DATATYPE_TIMELOCK, FT_DATATYPE_FREQ,
% FT_DATATYPE_SPIKE, FT_DATATYPE_SENS

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

% get the optional input arguments, which should be specified as key-value pairs
version       = ft_getopt(varargin, 'version', 'latest');
hassampleinfo = ft_getopt(varargin, 'hassampleinfo', 'ifmakessense'); % can be yes/no/ifmakessense
hastrialinfo  = ft_getopt(varargin, 'hastrialinfo',  'ifmakessense'); % can be yes/no/ifmakessense

% do some sanity checks
assert(isfield(data, 'trial') && isfield(data, 'time') && isfield(data, 'label'), 'inconsistent raw data structure, some field is missing');
assert(length(data.trial)==length(data.time), 'inconsistent number of trials in raw data structure');
for i=1:length(data.trial)
  assert(size(data.trial{i},2)==length(data.time{i}), 'inconsistent number of samples in trial %d', i);
  assert(size(data.trial{i},1)==length(data.label), 'inconsistent number of channels in trial %d', i);
end

% convert it into true/false
if isequal(hassampleinfo, 'ifmakessense')
  hassampleinfo = makessense(data, 'sampleinfo');
else
  hassampleinfo = istrue(hassampleinfo);
end
if isequal(hastrialinfo, 'ifmakessense')
  hastrialinfo = makessense(data, 'trialinfo');
else
  hastrialinfo = istrue(hastrialinfo);
end

if strcmp(version, 'latest')
  version = '2011';
end

if isempty(data)
  return;
end

switch version
  case '2011'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % ensure that the sensor structures are up to date
    if isfield(data, 'grad')
      data.grad = ft_datatype_sens(data.grad);
    end
    if isfield(data, 'elec')
      data.elec = ft_datatype_sens(data.elec);
    end
    if isfield(data, 'opto')
      data.opto = ft_datatype_sens(data.opto);
    end

    if ~isfield(data, 'fsample')
      for i=1:length(data.time)
        if length(data.time{i})>1
          data.fsample = 1/mean(diff(data.time{i}));
          break
        else
          data.fsample = nan;
        end
      end
      if isnan(data.fsample)
        ft_warning('cannot determine sampling frequency');
      end
    end

    if isfield(data, 'offset')
      data = rmfield(data, 'offset');
    end

    % the trialdef field should be renamed into sampleinfo
    if isfield(data, 'trialdef')
      data.sampleinfo = data.trialdef;
      data = rmfield(data, 'trialdef');
    end

    if (hassampleinfo && ~isfield(data, 'sampleinfo')) || (hastrialinfo && ~isfield(data, 'trialinfo'))
      % try to reconstruct the sampleinfo and trialinfo
      data = fixsampleinfo(data);
    end

    if ~hassampleinfo && isfield(data, 'sampleinfo')
      data = rmfield(data, 'sampleinfo');
    end

    if ~hastrialinfo && isfield(data, 'trialinfo')
      data = rmfield(data, 'trialinfo');
    end

  case '2010v2'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ~isfield(data, 'fsample')
      data.fsample = 1/mean(diff(data.time{1}));
    end

    if isfield(data, 'offset')
      data = rmfield(data, 'offset');
    end

    % the trialdef field should be renamed into sampleinfo
    if isfield(data, 'trialdef')
      data.sampleinfo = data.trialdef;
      data = rmfield(data, 'trialdef');
    end

    if (hassampleinfo && ~isfield(data, 'sampleinfo')) || (hastrialinfo && ~isfield(data, 'trialinfo'))
      % try to reconstruct the sampleinfo and trialinfo
      data = fixsampleinfo(data);
    end

    if ~hassampleinfo && isfield(data, 'sampleinfo')
      data = rmfield(data, 'sampleinfo');
    end

    if ~hastrialinfo && isfield(data, 'trialinfo')
      data = rmfield(data, 'trialinfo');
    end

  case {'2010v1' '2010'}
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ~isfield(data, 'fsample')
      data.fsample = 1/mean(diff(data.time{1}));
    end

    if isfield(data, 'offset')
      data = rmfield(data, 'offset');
    end

    if ~isfield(data, 'trialdef') && hascfg
      % try to find it in the nested configuration history
      data.trialdef = ft_findcfg(data.cfg, 'trl');
    end

  case '2007'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ~isfield(data, 'fsample')
      data.fsample = 1/mean(diff(data.time{1}));
    end

    if isfield(data, 'offset')
      data = rmfield(data, 'offset');
    end

  case '2003'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ~isfield(data, 'fsample')
      data.fsample = 1/mean(diff(data.time{1}));
    end

    if ~isfield(data, 'offset')
      data.offset = zeros(length(data.time),1);
      for i=1:length(data.time);
        data.offset(i) = round(data.time{i}(1)*data.fsample);
      end
    end

  otherwise
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ft_error('unsupported version "%s" for raw datatype', version);
end

% Numerical inaccuracies in the binary representations of floating point
% values may accumulate. The following code corrects for small inaccuracies
% in the time axes of the trials. See http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1390
data = fixtimeaxes(data);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function data = fixtimeaxes(data)

if ~isfield(data, 'fsample')
  fsample = 1/mean(diff(data.time{1}));
else
  fsample = data.fsample;
end

begtime   = nan(1, length(data.time));
endtime   = nan(1, length(data.time));
numsample = zeros(1, length(data.time));
for i=1:length(data.time)
  if ~isempty(data.time{i})
    begtime(i) = data.time{i}(1);
    endtime(i) = data.time{i}(end);
  end
  numsample(i) = length(data.time{i});
end

% compute the differences over trials and the tolerance
tolerance     = 0.01*(1/fsample);
begdifference = abs(begtime-begtime(1));
enddifference = abs(endtime-endtime(1));

% check whether begin and/or end are identical, or close to identical
begidentical  = all(begdifference==0);
endidentical  = all(enddifference==0);
begsimilar    = all(begdifference < tolerance);
endsimilar    = all(enddifference < tolerance);

% Compute the offset of each trial relative to the first trial, and express
% that in samples. Non-integer numbers indicate that there is a slight skew
% in the time over trials. This works in case of variable length trials.
offset = fsample * (begtime-begtime(1));
skew   = abs(offset - round(offset));

% try to determine all cases where a correction is needed
% note that this does not yet address all possible cases where a fix might be needed
needfix = false;
needfix = needfix || ~begidentical && begsimilar;
needfix = needfix || ~endidentical && endsimilar;
needfix = needfix || ~all(skew==0) && all(skew<0.01);

% if the skew is less than 1% it will be corrected
if needfix
  ft_warning('correcting numerical inaccuracy in the time axes');
  for i=1:length(data.time)
    % reconstruct the time axis of each trial, using the begin latency of
    % the first trial and the integer offset in samples of each trial
    data.time{i} = begtime(1) + ((1:numsample(i)) - 1 + round(offset(i)))/fsample;
  end
end

author: spm


text:
function P = polyhedron(C)
%POLYHEDRON (Overloaded)

P = Polyhedron(lmi(C));

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% test for active contours
%
%   Copyright (c) 2007 Gabriel Peyre

options.bound = 'per';
n = 32;
d = ones(n);
d = rand(n);
d = repmat(d,[1 1 2]);
D = zeros(n^2);
for i=1:n^2
    a = zeros(n); a(i)=1;
    a = divgrad(divgrad(a,options)./d,options);
    D(:,i) = a(:);
end

n = 150;

path(path, 'data/');
path(path, 'toolbox/');

motion = 'affine';
motion = 'errosion';
motion = 'chan-vese';
motion = 'mean';
motion = 'snake';

% load image for active contour
name = 'chan-vese';
name = 'disk';
name = 'brain';
M = rescale( sum( load_image(name, n), 3) );

options.solver = 'grad';
options.solver = 'cg';

%% load original shape
options.null = 0;
if strcmp(motion, 'snake')
    namec = 'square';
    options.width = 0.95*n;
elseif strcmp(motion, 'chan-vese')
    namec = 'small-disks';
else
    namec = 'circlerect1';
    namec = 'square';
    namec = 'circlerect2';
end
D0 = compute_levelset_shape(namec, n);

options.center = [0.15 0.15]*n;
options.radius = 0.1*n;
D1 = compute_levelset_shape('circle', n,options);

%% set up the parameters
options.Tmax = 1000;
options.redistance_freq = 150;
options.M = M;
options.dt = 0.1;
options.dt = 3;
options.display_freq = 5;
options.nb_svg = 100;
if strcmp(motion, 'snake')
    options.Tmax = 2000;
    options.redistance_freq = 30;
    options.dt = 1;
    % compute edge-based energy
    sigma = 4; % blurring size
    G = divgrad( perform_blurring(M,sigma) );
    G = sum( G.^2, 3);
    eta = 0.01;
    E = 1 ./ (eta + G);
    options.E = rescale(E, 0.7, 1);
elseif strcmp(motion, 'chan-vese')
    options.redistance_freq = 30;
    options.E = M;
    options.lambda = 0.8;
    options.update_c = 0;
    options.c1 = 0;
    options.c2 = 0.5;
    options.dt = 0.2;
    options.dt = 2;
end

D = perform_active_contour(D0, motion, options);


text:
function [theta,e,hy,hp] = GLM_covComp(y,X,Qy,Qp,my,mp)
% wraps GLM with multiple covariance components estimation into VBA
% function [theta,e] = GLM_covComp(y,X,Qy,Qp,my,mp)
% IN:
%   - y: px1 data vector
%   - X: pxn_phi GLM design matrix
%   - Qy: nqyx1 cell array of data covariance components
%   - Qp: nqpx1 cell array of GLM coef covariance components
%   - my: prior mean on data precision hyperparameters in log-space
%   - mp: prior mean on GLM coef precision hyperparameters in log-space
% OUT:
%   - theta: n_phix1 vector of estimated GLM coef
%   - e: px1 vector of estimated residuals

nqy = length(Qy);
nqp = length(Qp);

[ny,np] = size(X);

dim.p = size(y,1);
dim.n_phi = (np+1)*nqp + (ny+1)*nqy;
dim.n_theta = 0;
dim.n = 0;
dim.n_t = 1;

priors.muPhi = zeros(dim.n_phi,1);
priors.SigmaPhi = zeros(dim.n_phi,dim.n_phi);

ind = cell(4,1);
last = 0;
for i=1:nqy
    ind{1}(i) = last+1;
    ind2i = ind{1}(i)+1:ind{1}(i)+ny;
    ind{2} = [ind{2};ind2i];
    last = ind2i(end);
    priors.muPhi(ind{1}(i)) = my(i);
    priors.SigmaPhi(ind{1}(i),ind{1}(i)) = 4;
    priors.SigmaPhi(ind2i,ind2i) = Qy{i};
end
for j=1:nqp
    ind{3}(j) = last+1;
    ind4j = ind{3}(j)+1:ind{3}(j)+np;
    ind{4} = [ind{4};ind4j];
    last = ind4j(end);
    priors.muPhi(ind{3}(j)) = mp(j);
    priors.SigmaPhi(ind{3}(j),ind{3}(j)) = 4;
    priors.SigmaPhi(ind4j,ind4j) = Qp{j};
end
priors.a_sigma = 1e8;
priors.b_sigma = 1e0;

inG.X = X;
inG.ind = ind;
inG.dim = struct('ny',ny,'np',np,'nqy',nqy,'nqp',nqp);

options = struct(...
    'priors',priors,...
    'inG',inG,...
    'updateHP',0,...
    'MaxIter',128,...
    'GnMaxIter',64);
g_fname = @g_wrapGLM;

[posterior,out] = VBA_NLStateSpaceModel(y,[],[],g_fname,dim,options);

[theta,e,hy,hp] = collapseGLM(posterior.muPhi,inG);

function [gx] = g_wrapGLM(x,P,u,in)
[theta,e] = collapseGLM(P,in);
gx = in.X*theta + e;

function [theta,e,hy,hp] = collapseGLM(P,in)
P1 = P(in.ind{1});
P2 = reshape(P(in.ind{2}),in.dim.ny,in.dim.nqy);
P3 = P(in.ind{3});
P4 = reshape(P(in.ind{4}),in.dim.np,in.dim.nqp);
hp = exp(P3);
theta = P4*hp;
hy = exp(P1);
e = P2*hy;

author: MBB-team

repo: VBA-toolbox


text:
function [ h ] = ULA_fun( phi ,N)
    h=exp(1j*pi*sin(phi).*(0:N-1)');
end

author: guohuayan


text:
classdef CEC2017_F5 < PROBLEM
% <single> <real> <constrained>
% CEC'2017 constrained optimization benchmark problem

%------------------------------- Reference --------------------------------
% G. Wu, R. Mallipeddi, and P. N. Suganthan, Problem definitions and
% evaluation criteria for the CEC 2017 competition on constrained real-
% parameter optimization, National University of Defense Technology, China,
% 2016.
%------------------------------- 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".
%--------------------------------------------------------------------------

    properties
        O;      % Optimal decision vector
        Mat1;	% Rotation matrices
        Mat2;
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2017.mat'),'Data');
            obj.O = Data{5}.o;
            obj.M = 1;
            if isempty(obj.D) || obj.D < 30
                obj.D    = 10;
                obj.Mat1 = Data{5}.M1_10;
                obj.Mat2 = Data{5}.M2_10;
            elseif obj.D < 50
                obj.D    = 30;
                obj.Mat1 = Data{5}.M1_30;
                obj.Mat2 = Data{5}.M2_30;
            elseif obj.D < 100
                obj.D    = 50;
                obj.Mat1 = Data{5}.M1_50;
                obj.Mat2 = Data{5}.M2_50;
            else
                obj.D    = 100;
                obj.Mat1 = Data{5}.M1_100;
                obj.Mat2 = Data{5}.M2_100;
            end
            obj.lower    = zeros(1,obj.D) - 10;
            obj.upper    = zeros(1,obj.D) + 10;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopObj = sum(100*(Z(:,1:end-1).^2-Z(:,2:end)).^2+(Z(:,1:end-1)-1).^2,2);
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            Y = Z*obj.Mat1';
            W = Z*obj.Mat2';
            PopCon(:,1) = sum(Y.^2-50*cos(2*pi*Y)-40,2);
            PopCon(:,2) = sum(W.^2-50*cos(2*pi*W)-40,2);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%% Examples 10.17, 10.26, and 10.38: Benes-Daum and EKF/ERTS examples
%
% Copyright:
%   2018 - Simo Särkkä and Arno Solin
%
% License:
%   This software is provided under the MIT License. See the accompanying
%   LICENSE file for details.

%%
% Benes and Benes-Daum filter tests
%
    % Check that we have the correct normalization constants and
    % moments. We should have (for x = x(t) with t >= t_k)
    %
    % p(x | Y_k) =
    % 1/sqrt{2 pi P} exp(-(x-m)^2/(2P)) exp(-1/2 P) cosh(x) / cosh(m)
    %
    % and
    %
    %    E[x] = m + tanh(m) P
    %  E[x^2] = m^2 + P^2 + 2 m tanh(m) P + P
    %    V[x] = P + [1 - tanh(m)^2] P^2

    dx = 0.001;
    xx = -15:dx:15;
    dens = @(x,m,P) 1/sqrt(2*pi*P)*exp(-P/2)/cosh(m) * exp(-(x-m).^2/(2*P)) .* cosh(x);

    m = 1;
    P = 5;
    pp = dens(xx,m,P);
    figure(1); clf; plot(xx,pp);

    sum(pp .* dx)
    [sum(xx .* pp * dx)        m+tanh(m)*P]
    [sum(xx.^2 .* pp * dx)  m^2+P^2+2*m*tanh(m)*P+P]
    [sum((xx-m-tanh(m)*P).^2 .* pp * dx) P+(1-tanh(m)^2)*P^2]

%%
% Generate data
%
  % Lock random seed
  if exist('rng') % Octave doesn't have rng
    rng(1)
  else
    randn('state',9);
  end

    steps = 500;
    meas_mod = 20;
    dt = 0.01;
    R = 0.1;

    X = [];
    Y = [];
    T = [];
    x0 = 0;
    x = x0;
    t = 0;
    for k=1:steps
        x = x + tanh(x) * dt + sqrt(dt)*randn;
        if rem(k,meas_mod) == 0
            y = x + sqrt(R)*randn;
        else
            y = NaN;
        end
        t = t + dt;
        X(k) = x;
        Y(k) = y;
        T(k) = t;
    end

    figure(1); clf
    plot(T,X,T,Y,'o');

%%
% Benes-Daum filter
%
    EB = zeros(1,length(Y));
    VB = zeros(1,length(Y));
    MMB = zeros(1,length(Y));
    PPB = zeros(1,length(Y));

    mb0 = 0.1;
    Pb0 = 1;

    mb = mb0;
    Pb = Pb0;

    m0 = mb + Pb * tanh(mb);
    P0 = Pb + (1-tanh(mb)^2) * Pb^2;

    for k=1:length(Y)
        %
        % The Benes-Daum filter
        %
        Pb = Pb + dt;
        if ~isnan(Y(k))
            mb = mb + Pb / (R + Pb)*(Y(k) - mb);
            Pb = Pb - Pb^2/(R + Pb);
        end
        MMB(:,k) = mb;
        PPB(:,k) = Pb;
        EB(:,k) = mb + Pb * tanh(mb);
        VB(:,k) = Pb + (1-tanh(mb)^2) * Pb^2;
    end

    figure(1); clf;
    subplot(2,1,1);
      plot(T,MMB,T,EB);
    subplot(2,1,2);
      plot(T,PPB,T,VB);

    rmse_b = sqrt(mean((EB - X).^2))

%%
% Determine the quantiles of Benes solution
%
    QQ1 = zeros(size(MMB));
    QQ2 = zeros(size(MMB));
    for k=1:length(MMB)
        pp = dens(xx,MMB(:,k),PPB(:,k));
        cu = cumsum(pp ./ sum(pp));
        ind1 = find(cu >= 2.5/100,1);
        ind2 = find(cu >= 97.5/100,1);
        QQ1(k) = xx(ind1);
        QQ2(k) = xx(ind2);
    end

    figure(1); clf;
    plot(T,EB,T,QQ1,'--',T,QQ2,'--');

%%
% EKF solution to the same problem
%
    me = m0;
    Pe = P0;
    MMe = zeros(1,length(Y));
    PPe = zeros(1,length(Y));
    ekf_steps = 5;

    for k=1:length(Y)
        dt2 = dt/ekf_steps;
        for i=1:ekf_steps
            Pe = Pe + dt2 * (2 * (1 - tanh(me)^2) * Pe + 1);
            me = me + dt2 * tanh(me);
        end
        if ~isnan(Y(k))
            me = me + Pe / (R + Pe) * (Y(k) - me);
            Pe = Pe - Pe^2/(R + Pe);
        end
        MME(:,k) = me;
        PPE(:,k) = Pe;
    end

    rmse_e = sqrt(mean((MME - X).^2))

    figure(1); clf;
    subplot(2,2,1);
    plot(T,EB,T,QQ1,'--',T,QQ2,'--');
    subplot(2,2,2);
    plot(T,MME,T,MME-1.96*sqrt(PPE),'--',T,MME+1.96*sqrt(PPE),'--');

    subplot(2,2,3);
    plot(T,EB,T,MME)

    subplot(2,2,4);
    plot(T,VB,T,PPE)

%%
% ERTS Smoother
%
    MMS = MME;
    PPS = PPE;
    ms = MME(end);
    Ps = PPE(end);
    for k=length(Y)-1:-1:1
        m = MME(k);
        P = PPE(k);
        dms = dt * tanh(m) + dt * ((1 - tanh(m)^2) * P + 1)/P * (ms - m);
        dPs = dt * 2 * ((1 - tanh(m)^2) * P + 1)/P * Ps - dt;
        ms = ms - dms;
        Ps = Ps - dPs;
        MMS(k) = ms;
        PPS(k) = Ps;
    end

    rmse_s = sqrt(mean((MMS - X).^2))

    figure(1); clf;
    plot(T,MMS,T,MMS-1.96*sqrt(PPS),'--',T,MMS+1.96*sqrt(PPS),'--');

%%
% Plot the final figures
%

    %
    % Benes-Daum
    %
    figure(1); clf; hold on

    % Plot mean
    h0 = plot(T,EB,'-k');
    set(h0,'LineWidth',1);
    set(h0,'Color',[0.5 0.5 0.5]);

    % Plot uncetainty
    h1 = fill([T'; flipud(T')], [QQ1'; flipud(QQ2')],1);
    set(h1,'EdgeColor',[.7 .7 .7],'FaceColor',[.7 .7 .7])

    % Plot signal
    h2 = plot(T,X,'-k');
    set(h2,'LineWidth',1);

    % Plot mean
    h3 = plot(T,EB,'-k');
    set(h3,'LineWidth',1);
    set(h3,'Color',[0.5 0.5 0.5]);

    % Plot observations
    h4 = plot(T,Y,'+k');

    box on
    xlabel('Time, $t$'), ylabel('Signal, $x(t)$')

    legend([h0, h1, h2, h4],'BD mean','BD 95\% quantiles','Signal','Observations')

    %
    % Benes-Daum vs EKF
    %
    figure(2); clf; hold on

    % Plot BD mean
    h1 = plot(T,EB,'-k');
    set(h1,'LineWidth',1);
    set(h1,'Color',[0.5 0.5 0.5]);

    % Plot EKF mean
    h2 = plot(T,MME,'--k');
    set(h2,'LineWidth',1);
    %set(h2,'Color',[0.0 0.0 0.0]);

    % Plot BD uncertainty
    h3 = fill([T'; flipud(T')], [QQ1'; flipud(QQ2')],1);
    set(h3,'EdgeColor',[.7 .7 .7],'FaceColor',[.7 .7 .7])

    % Plot BD mean
    h4 = plot(T,EB,'-k');
    set(h4,'LineWidth',1);
    set(h4,'Color',[0.5 0.5 0.5]);

    % Plot EKF mean
    h5 = plot(T,MME,'--k');
    set(h5,'LineWidth',1);
    set(h5,'Color',[0.0 0.0 0.0]);

    % Plot EKF uncertainty
    h6 = plot(T,MME-1.96*sqrt(PPE),'-.k',T,MME+1.96*sqrt(PPE),'-.k');
    set(h6,'LineWidth',.5);

    box on
    xlabel('Time, $t$'), ylabel('Signal, $x(t)$')
    legend([h1, h2, h3, h6(1)],'BD mean','EKF mean','BD 95\% quantiles','EKF 95\% quantiles')

    %
    % ERTS
    %
    figure(3); clf; hold on

    % Plot mean
    h0 = plot(T,MMS,'-k');
    set(h0,'LineWidth',1);
    set(h0,'Color',[0.5 0.5 0.5]);

    QQ1S = MMS-1.96*sqrt(PPS);
    QQ2S = MMS+1.96*sqrt(PPS);

    % Plot uncetainty
    h1 = fill([T'; flipud(T')], [QQ1S'; flipud(QQ2S')],1);
    set(h1,'EdgeColor',[.7 .7 .7],'FaceColor',[.7 .7 .7])

    % Plot signal
    h2 = plot(T,X,'-k');
    set(h2,'LineWidth',1);

    % Plot mean
    h3 = plot(T,MMS,'-k');
    set(h3,'LineWidth',1);
    set(h3,'Color',[0.5 0.5 0.5]);

    % Plot observations
    h4 = plot(T,Y,'+k');

    box on
    xlabel('Time, $t$'), ylabel('Signal, $x(t)$')

    legend([h0, h1, h2, h4],'ERTS mean','ERTS 95\% quantiles','Signal','Observations')

author: AaltoML

repo: SDE

save_path: github-repos/MATLAB/AaltoML-SDE


text:
% test_all_mex
% make sure all mex files can execute
% by running the internal "check" of each.

list_spline = {
	'BsplCo2GdXMirr', ...
	'BsplCo2GdXTranMirr', ...
	'BsplCo2GdXTranZero', ...
	'BsplCo2GdXZero', ...
	'BsplCo2GdYMirr', ...
	'BsplCo2GdYTranMirr', ...
	'BsplCo2GdYTranZero', ...
	'BsplCo2GdYZero', ...
	'BsplCo2GdZMirr', ...
	'BsplCo2GdZTranMirr', ...
	'BsplCo2GdZTranZero', ...
	'BsplCo2GdZZero', ...
	'BsplCo2ValMirr', ...
	'BsplCo2ValTranMirr', ...
	'BsplCo2ValTranZero', ...
	'BsplCo2ValTranZeroFilt', ...
	'BsplCo2ValZero', ...
	'BsplCo2ValZeroFilt', ...
	'BsplExpand', ...
	'BsplReduce', ...
	'BsplVal2CoMirr', ...
	'BsplVal2CoZero'
};

list = {
	'jf_mex', ...
	'dtft_mex', ...
	'exp_xform_mex', ...
	'mri_exp_mult_mex', ...
...
	'interp1_table_adj_mex', ...
	'interp1_table_mex', ...
	'interp2_table_adj_mex', ...
	'interp2_table_mex', ...
	'interp3_table_adj_mex', ...
	'interp3_table_mex', ...
...
	'penalty_mex', ...
	'rotmex', ...
	'wtfmex', ...
	'f3d_mex', ...
...
	'delaysum1_mex', ...
	'ir_shrink1_mex', ...
	'ir_tridiag_inv_mex', ...
...
	list_spline{:}
};

% check for UM-only mex files
if exist('dd_ge1_mex') == 3
	list{end+1} = 'dd_ge1_mex';
end
if exist('dd_ge2_mex') == 3
	list{end+1} = 'dd_ge2_mex';
end

is_missing = false(1, numel(list));
for ii=1:numel(list)
	mex = list{ii};
%	pr mex
	if exist(mex) ~= 3
		is_missing(ii) = true;
	end
end

if any(is_missing)
	printf(' ')
	printm('These mex files are missing:')
%	pr find(is_missing)
%	list{is_missing}
	pr list(is_missing)
prompt
	list = list(~is_missing);
end

passed = '';
failed = '';
missing = '';
for ii=1:numel(list)
	mex = list{ii};
	pr mex
	if exist(mex) ~= 3
		missing = [missing ' ' mex];
		continue;
	end
	try
		fun = str2func(mex);
		fun('check')
		passed = [passed ' ' mex];
	catch
		failed = [failed ' ' mex];
	end
end

if ~isempty(missing) || ~isempty(failed)

	if ~isempty(missing)
		printm(['\nThese mex files are missing: ' missing])
	end

	if ~isempty(failed)
		printm(['\nThese mex files failed: ' failed])
	end

	if ~isempty(passed)
		printm(['\nThese mex files passed: ' passed])
		printm 'So perhaps some things will still work.'
	end

	printm 'Sorry, you seem to have mex problems. :-('
	printm 'Probably you are a PC user and Windoze is not supported.'
	printm 'Or (in linux) there may be a gcc library version issue?'
	printm 'Or you may have a path problem.'

else
	printm '\n----------------------------------------------------------'
	printm(['All mex files present and passed:\n' passed])
	printm '----------------------------------------------------------\n'
end


text:
function f = fracInt(f, mu)
%FRACINT  Fractional integral of a CHEBFUN.
%   FRACINT(F, MU) gives the order MU fractional integral of a CHEBFUN object F.
%
%   Currently this only supports the situation where F is smooth (i.e., it has
%   no breakpoints or endpoint singularities) and on a finite domain.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% No piecewise support yet:
if ( numel(domain(f)) > 2 )
    error('CHEBFUN:CHEBFUN:fracInt:breakpoints', ...
        'FRACINT does not currently support piecewise functions.');
end

% Extract the fractional part of mu:
mu_int = floor(mu);
mu_frac = mu - mu_int;

% Computer the integer CUMSUMs:
if ( mu_int > 0 )
    % Ensure mu is in [0, 1):
    f = cumsum(f, mu_int);
    f = fracInt(f, mu_frac);
    return
elseif ( mu_frac == 0 )
    % Nothing more to do (non-fractional integral).
    return
end

% Quasimatrix / Array-valued support:
% TODO: This should be overhauled once SINGFUN supports array-valuedness.
m = numColumns(f);
if ( m > 1 )
    % Convert to cell-array of scalar CHEBFUNs.
    f = mat2cell(f);
    % Loop over columns:
    for k = 1:m
        f{k} = fracInt(f{k}, mu_frac);
    end
    % Convert from cell array to quasimatrix.
    f = cat(2-f{1}.isTransposed, f{:});
    return
end

% From here on f is scalar-valued:

% Call BNDFUN/FRACINT():
f.funs{1} = fracInt(f.funs{1}, mu_frac);

% Update the pointValues:
f.pointValues = chebfun.getValuesAtBreakpoints(f.funs);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function varargout = plot(varargin)
% plot method for GIfTI objects
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: plot.m 7513 2019-01-15 15:17:29Z guillaume $

% if ishandle(varargin{1})
%     h = figure(varargin{1});
% else
%     h = figure;
%     %axis equal;
%     %axis off;
%     %camlight;
%     %camlight(-80,-10);
%     %lighting gouraud;
% end
% cameratoolbar;

cdata = [];
ax = [];
if nargin == 1
    this = varargin{1};
    h = gcf;
else
    if ishandle(varargin{1})
        ax = varargin{1};
        h = figure(get(ax,'parent'));
        this = varargin{2};
    else
        this = varargin{1};
        h = gcf;
        cdata = subsref(varargin{2},struct('type','.','subs','cdata'));
    end
    if nargin > 2
        indc = varargin{3};
    else
        indc = 1;
    end
end

if isempty(ax), ax = axes('Parent',h); end
axis(ax,'equal');
axis(ax,'off');
hp = patch(struct(...
    'vertices',  subsref(this,struct('type','.','subs','vertices')),...
    'faces',     subsref(this,struct('type','.','subs','faces'))),...
    'FaceColor', 'b',...
    'EdgeColor', 'none',...
    'Parent',ax);

if ~isempty(cdata)
    set(hp,'FaceVertexCData',cdata(:,indc), 'FaceColor','interp')
end

axes(ax);
camlight;
camlight(-80,-10);
lighting(ax,'gouraud');
if strcmp(spm_check_version,'matlab')
    cameratoolbar(h);
end

if nargout
    varargout{1} = hp;
end


text:
function img = rendervol(imgpath, volpath, texture)
% Renders a volume which has been generated by the VRN lua script.

vol = readvol(volpath, [192 192 200]);

im = imread(imgpath);

Q = isosurface(vol);
Q.vertices(:,3) = Q.vertices(:,3) * 0.5;
Q.EdgeColor = 'none';
Q.FaceColor = [120 154 194]/255;

if texture
    red = double(im(:,:,1));
    green = double(im(:,:,2));
    blue = double(im(:,:,3));
    [vcx,vcy] = meshgrid(1:192,1:192);
    vc = [vcx(:) vcy(:) red(:) green(:) blue(:)];
    cols = knnsearch(vc(:,1:2), Q.vertices(:,1:2));
    Q.FaceVertexCData = uint8(vc(cols,3:5));
    Q.FaceColor = 'interp';
end

figure
imshow(im, 'Border', 'tight');
hold on
patch(Q);
hold on
h = camlight;
lightangle(h, -10, 80)
lighting gouraud
material dull
view(0,90)
cap = getframe;
img = cap.cdata;

author: AaronJackson

repo: vrn

save_path: github-repos/MATLAB/AaronJackson-vrn


text:
function a = antihadamard ( n )

%*****************************************************************************80
%
%% ANTIHADAMARD returns the ANTIHADAMARD matrix.
%
%  Discussion:
%
%    An Anti-Hadamard matrix is one whose elements are all 0 or 1,
%    and for which the Frobenius norm of the inverse is as large as
%    possible.  This routine returns a matrix for which the Frobenius norm
%    of the inverse is large, though not necessarily maximal.
%
%  Formula:
%
%    if ( I = J )
%      A(I,J) = 1
%    else if ( I < J and mod ( I+J, 2 ) = 1 )
%      A(I,J) = 1
%    else
%      A(I,J) = 0
%
%  Example:
%
%    N = 5
%
%    1 1 0 1 0
%    0 1 1 0 1
%    0 0 1 1 0
%    0 0 0 1 1
%    0 0 0 0 1
%
%  Properties:
%
%    A is generally not symmetric: A' /= A.
%
%    A is Toeplitz: constant along diagonals.
%
%    A is upper triangular.
%
%    A is integral: int ( A ) = A.
%
%    A is a zero-one matrix.
%
%    The inverse of A is integral.
%
%    det ( A ) = 1.
%
%    A is unimodular.
%
%    LAMBDA(1:N) = 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ronald Graham, Neal Sloane,
%    Anti-Hadamard Matrices,
%    Linear Algebra and Applications,
%    Volume 62, 1984, pages 113-137.
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    for j = 1 : n

      if ( j < i )
        a(i,j) = 0.0;
      elseif ( i == j )
        a(i,j) = 1.0;
      elseif ( mod ( i+j, 2 ) == 1 )
        a(i,j) = 1.0;
      else
        a(i,j) = 0.0;
      end

    end
  end

  return
end


load uav1u
a1 = A;
load uav2u
a2 = A;
load uav3u
a3 = A;

n1 = 3890;
plot3(a1(1:n1,2), a1(1:n1,3),a1(1:n1,4)); hold on;
plot3(a2(1:n1,2), a2(1:n1,3),a2(1:n1,4)); hold on;
plot3(a3(1:n1,2), a3(1:n1,3),a3(1:n1,4)); hold on;
axis equal


text:
function [ t, w ] = t_quadrature_rule ( n )

%*****************************************************************************80
%
%% T_QUADRATURE_RULE: quadrature rule for T(n,x).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the rule.
%
%    Output, real T(N,1), W(N,1), the points and weights of the rule.
%
  aj = zeros ( n, 1 );

  bj = 0.5 * ones ( n, 1 );
  bj(1) = sqrt ( 0.5 );

  w = zeros ( n, 1 );
  w(1,1) = sqrt ( pi );

  [ t, w ] = imtqlx ( n, aj, bj, w );

  w(1:n,1) = w(1:n,1).^2;

  return
end


text:
function test_suite = test_linePosition
%TESTLINEPOSITION  One-line description here, please.
%
%   output = testLinePosition(input)
%
%   Example
%   testLinePosition
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-06-15,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testBasic(testCase) %#ok<*DEFNU>

point = [20 60];
line = createLine([10 30], [30 90]);

exp = .5;
pos = linePosition(point, line);

testCase.assertEqual(exp, pos);

function testPointArray(testCase)

point = [20 60;10 30;25 75];
line = createLine([10 30], [30 90]);

exp = [.5; 0; .75];
pos = linePosition(point, line);

testCase.assertEqual(exp, pos);

function testLineArray(testCase)

point = [20 60];
line1 = createLine([10 30], [30 90]);
line2 = createLine([0 0], [20 60]);
line3 = createLine([20 60], [40 120]);
lines = [line1;line2;line3];

exp = [.5  1  0];
pos = linePosition(point, lines);

testCase.assertEqual(exp, pos);

function testArrayArray(testCase)

lines = [10 20 1 0;10 20 0 1;10 20 1 1];
points = [30 20;10 40; 30 40];

exp = [20 0 10; 0 20 10; 20 20 20];
pos = linePosition(points, lines);

testCase.assertEqual(exp, pos, 'AbsTol', .01);

function testArrayArrayDiag(testCase)

lines = [10 20 1 0;10 20 0 1;10 20 1 1];
points = [30 20;10 40; 30 40];

exp = [20; 20; 20];
pos = linePosition(points, lines, 'diag');

testCase.assertEqual(exp, pos, 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
classdef FeatureExtractionConfiguratorUIElements < handle
    properties
        defaultFeaturesList;
        addFeatureButton;
        removeFeatureButton;
        selectedFeaturesList;
        featureStartRangeEditText;
        featureEndRangeEditText;
        featureFullSegmentCheckBox;
        featureAxisEditText;
        loadedFeatureExtractorsList;
        manualFeatureExtractionPanel;
        featuresSourceButtonGroup;
        featuresSourceFromFileRadio;
        featuresSourceManuallyRadio;
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function linpack_c_test17 ( )

%*****************************************************************************80
%
%% TEST17 tests CPBCO.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  m = 1;
  n = 3;
  lda = m + 1;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST17\n' );
  fprintf ( 1, '  For a single precision complex (C)\n' );
  fprintf ( 1, '  positive definite hermitian band matrix (PB),\n' );
  fprintf ( 1, '  CPBCO estimates the reciprocal condition number.\n' );
  fprintf ( 1, '  The matrix size is N = %d\n', n );
%
%  Set the value of the superdiagonal and diagonal.
%
  a(1,1) = complex ( 0.0000,  0.0000 );
  a(1,2) = complex ( 2.1341, -0.2147 );
  a(1,3) = complex ( 2.0905,  1.1505 );

  a(2,1) = complex ( 4.5281,  0.0000 );
  a(2,2) = complex ( 5.0371,  0.0000 );
  a(2,3) = complex ( 4.7638,  0.0000 );
%
%  Estimate the condition.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Estimate the condition.\n' );

  [ a, rcond, info ] = cpbco ( a, lda, n, m );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  CPBCO returned INFO = %d\n', info );
    fprintf ( 1, '  The factorization was not completed.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Reciprocal condition  = %f\n', rcond );

  return
end


text:
function calpak_test327 ( )

%*****************************************************************************80
%
%% CALPAK_TEST327 tests MONTH_CAL_STORE_COMMON.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 March 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST327\n' );
  fprintf ( 1, '  For the Common calendar:\n' );
  fprintf ( 1, '  MONTH_CAL_STORE_COMMON writes the day numbers for\n' );
  fprintf ( 1, '  a monthly calendar into a data structure.\n' );

  y = 1984;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Calendar:\n' );
  fprintf ( 1, '  Year = %d\n', y );

  for m = 1 : 12

    lines = month_cal_store_common ( y, m );
    [ r, ~ ] = size ( lines );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Month = %d\n', m );
    fprintf ( 1, '\n' );
    for i = 1 : r
      fprintf ( 1, '  %d    %s\n', i, lines(i,:) );
    end

  end

  return
end


text:
function scanStruct = setScanXYZVals(scanStruct, xVals, yVals, zVals);
%"setScanXYZVals"
%   Sets the x, y, and z values of the cols, rows, and slices of the
%   passed scanStruct.
%
%   xVals must be EVENLY spaced and INCREASING.
%   yVals must be EVENLY spaced and DECREASING.
%   zVals must be simply INCREASING.
%
%   By JRA 1/26/05
%
%   scanStruct        : a planC{indexS.scan} struct.
%   xVals yVals zVals : new x,y,z values for scan.
%
% Usage:
%   function scanStruct = setScanXYZVals(scanStruct, xVals, yVals, zVals);
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

scanSize = size(getScanArray(scanStruct));
if length(xVals) ~= scanSize(2)
    error('Invalid number of xVals, must have the same number as columns in scanArray.');
end
if length(yVals) ~= scanSize(1)
    error('Invalid number of yVals, must have the same number as rows in scanArray.');
end
if length(zVals) ~= scanSize(3)
    error('Invalid number of zVals, must have the same number as slices in scanArray.');
end

sizeDim1 = scanSize(1);
sizeDim2 = scanSize(2);
sizeDim3 = scanSize(3);

grid2Units = xVals(2) - xVals(1);
grid1Units = -(yVals(2) - yVals(1));

for i=1:scanSize(3)
    scanStruct.scanInfo(i).grid1Units = grid1Units;
    scanStruct.scanInfo(i).grid2Units = grid2Units;
    scanStruct.scanInfo(i).sizeOfDimension1 = sizeDim1;
    scanStruct.scanInfo(i).sizeOfDimension2 = sizeDim2;
    scanStruct.scanInfo(i).zValue = zVals(i);
    scanStruct.scanInfo(i).xOffset = mean(xVals);
    scanStruct.scanInfo(i).yOffset = mean(yVals);
end

scanStruct.uniformScanInfo.grid1Units = grid1Units;
scanStruct.uniformScanInfo.grid2Units = grid2Units;
scanStruct.uniformScanInfo.sizeOfDimension1 = sizeDim1;
scanStruct.uniformScanInfo.sizeOfDimension2 = sizeDim2;
scanStruct.uniformScanInfo.zValue = zVals(i);
scanStruct.uniformScanInfo.xOffset = mean(xVals);
scanStruct.uniformScanInfo.yOffset = mean(yVals);
scanStruct.uniformScanInfo.firstZValue = scanStruct.scanInfo(1).zValue;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [data] = tapas_sem_multiv_prepare_data(data, model, inference)
%% Return a structure ready to be used.
%
% Input
%
% Output
%

% aponteeduardo@gmail.com
% copyright (C) 2016
%

data = data;

end

author: translationalneuromodeling

repo: tapas


text:
% This script demonstrate that even if several sparse linear arrays
% share the same virtual ULA, they exhibit different performances under
% different SNR settings and different numbers of sources.
% This script will produce results similar to Fig. 1 in the following
% paper:
% * M. Wang, Z. Zhang, and A. Nehorai, "Performance analysis of
%   coarray-based MUSIC and the Cramér-Rao bound," in 2017 IEEE
%   International Conference on Acoustics, Speech and Signal Processing
%   (ICASSP), 2017, pp. 3061-3065.

clear(); close all;

wavelength = 1; % normalized wavelength
d_0 = wavelength / 2;
designs = {...
    design_array_1d('nested', [5  6], d_0, 'Nested (5, 6)') ...
    design_array_1d('nested', [2 12], d_0, 'Nested (2, 12)') ...
    design_array_1d('nested', [3  9], d_0, 'Nested (3, 9)') ...
    design_array_1d('nested', [1 18], d_0, 'Nested (1, 18)') ...
};
n_designs = length(designs);

power_source = 1;
n_snaphots = 1000;

n_grid = 20;
SNRs = linspace(-10, 20, n_grid);

doas1 = linspace(-pi/3, pi/3, 8);

MSEs_SNR_ana1 = zeros(n_designs, n_grid);
for dd = 1:n_designs
    design = designs{dd};
    A = steering_matrix(design, wavelength, doas1);
    for ii = 1:n_grid
        power_noise = power_source*10^(-SNRs(ii)/10);
        MSEs_SNR_ana1(dd, ii) = mean(ecov_coarray_music_1d(design, wavelength, ...
                doas1, power_source, power_noise, n_snaphots, 'DiagonalsOnly'));
    end
end

doas2 = linspace(-pi/3, pi/3, 20);

MSEs_SNR_ana2 = zeros(n_designs, n_grid);
for dd = 1:n_designs
    design = designs{dd};
    A = steering_matrix(design, wavelength, doas2);
    for ii = 1:n_grid
        power_noise = power_source*10^(-SNRs(ii)/10);
        MSEs_SNR_ana2(dd, ii) = mean(ecov_coarray_music_1d(design, wavelength, ...
                doas2, power_source, power_noise, n_snaphots, 'DiagonalsOnly'));
    end
end

figure;
subplot(1,2,1);
semilogy(SNRs, rad2deg(sqrt(MSEs_SNR_ana1)));
xlabel('SNR (dB)'); ylabel('RMSE (deg)'); grid on;
legend(arrayfun(@(x) x{1}.name, designs, 'UniformOutput', false));
title(sprintf('K = %d', length(doas1)));
subplot(1,2,2);
semilogy(SNRs, rad2deg(sqrt(MSEs_SNR_ana2)));
xlabel('SNR (dB)'); ylabel('RMSE (deg)'); grid on;
legend(arrayfun(@(x) x{1}.name, designs, 'UniformOutput', false));
title(sprintf('K = %d', length(doas2)));

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
function ok = validate(qual)
% Checks that an object of type Seg_Quality obeys certain rules
% that must always be true.
% Inputs:
%   qual: the object to be checked.
% Outputs:
%   ok: a boolean value indicating whether the object is valid.

assert(nargin==1)
assert(isa(qual,'Seg_Quality'))

ok = iscell(qual.ids);
ok = ok && (size(qual.values,2)==length(qual.ids));

author: nesl


text:
 function [fw, hr, hl] = fwhm(psfs, ii, doplot)
%function [fw, hr, hl] = fwhm(psfs, ii, doplot)
% compute fwhm of point-spread function centered at pixel ii
% and half-right width and half-left width (fw = hr + hl)

if ~nargin, ir_usage, end

warning 'fwhm is obsolete.  use fwhm1 for 1D fwhm'

[np nc] = size(psfs);
if (np == 1)
	psfs = psfs';
	[np, nc] = size(psfs);
end

warned = false;
for ic = 1:nc
	psf = psfs(:,ic);
	if (nargin < 2)
		ii = imax(psf);
	end

	% normalize
	psf = psf / psf(ii);
	if ~warned && (1 ~= max(psf))
		warning('peak not at center')
		warned = true;
	end

	% right
	ir = sum(cumprod(double6(psf((ii+1):np) >= 0.5)));
	high	= psf(ii + ir);
	low	= psf(ii + ir + 1);
	hr(ic,1) = ir + (high - 1/2) / (high-low);

	% left
	il = sum(cumprod(double6(psf((ii-1):-1:1) >= 0.5)));
	high	= psf(ii - il);
	low	= psf(ii - il - 1);
	hl(ic,1) = il + (high - 1/2) / (high-low);
end

fw = hr + hl;

if nargin > 2
	plot(1:np, psf, 'o', ii+[-left right], [0.5 0.5], '-')
end

author: JeffFessler

repo: mirt


text:
function [Model,probability] = LocalPCA(PopDec,M,K)
% Partitioning the population by Local PCA algorithm

%------------------------------- 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".
%--------------------------------------------------------------------------

% This function is modified from the code in
% http://dces.essex.ac.uk/staff/zhang/IntrotoResearch/RegEDA.htm

    [N,D] = size(PopDec);
    Model = struct('mean',   num2cell(PopDec(1:K,:),2),...  % The mean of the model
                   'PI',     eye(D),...                     % The matrix PI
                   'eVector',[],...                         % The eigenvectors
                   'eValue', [],...                         % The eigenvalues
                   'a',      [],...                         % The lower bound of the projections
                   'b',      []);                           % The upper bound of the projections

    %% Modeling
    for iter = 1 : 50
        % Calculte the distance between each solution and its projection in
        % affine principal subspace of each cluster
        distance = zeros(N,K);
        for k = 1 : K
            distance(:,k) = sum((PopDec-repmat(Model(k).mean,N,1))*Model(k).PI.*(PopDec-repmat(Model(k).mean,N,1)),2);
        end
        % Partition
        [~,partition] = min(distance,[],2);
        % Update the model of each cluster
        updated = false(1,K);
        for k = 1 : K
            oldMean = Model(k).mean;
            current = partition == k;
            if sum(current) < 2
                if ~any(current)
                    current = randi(N);
                end
                Model(k).mean    = PopDec(current,:);
                Model(k).PI      = eye(D);
                Model(k).eVector = [];
                Model(k).eValue  = [];
            else
                Model(k).mean    = mean(PopDec(current,:),1);
                [eVector,eValue] = eig(cov(PopDec(current,:)-repmat(Model(k).mean,sum(current),1)));
                [eValue,rank]    = sort(diag(eValue),'descend');
                Model(k).eValue  = real(eValue);
                Model(k).eVector = real(eVector(:,rank));
                Model(k).PI      = Model(k).eVector(:,M:end)*Model(k).eVector(:,M:end)';
            end
            updated(k) = ~any(current) || sqrt(sum((oldMean-Model(k).mean).^2)) > 1e-5;
        end
        % Break if no change is made
        if ~any(updated)
            break;
        end
    end

	%% Calculate the smallest hyper-rectangle of each model
    for k = 1 : K
        if ~isempty(Model(k).eVector)
            hyperRectangle = (PopDec(partition==k,:)-repmat(Model(k).mean,sum(partition==k),1))*Model(k).eVector(:,1:M-1);
            Model(k).a     = min(hyperRectangle,[],1);
            Model(k).b     = max(hyperRectangle,[],1);
        else
            Model(k).a = zeros(1,M-1);
            Model(k).b = zeros(1,M-1);
        end
    end

    %% Calculate the probability of each cluster for reproduction
    % Calculate the volume of each cluster
    volume = prod(cat(1,Model.b)-cat(1,Model.a),2);
    % Calculate the cumulative probability of each cluster
    probability = cumsum(volume/sum(volume));
end


text:
function [dmeg, dspc] = block_delay(bmeg, bspc, d)

Nblock = length(bmeg);
dmeg = cell(Nblock,1);
dspc = cell(Nblock,1);
for bi=1:Nblock
    common = min(size(bmeg{bi},1),length(bspc{bi}));
    dmeg{bi} = bmeg{bi}((1+d):common,:);
    dspc{bi} = bspc{bi}(1:(common-d));
end
dmeg = cell2mat(dmeg);
dspc = cell2mat(dspc);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function out = DN_Spread(y,spreadMeasure)
% DN_Spread     Measure of spread of the input time series.
%
% Returns the spread of the raw data vector, as the standard deviation,
% inter-quartile range, mean absolute deviation, or median absolute deviation.
%
%---INPUTS:
% y, the input data vector
%
% spreadMeasure, the spead measure:
%               (i) 'std': standard deviation
%               (ii) 'iqr': interquartile range
%               (iii) 'mad': mean absolute deviation
%               (iv) 'mead': median absolute deviation

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
% Check Inputs
% ------------------------------------------------------------------------------
if nargin < 2 || isempty(spreadMeasure)
    spreadMeasure = 'std'; % return std by default
end

% ------------------------------------------------------------------------------
% Evaluate the spread measure
% ------------------------------------------------------------------------------
switch spreadMeasure
	case 'std'
        % Standard deviation
		out = std(y);

	case 'iqr'
        % Interquartile range
		out = iqr(y);

	case 'mad'
        % Mean absolute deviation
		out = mad(y,0);

    case 'mead'
        % Median absolute deviation
        out = mad(y,1);

    otherwise
        error('Unknown spread measure ''%s''',spreadMeasure)
end

end

author: benfulcher

repo: hctsa


text:
function [flow_out, time_shifted_points_out] = em1_flow(...
    events, ...
    feature_pos, ...
    flow_init, ...
    params, ...
    fig)
%EM1_FLOW Estimates optical flow for an event feature.
%
% EM1_FLOW estimates the optical flow of a set of events % combined to form
% a 'feature'. This EM step runs from the events alone, as in:
% Alex Zihao Zhu, Nikolay Atanasov and Kostas Daniilidis.
% "Event-based Feature Tracking with Probabilistic Data Associations", ICRA 2017.
%
% Syntax:  [flow_out, time_shifted_points_out] = EM1_FLOW(...
%            events, ...
%            feature_pos, ...
%            flow_init, ...
%            params, ...
%            fig)
%
% Inputs:
%    events               - 4xN, each column is (x,y,t,p).
%    feature_pos          - 2x1, pixel position of the feature.
%    flow_init            - 2x1, initialization for the flow.
%    params               - parameters, defined in get_params().
%    fig                  - figure handle for plotting.
%
% Outputs:
%    flow_out                - 2x1, estimated flow.
%    time_shifted_points_out - 2xN, input events in the feature window
%                              shifted by the flow [x;y] + dt * flow.
%
% See also GET_PARAMS
%
% Author: Alex Zihao Zhu, University of Pennsylvania
% Email: alexzhu(at)seas.upenn.edu
% Copyright 2018 University of Pennsylvania
% Alex Zihao Zhu, Nikolay Atanasov, Kostas Daniilidis
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS, CONTRIBUTORS, AND THE
% TRUSTEES OF THE UNIVERSITY OF PENNSYLVANIA "AS IS" AND ANY EXPRESS OR
% IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
% OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
% IN NO EVENT SHALL THE COPYRIGHT OWNER, CONTRIBUTORS OR THE TRUSTEES OF
% THE UNIVERSITY OF PENNSYLVANIA BE LIABLE FOR ANY DIRECT, INDIRECT,
% INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
% NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
% THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
% (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
% THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

%% Initialization
flow = flow_init;
% Estimated flow
flow_out = [nan; nan];
% Final time shifted events.
time_shifted_points_out = [];
% Store change in flow norm at each iteration for plotting.
delta_flows = zeros(params.em1_params.max_iters, 1);
num_iter = 0;
% Plot handle.
scatter_plot_handle = [];
% Boolean array of size equal to the number of events. Is true if an event
% is in the spatial window of the feature, and thus needs to be processed.
event_window = [];
n_in_window = 0;

target_time = events(3, 1);

% Translate the events so that they are centered at the last position of
% the feature.
centered_events = events;
centered_events(1:2, :) = bsxfun(@minus, events(1:2, :), feature_pos);
prev_flow = flow;

%% Main EM Loop
while true
    if num_iter > params.em1_params.max_iters
        return
    end
    % Shift the centered events in time by the current estimate of the
    % flow.
    time_shifted_points = centered_events(1:2,:) + ...
        bsxfun(@times, flow,(target_time - centered_events(3,:)));

    % Only compute the event window once.
    if isempty(event_window)
        event_window = time_shifted_points(1, :) >= -params.window_size/2 & ...
            time_shifted_points(2, :) >= -params.window_size/2 & ...
            time_shifted_points(1, :) <= params.window_size/2 & ...
            time_shifted_points(2, :) <= params.window_size/2;
        n_in_window = sum(event_window);

        % Don't bother optimizing if there aren't enough events.
        if n_in_window < params.min_events_for_em
            return
        end

        centered_events = centered_events(:, event_window);
        time_shifted_points = time_shifted_points(:, event_window);
        %         params.max_distance = n_in_window / 400;
    end

    % Scale the events so that the computed distances are equal to
    % (x1-x2)/(2*sigma^2). Saves us having to divide all the
    % correspondences later.
    normalized_events = time_shifted_points' / (sqrt(2) * params.em1_params.sigma);

    % The KD tree allows us to compute distances between neighboring
    % events, while also performing outlier rejection. Unfortunately, as
    % the time shifted events change every iteration, it must also be
    % reconstructed at every iteration.
    kdtree = KDTreeSearcher(...
        normalized_events, ...
        'Distance', ...
        'euclidean');

    % Threshold distances by the Malhalanobis distance.
    [neighbors_cell, distances_cell] = rangesearch(...
        kdtree, ...
        normalized_events, ...
        params.em1_params.max_distance);

    distancesstacked = cell2mat(cellfun(...
        @transpose, distances_cell, 'UniformOutput', false))';

    % Can't solve for flow without any correspondences :(
    if isempty(distancesstacked)
        return
    end

    % Number of neighbors for each event.
    % NOTE: 'length' is not the same as @length
    num_neighbors_per_event = cellfun('length', neighbors_cell);

    neighbor_inds = cell2mat(cellfun(...
        @transpose, neighbors_cell, 'UniformOutput', false))';

    event_inds = repelem(1:n_in_window, num_neighbors_per_event);

    % The event to neighbor graph is undirected, so no need to double count
    % the event-neighbor correspondences.
    valid_correspondences = neighbor_inds > event_inds;

    neighbor_inds = neighbor_inds(valid_correspondences);
    event_inds = event_inds(valid_correspondences);
    distancesstacked = distancesstacked(valid_correspondences);

    % Distances are already scaled by the variance. Note that the original
    % equation in the paper is the sum of the product of two weights. Here
    % we simplify it with a single weight for speed.
    weights = exp(-distancesstacked);

    %% Simultaneously minimize over flow and translation.
    neighbor_events = centered_events(:, neighbor_inds);
    original_events = centered_events(:, event_inds);

    % These are the X and D matrices specified in the original paper,
    % except the weighting is handled by multiplying the D with the full
    % weight.
    X = original_events(1:2, :) - neighbor_events(1:2, :);
    D = original_events(3, :)-neighbor_events(3, :);
    weighted_D = bsxfun(@times, D, weights);
    DDT = weighted_D*D';
    XDT = X*weighted_D';

    flow = XDT/DDT;

    %% Calculate change in flow, plot debug information.
    if (norm(flow - prev_flow) < params.em1_params.min_err)
        break;
    end

    delta_flows(num_iter+1) = norm(flow - prev_flow);
    prev_flow = flow;

    if params.debug
        set(0, 'CurrentFigure', fig)
        subplot(2,1,1)
        plot(delta_flows(delta_flows > 0),'b')
        title('EM1 change in flow (convergence criterion)')
        xlim([0 params.em1_params.max_iters])

        subplot(2,1,2)
        if (~isempty(scatter_plot_handle))
            delete(scatter_plot_handle)
        end
        scatter_plot_handle = scatter(...
            time_shifted_points(1, :), ...
            time_shifted_points(2, :), 'r.');

        axis equal
        axis([-params.window_size/2-5 params.window_size/2+5 -params.window_size/2-5 params.window_size/2+5])
        axis ij
        title('EM1 Time shifted events')
        pause(0.01)
    end

    num_iter = num_iter + 1;
end

if params.debug
    pause(0.5)
end

flow_out = flow;

% Calculate the final shifted events to be used in EM2 later.
dt = events(3, end) - events(3, 1);
centered_events = events;
centered_events(1:2, :) = bsxfun(@minus, events(1:2, :), feature_pos + flow * dt);
time_shifted_points = centered_events(1:2, :) + ...
    bsxfun(@times, flow,(events(3, end)-centered_events(3, :)));

% Make the window a little bigger.
window_size = round(params.window_size * 1.5);

event_window = time_shifted_points(1, :) >= -window_size/2 & ...
    time_shifted_points(2, :) >= -window_size/2 & ...
    time_shifted_points(1, :) <= window_size/2 & ...
    time_shifted_points(2, :) <= window_size/2;

time_shifted_points_out = time_shifted_points(:, event_window);
end

author: daniilidis-group


temp4 = fact(10) / (fact(4) * fact(10-4)) * 0.5^4 * (1-0.5)^(10-4)
temp5 = fact(10) / (fact(5) * fact(10-5)) * 0.5^5 * (1-0.5)^(10-5)


text:
function bf = cubic_spline_fourier(f, a)

% The continuous Fourier transform of a cubic spline kernel.

bf = -(- 12*a + 12*exp(-pi*f*2i) + 12*exp(pi*f*2i) + 6*a*exp(-pi*f*4i) + ...
    6*a*exp(pi*f*4i) + f.*(pi*exp(-pi*f*2i)*12i) - f.*(pi*exp(pi*f*2i)*12i) + ...
    a*f.*(pi*exp(-pi*f*2i)*16i) - a*f.*(pi*exp(pi*f*2i)*16i) + ...
    a*f.*(pi*exp(-pi*f*4i)*4i) - a*f.*(pi*exp(pi*f*4i)*4i) - 24)./(16*f.^4*pi^4);

bf(f == 0) = 1;

author: martin-danelljan

repo: ECO