text:
function lik = lik_gaussiansmt(varargin)
%LIK_GAUSSIANSMT  Create a Gaussian scale mixture likelihood structure
%                 with priors producing approximation of the Student's t
%
%  Description
%    LIK = LIK_GAUSSIANSMT('ndata',N,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates a scale mixture noise covariance function structure
%    (with priors producing approximation of the Student's t) in
%    which the named parameters have the specified values. Any
%    unspecified parameters are set to default values. Obligatory
%    parameter is 'ndata', which tells the number of data points,
%    that is, number of mixture components.
%
%    LIK = LIK_GAUSSIANSMT(LIK,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a covariance function structure with the named
%    parameters altered with the specified values.
%
%    Parameters for the Gaussian scale mixture approximation of the
%    Student's t
%      sigma2    - Variances of the mixture components.
%                  The default is 1 x ndata vector of 0.1s.
%      U         - Part of the parameter expansion, see below.
%                  The default is 1 x ndata vector of 1s.
%      tau2      - Part of the parameter expansion, see below.
%                  The default is 0.1.
%      alpha     - Part of the parameter expansion, see below.
%                  The default is 0.5.
%      nu        - Degrees of freedom. The default is 4.
%      nu_prior  - Prior for nu. The default is prior_fixed().
%      gibbs     - Whether Gibbs sampling is 'on' (default) or 'off'.
%
%    Parametrisation and non-informative priors for alpha and tau
%    are same as in Gelman et. al. (2004) page 304-305:
%      y-E[y] ~ N(0, alpha^2 * U),
%      where U = diag(u_1, u_2, ..., u_n)
%          u_i ~ Inv-Chi^2(nu, tau^2)
%
%    The parameters of this likelihood can be inferred only by
%    Gibbs sampling by calling GP_MC.
%
%    If degrees of freedom nu is given a prior (other than
%    prior_fixed), it is sampled using slice sampling within Gibbs
%    sampling with limits [0,128].
%
%  See also
%    GP_SET, PRIOR_*, LIK_*

% Copyright (c) 1998,1999,2010 Aki Vehtari
% Copyright (c) 2007-2010 Jarno Vanhatalo

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'LIK_GAUSSIANSMT';
  ip.addOptional('lik', [], @isstruct);
  ip.addParamValue('ndata',[], @(x) isscalar(x) && x>0 && mod(x,1)==0);
  ip.addParamValue('sigma2',[], @(x) isvector(x) && all(x>0));
  ip.addParamValue('U',[], @isvector);
  ip.addParamValue('tau2',0.1, @isscalar);
  ip.addParamValue('alpha',0.5, @isscalar);
  ip.addParamValue('nu',4, @isscalar);
  ip.addParamValue('nu_prior',[], @(x) isstruct(x) || isempty(x));
  ip.addParamValue('censored',[], @(x) isstruct);
  ip.addParamValue('gibbs','on', @(x) ismember(x,{'on' 'off'}));
  ip.parse(varargin{:});
  lik=ip.Results.lik;

  if isempty(lik)
    init=true;
    lik.type = 'Gaussian-smt';
  else
    if ~isfield(lik,'type') || ~isequal(lik.type,'Gaussian-smt')
      error('First argument does not seem to be a valid likelihood function structure')
    end
    init=false;
  end

  % Initialize parameters
  if init || ~ismember('ndata',ip.UsingDefaults)
    ndata = ip.Results.ndata;
    lik.ndata=ndata;
    lik.r = zeros(ndata,1);
  end
  if isempty(ndata)
    error('NDATA has to be defined')
  end
  if init || ~ismember('sigma2',ip.UsingDefaults)
    sigma2=ip.Results.sigma2;
    if isempty(sigma2)
      lik.sigma2 = repmat(0.1,ndata,1);
    else
      if (size(sigma2,1) == lik.ndata && size(sigma2,2) == 1)
        lik.sigma2 = sigma2;
      else
        error('The size of sigma2 has to be NDATAx1')
      end
    end
    lik.sigma2 = sigma2;
  end
  if init || ~ismember('U',ip.UsingDefaults)
    U=ip.Results.U;
    if isempty(U)
      lik.U = ones(ndata,1);
    else
      if size(U,1) == lik.ndata
        lik.U = U;
      else
        error('the size of U has to be NDATAx1')
      end
    end
  end
  if init || ~ismember('tau2',ip.UsingDefaults)
    lik.tau2=ip.Results.tau2;
  end
  if init || ~ismember('alpha',ip.UsingDefaults)
    lik.alpha=ip.Results.alpha;
  end
  if init || ~ismember('nu',ip.UsingDefaults)
    lik.nu=ip.Results.nu;
  end
  if init || ~ismember('censored',ip.UsingDefaults)
    censored=ip.Results.censored;
    if ~isempty(censored)
      lik.censored = censored{1};
      yy = censored{2};
      if lik.censored(1) >= lik.censored(2)
        error('lik_gaussiansmt -> if censored model is used, the limits must be given in increasing order.')
      end

      imis1 = [];
      imis2 = [];
      if lik.censored(1) > -inf
        imis1 = find(yy<=lik.censored(1));
      end
      if lik.censored(1) < inf
        imis2 = find(yy>=lik.censored(2));
      end
      lik.cy = yy([imis1 ; imis2])';
      lik.imis = [imis1 ; imis2];
    end
  end
  % Initialize prior structure
  lik.p=[];
  lik.p.sigma=[];
  if init || ~ismember('nu_prior',ip.UsingDefaults)
    lik.p.nu=ip.Results.nu_prior;
  end
  % using Gibbs or not
  if init || ~ismember('gibbs',ip.UsingDefaults)
    lik.gibbs = ip.Results.gibbs;
  end
  if init
    % Set the function handles to the subfunctions
    lik.fh.pak = @lik_gaussiansmt_pak;
    lik.fh.unpak = @lik_gaussiansmt_unpak;
    lik.fh.lp = @lik_gaussiansmt_lp;
    lik.fh.lpg = @lik_gaussiansmt_lpg;
    lik.fh.cfg = @lik_gaussiansmt_cfg;
    lik.fh.trcov  = @lik_gaussiansmt_trcov;
    lik.fh.trvar  = @lik_gaussiansmt_trvar;
    lik.fh.gibbs = @lik_gaussiansmt_gibbs;
    lik.fh.recappend = @lik_gaussiansmt_recappend;
  end

end

function [w,s] = lik_gaussiansmt_pak(lik)
  w = []; s = {};
end

function [lik, w] = lik_gaussiansmt_unpak(lik, w)

end

function lp =lik_gaussiansmt_lp(lik)
  lp = 0;
end

function lpg  = lik_gaussiansmt_lpg(lik)
  lpg = [];
end

function DKff  = lik_gaussiansmt_cfg(lik, x, x2)
  DKff = [];
end

function C = lik_gaussiansmt_trcov(lik, x)
%LIK_GAUSSIANSMT_TRCOV  Evaluate training covariance matrix
%                    corresponding to Gaussian noise
%  Description
%    C = LIK_GAUSSIANSMT_TRCOV(GP, TX) takes in covariance function
%    of a Gaussian process GP and matrix TX that contains
%    training input vectors. Returns covariance matrix C. Every
%    element ij of C contains covariance between inputs i and j
%    in TX. This subfunction is needed only in Gaussian likelihoods.
%
%  See also
%    LIK_GAUSSIANSMT_COV, LIK_GAUSSIANSMT_TRVAR, GP_COV, GP_TRCOV

  [n, m] =size(x);
  n1=n+1;

  if n ~= lik.ndata
    error(['lik_gaussiansmt -> _trvar: The training variance can be evaluated'...
           '      only for training data.                                 '])
  end

  C = sparse(1:n, 1:n, lik.sigma2, n, n);
end

function C = lik_gaussiansmt_trvar(lik, x)
%LIK_GAUSSIANSMT_TRVAR  Evaluate training variance vector
%                    corresponding to Gaussian noise
%
%  Description
%    C = LIK_GAUSSIANSMT_TRVAR(LIK, TX) takes in covariance function
%    of a Gaussian process LIK and matrix TX that contains
%    training inputs. Returns variance vector C. Every
%    element i of C contains variance of input i in TX. This
%    subfunction is needed  only in Gaussian likelihoods.
%
%
%  See also
%    LIK_GAUSSIANSMT_COV, GP_COV, GP_TRCOV

  [n, m] =size(x);
  if n ~= lik.ndata
    error(['lik_gaussiansmt -> _trvar: The training variance can be evaluated'...
           '      only for training data.                                 '])
  end
  C = lik.sigma2;

end

function [lik, y] = lik_gaussiansmt_gibbs(gp, lik, x, y)
%LIK_GAUSSIANSMT_GIBBS  Function for sampling the sigma2's
%
%  Description
%    Perform Gibbs sampling for the scale mixture variances. This
%    function is likelihood specific.

  [n,m] = size(x);

  % Draw a sample of the mean of y. Its distribution is
  % f ~ N(K*inv(C)*y, K - K*inv(C)*K')
  switch gp.type
    case 'FULL'
      sampy = gp_rnd(gp, x, y, x);
    case 'FIC'
      sampy = gp_rnd(gp, x, y, x, 'tstind', 1:n);
    case {'PIC' 'PIC_BLOCK'}
      sampy = gp_rnd(gp, x, y, x, 'tstind', gp.tr_index);
  end
  % Calculate the residual
  r = y-sampy;

  U = lik.U;
  t2 = lik.tau2;
  alpha = lik.alpha;
  nu = lik.nu;
  rss2=alpha.^2.*U;

  % Perform the gibbs sampling (Gelman et.al. (2004) page 304-305)
  % Notice that 'sinvchi2rand' is parameterized as in Gelman et. al.
  U=sinvchi2rand(nu+1, (nu.*t2+(r./alpha).^2)./(nu+1));
  shape = n*nu./2;                               % These are parameters...
  invscale = nu.*sum(1./U)./2;                   % used in Gelman et al
  t2=gamrnd(shape, 1./invscale);                 % Notice! The matlab parameterization is different
  alpha2=sinvchi2rand(n,mean(r.^2./U));
  rss2=alpha2.*U;
  if ~isempty(lik.p.nu)
    % Sample nu using Gibbs sampling
    pp = lik.p.nu;
    opt=struct('nomit',4,'display',0,'method','doubling', ...
               'wsize',4,'plimit',5,'unimodal',1,'mmlimits',[0; 128]);
    nu=sls(@(nu) (-sum(sinvchi2_lpdf(U,nu,t2))-pp.fh.lp(nu, pp)),nu,opt);
  end
  lik.sigma2 = rss2;
  lik.U = U;
  lik.tau2 = t2;
  lik.alpha = sqrt(alpha2);
  lik.nu = nu;
  lik.r = r;
  if isfield(lik, 'censored')
    imis1 = [];
    imis2 = [];
    if lik.censored(1) > -inf
      imis1 = find(y<=lik.censored(1));
      y(imis1)=normrtrand(sampy(imis1),alpha2*U(imis1),lik.censored(1));
    end

    if lik.censored(1) < inf
      imis2 = find(y>=lik.censored(2));
      y(imis2)=normltrand(sampy(imis2),alpha2*U(imis2),lik.censored(2));
    end
    lik.cy = y([imis1 ; imis2]);
  end
end

function reccf = lik_gaussiansmt_recappend(reccf, ri, lik)
%RECAPPEND  Record append
%
%  Description
%    RECCF = LIK_GAUSSIANSMT_RECAPPEND(RECCF, RI, LIK)
%    takes a likelihood record structure RECCF, record
%    index RI and likelihood structure LIK with the
%    current MCMC samples of the parameters. Returns
%    RECCF which contains all the old samples and the
%    current samples from LIK . This subfunction is
%    needed when using MCMC sampling (gp_mc).
%
%  See also
%    GP_MC and GP_MC -> RECAPPEND

  if nargin == 2
    % Initialize the record
    reccf.type = 'Gaussian-smt';
    lik.ndata = [];

    % Initialize parameters
    reccf.sigma2 = [];

    % Set the function handles
    reccf.fh.pak = @lik_gaussiansmt_pak;
    reccf.fh.unpak = @lik_gaussiansmt_unpak;
    reccf.fh.lp = @lik_gaussiansmt_lp;
    reccf.fh.lpg = @lik_gaussiansmt_lpg;
    reccf.fh.cfg = @lik_gaussiansmt_cfg;
    reccf.fh.cov = @lik_gaussiansmt_cov;
    reccf.fh.trcov  = @lik_gaussiansmt_trcov;
    reccf.fh.trvar  = @lik_gaussiansmt_trvar;
    reccf.fh.gibbs = @lik_gaussiansmt_gibbs;
    reccf.fh.recappend = @lik_gaussiansmt_recappend;
  else
    % Append to the record
    reccf.ndata = lik.ndata;
    gpp = lik.p;

    % record noiseSigma
    reccf.sigma2(ri,:)=lik.sigma2;
    if ~isempty(lik.nu)
      reccf.nu(ri,:)=lik.nu;
      reccf.U(ri,:) = lik.U;
      reccf.tau2(ri,:) = lik.tau2;
      reccf.alpha(ri,:) = lik.alpha;
      reccf.r(ri,:) = lik.r;
    end
    if isfield(lik, 'censored')
      reccf.cy(ri,:) = lik.cy';
    end
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%%%  HeightControl_HIL
clear
path( './icon/',path);
load vehicle_local_position.mat;
% constant value
RAD2DEG = 57.2957795;
DEG2RAD = 0.0174533;
% throttle when UAV is hovering
THR_HOVER = 0.609;

%% control parameter
% attitude PID parameters
Kp_PITCH_ANGLE = 6.5;
Kp_PITCH_AngleRate = 0.1;
Ki_PITCH_AngleRate = 0.02;
Kd_PITCH_AngleRate = 0.001;
Kp_ROLL_ANGLE = 6.5;
Kp_ROLL_AngleRate = 0.1;
Ki_ROLL_AngleRate = 0.02;
Kd_ROLL_AngleRate = 0.001;

Kp_YAW_AngleRate = 0.5;
Ki_YAW_AngleRate = 0.01;
Kd_YAW_AngleRate = 0.00;
% position PID parameters
Kpxp = 1.0;
Kpyp = 1.0;
Kpzp = 4.0;
Kvxp = 2.5; Kvxi = 0.4; Kvxd = 0.01;
Kvyp = 2.5; Kvyi = 0.4; Kvyd = 0.01;
Kvzp = 0.45; Kvzi = 0.01; Kvzd = 0.005;
% integral saturation
Saturation_I_RP_Max = 0.3;
Saturation_I_RP_Min = -0.3;
Saturation_I_Y_Max = 0.2;
Saturation_I_Y_Min = -0.2;
Saturation_I_ah = 3.43;
Saturation_I_az = 5;

% max control angle,default 35deg
MAX_CONTROL_ANGLE_ROLL = 35;
MAX_CONTROL_ANGLE_PITCH  = 35;
% max control angle rate,rad/s
MAX_CONTROL_ANGLE_RATE_PITCH = 220;
MAX_CONTROL_ANGLE_RATE_ROLL = 220;
MAX_CONTROL_ANGLE_RATE_Y = 200;
% max control speed, m/s
MAX_CONTROL_VELOCITY_XY = 5;
MAX_CONTROL_VELOCITY_Z = 3;
% throttle amplitude
MAX_MAN_THR = 0.9;
MIN_MAN_THR = 0.05;
%% run simulink model
HeightControl_HIL

author: RflySim

repo: RflyExpCode

save_path: github-repos/MATLAB/RflySim-RflyExpCode


text:
function [ Detect ] = Truelabel2Trueframe( TrueLabel_bin,wsize,wstep )

iidx=1;
Frame_iidx=1;
Frame_len=Frame_Length( TrueLabel_bin,wstep,wsize );
Detect=zeros(Frame_len,1);

while(1)

    %% Windowing
    if(iidx+wsize<length(TrueLabel_bin))
        TrueLabel_frame=TrueLabel_bin(iidx:iidx+wsize-1)*10; % Because the size of TrueLabel_bin is 0.1 in TrueLabel function
    else
        TrueLabel_frame=TrueLabel_bin(iidx:end)*10; % Because the size of TrueLabel_bin is 0.1 in TrueLabel function
    end

    if(sum(TrueLabel_frame)>=wsize/2)
        TrueLabel_frame=1;
    else
        TrueLabel_frame=0;
    end
    if(Frame_iidx>length(Detect))
        break;
    end
    Detect(Frame_iidx)=TrueLabel_frame;
    iidx=iidx+wstep;
    Frame_iidx=Frame_iidx+1;
    if(iidx>length(TrueLabel_bin))
        break;
    end

end

end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function aux_movingwindow(params, hParentFigure)
% function aux_FMD(params, hParentFigure);
%-------------------------------------------
%
% Incoming variables:
% params        : all variables
% hParentFigure : Handle of the parent figure
%
% Thomas van Stiphout, thomas@sed.ethz.ch
% last update: 7.9.2005

report_this_filefun(mfilename('fullpath'));

fStartTime_ = params(1).fStartTime;
fTimeCut_   = params(1).fTimeCut;
vTimeCut_=[params(1).fStartTime+params(1).fTwLength:1:params(1).fTimeCut]
for i = 1:1:length(vTimeCut_)
    params(i).fTimeCut = vTimeCut_(i);
    params(i)=sr_calcZ(params(i))
    % Add parameter to params.sComment
     params(i).sComment = [params(i).rContainer ' Time Cut ' num2str(params(i).fTimeCut) ...
                   ', Time period ' num2str(params(i).fTwLength) ' yrs'];
    if i<length(vTimeCut_)
        params=[params;params(1)]
    end
end
gui_result(params(:))
clear fStartTime_ fTimeCut_ vTimeCut_ i

author: CelsoReyes


text:
function d = edge_degree ( n_node, n_edge, t )

%*****************************************************************************80
%
%% EDGE_DEGREE returns the degree of the nodes of a graph stored by edges.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer N_NODE, the number of nodes in the graph.
%    N_NODE must be positive.
%
%    Input, integer N_EDGE, the number of edges in the graph.
%    N_EDGE must be positive.
%
%    Input, integer T(2,N_EDGE), describes the edges of the tree
%    as pairs of nodes.
%
%    Output, integer D(N_NODE), the degree of each node.
%

%
%  Check.
%
  ierror = edge_check ( n_node, n_edge, t );

  if ( ierror ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'EDGE_DEGREE - Fatal error!\n' );
    fprintf ( 1, '  The input array is illegal.\n' );
    fprintf ( 1, '  IERROR = %d\n', ierror );
    error ( 'EDGE_DEGREE - Fatal error!' );
  end
%
%  Compute the degree of each node.
%
  d(1:n_node) = 0;

  for j = 1 : n_edge
    d(t(1,j)) = d(t(1,j)) + 1;
    d(t(2,j)) = d(t(2,j)) + 1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [feat boxes] = lbp_feature(im)

if(size(im,3) < 3)
    %black and white image
    im = cat(3, im, im, im); %make it a trivial color image
end
boxes = [1;1;size(im,1);size(im,2)];

[lbp_L0 lbp_L1 lbp_L2] = lbp_original_4x4(im);

feat.hists.L0 = lbp_L0;
feat.hists.L1 = lbp_L1;
feat.hists.L2 = lbp_L2;

author: adikhosla

repo: feature-extraction


text:
function msm_to_mm_test18 ( )

%*****************************************************************************80
%
%% MSM_TO_MM_TEST18 tests MSM_TO_MM_COORDINATE_REAL_GENERAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MSM_TO_MM_TEST18\n' );
  fprintf ( 1, '  Convert an MSM to MM coordinate real general format.\n' );

  output_filename = 'msm_to_mm_test18.mm';

  a = r8mat_indicator ( 5, 3 );
%
%  Have MSM_TO_MM write the matrix to a file.
%
  msm_to_mm ( output_filename, a, 'coordinate', 'real', 'general' );

  return
end


text:
function JTotal=polAngGradient(xG,systemType,lRx)
%%POLANGGRADIENT Determine the gradient of the angle of a 2D polar
%           measurement with respect to position (gradient components for
%           velocity etc. are zero and are not provided).Atmospheric and
%           other propagation effects are not taken into account. The
%           measurement can be monostatic or bistatic, but the location of
%           the transmitter does not matter.
%
%INPUTS: xG A 2XN set of Cartesian locations in the format [x;y] where the
%          gradient is desired.
% systemType An optional parameter specifying the axis from which the
%          azimuth angle is measured. It is assumed that the azimuth
%          angle is given in radians. Possible values are
%          0 (The default if omitted) The azimuth angle is
%             counterclockwise from the x axis.
%          1 The azimuth angle is measured clockwise from the y axis.
%      lRx The 2X1 [x;y] location vector of the receiver in Cartesian
%          coordinates. If this parameter is omitted or an empty matrix is
%          passed, then the receiver is assumed to be at the origin.
%
%OUTPUTS: JTotal A 1X2XN set of gradient matrices of the polar angle taken
%           with respect to the components [x,y] in 2D in that order for
%           each of the N points in xG.
%
%The conversion utilizing bistatic polar measurements in 2D is similar to
%that using bistatic r-u-v measurements in 3D, which is discussed in [1].
%Basic differentiation was analytically performed to obtain the expressions
%used in this file.
%
%Note that normally we would convert the point to the local coordinate
%system as
%xLocal=M*(xG(1:2,curPoint)-lRx(1:2))
%then compute the gradient J, and then undo the effects of the rotation
%matrix M in the end to the global coordinates as JTotal=J*M.
%However, that is not necessary with polar angles, due to cancellations.
%Say that M=[m11,m12;m21,m22] is a rotation matrix. Then, the partial
%derivative with respect to x is explicitly:
%((m12*m21-m11*m22)*y)/((m11^2+m21^2)*x^2+2*(m11*m12+m21*m22)*x*y+(m12^2+m22^2)*y^2)
%However, one will see that because M is a rotation matrix and M*M'=eye(2)
%then m12*m21-m11*m22=-1, m11^2+m21^2=m12^2+m22^2=1, and m11*m12+m21*m22=0,
%which is the same as if no rotation matrix were present. Basically, all M
%does is add a constant to the direction angle; it does not change any of
%the derivatives. Thus, M does not matter for the Jacobian.
%
%EXAMPLE:
%Here, we verify that a numerically differentiated Jacobian is consistent
%with the analytic one produced by this function.
% x=[100;-1000];
% lRx=[500;20];
% epsVal=1e-5;
% systemType=0;
% M=randRotMat(2);
% J=polAngGradient(x,systemType,lRx);
%
% theta=getPolAngle(x,systemType,lRx,M);
% thetadX=getPolAngle(x+[epsVal;0],systemType,lRx,M);
% thetadY=getPolAngle(x+[0;epsVal],systemType,lRx,M);
%
% JNumDiff=[(thetadX-theta)/epsVal;(thetadY-theta)/epsVal];
% max(abs((J(:)-JNumDiff(:))./J(:)))
%The relative error will be on the order of 1e-8, indicating good agreement
%between the numerical Jacobian vector and the actual Jacobian vector.
%
%REFERENCES:
%[1] David F. Crouse , "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems
%    Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%
%February 2017 David F.Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

N=size(xG,2);

%M does not matter, so we do not need to check whether it is given.

if(nargin<3||isempty(lRx))
    lRx=zeros(2,1);
end

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

JTotal=zeros(1,2,N);
for curPoint=1:N
    %Convert the point into the local coordinate system of the receiver.
    xLocal=xG(1:2,curPoint)-lRx(1:2);

    xL=xLocal(1);
    yL=xLocal(2);

    J=zeros(1,2);
    switch(systemType)
        case 0
            %Derivative with respect to x.
            J(1,1)=-yL/(xL^2+yL^2);

            %Derivative with respect to y.
            J(1,2)=xL/(xL^2+yL^2);
        case 1
            %Derivative with respect to x.
            J(1,1)=yL/(xL^2+yL^2);

            %Derivative with respect to y.
            J(1,2)=-xL/(xL^2+yL^2);
        otherwise
            error('Invalid system type specified.')
    end

    %Undo the rotation to move the gradient into the global coordinate
    %system.
    JTotal(:,:,curPoint)=J;
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 [Yi]=interp1_ND(X,Y,Xi,interpDim,interpMethod)

% function [Yi]=interp1_ND(X,Y,Xi,interpDim,interpMethod)
% ------------------------------------------------------------------------
% This function is similar to interp1. However it can perform 1D
% interpolation for multidimensional arrays. E.g. Time interpolation for 3D
% image data varying in time. The direction of interpolation is specified
% by interpDim, and the method by interpMethod
% ('nearest','linear','cubic','pchip').
%
%
% Reference for pchip method:
%   F. N. Fritsch and R. E. Carlson, "Monotone Piecewise Cubic
%   Interpolation", SIAM J. Numerical Analysis 17, 1980, 238-246.
%   David Kahaner, Cleve Moler and Stephen Nash, Numerical Methods
%   and Software, Prentice Hall, 1988.
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2015/03/04
%------------------------------------------------------------------------

%% Check input

%Check number of dimensions
nDimsX = ndims(X);
nDimsXi = ndims(Xi);

if nDimsX~=nDimsXi
    if size(Xi,nDimsX)~=1
        error('Input arrays have different number of dimensions');
    end
end

%Check column/row form
if iscolumn(X)
    X=X';
end

if iscolumn(Xi)
    Xi=Xi';
end

if isvector(X)
    interpDim=2;
end

%% Cast in 2D array format

siz=size(X);
siz_i=size(Xi);

dimInd=1:nDimsX;
logicNotInterpDim=dimInd~=interpDim;
permOrder=[sort(dimInd(logicNotInterpDim)) interpDim];

%Permuting array order so that interpolation dimension is last
X = permute(X,permOrder);
Y = permute(Y,permOrder);
Xi = permute(Xi,permOrder);

%Convert to 2D array with interpolation performed allong column direction
siz1=prod(siz(logicNotInterpDim));
X=reshape(X(:),[siz1, numel(X)/siz1]);
Y=reshape(Y(:),[siz1, numel(Y)/siz1]);
Xi=reshape(Xi(:),[siz1, numel(Xi)/siz1]);

%%

siz_2D=size(X);
siz_i_2D=size(Xi);

switch interpMethod
    case {'cubic','pchip'}
        diff_Y=diff(Y,1,2);
        diff_X=diff(X,1,2);
        K=diff_Y./diff_X;
end

switch interpMethod
    case 'cubic'
        M=zeros(siz_2D);
        M(:,1)=K(:,1);
        M(:,end)=K(:,end);
        M(:,2:end-1)=0.5*(K(:,1:end-1)+K(:,2:end));
    case 'pchip'
        M = pchipSlopes(X,Y,K);
end

%%

I=(1:siz_2D(1))';

%Initialize
Yi=nan(siz_i_2D);
p0=zeros(size(X,1),1); p1=p0;
switch interpMethod
    case 'cubic'
        m0=p0; m1=p0;
    case 'pchip'
        m0=p0; m1=p0;
end

for q=1:1:siz_i_2D(end)

    %Find close point pair to define interval
    D=abs(X-Xi(:,q*ones(1,size(X,2))));

    [~,J1]=min(D,[],2);
    indMin1=sub2ind(siz_2D,I,J1);
    D(indMin1)=inf;
    [~,J2]=min(D,[],2);
    indSort=[J1 J2];
    indSort=sort(indSort(:,1:2),2); %Sorted indices of closest two

    ind1=sub2ind(siz_2D,I,indSort(:,1));
    ind2=sub2ind(siz_2D,I,indSort(:,2));

    w=X(ind2)-X(ind1); %interval width
    t=(Xi(:,q)-X(ind1))./w; %t coordinate in domain
    p0(1:end)=Y(ind1); %First point
    p1(1:end)=Y(ind2); %End point
    switch interpMethod
        case {'cubic','pchip'}
            m0(1:end)=M(ind1); %First derivative
            m1(1:end)=M(ind2); %Final derivative

            h00=2.*t.^3-3.*t.^2+1;
            h10=t.^3-2*t.^2+t;
            h01=-2*t.^3+3*t.^2;
            h11=t.^3-t.^2;
            pt=h00.*p0 + w.*h10.*m0 + h01.*p1 + w.*h11.*m1;
        case 'linear'
            pt=t.*p1+(1-t).*p0;
        case 'nearest'
            logicPrev=t<=0.5;
            pt(logicPrev)=p0(logicPrev);
            pt(~logicPrev)=p1(~logicPrev);
    end

    Yi(:,q)=pt; %Store in array

end

%% Cast back in original form

% Xi=reshape(Xi(:),siz_i(permOrder));
% Xi=ipermute(Xi,permOrder);

Yi=reshape(Yi(:),siz_i(permOrder));
Yi=ipermute(Yi,permOrder);

end
%%
function M = pchipSlopes(X,Y,dYdX)
%PCHIPSLOPES  Derivative values for shape-preserving Piecewise Cubic Hermite
% Interpolation.
% d = pchipslopes(x,y,del) computes the first derivatives, d(k) = P'(x(k)).

n = size(X,2);

%  Slopes at interior points.
%  d(k) = weighted average of del(k-1) and del(k) when they have the same sign.
%  d(k) = 0 when del(k-1) and del(k) have opposites signs or either is zero.

M = zeros(size(Y));

[I,J]= find(sign(dYdX(:,1:n-2)).*sign(dYdX(:,2:n-1)) > 0);
ind=sub2ind(size(X),I,J);
ind1=sub2ind(size(X),I,J+1);

h = diff(X,1,2);
hs = h(ind)+h(ind1);
w1 = (h(ind)+hs)./(3*hs);
w2 = (hs+h(ind1))./(3*hs);
dmax = max(abs(dYdX(ind)), abs(dYdX(ind1)));
dmin = min(abs(dYdX(ind)), abs(dYdX(ind1)));
M(ind1) = dmin./conj(w1.*(dYdX(ind)./dmax) + w2.*(dYdX(ind1)./dmax));

%  Slopes at end points.
%  Set d(1) and d(n) via non-centered, shape-preserving three-point formulae.

M(:,1) = ((2*h(:,1)+h(:,2)).*dYdX(:,1) - h(:,1).*dYdX(:,2))./(h(:,1)+h(:,2));
L = (sign(M(:,1)) ~= sign(dYdX(:,1)));
M(L,1) = 0;

L = ((sign(dYdX(:,1)) ~= sign(dYdX(:,2))) & (abs(M(:,1)) > abs(3*dYdX(:,1))));
M(L,1) = 3*dYdX(L,1);

M(:,n) = ((2*h(:,n-1)+h(:,n-2)).*dYdX(:,n-1) - h(:,n-1).*dYdX(:,n-2))./(h(:,n-1)+h(:,n-2));

L= sign(M(:,n)) ~= sign(dYdX(:,n-1));
M(L,n) = 0;

L= (sign(dYdX(:,n-1)) ~= sign(dYdX(:,n-2))) & (abs(M(:,n)) > abs(3*dYdX(:,n-1)));
M(L,n) = 3*dYdX(L,n-1);

end

%%
% _*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 sr_parameterstudy

disp('~/zmap/src/thomas/seismicrates/sr_parameterstudy.m')

% reset random number generator
rand('state',sum(100*clock));
% set default sample size
nN=[50];

% create random number generator

for iCat=1:200
    vRand=rand(nN(1),1);
    iCat
    % create different time intervals from 2 to 50 year
    vYY=[5:5:40]';
    for iYY=1:size(vYY,1)
        fYr1=1980;
        fYr2=1980+vYY(iYY);
        fTw=(fYr2-fYr1)/2;

        % create range of time bins
        %     nTbin1=[1/10000 1/2000 1/1000 1/500 1/100 1/50 1/25 1/10];
        %     nTbin1=1./([10 20 30 40 50 60 70 80 90 100 125 150 175 200 250 300 400 500 600 700 800 900 1000]);
        nTbin1=1./([20 40 50 75 100 125 150 175 200 300 400 500]);

        % round time bins
        nTbin=ceil((fYr2-fYr1)*nTbin1*365)';

        % set default catalog size
        nCatSize=1000;
        for j=1:size(nN)     % for each nN
            %             mCat00=rand(nN(j),1)*(fYr2-fYr1)+fYr1;

            mCat00=vRand*(fYr2-fYr1)+fYr1;
            mCat=rand(nCatSize,1)*(fYr2-fYr1)+fYr1;
            mCat20=mCat00(mCat00>(fYr2-fTw));
            for i=1:size(nTbin)
                [mLTA(i,iYY,iCat), mLTAprob(i,iYY,iCat)] =calc_zlta(mCat,mCat00,...
                    mCat20,fYr1, fYr2,fTw,nTbin(i),nN(1));
                [mBeta(i,iYY,iCat), mBetaprob(i,iYY,iCat)] =calc_beta(mCat,mCat00,...
                    mCat20,fYr1, fYr2,fTw,nTbin(i),nN(1));
            end
        end
        %         disp(n);
        %     subplot(2,2,1)
        %     hold on;plot(1./nTbin1,mLTA);
        %     subplot(2,2,2)
        %     hold on;plot(1./nTbin1,mLTAprob);
        %     subplot(2,2,3)
        %     hold on;plot(1./nTbin1,mBeta);
        %     subplot(2,2,4)
        %     hold on;plot(1./nTbin1,mBetaprob);

    end
end

for i=1:size(mLTA,1)
    mLTA1(:,i)=reshape(mLTA(i,:,:),size(mLTA,2)*size(mLTA,3),1);
    mLTAprob1(:,i)=reshape(mLTAprob(i,:,:),size(mLTAprob,2)*size(mLTAprob,3),1);
    mBeta1(:,i)=reshape(mBeta(i,:,:),size(mBeta,2)*size(mBeta,3),1);
    mBetaprob1(:,i)=reshape(mBetaprob(i,:,:),size(mBetaprob,2)*size(mBetaprob,3),1);
end

save test.mat
figure;
subplot(2,1,1)
errorbar(1./nTbin1,mean(mLTA1),std(mLTA1),'b--x','LineWidth',3);
hold on;
errorbar(1./nTbin1,mean(mBeta1),std(mBeta1),'r:o','LineWidth',2);
xlabel('Bin Size [1/X]');
ylabel('Z and \beta [ ]');
legend('z','\beta');
subplot(2,1,2);
errorbar(1./nTbin1,mean(mLTAprob1),std(mLTAprob1),'b--x','LineWidth',3);
YLim([0,1]);
hold on;
errorbar(1./nTbin1,mean(mBetaprob1),std(mBetaprob1),'r:o','LineWidth',2);
xlabel('Bin Size [1/X]');
ylabel('Probability of Z and \beta [ ]');
legend('p(z)','p(\beta)');

sPrint=sprintf('print -dpng -r300 Tbin_N-%04.0f_Sim-%04.0f.png',nN,iCat);
eval(sPrint);
 disp(sPrint);


text:
function XYZ = rspd2xyz(varargin)
%RSPD2XYZ Convert relative spectral power density to XYZ.
%   RSPD2XYZ(lambda,phi) converts a relative spectral power density
%   function to 1931 CIE XYZ tristimulus values. The inputs, lambda and
%   phi, are both vectors with the same length. The input lambda is the
%   wavelength (in nm). The input phi is the relative spectral power
%   density.
%
%   RSPD2XYZ(S) converts a relative spectral power density function
%   represented as either a Px2 numeric matrix or a table. If S is a
%   numeric matrix, then S(:,1) is the wavelength, lambda, and S(:,2) is
%   phi. If S is a table with a variable called "lambda", then that
%   variable is used as lambda and the table's other variable is used as
%   phi. If the table does not have a variable named "lambda", then the
%   table's first variable is used as lambda, and the second variable is
%   used as phi.
%
%   RSPD2XYZ(I,R) converts an illumination spectral power density curve
%   (I) and a relative reflectance curve (R) to tristimulus values. I
%   and R can each be either a Px2 matrix or a table.

%   Written by Steve Eddins to accompany Digital Image Processing Using
%   MATLAB, 3rd edition, Gatesmark Press, 2020,
%   http://imageprocessingplace.com.
%
%   Copyright 2019 The MathWorks, Inc.
%   License: https://github.com/mathworks/matlab-color-tools/blob/master/license.txt

narginchk(1,2);

if nargin == 2
   if isvector(varargin{1}) && isvector(varargin{2})
      % rspd2xyz(lambda,phi)
      % Treat the input as I (the illumination curve). Construct the
      % reflectance curve as a constant from lambda = 300 to lambda =
      % 830.
      lambda = varargin{1};
      phi = varargin{2};
      I = [lambda(:) phi(:)];
      R = [300 1; 830 1];
   else
      % rspd2xyz(I,R)
      I = varargin{1};
      R = varargin{2};
   end
else
   % rspd2xyz(S)
   % Treat the input as I (the illumination curve). Construct the
   % reflectance curve as a constant from lambda = 300 to lambda = 830.
   I = varargin{1};
   R = [300 1; 830 1];
end

% Convert the inputs, in either Px2 matrix or table form, to function
% handles that return phi(lambda) for any lambda using spline
% interpolation and constant extrapolation.
fI = phiFunction(I);
fR = phiFunction(R);

matching_fcns = colorMatchingFunctions;

L = matching_fcns.lambda;
x = matching_fcns.x;
y = matching_fcns.y;
z = matching_fcns.z;
x = x/(100*sum(x));
y = y/(100*sum(y));
z = z/(100*sum(z));

% Perform a simple numerical integration.
XYZ = sum(fI(L) .* fR(L) .* [x y z],1);

%----------------------------------------------------------------------%
function f = phiFunction(S)
%   Create a function handle that returns phi(lambda). The input can be
%   a Px2 numeric matrix or a table as described in the documentation
%   above for rspd2xyz.

if istable(S)
   % Look for table variable called lambda.
   idx = find(S.Properties.VariableNames == "lambda",1);
   if isempty(idx)
      % The table does not have a lambda variable. Assume the first
      % variable is lambda and the second variable is phi.
      lambda = S{:,1};
      phi = S{:,2};
   else
      lambda = S{:,idx};

      % Use the first non-lambda table variable as phi.
      kk = 1:width(S);
      kk(idx) = [];
      phi = S{:,kk(1)};
   end
else
   % S is a matrix. Use the first column as lambda and the second
   % column as phi.
   S = double(S);
   lambda = S(:,1);
   phi = S(:,2);
end

% Add zeros to the beginning and end of phi. When using the spline
% function for cubic spline interpolation, this has the effect of
% enabling constant extrapolation using the first and last values of
% phi. See section 7.2.2.1, Extrapolation, CIE 15:2004, Technical Report
% on Colorimetry.
phi = [0 ; phi ; 0];

% Create a function handle that interpolates phi using cubic spline
% interpolation.  Cubic spline is one of the recommended interpolation
% methods. See section 7.2.1.1, Interpolation, CIE 15:2004, Technical
% Report on Colorimetry.
f = @(lambda_q) spline(lambda,phi,lambda_q);

author: dipum

repo: dipum-toolbox

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


text:
function [status] = isplottingfunction(function_name)

% ISPLOTTINGFUNCTION is a helper function for reproducescript, and
% is used for the cfg.reproducescript functionality. It compares the input
% function name with the list of known FieldTrip plotting functions and
% returns 1 if it is a plotting function, and 0 otherwise.

% determine FieldTrip path
[ftver, ftpath] = ft_version;

functionlist = [];
% browse the main FieldTrip directory
functionlist = cat(1, functionlist, dir(fullfile(ftpath, '*plot*.m')));
functionlist = cat(1, functionlist, dir(fullfile(ftpath, '*movie*.m')));
% browse the plotting directory
functionlist = cat(1, functionlist, dir(fullfile(ftpath, 'plotting', '*plot*.m')));
functionlist = cat(1, functionlist, dir(fullfile(ftpath, 'plotting', '*movie*.m')));

% concatenate all function names
functionlist = struct2cell(functionlist);
functionlist(2:end,:) = [];

% compare the function with the list of plotting functions
isinlist = strfind(functionlist, function_name);
status = any(cat(1,isinlist{:}));


text:
clear, clc;

% This is an example for running the function epp
%
%  Problem:
%
%  min  1/2 || x - y||^2 + rho * \|x\|_q
%
%  which is a subproblem of the Euclidean projection:
%
%  min \sum_i 1/2 || X_i - Y_i||^2 + rho * \sum_i \|X_i\|_q
%
% For detailed description of the function, please refer to the Manual.
%
%% ------------   History --------------------
% First version on August 10, 2009.
%
%
% For any problem, please contact Jun Liu (j.liu@asu.edu)

cd ..
cd ..

root=cd;
addpath(genpath([root '/SLEP']));
                     % add the functions in the folder SLEP to the path

% change to the original folder
cd Examples/L1Lq;

n=100;
v=randn(n, 1);
lambda=1;
q=3;
c0=0.1;

% run the function epp
[x, c, iter_step]=epp(v, n, lambda, q, c0);

author: zzstefan

repo: BrainNetClass


text:
function webKeggMapper(id,str)
% Small script that links to kegg mapper.
%
% USAGE:
%
%    webKeggMapper(id,str)
%
% INPUTS:
%    id:    Either 'ec' or 'kegg'
%    str:   ID of themetabolite / reaction
%
% .. Author: - Stefan G. Thorleifsson March 2011
% .. In case of changes of the website  I keep this script.
%
% .. rBioNet is published under GNU GENERAL PUBLIC LICENSE 3.0+
% .. Thorleifsson, S. G., Thiele, I., rBioNet: A COBRA toolbox extension for
% .. reconstructing high-quality biochemical networks, Bioinformatics, Accepted.
% .. rbionet@systemsbiology.is

if isempty(str) %No id
    msgbox(['There is no ' id ' ID available for this metabolite/reaction.'],'Help','help');
    return;
end

switch id
    case 'ec'
        web(['http://www.genome.jp/kegg-bin/show_pathway?ec01100+' str],'-browser');
    case 'kegg'
        web(['http://www.genome.jp/kegg-bin/show_pathway?rn01100+' str],'-browser');
    otherwise
        msgbox('unknown source');
end

author: opencobra

repo: cobratoolbox


text:
function [C, sigma] = dataset3Params(X, y, Xval, yval)
%EX6PARAMS returns your choice of C and sigma for Part 3 of the exercise
%where you select the optimal (C, sigma) learning parameters to use for SVM
%with RBF kernel
%   [C, sigma] = EX6PARAMS(X, y, Xval, yval) returns your choice of C and
%   sigma. You should complete this function to return the optimal C and
%   sigma based on a cross-validation set.
%

% You need to return the following variables correctly.
C = 1;
sigma = 0.3;

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return the optimal C and sigma
%               learning parameters found using the cross validation set.
%               You can use svmPredict to predict the labels on the cross
%               validation set. For example,
%                   predictions = svmPredict(model, Xval);
%               will return the predictions on the cross validation set.
%
%  Note: You can compute the prediction error using
%        mean(double(predictions ~= yval))
%

possibles = [0.01 0.03 0.1 0.3 1 3 10 30];
error = realmax();
for i=1:size(possibles, 2)
	for j=1:size(possibles, 2)
		C_pos = possibles(i);
		sigma_pos = possibles(j);

		model = svmTrain(X, y, C_pos, @(x1, x2) gaussianKernel(x1, x2, sigma_pos));
		predictions = svmPredict(model, Xval);
		new_error = mean(double(predictions ~= yval));

		if new_error < error
			error = new_error;
			C = C_pos;
			sigma = sigma_pos;
		end
	end
end

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

end

author: gopaczewski

repo: coursera-ml


text:
function h = cmpndNoisePointPlot(noise, X, y, fontName, fontSize, markerSize, lineWidth)

% CMPNDNOISEPOINTPLOT Plot the data-points for the CMPND noise model.
% FORMAT
% DESC plots the data point locations for the compound
%  noise structure.
% ARG noise : the noise structure which is to be plotted.
% ARG X : input locations to be plotted.
% ARG y : target locations to be plotted.
% ARG fontName : name of any fonts to be used.
% ARG fontSize : size of any fonts to be used.
% ARG markerSize : size of the makers to be used.
% ARG lineWidth : size of any lines to be plotted.
%
% SEEALSO : cmpndNoiseParamInit, noisePointPlot
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005

% NOISE

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ y, m, d ] = day_carry_roman ( y, m, d )

%*****************************************************************************80
%
%% DAY_CARRY_ROMAN carries days to months in a Roman date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, M, D, the YMD date.
%    On output, D is between 1 and the number of days in M.
%
  days = month_length_roman ( y, m );

  while ( days < d )

    d = d - days;
    m = m + 1;
    days = month_length_roman ( y, m );
%
%  Make sure the month isn't too big.
%
    [ y, m ] = month_carry_roman ( y, m );

  end

  return
end


text:
function z = log_sigmoid(x)
% computes and returns the logistic sigmoid of x. Operation is performed
% element-wise if x is a vector or matrix.
%
%       z = 1 ./ (1-exp(-x))

z = (1+exp(-x)).^-1;

end

author: OHBA-analysis

repo: HMM-MAR


text:
function M = fixedrankfactory_2factors_preconditioned(m, n, k)
% Manifold of m-by-n matrices of rank k with two factor quotient geometry.
%
% function M = fixedrankfactory_2factors_preconditioned(m, n, k)
%
% This geometry is tuned to least-squares problems such as low-rank matrix
% completion with ell-2 loss.
%
% A point X on the manifold is represented as a structure with two
% fields: L and R. The matrices L (m-by-k) and R (n-by-k) are
% full column-rank matrices such that X = L*R'.
%
% Tangent vectors are represented as a structure with two fields: L, R.
%
% Please cite the Manopt paper as well as the research paper:
%     @Techreport{mishra2012optimized,
%       Title   = {A {R}iemannian geometry for low-rank matrix completion},
%       Author  = {Mishra, B. and Adithya Apuroop, K. and Sepulchre, R.},
%       Journal = {Arxiv preprint arXiv:1211.1550},
%       Year    = {2012}
%     }
%
%
% See also: fixedrankembeddedfactory fixedrankfactory_2factors fixedrankfactory_3factors_preconditioned

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, Dec. 30, 2012.
% Contributors:
% Change log:
%
%	April 04, 2015 (BM):
%      Cosmetic changes including avoiding storing the inverse of a
%       k-by-k matrix.

    M.name = @() sprintf('LR''(tuned to least square problems) quotient manifold of %dx%d matrices of rank %d', m, n, k);

    M.dim = @() (m+n-k)*k;

    % Some precomputations at the point X to be used in the inner product
    % (and pretty much everywhere else).
    function X = prepare(X)
        if ~all(isfield(X,{'LtL','RtR'}))
            L = X.L;
            R = X.R;
            X.LtL = L'*L;
            X.RtR = R'*R;
        end
    end

    % The choice of metric is motivated by symmetry and
    % tuned to least-squares cost function.
    M.inner = @iproduct;
    function ip = iproduct(X, eta, zeta)
        X = prepare(X);
        ip = trace(X.RtR*(eta.L'*zeta.L)) + trace(X.LtL*(eta.R'*zeta.R)); % Scaled metric
    end

    M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));

    M.dist = @(x, y) error('fixedrankfactory_2factors_preconditioned.dist not implemented yet.');

    M.typicaldist = @() 10*k;

    M.egrad2rgrad = @egrad2rgrad;
    function rgrad = egrad2rgrad(X, egrad)
        X = prepare(X);

        % Riemannian gradient
        rgrad.L = egrad.L/X.RtR;
        rgrad.R = egrad.R/X.LtL;
    end

    M.ehess2rhess = @ehess2rhess;
    function Hess = ehess2rhess(X, egrad, ehess, eta)
        X = prepare(X);

        % Riemannian gradient.
        rgrad = egrad2rgrad(X, egrad);

        % Directional derivative of the Riemannian gradient.
        Hess.L = ehess.L/X.RtR - 2*egrad.L*(X.RtR \ (symm(eta.R'*X.R) / X.RtR));
        Hess.R = ehess.R/X.LtL - 2*egrad.R*(X.LtL \ (symm(eta.L'*X.L) / X.LtL));

        % We still need a correction factor for the non-constant metric.
        Hess.L = Hess.L + rgrad.L*(symm(eta.R'*X.R)/X.RtR) + eta.L*(symm(rgrad.R'*X.R)/X.RtR) - X.L*(symm(eta.R'*rgrad.R)/X.RtR);
        Hess.R = Hess.R + rgrad.R*(symm(eta.L'*X.L)/X.LtL) + eta.R*(symm(rgrad.L'*X.L)/X.LtL) - X.R*(symm(eta.L'*rgrad.L)/X.LtL);

        % Project on the horizontal space.
        Hess = M.proj(X, Hess);
    end

    M.proj = @projection;
    function etaproj = projection(X, eta)
        X = prepare(X);

        % Projection onto the horizontal space.
        Lambda = 0.5*((eta.R'*X.R)/X.RtR  -   X.LtL\(X.L'*eta.L));
        etaproj.L = eta.L + X.L*Lambda;
        etaproj.R = eta.R - X.R*Lambda';
    end

    M.tangent = M.proj;

    M.tangent2ambient = @(X, eta) eta;

    M.retr = @retraction;
    function Y = retraction(X, eta, t)
        if nargin < 3
            t = 1.0;
        end
        Y.L = X.L + t*eta.L;
        Y.R = X.R + t*eta.R;

        % Numerical conditioning step: a simpler version.
        % We need to ensure that L and R are do not have very relative
        % skewed norms.

        scaling = norm(X.L, 'fro')/norm(X.R, 'fro');
        scaling = sqrt(scaling);
        Y.L = Y.L / scaling;
        Y.R = Y.R * scaling;

        % These are reused in the computations of gradient and Hessian.
        Y = prepare(Y);
    end

    M.exp = @exponential;
    function Y = exponential(X, eta, t)
        if nargin < 3
            t = 1.0;
        end

        Y = retraction(X, eta, t);
        warning('manopt:fixedrankfactory_2factors_preconditioned:exp', ...
            ['Exponential for fixed rank ' ...
            'manifold not implemented yet. Used retraction instead.']);
    end

    M.hash = @(X) ['z' hashmd5([X.L(:) ; X.R(:)])];

    M.rand = @random;

    function X = random()
        X.L = randn(m, k);
        X.R = randn(n, k);
    end

    M.randvec = @randomvec;
    function eta = randomvec(X)
        eta.L = randn(m, k);
        eta.R = randn(n, k);
        eta = projection(X, eta);
        nrm = M.norm(X, eta);
        eta.L = eta.L / nrm;
        eta.R = eta.R / nrm;
    end

    M.lincomb = @lincomb;

    M.zerovec = @(X) struct('L', zeros(m, k),'R', zeros(n, k));

    M.transp = @(x1, x2, d) projection(x2, d);

    % vec and mat are not isometries, because of the scaled inner metric.
    M.vec = @(X, U) [U.L(:) ; U.R(:)];

    M.mat = @(X, u) struct('L', reshape(u(1:(m*k)), m, k), ...
        'R', reshape(u((m*k+1):end), n, k));

    M.vecmatareisometries = @() false;

    % Auxiliary functions
    symm = @(M) .5*(M+M');
end

% Linear combination of tangent vectors.
function d = lincomb(x, a1, d1, a2, d2) %#ok<INUSL>

    if nargin == 3
        d.L = a1*d1.L;
        d.R = a1*d1.R;
    elseif nargin == 5
        d.L = a1*d1.L + a2*d2.L;
        d.R = a1*d1.R + a2*d2.R;
    else
        error('Bad use of fixedrankfactory_2factors_preconditioned.lincomb.');
    end

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function out = ggrnd(n,p,g)
% GG - random numbers from generalized Gaussian distribution (requires
% Statistics toolbox)
%   R = GGRND(N,P,G) returns a vector of N random numbers having a
%   generalized gaussian distribution with shape parameter P, and scale
%   parameter G, i.e. p(x) \propto exp(-G*|x|^P).
%
%   If P=2 it GGRND returns normally distributed random numbers
%
%   Reference: Johnson, M. E., Computer generation of the exponential power
%   distributions", Journal of Statistical Computation and Simulation, 9,
%   pp. 239--240, 1979
%
%   See also: RANDN, GAMRND

%   Author: G. Gomez-Herrero, german.gomezherrero@ieee.org
%   $Revision: 1.0 $Date: 21/06/2007

if nargin < 3, g=1; end
if abs(p-2)<1e-3,
    x = randn(1,n);
else
    y = gamrnd(1/p,1,1,n);
    x = abs(g.*y).^(1/p);
    pos = find(rand(1,n)<0.5);
    x(pos) = -x(pos);
end
% normalize variance and return
out = x'./sqrt(var(x));
return;

repo: REST


text:
function GUI = sessionGUI_plotVoxelData(rois, scans, viewer);
%
%  GUI = sessionGUI_plotVoxelData(<rois>, <scans=scan group>, <viewer>);
%
% Open a Multi Voxel UI(s) for the selected ROIs.
%
% If a viewer is open, will open an MVUI for the currently-selected ROI
% in the selected viewer. If no viewers are open, will open it for the
% selected ROI in the ROI listbox in the session GUI. These are both
% over-ridden if an roi list and/or viewer are manually specified.
%
%
%
% ras, 07/06.
mrGlobals2;

if notDefined('viewer')
    if ~isempty(GUI.viewers)
        viewer = GUI.viewers(GUI.settings.viewer);
    else
        viewer = 0;
    end
end

if notDefined('rois')
    if viewer==0 % no viewer, get from saved list
        rois = guiGet('selectedRois');
    else
        rois = {mrViewGet(viewer, 'selectedRoi')};
    end
end

if ~checkfields(GUI, 'controls', 'mvui')
    GUI.controls.mvui = [];
end

% check if an event-related scan group is assigned.
% If so, use that; otherwise use the GUI selected scans.
if notDefined('scans')
	scans = GUI.settings.scan;
	dt = GUI.settings.dataType;
	if isfield(dataTYPES(dt).scanParams(scans(1)), 'scanGroup') & ...
			~isempty(dataTYPES(dt).scanParams(scans(1)).scanGroup)
		[scans dt] = er_getScanGroup(INPLANE{1});
	end
end

if notDefined('dt')
	dt = guiGet('dataType');
end

for i = 1:length(rois)
    mv = multiVoxelUI(INPLANE{1}, rois{i}, scans, dt);
    GUI.controls.mvui(end+1) = mv.ui.fig;
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% LAPPLOT -  Compute the discrete laplacian of EEG scalp distribution(s)
%
% Usage:
%   >> laplace = lapplot(map,eloc_file,draw)
%
% Inputs:
%    map        - Activity levels, size (nelectrodes,nmaps)
%    eloc_file	- Electrode location filename (.loc file)
%                 For format, see  >> topoplot example
%    draw       - If defined, draw the map(s) {default: no}
%
% Output:
%    laplace    - Laplacian map, size (nelectrodes,nmaps)
%
% Note: uses DEL2
%
% Author: Scott Makeig, SCCN/INC/UCSD, La Jolla, 1998
%
% See also: TOPOPLOT, GRADPLOT

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

% 01-25-02 reformated help & license, added links -ad

function [laplac] = lapplot(map,filename,draw)

if nargin < 2
	help lapplot;
	return;
end

MAXCHANS = size(map,1);
GRID_SCALE = 2*MAXCHANS+5;
MAX_RADIUS = 0.5;

% ---------------------
% Read the channel file
% ---------------------
if ischar( filename )
	fid = fopen(filename);
	locations = fscanf(fid,'%d %f %f %s',[7 MAXCHANS]);
	fclose(fid);
	locations = locations';
	Th = pi/180*locations(:,2);   % convert degrees to rads
	Rd = locations(:,3);
	ii = find(Rd <= MAX_RADIUS); % interpolate on-scalp channels only
	Th = Th(ii);
	Rd = Rd(ii);
	[x,y] = pol2cart(Th,Rd);
else
	x = real(filename);
	y = imag(filename);
end;

% ---------------------------------------------------
% Locate nearest position of an electrode in the grid
% ---------------------------------------------------
xi = linspace(-0.5,0.5,GRID_SCALE);   % x-axis description (row vector)
yi = linspace(-0.5,0.5,GRID_SCALE);   % y-axis description (row vector)
for i=1:MAXCHANS
   [useless_var horizidx(i)] = min(abs(y(i) - xi));    % find pointers to electrode
   [useless_var vertidx(i)] = min(abs(x(i) - yi));     % positions in Zi
end

% -----------------
% Compute laplacian
% -----------------
for i=1:size(map,2)
   	[Xi,Yi,Zi] = griddata(y,x,map(:,i),yi',xi, 'v4');   % interpolate data

   	laplac2D = del2(Zi);
	positions = horizidx + (vertidx-1)*GRID_SCALE;
	laplac(:,i) = laplac2D(positions(:));

	% ------------------
	% Draw laplacian map
	% ------------------
	if exist('draw');
        mask = (sqrt(Xi.^2+Yi.^2) <= MAX_RADIUS);
        laplac2D(find(mask==0)) = NaN;

		subplot(ceil(sqrt(size(map,2))), ceil(sqrt(size(map,2))), i);
		contour(laplac2D);
		title( int2str(i) );

% %%% Draw Head %%%%
ax = axis;
width = ax(2)-ax(1);
axis([ax(1)-width/3 ax(2)+width/3 ax(3)-width/3 ax(4)+width/3])
steps = 0:2*pi/100:2*pi;
basex = .18*MAX_RADIUS;
tip = MAX_RADIUS*1.15;
base = MAX_RADIUS-.004;
EarX = [.497 .510 .518 .5299 .5419 .54 .547 .532 .510 .489];
EarY = [.0555 .0775 .0783 .0746 .0555 -.0055 -.0932 -.1313 -.1384 -.1199];

HCOLOR = 'k';
HLINEWIDTH = 1.8;

% Plot Head, Ears, Nose
hold on
plot(1+width/2+cos(steps).*MAX_RADIUS*width,...
     1+width/2+sin(steps).*MAX_RADIUS*width,...
    'color',HCOLOR,'Linestyle','-','LineWidth',HLINEWIDTH); % head

plot(1+width/2+[.18*MAX_RADIUS*width;0;-.18*MAX_RADIUS*width],...
     1+width/2+[base;tip;base]*width,...
    'Color',HCOLOR,'LineWidth',HLINEWIDTH);                 % nose

plot(1+width/2+EarX*width,...
     1+width/2+EarY*width,...
           'color',HCOLOR,'LineWidth',HLINEWIDTH)           % l ear
plot(1+width/2-EarX*width,...
     1+width/2+EarY*width,...
           'color',HCOLOR,'LineWidth',HLINEWIDTH)           % r ear

hold off
axis off

	end
end;

return;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function planC = copyStructsFromFilesToPlanC(scanFileNames, planC)
% function planC = copyStructsFromFilesToPlanC(scanFileNames, planC)
%
% APA 07/17/2013

if ~exist('planC','var')
    global planC
end

indexS = planC{end};

[xValsBase, yValsBase, zValsBase] = getScanXYZVals(planC{indexS.scan});

for scanFileIndex = length(scanFileNames)
    scanBasePlanC = loadPlanC(scanFileNames{scanFileIndex},tempdir);
    indexSbaseScan = scanBasePlanC{end};
    annotROIIndV = strcmpi('Annotation ROI',{scanBasePlanC{indexSbaseScan.structures}.structureName});
    annotStrV = find(annotROIIndV);
    [xValsBase1, yValsBase1, zValsBase1] = getScanXYZVals(scanBasePlanC{indexSbaseScan.scan});
    for structIndex = 1:length(annotStrV)
        structNum = annotStrV(structIndex);
        sliceNumsV = [];
        clear pointsC
        for sliceNum = 1:length(scanBasePlanC{indexSbaseScan.structures}(structNum).contour)
            for segNum = 1:length(scanBasePlanC{indexSbaseScan.structures}(structNum).contour(sliceNum).segments)
                points = scanBasePlanC{indexSbaseScan.structures}(structNum).contour(sliceNum).segments(segNum).points;
                if ~isempty(points)
                    zValue = points(1,3);
                    newSliceNum = findnearest(zValsBase,zValue);
                    sliceNumsV = [sliceNumsV newSliceNum];
                    points(:,3) = zValsBase(newSliceNum);
                    pointsC{newSliceNum} = points;
                end
            end
        end

        [sliceNumsV, indUniq] = unique(sliceNumsV);

        % Create Structures segments on sacn slices
        newStructS = newCERRStructure(1, planC);
        for slcNum = 1:length(zValsBase)
            if ismember(slcNum,sliceNumsV)
                newStructS.contour(slcNum).segments(1).points = pointsC{slcNum};
            else
                newStructS.contour(slcNum).segments.points = [];
            end
        end

        newStructNum = length(planC{indexS.structures}) + 1;
        newStructS.structureName = 'Annotation ROI';

        planC{indexS.structures} = dissimilarInsert(planC{indexS.structures}, newStructS, newStructNum);
        planC = getRasterSegs(planC, newStructNum);
        planC = updateStructureMatrices(planC, newStructNum, sliceNumsV);
    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [finalVectors] = detMinSpan(model, params, vectors)
% Calculates the MinSpan vectors for a COBRA model. The
% algorithm determines a set of linearly independent basis vectors that
% span the nullspace that meet the criteria of the flux bounds of the
% network and minimizes the number of reactions used. The algorithm
% operates in an interative manner and checkes for convergence after each
% iteration. Parameters may be provided to skip convergence and terminate
% the problem when models are large. See Bordbar et al. Mol Syst Biol 2014
% for more details. This algorithm has only been tested and requires Gurobi
% as MILP solver.
%
% INPUTS:
%    model:             COBRA model structure (requires S, lb, ub)
%                       Note: MinSpan calculations are done w/o biomass
%                       reaction and all reactions must be able to carry a
%                       flux = 0 for the trivial solution to be feasible.
%                       model is auto corrected by bounds but biomass
%                       must be removed manually
%
% OPTIONAL INPUTS:
%    params:            Optional parameters to calculate MinSpan.
%                       Determining the MinSpan is a NP hard calculation.
%                       For large models, the algorithm may not converge
%                       and parameters must be provided to stop the
%                       algorithm to provide an approximate solution.
%
%                       * .coverage - Number of iterations to run the algorithm, if not
%                       converge (Default = 10)
%                       * .timeLimit - Time to spend on each MinSpan calculation (sec),
%                       (Default = 30)
%                       * .saveIntV - Save intermediate vectors in order to restart from
%                       latest iteration (Default = 0)
%                       * .cores - Number of cores to use (Default = 1)
%
%    vectors:           Set of intermediate MinSpan vectors that may
%                       have not yet reached convergence, allowing to
%                       pickup calculation from last spot.
%
% OUTPUTS:
%    finalVectors:      MinSpan vectors for COBRA model
%
% .. Author: Aarash Bordbar 05/15/2017
%            Ronan Fleming, nullspace computation with LUSOL

global CBT_MILP_SOLVER
global CBT_LP_SOLVER
global CBT_QP_SOLVER

if ~strcmp(CBT_MILP_SOLVER, 'gurobi') || ~strcmp(CBT_LP_SOLVER, 'gurobi') || ~strcmp(CBT_QP_SOLVER, 'gurobi')
    error('detMinSpan only runs with Gurobi.\nTry to run `changeCobraSolver(''gurobi'', ''MILP'')` to use Gurobi.');
end

% Ensure model has correct bounds
origModel = model;
model.lb(model.lb < 0) = -1000;
model.ub(model.ub > 0) = 1000;
model.lb(model.lb > 0) = 0;
model.ub(model.ub < 0) = 0;

% Reduce model to just reactions that can carry a flux
[minF, maxF] = fluxVariability(model, 0);
minF(abs(minF) < 1e-8) = 0;
maxF(abs(maxF) < 1e-8) = 0;
remRxns = model.rxns(minF == 0 & maxF == 0);
model = removeRxns(model, remRxns);

% Setup MILP and solving parameters
if ~exist('params', 'var')
    params.coverage = 10;
    params.timeLimit = 30;
    params.saveIntV = 0;
    params.cores = 1;
    params.nullAlg = 'matlab';
else
    if ~isfield(params, 'coverage')
        params.coverage = 10;
    end
    if ~isfield(params, 'timeLimit')
        params.timeLimit = 30;
    end
    if ~isfield(params, 'saveIntV')
        params.saveIntV = 1;
    end
    if ~isfield(params, 'cores')
        params.cores = 1;
    end
    if ~isfield(params, 'nullAlg')
        params.nullAlg = 'matlab';
    end
end

MILPparams.TimeLimit = params.timeLimit;
MILPparams.outputFlag = 1;
MILPparams.Presolve = 2;
MILPparams.Threads = params.cores;
MILPparams.DisplayInterval = 10;

if ~exist('vectors', 'var')
    if strcmp('params.nullAlg','lusol')
        error('nearly but not quite implemented yet')
        [vectors, ~] = getNullSpace(model.S, 0);
    else
        vectors = null(full(model.S));
    end
end

% Prepratory steps for MinSpan determination
rng('shuffle');

if strcmp('params.nullAlg','lusol')
    [N, ~] = getNullSpace(S, 0);
else
    N = null(full(model.S));
end

[m, n] = size(model.S);

locRev = find(model.lb < 0);
lengthRev = length(locRev);
revConstraintMat = zeros(lengthRev, n);
for i = 1:lengthRev
    revConstraintMat(i, locRev(i)) = 1;
end

% Run MinSpan
tmpNvProd = [];
totalNNZ = [];
for k = 1:params.coverage

    numToCheck = randperm(size(N, 2));
    prevNNZ = nnz(vectors);

    for i = 1:length(numToCheck)
        tic;
        oldPath = vectors(:, numToCheck(i));
        pathLength = nnz(vectors(:, numToCheck(i)));
        vectors(:, numToCheck(i)) = zeros(n, 1);

        sizeN = 1;
        theta = N \ vectors;

        if strcmp('params.nullAlg','lusol')
            [Z, ~] = getNullSpace(theta', 0);
            tmpN = sparse(N * Z);
        else
            tmpN = sparse(N * null(theta'));
        end

        tmptmpNprod = tmpN' * oldPath;
        tmpN = tmpN * (1 / tmptmpNprod);

        % Model Formulation: A, b, csense, lb, ub, vartype, c
        MILPproblem.A = [model.S, sparse(m, n + 2 * sizeN);  % S. v = 0
                         revConstraintMat, 1e4 * revConstraintMat, sparse(lengthRev, 2 * sizeN);  % v >= -10000*b
                         speye(n), -1e4 * speye(n), sparse(n, 2 * sizeN);  % v <= 10000*b
                         tmpN', sparse(sizeN, n), (-1001) * speye(sizeN), sparse(sizeN, sizeN);  % N.v >= 1001*fi+ - 1000
                         -tmpN', sparse(sizeN, n), sparse(sizeN, sizeN), (-1001) * speye(sizeN);  % -N.v >= 1001*fi- - 1000
                         sparse(1, 2 * n), ones(1, 2 * sizeN);  % sum(fi+, fi-) >= 1
                         ];

        MILPproblem.b = [zeros(m, 1);  % S. v = 0
                         zeros(lengthRev, 1);  % v >= -10000*b
                         zeros(n, 1);  % v <= 10000*b
                         -1000 * ones(sizeN, 1);  % N.v >= 1000*fi+ - 1000
                         -1000 * ones(sizeN, 1);  % -N.v >= 1000*fi+ - 1000
                         1  % sum(vi+, vi-) >= 1
                         ];

        MILPproblem.csense = '';
        for l = 1:m
            MILPproblem.csense(end + 1, 1) = 'E';  % S.v = 0
        end
        for l = 1:lengthRev
            MILPproblem.csense(end + 1, 1) = 'G';  % v >= -10000*b
        end
        for l = 1:n
            MILPproblem.csense(end + 1, 1) = 'L';  % v <= 10000*b
        end
        for l = 1:sizeN
            MILPproblem.csense(end + 1, 1) = 'G';  % N.v >= 1000*fi+ - 1000
        end
        for l = 1:sizeN
            MILPproblem.csense(end + 1, 1) = 'G';  % -N.v >= 1000*fi+ - 1000
        end
        MILPproblem.csense(end + 1, 1) = 'G';  % sum(vi+, vi-) > 1

        MILPproblem.lb = [model.lb;  % v
                          zeros(n, 1);  % a
                          zeros(2 * sizeN, 1)];  % fi+, fi-

        MILPproblem.ub = [model.ub;  % v
                          ones(n, 1);  % a
                          ones(2 * sizeN, 1)];  % k+, k-

        MILPproblem.vartype = '';
        for l = 1:n
            MILPproblem.vartype(end + 1, 1) = 'C';  % v
        end
        for l = 1:n
            MILPproblem.vartype(end + 1, 1) = 'B';  % b
        end
        for l = 1:2 * sizeN
            MILPproblem.vartype(end + 1, 1) = 'B';  % fi+, fi-
        end

        MILPproblem.c = [zeros(n, 1);  % v
                         ones(n, 1);  % b
                         zeros(2 * sizeN, 1);  % fi+, fi-
                         ];

        MILPproblem.osense = 1;  % minimize

        % Setup initial solution
        binOldPath = zeros(length(oldPath), 1);
        binOldPath(find(oldPath)) = 1;
        MILPproblem.x0 = [oldPath; binOldPath; 1e101; 1e101];

        MILPsolution = solveCobraMILP(MILPproblem, MILPparams);

        % Check solution
        % If unable to find solution, break iteration
        if length(MILPsolution.cont) < n
            break
        end

        % If solution found, normalize and replace vector in intermediate
        % matrix
        vector = MILPsolution.full(1:n);
        vector(abs(vector) < 1e-6) = 0;

        if strcmp('params.nullAlg','lusol')
            [Z, ~] = getNullSpace((N \ [vectors, vector])', 0);
            tmpNullCheck = N * Z;
        else
            tmpNullCheck = N * null((N \ [vectors, vector])');
        end

        if nnz(vector) > 0 && isempty(tmpNullCheck)
            vector = vector / norm(vector);
            vectors(:, numToCheck(i)) = vector;
        else
            vectors(:, numToCheck(i)) = oldPath;
            vector = oldPath;
        end
        tmpNvProd = [tmpNvProd; tmpN' * vector];
        totalNNZ = [totalNNZ; nnz(vectors)];

        time(i, 1) = toc;

        % Save intermediate matrices (within iteration)
        if params.saveIntV == 1
            filename = strcat('save_', num2str(k), '_', num2str(i));
            save(filename, 'MILPproblem', 'MILPsolution', 'vectors', ...
                 'time', 'tmpNvProd', 'totalNNZ', 'numToCheck');
        end

        clear mex
    end

    newNum = nnz(vectors);

    % Save intermediate matrices (after a completed iteration)
    if params.saveIntV == 1
        filename = strcat('save_finalround_', num2str(k));
        save(filename, 'vectors', 'time');
    end

    % If MinSpan solution has converged (same NNZ as previous iteration)
    if newNum == prevNNZ
        break
    end
end

clear x
for i = 1:size(vectors, 2)
    x(i, 1) = nnz(vectors(:, i));
end

completedPaths = find(x < n & x > 0);
vectors = vectors(:, completedPaths);

% Normalize vectors such that smallest flux value in vector is 1
for i = 1:size(vectors, 2)
    loc = find(vectors(:, i));
    tmp = min(abs(vectors(loc, i)));
    vectors(:, i) = vectors(:, i) / tmp;
end

% Cast vectors from reduce model size to full model size
finalVectors = zeros(length(origModel.rxns), size(vectors, 2));
loc = find(ismember(origModel.rxns, model.rxns));
for i = 1:size(vectors, 2)
    finalVectors(loc, i) = vectors(:, i);
end


text:
function w = waveform(TC)

%WAVEFORM Get waveforms from a threecomp object
%   W = WAVEFORM(TC) Extracts waveforms from a threecomp object where W is
%   a waveform object and TC is a threecomp object.

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

w = TC.traces;

author: geoscience-community-codes

repo: GISMO


text:
function [sma] = computeSMAFromPeriod(period, gmu)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

    sma = (((period./(2.*pi)).^2).*gmu).^(1/3);

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function s = weekday_to_name_republican ( w )

%*****************************************************************************80
%
%% WEEKDAY_TO_NAME_REPUBLICAN returns the name of a Republican weekday.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer W, the weekday index.
%
%    Output, string S, the weekday name.
%
  weekday_name = [ ...
    'Primedi '; ...
    'Duodi   '; ...
    'Tridi   '; ...
    'Quartidi'; ...
    'Quintidi'; ...
    'Sextidi '; ...
    'Septidi '; ...
    'Octidi  '; ...
    'Nonidi  '; ...
    'Decadi  ' ];

  if ( w < 1 || 10 < w )
    s = '?????';
  else
    s = weekday_name(w,:);
  end

  return
end

repo: jburkardt-m


text:
classdef CEC2017_F8 < 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
    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{8};
            obj.M = 1;
            if isempty(obj.D); obj.D = 10; end
            obj.D = min(obj.D,length(obj.O));
            obj.lower    = zeros(1,obj.D) - 100;
            obj.upper    = zeros(1,obj.D) + 100;
            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 = max(Z,[],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(:,1:2:end);
            W = Z(:,2:2:end);
            PopCon(:,1) = abs(sum(cumsum(Y,2).^2,2)) - 1e-4;
            PopCon(:,2) = abs(sum(cumsum(W,2).^2,2)) - 1e-4;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
 function ftab = de_ftab(xrs, mas, varargin)
%function ftab = de_ftab(xrs, mas, [options])
%|
%| For multiple-kVp X-ray imaging, we must evaluate functions f_m(s_1, ..., s_L)
%| for m=1,..,M, where M is the number of kVp settings,
%| and L is the number of material components.
%| f_m(s1, s2) = -log( \int exp(- (m1(E) s1 + m2(E) s2)) dI_m(E) / I_m(E) )
%| (For "dual-energy" imaging, M = L = 2.)
%|
%| This routine builds tables/models of f_m, its inverse, and its derivatives.
%| These are needed for (dual-energy) X-ray CT reconstruction.
%|
%| in
%|	xrs	strum	X-ray spectra; see xray_read_spectra.m
%|	mas	strum	mass attenuation coefficients; see xray_read_mac.m
%|
%| option
%|	'sls'	strum	see de_ftab_sls.m
%|	'ftype'	char	fit type for de_ftab_fit() (default: '')
%|	'fit_args' {}	options for de_ftab_fit() (default: {})
%|			e.g.: {'wt', num2cell(ones(M,1))} weighting for fitting
%|	'ctype'	char	curv type for de_ftab_curv() (default: '')
%|	'show'	0|1	visualize
%|
%| out
%|	ftab	strum
%|		data:
%|			.xrs,.mas,.sls,.mac	based on input
%|			.fit	strum with fm approximation methods
%|				see de_ftab_fit.m
%|			.inv1	strum to invert 1 component BH, e.g., water
%|				see de_ftab_inv1.m
%|			.inv2	strum to invert BH by polynomials,
%|				see de_ftab_inv1.m
%|		methods:
%|			.fm_fun(s1, s2, ...) (probably not needed by user)
%|			.plot_fm()
%|			.plot_mac()
%|			.plot_jac()
%|			.plot_inv()
%|
%| Copyright 2008-6-15, Jeff Fessler, University of Michigan

if nargin == 1 && streq(xrs, 'test'), de_ftab_test, return, end
if nargin < 2, ir_usage, end

ftab.show = false;
ftab.sls = struct([]);
ftab.ftype = ''; % defer to de_ftab_fit() default
ftab.fit_args = {}; % defer to de_ftab_fit() default(s)
ftab.ctype = ''; % defer to de_ftab_curv() default ('newt' is fastest)
ftab = vararg_pair(ftab, varargin);

% X-ray spectra
if isempty(xrs)
%	xrs = {'mono,60,100'}; % monoenergetic for testing
	xrs = {'ps1'};
end
if iscell(xrs)
	xrs = xray_read_spectra(xrs{:});
end
if ftab.show
	clf, xrs.plot
prompt
end
ftab.xrs = xrs;

ftab.MM = size(xrs.sp, 2); % # of spectra

% mas: mass attenuation coefficient (strum), at tabulated energies
if isempty(mas)
	mas = {'soft', 'bone'};	% default materials
end
if iscell(mas)
	mas = xray_read_mac(mas);
end
if ftab.show
	clf, mas.plot('kev', xrs.en)
prompt
end
ftab.mas = mas;

% mac: mass attenuation coefficient, at spectrum energies
ftab.mac = xray_make_mac(xrs, mas);
ftab.LL = ncol(ftab.mac.mac);
if ftab.show
	pr ftab.mac.bar % 'bmassml'
%	pr cond(ftab.mac.bar) % examine condition number
	pr sqrt(cond(ftab.mac.bar' * ftab.mac.bar)) % examine condition number
	de_ftab_plot_mac(ftab);
prompt
end

% material integral sampling:
% s_l is samples of the integral of the density of the lth material type.
if isempty(ftab.sls)
	ftab.sls = de_ftab_sls;
end

% build tables of f_m(s1, s2)
fm_fun = @de_ftab_make_fm;
ftab.fm = fm_fun(ftab);
if ftab.show % figure showing f1, f2.  it looks linear, but it not quite!
	de_ftab_plot_fm(ftab, 'down', 4);
prompt
end

% parametric fit to each fm to form a continuous function!
printm 'fit'
if streq(ftab.ftype, 'exp3')
	ftab.fit = de_ftab_fit_exp3(ftab.xrs, ftab.mac.mac, ftab.fit_args{:});
else
	ftab.fit = de_ftab_fit(ftab.sls.sl, ftab.fm, 'type', ftab.ftype, ...
		'mtype', mas.type, ...
		'kev', xrs.en, 'mac', ftab.mac.mac, ... % needed for exp
		'macbar', ftab.mac.bar, ftab.fit_args{:});
% todo: may be bad (and slow) to use finely sampled energies (kev) for exp fit
%		'type', 'exp');
%ftab.fit = de_ftab_fit(ftab.sl, ftab.fm, 'type', 'poly', 'order', 3, 'dc', 1);
end

ftab.fit = de_ftab_curv(ftab.fit, 'ctype', ftab.ctype);

% Jacobian transformation
ftab.T = pinv(ftab.mac.bar); % [L M] transformation, fs = T f approx s !

% Build 1D inverse of 1st material component (usually water)
% for conventional "water only" beam-hardening correction.
ftab.inv1 = de_ftab_inv1(ftab.fit, ftab.sls.sl{1});

% build polynomial inverse approximation
ftab.inv2 = de_ftab_inv2(ftab.fit, ftab.sls.sl, 'T', ftab.T); % todo: options

%ftab.inv.eval = @(ftab,fhat) de_ftab_invert(ftab,fhat);
%ftab.inv = de_ftab_inv_setup(ftab);

meth = {'fm_fun', fm_fun, '(cell: sl)';
	'show_fm', @de_ftab_plot_fm, '()'; % backward compat
	'plot_fm', @de_ftab_plot_fm, '()';
	'plot_mac', @de_ftab_plot_mac, '()';
	'plot_jac', @de_ftab_plot_jac, '()';
	};

ftab = strum(ftab, meth);

ftab = de_ftab_inv_setup(ftab);

if ftab.show % compare fit to sampled fm, etc.
	ftab.inv1.plot(ftab.fit);
	if ~isempty(ftab.inv2)
		ftab.inv2.plot(ftab.fit, ftab.sls.sl);
	end
%	ftab.fit.show_fm(ftab.sls.sl, ftab.fm)
%	ftab.fit.show_err(ftab.sls.sl, ftab.fm)
	ftab.plot_jac;
	ftab.plot_inv;
%keyboard
end

% for safefty, remove table, forcing user to use ftab.fm_fun or ftab.fit.fmfun !
% ftab = rmfield(ftab, 'fm');

end % de_ftab()

%
% de_ftab_inv_setup()
% add methods related to inverse to ftab
%
function ftab = de_ftab_inv_setup(ftab)

printm 'inv: f -> s'

arg.inv = struct; % place holder
fun = @(ftab, fhat) de_ftab_invert(ftab, fhat);
meth = {'inv_fun', fun, '(fhat)';
	'plot_inv', @de_ftab_plot_inv, '()'};
ftab = strum(arg, meth, 'base', ftab);

end % de_ftab_inv_setup()

%
% de_ftab_plot_inv()
% plot inverse, testing with a different set of samples to see errors
%
function dummy = de_ftab_plot_inv(ftab)
dummy = [];
LL = ftab.LL;
MM = ftab.MM;
if LL ~= 2, warn 'L=2 done only', return, end
if MM < LL, warn('plot_inv with MM=%d < LL=%d skipped', MM, LL), return, end

if 1
	sls = ftab.sls;
	s1 = linspace(0, sls.max(1)+1, length(sls.sl{1}) + 8);
	s2 = linspace(0, sls.max(2)+1, length(sls.sl{2}) + 8);
	ss = ndgrid_jf('mat', s1, s2);
	ftmp = ftab.fit.fmfun(ss);
	stmp = ftab.inv_fun(ftmp);

	err = stmp - ss;
	err = abs(err);
	for ll=1:LL
		printm('worst inverse error l=%d: %g of %g', ll, ...
			max(col(stackpick(err,ll))), max(col(stackpick(ss,ll))))
	end
end

if im && usejava('jvm')
	ss1 = ss(:,:,1);
	ss2 = ss(:,:,2);
	clf, pl = @(n) subplot(220 + n);
	pl(1), mesh(s1, s2, ss1')
	axis tight, xtick, ytick, zwhite, grid
	hold on, plot3(ss1, ss2, stmp(:,:,1), 'y.'), hold off
	xlabel 's1', ylabel 's2', title 's1'

	pl(2), mesh(s1, s2, ss2')
	axis tight, xtick, ytick, zwhite, grid
	hold on, plot3(ss1, ss2, stmp(:,:,2), 'y.'), hold off
	xlabel 's1', ylabel 's2', title 's2'

	pl(3), mesh(s1, s2, err(:,:,1)')
	colormap hsv, caxis([0 max(err(:))]), cbar
	axis tight, xtick, ytick, zwhite, grid
	xlabel 's1', ylabel 's2', title 'error 1'

	pl(4), mesh(s1, s2, err(:,:,2)')
	axis tight, xtick, ytick, zwhite, grid
	xlabel 's1', ylabel 's2', title 'error 2'
prompt
end

end % de_ftab_plot_inv()

%
% de_ftab_plot_mac()
%
function dummy = de_ftab_plot_mac(ftab, varargin)
dummy = [];
xrs = ftab.xrs;
mas = ftab.mas;
mac = ftab.mac;

ltype = {'c:', 'y:', 'g:', 'r:'};
mtype = {'g+', 'r^', 'm*', 'yo'};
arg = {};
for ll=1:ftab.LL
	arg = {arg{:}, xrs.en, mac.mac(:,ll), ltype{ll}};
end
for mm=1:ftab.MM
	arg = {arg{:}, xrs.eff(mm), mac.bar(mm,:), mtype{mm}};
end

plot(arg{:})
axisy(0.1, max(0.7, 1.1*max(mac.bar(:))))
xlabel 'E', ylabel 'mac(E)'
mtype = mas.type;
legend(mtype{:})

prompt
end % de_ftab_plot_mac()

%
% de_ftab_make_fm()
% build tables of f_m(s1, s2, ...)
%
function fm = de_ftab_make_fm(ftab, varargin)
if length(varargin) == 0
	varargin = ftab.sls.sl; % default
end
if length(varargin) == 1 && iscell(varargin{1})
	sll = varargin{1};
else
	sll = ndgrid_jf('cell', varargin{:});
end
fm = de_ftab_fm(sll, ftab.mac.mac, ftab.xrs.Ide);

end % de_ftab_make_fm()

%
% de_ftab_plot_fm()
%
function dummy = de_ftab_plot_fm(ftab, varargin)
dummy = [];
arg.down = 1;
arg = vararg_pair(arg, varargin);

fm = ftab.fm;
sl = ftab.sls.sl;

if ftab.LL == 1
	plot(sl{1}, fm, '-', sl{1}, sl{1} * ftab.mac.bar(:)', ':')
	xlabel 's1'
	ylabel 'f1(s1)'
	name = ftab.xrs.name;
	legend(name{:}, 2)

elseif ftab.LL == 2

	fmax = 1.01 * max(fm(:));
	fmax = ceil(fmax);
	s_max(1) = max(sl{1}(:));
	s_max(2) = max(sl{2}(:));

	i1 = 1:arg.down:length(sl{1});
	i2 = 1:arg.down:length(sl{2});
	s1 = sl{1}(i1);
	s2 = sl{2}(i2);

	pl = @(mm) subplot(1,ftab.MM,mm);

	clf
	for mm=1:ftab.MM
		f1 = fm(i1,i2,mm);
		pl(mm)
		plot(s1, f1', '.-')
		axisx([0 s_max(1)]), xlabel 's_1'
		titlef('f_%d(s)', mm)
		text(-60, 10, '[cm^2/g]')
	end

	if usejava('jvm')
		prompt
	else
		warn 'need jvm for plot_fm mesh'
	return
	end

	clf
	for mm=1:ftab.MM
		f1 = fm(i1,i2,mm);
		pl(mm)
		mesh(s1, s2, f1')
		colormap hsv, caxis([0 fmax]), cbar
		axis([0 s_max(1) 0 s_max(2) 0 fmax])
		xtick, ytick, ztick, zwhite, xlabel s_1, ylabel s_2
		titlef('f_%d(s)', mm)
		text(-60, 10, '[cm^2/g]')
	end
end

end % de_ftab_plot_fm()

%
% de_ftab_test
%
function de_ftab_test
xrs = xray_read_spectra('ps1'); % M=2
%xrs = xray_read_spectra('poly1,60,90,120'); % M=3 test
%xrs = xray_read_spectra('poly1,160', 'filters', {{'copper', 0.1}}); % M < L test
xrs.plot
prompt
mas = xray_read_mac({'water', 'bone'});
if 0
	mas = xray_read_mac({'soft', 'iodine'});
	tmp = mas.mac_raw;
	tmp{2}(:) = tmp{2}(:) / 10; % todo: explore scaling to improve cond #
	mas.mac_raw = tmp;
end
ftab = de_ftab(xrs, mas, 'show', im, ...
	 'ftype','exp', 'fit_args', {'kev', 10:5:160, 'mac', []});
if im
	ftab.plot_fm;
	prompt
	ftab.inv1.plot(ftab.fit);
	if ~isempty(ftab.inv2)
		ftab.inv2.plot(ftab.fit, ftab.sls.sl);
	end
	clf, ftab.plot_mac;
	ftab.plot_jac;
	ftab.plot_inv;
end
end % de_ftab_test

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% Stability domain of Adams-Bashforth_Euler scheme

% Copyright 2001 P. Wesseling
% This program and its subprograms may be freely used, modified and distributed
% under the GNU General Public License: http://www.gnu.org/copyleft/gpl.html

% Theory is given in Section 5.8 of
% 	P. Wesseling: Principles of Computational Fluid Dynamics
% 	Springer, Heidelberg, 2000 ISBN 3-5453-0. XII, 642 pp.
% See http://ta.twi.tudelft.nl/nw/users/wesseling/cfdbook.html

% This program generates Fig. 5.5  in the book

figure(1), clf, hold on

		% Plot of stability domain
theta = 0:0.01:1; theta = theta*pi;
z = - 2*(1-cos(theta)).^2./(3-cos(theta)) + i*2*sin(theta)./(3-cos(theta));
plot(z)

		% Stability domain of second order Adams-Bashforth scheme
z = exp(2*i*theta) - exp(i*theta); z = z./(1.5*exp(i*theta) - 0.5); plot(z,'-.')

		% Plot of oval
b = 1; a = 1; z =  - a*(1-cos(theta)) + i*b*(sin(theta)).^0.5; plot(z,'--')

		% Plot of half ellipse
a = 2.0; b = sqrt(1/3); theta = -pi/2 +0.5*theta;
z = - a*cos(theta) - i*b*sin(theta); plot(z,'.')

repo: Aero-matlab


text:
function title = p31_title ( )

%*****************************************************************************80
%
%% P31_TITLE returns a title for problem 31.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = 'The Shekel SQRN5 Function';

  return
end

author: johannesgerer


text:
function test8 (nmat)
%TEST8 order a large range of sparse matrices, test symbfact2
% compare AMD and METIS
% Example:
%   test8(nmat)
% See also cholmod_test

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

fprintf ('=================================================================\n');
fprintf ('test8: factorize a large range of sparse matrices\n') ;

% get list of test matrices

index = UFget ;

% GHS posdef test set (more or less)
f = find (...
    ((index.numerical_symmetry == 1 & index.isBinary) | (index.posdef)) ...
    & (index.nnzdiag == index.nrows) ...
    & (index.nrows > 10000 | index.nrows == 9000) ...
    & (index.nrows < 600000) & (index.nnz > index.nrows)) ;		    %#ok

% include small matrices
f = find (...
    ((index.numerical_symmetry == 1 & index.isBinary) | (index.posdef)) ...
    & (index.nnzdiag == index.nrows) ...
    & (index.nrows < 600000) & (index.nnz > index.nrows)) ;

for k = 1:length (f)
    names {k} = index.Name {f(k)} ;			%#ok
end

[ignore i] = sort (names) ;

f = f (i) ;

% fprintf ('test matrices sorted by name:\n') ;
% for i = f
%     fprintf ('%4d: %-20s %-20s %12d %d\n', i,  ...
% 	index.Group {i}, index.Name {i}, index.nrows (i), index.posdef (i)) ;
% end

[ignore i] = sort (index.nrows (f)) ;
f = f (i) ;

if (nargin > 0)
    nmat = max (0,nmat) ;
    nmat = min (nmat, length (f)) ;
    f = f (1:nmat) ;
end

fprintf ('test matrices sorted by dimension:\n') ;
for i = f
    fprintf ('%4d: %-20s %-20s %12d %d\n', i,  ...
	index.Group {i}, index.Name {i}, index.nrows (i), index.posdef (i)) ;
end

junk = sparse (1) ;

% input ('hit enter to continue: ') ;

for k = 1:length (f)

    Problem = UFget (f(k)) ;
    A = Problem.A ;
    fprintf ('\n================== Problem: %s  n: %d nnz: %d\n', ...
	Problem.name, size (A,1), nnz (A)) ;
    fprintf ('title: %s\n\n', Problem.title) ;
    clear Problem
    n = size (A,1) ;							    %#ok

    amd2 (junk) ;
    metis (junk) ;

    tic ;
    [p1,info] = amd2 (A) ;						    %#ok
    t1 = toc ;
    S1 = A (p1,p1) ;
    tic ;
    c1 = symbfact (S1) ;
    ts1 = toc ;
    tic ;
    d1 = symbfact (S1) ;
    ts2 = toc ;
    if (any (c1 ~= d1))
	error ('!')
    end
    fprintf ('symbfact time: MATLAB %9.4f  CHOLMOD %9.4f  speedup %8.2f\n', ...
	ts1, ts2, ts1/ts2) ;

    lnz1 = sum (c1) ;
    fl1 = sum (c1.^2) ;
    fprintf ('time: amd     %10.4f mnnz(L) %8.1f mfl %8.0f  fl/nnz(L) %8.1f\n', ...
	t1, lnz1/1e6, fl1 /1e6, fl1/lnz1) ;

    tic ;
    p2 = metis (A) ;
    t2 = toc ;
    S2 = A (p2,p2) ;
    c2 = symbfact (S2) ;
    lnz2 = sum (c2) ;
    fl2 = sum (c2.^2) ;

    fprintf ('time: metis   %10.4f mnnz(L) %8.1f mfl %8.0f  fl/nnz(L) %8.1f\n', ...
	t2, lnz2/1e6, fl2/1e6, fl2/lnz2) ;

    r = lnz2 / lnz1 ;							    %#ok
    fprintf ('\nmetis/amd time: %8.4f nnz(L): %8.4f\n', t2/t1, lnz2/lnz1) ;

    % save results
    lnz (k,1) = lnz1 ;			    %#ok
    lnz (k,2) = lnz2 ;			    %#ok
    fl1 (k,1) = fl1 ;			    %#ok
    fl2 (k,2) = fl2 ;			    %#ok
    t (k,1) = t1 ;			    %#ok
    t (k,2) = t2 ;			    %#ok

end

fprintf ('test8 passed\n') ;

author: fieldtrip

repo: fieldtrip