text:
function toms655_test ( )

%*****************************************************************************80
%
%% TOMS655_TEST tests the TOMS655 library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TOMS655_PRB\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the TOMS655 library.\n' );

  test01 ( );
  test02 ( );
  test03 ( );
  test04 ( );
  test05 ( );
  test06 ( );
  test07 ( );
  test08 ( );
  test09 ( );
%
%  Compute 15 points of an example of each rule.
%
  for kind = 1 : 9
    nt = 15;
    if ( kind == 8 )
      alpha = 1.0;
      beta = - alpha - 2 * nt - 2;
    else
      alpha = 0.0;
      beta = 0.0;
    end
    test10 ( nt, kind, alpha, beta );
  end
%
%  Compute 15 points of an example of each rule using nondefault A, B.
%
  for kind = 1 : 9

    nt = 15;

    if ( kind == 1 )
      alpha = 0.0;
      beta = 0.0;
      a = 0.0;
      b = 1.0;
    elseif ( kind == 2 )
      alpha = 0.0;
      beta = 0.0;
      a = 0.0;
      b = 1.0;
    elseif ( kind == 3 )
      alpha = 1.0;
      beta = 0.0;
      a = 0.0;
      b = 1.0;
    elseif ( kind == 4 )
      alpha = 1.5;
      beta = 0.5;
      a = 0.0;
      b = 1.0;
    elseif ( kind == 5 )
      alpha = 1.0;
      beta = 0.0;
      a = 1.0;
      b = 1.0;
    elseif ( kind == 6 )
      alpha = 1.0;
      beta = 0.0;
      a = 0.0;
      b = 0.5;
    elseif ( kind == 7 )
      alpha = 1.0;
      beta = 0.0;
      a = 0.0;
      b = 1.0;
    elseif ( kind == 8 )
      alpha = 1.0;
      beta = - alpha - 2 * nt - 2;
      a = 0.0;
      b = 1.0;
    elseif ( kind == 9 )
      alpha = 0.0;
      beta = 0.0;
      a = 0.0;
      b = 1.0;
    end

    test11 ( nt, kind, alpha, beta, a, b );

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TOMS655_PRB\n' );
  fprintf ( 1, '  Normal end of TOMS655 tests.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function test01 ( )

%*****************************************************************************80
%
%% TEST01 tests CIQFS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
  fprintf ( 1, '  ----------------------------------------\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  Test CIQFS.\n' );
%
%  Number of knots.
%
  nt = 5;
%
%  Set the knots in the default interval [-1,+1].
%
  t = zeros ( nt, 1 );
  for i = 1 : nt
    t(i) = cos ( ( 2 * i - 1 ) * pi / 2.0 / nt );
  end
%
%  Set the knot multiplicities.
%
  mlt(1:nt) = 2;
%
%  Set the size of the weights array.
%
  nwts = sum ( mlt(1:nt) );
%
%  Because KEY = 1, NDX will be set up for us.
%
  ndx = zeros(nt,1);
%
%  KEY = 1 indicates that the WTS array should hold the weights
%  in the usual order.
%
  key = 1;
%
%  Request Legendre weight function.
%
  kind = 1;
%
%  ALPHA, BETA not used in Legendre weight function but set anyway.
%
  alpha = 0.0;
  beta  = 0.0;
%
%  LU controls printing.
%  A positive value requests that we compute and print weights, and
%  conduct a moments check.
%
  lu = 6;

  [ wts, ndx ] = ciqfs ( nt, t, mlt, nwts, ndx, key, kind, alpha, beta, lu );

  return
end
function test02 ( )

%*****************************************************************************80
%
%% TEST02 tests CIQFS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
  fprintf ( 1, '  ----------------------------------------\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  Test CIQF, CIQFS, CGQF and CGQFS\n' );
  fprintf ( 1, '  with all classical weight functions.\n' );
%
%  Try all weight functions.
%
  for kind = 1 : 9
%
%  Number of knots.
%
    nt = 5;
%
%  Set parameters ALPHA and BETA.
%
    alpha = 0.5;
    if ( kind ~= 8 )
      beta  = 2.0;
    else
      beta = - 16.0;
    end
%
%  Set A and B.
%
    lo = 6;
    a = - 0.5;
    b = 2.0;
%
%  Have CGQF compute the knots and weights.
%
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Knots and weights of Gauss quadrature formula\n' );
    fprintf ( 1, '  computed by CGQF.\n' );
    [ t, wts ] = cgqf ( nt, kind, alpha, beta, lo, a, b );
%
%  Now compute the weights for the same knots by CIQF.
%
%  Set the knot multiplicities.
%
    mlt = zeros(nt,1);
    mlt(1:nt) = 2;
%
%  Set the size of the weights array.
%
    nwts = sum ( mlt(1:nt) );
%
%  Because KEY = 1, NDX will be set up for us.
%
    ndx = zeros(nt,1);
%
%  KEY = 1 indicates that the WTS array should hold the weights
%  in the usual order.
%
    key = 1;
%
%  LU controls printing.
%  A positive value requests that we compute and print weights, and
%  conduct a moments check.
%
    lu = 6;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Weights of Gauss quadrature formula computed from the\n' );
    fprintf ( 1, '  knots by CIQF.\n' );

    wts = ciqf ( nt, t, mlt, nwts, ndx, key, kind, alpha, beta, a, b, lu );

  end

  return
end
function test03 ( )

%*****************************************************************************80
%
%% TEST03 tests CEIQFS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
  fprintf ( 1, '  ----------------------------------------\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST03\n' );
  fprintf ( 1, '  Test CEIQFS.\n' );
%
%  Number of knots.
%
  nt = 5;
%
%  Set the knots in the default interval [-1,+1].
%
  t = zeros ( nt, 1 );

  for i = 1 : nt
    t(i) = cos ( ( 2 * i - 1 ) * pi / 2.0 / nt );
  end
%
%  Set the knot multiplicities.
%
  mlt = zeros ( nt, 1 );
  mlt(1:nt) = 2;
%
%  Set KIND to the Legendre weight function.
%
  kind = 1;
%
%  ALPHA, BETA not used in Legendre weight function but set anyway.
%
  alpha = 0.0;
  beta  = 0.0;
%
%  Call CEIQFS to set up the quadrature formula and evaluate it on F.
%
  qfsum = ceiqfs ( nt, t, mlt, kind, alpha, beta, @f );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Integral of sin(x) on -1, 1 by Fejer type rule\n' );
  fprintf ( 1, '  with %d points of multiplicity 2.\n', nt );
  fprintf ( 1, '  Quadrature formula: %24.16f\n', qfsum );

  qfsx = cos ( - 1.0 ) - cos ( 1.0 );
  fprintf ( 1, '  Exact value       : %24.16f\n', qfsx );
  fprintf ( 1, '  Error             : %e\n', abs ( qfsum - qfsx ) );

  return
end
function test04 ( )

%*****************************************************************************80
%
%% TEST04 tests CEIQF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
  fprintf ( 1, '  ----------------------------------------\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST04\n' );
  fprintf ( 1, '  Test CEIQF.\n' );
%
%  Number of knots.
%
  nt = 5;
%
%  Set the knots in the default interval [-1,+1].
%
  t = zeros ( nt, 1 );

  for i = 1 : nt
    t(i) = cos ( ( 2 * i - 1 ) * pi / 2.0 / nt );
  end
%
%  Set the knot multiplicities.
%
  mlt = zeros ( nt, 1 );
  mlt(1:nt) = 2;
%
%  Set KIND to the Legendre weight function.
%
  kind = 1;
%
%  ALPHA, BETA not used in Legendre weight function but set anyway.
%
  alpha = 0.0;
  beta  = 0.0;
%
%  Set nonstandard interval A, B.
%
  a = - 0.5;
  b = 2.0;
%
%  Shift knots from [-1,1] to [A,B].
%
  for i = 1 : nt
    t(i) = ( ( b - a ) * t(i) + ( a + b ) ) / 2.0;
  end
%
%  Call CEIQF to set up the quadrature formula and evaluate it on F.
%
  qfsum = ceiqf ( nt, t, mlt, kind, alpha, beta, a, b, @f );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Integral of sin(x) from %f to %f by Fejer type rule\n', a, b );
  fprintf ( 1, '  with %d points of multiplicity 2.\n', nt );
  fprintf ( 1, '  Quadrature formula: %24.16f\n', qfsum );

  qfsx = cos ( a ) - cos ( b );
  fprintf ( 1, '  Exact value       : %24.16f\n', qfsx );
  fprintf ( 1, '  Error             : %e\n', abs ( qfsum - qfsx ) );

  return
end
function test05 ( )

%*****************************************************************************80
%
%% TEST05 tests CLIQFS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
  fprintf ( 1, '  ----------------------------------------\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST05\n' );
  fprintf ( 1, '  Test CLIQFS.\n' );
%
%  Number of knots.
%
  nt = 5;
%
%  Set the knots in the default interval [-1,+1].
%
  t = zeros(nt,1);

  for i = 1 : nt
    t(i) = cos ( ( 2 * i - 1 ) * pi / ( 2 * nt ) );
  end
%
%  Request Legendre weight function.
%
  kind = 1;
%
%  ALPHA, BETA not used in Legendre weight function but set anyway.
%
  alpha = 0.0;
  beta  = 0.0;
%
%  LU controls printing.
%  A positive value requests that we compute and print weights, and
%  conduct a moments check.
%
  lu = 6;
%
%  This call returns the WTS array.
%
  wts = cliqfs ( nt, t, kind, alpha, beta, lu );

  return
end
function test06 ( )

%*****************************************************************************80
%
%% TEST06 tests CLIQF and EIQFS..
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
  fprintf ( 1, '  ----------------------------------------\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST06\n' );
  fprintf ( 1, '  Test CLIQF and EIQFS.\n' );
%
%  Number of knots.
%
  nt = 5;
%
%  Set the knots in the default interval [-1,+1].
%
  t = zeros(nt,1);

  for i = 1 : nt
    t(i) = cos ( ( 2 * i - 1 ) * pi / ( 2 * nt ) );
  end
%
%  Set KIND to the Legendre weight function.
%
  kind = 1;
%
%  ALPHA, BETA not used in Legendre weight function but set anyway.
%
  alpha = 0.0;
  beta  = 0.0;
%
%  Set nonstandard interval A, B.
%
  a = - 0.5;
  b = 2.0;
%
%  Shift knots from [-1,1] to [A,B].
%
  for i = 1 : nt
    t(i) = ( ( b - a ) * t(i) + ( a + b ) ) / 2.0;
  end
%
%  LU controls printout.
%
  lu = 6;
%
%  Call CLIQF to set up the quadrature formula.
%
  wts = cliqf ( nt, t, kind, alpha, beta, a, b, lu );
%
%  Call EIQFS to evaluate the quadrature formula.
%
  qfsum = eiqfs ( nt, t, wts, @f );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Integral of sin(x) from %f to %f\n', a, b );
  fprintf ( 1, '  by Fejer type rule with %d points\n', nt );
  fprintf ( 1, '  of multiplicity 1.\n' );
  fprintf ( 1, '  Quadrature formula: %24.16f\n', qfsum );

  qfsx = cos ( a ) - cos ( b );
  fprintf ( 1, '  Exact value       : %24.16f\n', qfsx );
  fprintf ( 1, '  Error             : %e\n', abs ( qfsum - qfsx ) );

  return
end
function test07 ( )

%*****************************************************************************80
%
%% TEST07 tests CEGQF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
  fprintf ( 1, '  ----------------------------------------\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST07\n' );
  fprintf ( 1, '  Test CEGQF.\n' );
%
%  Number of knots.
%
  nt = 12;
%
%  Request exponential weight function.
%
  kind = 7;
%
%  Set ALPHA and BETA.
%
  alpha = 1.0;
  beta  = 0.0;
%
%  Set interval [A,B].
%
  a = - 0.5;
  b = 2.0;
%
%  Call CEGQF to compute and evaluate the Gauss quadrature formula.
%
  qfsum = cegqf ( nt, kind, alpha, beta, a, b, @f );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Integral of x*sin(x) from %f to %f\n', a, b );
  fprintf ( 1, '  by Gauss-exponential rule with %d points\n', nt );
  fprintf ( 1, '  Quadrature formula: %24.16f\n', qfsum );

  qfsx = ( b - a ) * 0.5 * ( cos ( a ) - cos ( b ) ) ...
    + sin ( b ) + sin ( a ) - 2.0 * sin ( ( a + b ) / 2.0 );

  fprintf ( 1, '  Exact value       : %24.16f\n', qfsx );
  fprintf ( 1, '  Error             : %e\n', abs ( qfsum - qfsx ) );

  return
end
function test08 ( )

%*****************************************************************************80
%
%% TEST08 tests CEGQFS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
  fprintf ( 1,'  ----------------------------------------\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST08\n' );
  fprintf ( 1, '  Test CEGQFS.\n' );
%
%  Number of knots.
%
  nt = 12;
%
%  Request exponential weight function.
%
  kind = 7;
%
%  Set ALPHA and BETA.
%
  alpha = 1.0;
  beta  = 0.0;
%
%  Call CEGQFS to compute and evaluate the Gauss quadrature formula.
%
  qfsum = cegqfs ( nt, kind, alpha, beta, @f );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Integral of x*sin(x) from -1 to +1\n' );
  fprintf ( 1, '  by Gauss-exponential rule with %d points.\n', nt )
  fprintf ( 1, '  Quadrature formula: %24.16f\n', qfsum );

  qfsx = cos ( -1.0 ) - cos ( +1.0 );

  fprintf ( 1, '  Exact value       : %24.16f\n', qfsx );
  fprintf ( 1, '  Error             : %e\n', abs ( qfsum - qfsx ) );

  return
end
function test09 ( )

%*****************************************************************************80
%
%% TEST09 calls CGQFS to compute and print generalized Gauss-Hermite rules.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 January 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST09\n' );
  fprintf ( 1, '  Call CGQFS for a generalized Gauss Hermite rule.\n' );

  nt = 15;
  kind = 6;
  alpha = 1.0;
  beta = 0.0;
  io = - 6;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  NT = %d\n', nt );
  fprintf ( 1, '  ALPHA = %f\n', alpha );

  [ t, wts ] = cgqfs ( nt, kind, alpha, beta, io );

  return
end
function test10 ( nt, kind, alpha, beta )

%*****************************************************************************80
%
%% TEST10 calls CDGQF to compute a quadrature formula.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 January 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST10\n' );
  fprintf ( 1, '  Call CDGQF to compute a quadrature formula.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  KIND = %d\n', kind );
  fprintf ( 1, '  ALPHA = %f\n', alpha );
  fprintf ( 1, '  BETA  = %f\n', beta );

  [ t, wts ] = cdgqf ( nt, kind, alpha, beta );

  fprintf ( 1, '\n' );
  fprintf ( 1, ' Index     Abscissas                 Weights\n' );
  fprintf ( 1, '\n' );
  for i = 1 : nt
    fprintf ( 1, '  %4d  %24.16e  %24.16e\n', i, t(i), wts(i) );
  end

  return
end
function test11 ( nt, kind, alpha, beta, a, b )

%*****************************************************************************80
%
%% TEST11 calls CGQF to compute a quadrature formula.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST11\n' );
  fprintf ( 1, '  Call CGQF to compute a quadrature formula\n' );
  fprintf ( 1, '  with nondefault values of A and B.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  KIND = %d\n', kind );
  fprintf ( 1, '  ALPHA = %f\n', alpha );
  fprintf ( 1, '  BETA  = %f\n', beta );
  fprintf ( 1, '  A     = %f\n', a );
  fprintf ( 1, '  B     = %f\n', b );

  lo = 0;
  [ t, wts ] = cgqf ( nt, kind, alpha, beta, lo, a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, ' Index     Abscissas                 Weights\n' );
  fprintf ( 1, '\n' );
  for i = 1 : nt
    fprintf ( 1, '  %4d  %24.16e  %24.16e\n', i, t(i), wts(i) );
  end

  return
end
function value = f ( x, i )

%*****************************************************************************80
%
%% F returns values of the integrand or its derivatives.
%
%  Discussion:
%
%    This function is an example of an integrand function.
%
%    The package can generate quadrature formulas that use derivative
%    information as well as function values.  Therefore, this routine is
%    set up to provide derivatives of any order as well as the function
%    value.  In an actual application, the highest derivative needed
%    is of order one less than the highest knot multiplicity.
%
%    In other words, in the usual case where knots are not repeated,
%    this routine only needs to return function values, not any derivatives.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Input, integer I, the order of the derivative of F to
%    be evaluated.
%
%    Output, real VALUE, the value of the I-th derivative of F at X.
%
  l = mod ( i, 4 );

  if ( l == 0 )
    value = sin ( x );
  elseif ( l == 1 )
    value = cos ( x );
  elseif ( l == 2 )
    value = - sin ( x );
  elseif ( l == 3 )
    value = - cos ( x );
  end

  return
end


text:
function [varargout]=subEdge(varargin)

%%

switch nargin
    case 2
        F=varargin{1};
        V=varargin{2};
        n=1;
        uniqueOpt=1;
    case 3
        F=varargin{1};
        V=varargin{2};
        n=varargin{3};
        uniqueOpt=1;
    case 4
        F=varargin{1};
        V=varargin{2};
        n=varargin{3};
        uniqueOpt=varargin{4};
end

%%

%Check if n can be achieved through splitting
nSplitIterations=log2(n+1); %Check for integer solution
logicInteger=abs(round(nSplitIterations)-nSplitIterations)<eps(nSplitIterations);

if uniqueOpt==1 && logicInteger
    subMethod='split';
else
    subMethod='seed';
end

switch subMethod
    case 'split' %iteratively split edges (no double points created, unique operation avoided)
        Fs=F; Vs=V;
        for q=1:1:nSplitIterations
            F=Fs; V=Vs;

            Fs=zeros(size(Fs,1),size(Fs,2)*2);
            Fs(:,1:2:end)=F;

            [E]=patchEdges(F,1);

            numPoints = size(V,1);
            numEdges = size(E,1);

            % Get indices of the edges associated with each face
            A = sparse(E(:,1),E(:,2),(1:numEdges)+numPoints,numPoints,numPoints,numEdges);
            A = max(A,A'); %Copy symmetric

            indFs=2:2:size(Fs,2);
            for qp=1:1:size(F,2)
                if qp<size(F,2)
                    indA_now=F(:,qp)+(F(:,qp+1)-1)*numPoints;
                else
                    indA_now=F(:,qp)+(F(:,1)-1)*numPoints;
                end

                indV_now=full(A(indA_now));
                Fs(:,indFs(qp))=indV_now;
            end

            %Create vertex array
            Vn=0.5*(V(E(:,1),:)+V(E(:,2),:)); %new mid-edge points
            Vs = [V; Vn]; %Join point sets
        end
    case 'seed' %Seed edge points and remove doubles (more memory intensive)
        %Edges matrix
        [E]=patchEdges(F,0);

        Vs=zeros(size(E,1)*(n+1),size(V,2));
        for q=1:1:size(V,2)
            X=V(:,q);
            XE=X(E);
            X_add=linspacen(XE(:,1),XE(:,2),n+2);
            X_add=X_add(:,1:end-1)';
            Vs(:,q)=X_add(:);
        end

        ind=(1:1:size(Vs,1))';
        Fs=reshape(ind,[(n+1)*size(F,2),size(Vs,1)./((n+1)*size(F,2))])';

        %Merge non-unique nodes
        [Fs,Vs]=mergeVertices(Fs,Vs);

end

%% Collect output
varargout{1}=Fs;
varargout{2}=Vs;

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 [Forecasteval]=TVEmafeval(data_endo_a,data_endo_c,data_endo_c_lags,data_exo_c,data_exo_c_lags,stringdates3,Fstartdate,Fcenddate,Fcperiods,Fcomp,const,n,m,p,k1,k3,It,Bu,beta_gibbs,sigma_gibbs,forecast_record,forecast_estimates,names,endo,pref,theta_gibbs,TVEHfuture,ss_record,indH)

Forecasteval.RMSE=[];
Forecasteval.MAE=[];
Forecasteval.MAPE=[];
Forecasteval.Ustat=[];
Forecasteval.CRPS_estimates=[];
Forecasteval.S1_estimates=[];
Forecasteval.S2_estimates=[];

% first, note that forecast evaluation can only be conducted if there is some observable data after the beginning of the forecast
if Fcomp==1

   % preliminary task: obtain a matrix of forecasts over the common periods
   for ii=1:n
   forecast_c(:,ii)=forecast_estimates{ii,1}(2,1:Fcperiods)';
   end
   % then compute the matrix of forecast errors
   ferrors=data_endo_c-forecast_c;

   % compute first the sequential RMSE, defined in (a.8.11)

   % square the forecast error matrix entrywise
   sferrors=ferrors.^2;
   % sum entries sequentially
   sumsferrors=sferrors(1,:);
   for ii=2:Fcperiods
   sumsferrors(ii,:)=sumsferrors(ii-1,:)+sferrors(ii,:);
   end
   % divide by the number of forecast periods and take square roots to obtain RMSE
   for ii=1:Fcperiods
   Forecasteval.RMSE(ii,:)=((1/ii)*sumsferrors(ii,:)).^0.5;
   end

   % compute then the sequential MAE, defined in (a.8.12)

   % take the absolute value of the forecast error matrix
   absferrors=abs(ferrors);
   % sum entries sequentially
   sumabsferrors=absferrors(1,:);
   for ii=2:Fcperiods
   sumabsferrors(ii,:)=sumabsferrors(ii-1,:)+absferrors(ii,:);
   end
   % divide by the number of forecast periods to obtain MAE
   for ii=1:Fcperiods
   Forecasteval.MAE(ii,:)=(1/ii)*sumabsferrors(ii,:);
   end

   % compute the sequential MAPE, defined in (a.8.13)

   % divide entrywise by actual values and take absolute values
   absratioferrors=abs(ferrors./data_endo_c);
   % sum entries sequentially
   sumabsratioferrors=absratioferrors(1,:);
   for ii=2:Fcperiods
   sumabsratioferrors(ii,:)=sumabsratioferrors(ii-1,:)+absratioferrors(ii,:);
   end
   % divide by 100*(number of forecast periods) to obtain MAPE
   for ii=1:Fcperiods
   Forecasteval.MAPE(ii,:)=(100/ii)*sumabsratioferrors(ii,:);
   end

   % compute the Theil's inequality coefficient, defined in (a.8.14)

   % first compute the left term of the denominator
   % square entrywise the matrix of actual data
   sendo=data_endo_c.^2;
   % sum entries sequentially
   sumsendo=sendo(1,:);
   for ii=2:Fcperiods
   sumsendo(ii,:)=sumsendo(ii-1,:)+sendo(ii,:);
   end
   % divide by the number of forecast periods and take square roots
   for ii=1:Fcperiods
   leftterm(ii,:)=((1/ii)*sumsendo(ii,:)).^0.5;
   end
   % then compute the right term of the denominator
   % square entrywise the matrix of forecast values
   sforecasts=forecast_c.^2;
   % sum entries sequentially
   sumsforecasts=sforecasts(1,:);
   for ii=2:Fcperiods
   sumsforecasts(ii,:)=sumsforecasts(ii-1,:)+sforecasts(ii,:);
   end
   % divide by the number of forecast periods and take square roots
   for ii=1:Fcperiods
   rightterm(ii,:)=((1/ii)*sumsforecasts(ii,:)).^0.5;
   end
   % finally, compute the U stats
   Forecasteval.Ustat=Forecasteval.RMSE./(leftterm+rightterm);

   % now compute the continuous ranked probability score

   % create the cell storing the results
   CRPS=cell(n,1);
   % loop over endogenous variables
   for ii=1:n
      % loop over forecast periods on which actual data is known
      for jj=1:Fcperiods
      % compute the continuous ranked probability score
      score=bear.crps(forecast_record{ii,1}(:,jj),forecast_estimates{ii,1}(2,jj));
      CRPS{ii,1}(1,jj)=score;
      end
   end
   Forecasteval.CRPS_estimates=(cell2mat(CRPS))';

   % finally, consider the big piece: the computation of the log predictive score
   % this part implements algorithm a.8.1, described p 115 of technical guide
   % details of the procedure can be found p 111-115 of the technical guide

   % preliminary task: if there is a constant in the model, concatenate a column of ones to data_exo_c
   if const==1
   data_exo_c=[ones(Fcperiods,1) data_exo_c];
   data_exo_c_lags=[ones(p,1) data_exo_c_lags];
   end

   % create the cells storing the results
   S1=cell(n,1);
   S2=cell(n,1);

   for ll=1:It-Bu
      % step 2: draw beta, Delta and sigma
      beta=beta_gibbs(:,ll);
      sigma=reshape(sigma_gibbs(:,ll),n,n);
      % reshape beta to obtain B
      B=reshape(beta,k1,n);
      theta=theta_gibbs(:,ii);
      tempeq = [];
      for nnn = 1:n
        tempeq = [tempeq ss_record{nnn,1}(ii,end-p+1:end)'];
      end

      % step 3: obtain mu, the mean vector
      temp=bear.TVEmaforecastsim(data_endo_a,B,p,n,Fcperiods,theta,TVEHfuture,tempeq, indH);
      mu=reshape(temp',n*Fcperiods,1);

      % step 4: obtain upsilon, the covariance matrix
      % recover the A1,A2,... matrices by reshaping beta
      temp=B';
      Amatrices=cell(1,p);
      for ii=1:p
      Amatrices{1,ii}=temp(:,n*(ii-1)+1:n*ii);
      end
      % initiate upsilon with upsilon_1,1
      upsilon=cell(Fcperiods,Fcperiods);
      upsilon{1,1}=sigma;
      % complete the first column  using (a.8.23)
      for ii=2:Fcperiods
      upsilon_ij=Amatrices{1,1}*upsilon{ii-1,1};
      summ=min(ii-1,p);
         for jj=2:summ
         upsilon_ij=upsilon_ij+Amatrices{1,jj}*upsilon{ii-jj,1};
         end
      upsilon{ii,1}=upsilon_ij;
      end
      % complete the first row, using symmetry
      for ii=2:Fcperiods
      upsilon{1,ii}=upsilon{ii,1}';
      end
      % now take care of the other columns
      for jj=2:Fcperiods
      % first complete the variance matrix (the diagonal one), using (a.8.24)
      upsilon_jj=sigma;
      summ=min(jj-1,p);
         for kk=1:summ
         upsilon_jj=upsilon_jj+Amatrices{1,kk}*upsilon{jj-kk,jj};
         end
      upsilon{jj,jj}=upsilon_jj;
      % then complete the rest of the column, using (a.8.23)
         for ii=jj+1:Fcperiods
         % initiate upsilon_ij
         upsilon_ij=Amatrices{1,1}*upsilon{ii-1,jj};
         % continue the summation
         summ=min(ii-1,p);
            for kk=2:summ
            upsilon_ij=upsilon_ij+Amatrices{1,kk}*upsilon{ii-kk,jj};
            end
         upsilon{ii,jj}=upsilon_ij;
         end
         % complete the corresponding row, using symmetry
         for kk=jj+1:Fcperiods
         upsilon{jj,kk}=upsilon{kk,jj}';
         end
      end
      upsilon=cell2mat(upsilon);

      % step 5: obtain the predictive density

      % first consider scenario 1 (forecast at period T+i)
      % loop over variables
      for ii=1:n
         % loop over periods
         for jj=1:Fcperiods
         % define R
         R=zeros(1,n*Fcperiods);
         R(1,n*(jj-1)+ii)=1;
         % define the mean vector
         mean=R*mu;
         % define the covariance matrix
         covar=R*upsilon*(R');
         % define the vector of actual values
         actual=data_endo_c(jj,ii);
         % determine the density
         [~,density]=bear.mndensity(actual,mean,covar,1);
         % record the result
         S1{ii,1}(ll,jj)=density;
         end
      end
      % then consider scenario 2 (forecast from beginning until period T+i)
      % loop over variables
      for ii=1:n
      R=[];
         % loop over periods
         for jj=1:Fcperiods
         % define R
         R=[R;zeros(1,n*Fcperiods)];
         R(jj,n*(jj-1)+ii)=1;
         % define the mean vector
         mean=R*mu;
         % define the covariance matrix
         covar=R*upsilon*(R');
         % define the vector of actual values
         actual=data_endo_c(1:jj,ii);
         % determine the density
         [~,density]=bear.mndensity(actual,mean,covar,jj);
         % record the result
         S2{ii,1}(ll,jj)=density;
         end
      end
   end

   % step 6: compute the log predictive score from (XXX.11)
   Forecasteval.S1_estimates=cell(n,1);
   Forecasteval.S2_estimates=cell(n,1);
   for ii=1:n
      for jj=1:Fcperiods
      Forecasteval.S1_estimates{ii,1}(1,jj)=log((1/(It-Bu))*sum(S1{ii,1}(:,jj)));
      Forecasteval.S2_estimates{ii,1}(1,jj)=log((1/(It-Bu))*sum(S2{ii,1}(:,jj)));
      end
   end
   Forecasteval.S1_estimates=(cell2mat(Forecasteval.S1_estimates))';
   Forecasteval.S2_estimates=(cell2mat(Forecasteval.S2_estimates))';

% if forecast evaluation is not possible, do not do anything
elseif Fcomp==0
end

% now, print the results and display them

filelocation=fullfile(pref.results_path, [pref.results_sub '.txt']);
fid=fopen(filelocation,'at');

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

Fevalinfo='Forecast evaluation:';
fprintf('%s\n',Fevalinfo);
fprintf(fid,'%s\n',Fevalinfo);

fprintf('%s\n','');
fprintf(fid,'%s\n','');

% if forecast evaluation is not possible, return a message to signal it

if Fcomp==0

finfo1=['Forecast evaluation cannot be conducted.'];
fprintf('%s\n',finfo1);
fprintf(fid,'%s\n',finfo1);
finfo2=['Forecasts start in ' Fstartdate ', while observable data is available only until ' names{end,1} '.'];
fprintf('%s\n',finfo2);
fprintf(fid,'%s\n',finfo2);
finfo3=['To obtain forecast evaluation, the forecast start date must be anterior to the end of the data set.'];
fprintf('%s\n',finfo3);
fprintf(fid,'%s\n',finfo3);

% if forecast evaluation is possible, display the results
elseif Fcomp==1

finfo1=['Evaluation conducted over ' num2str(Fcperiods) ' periods (from ' Fstartdate ' to ' Fcenddate ').'];
fprintf('%s\n',finfo1);
fprintf(fid,'%s\n',finfo1);

   % loop over endogenous variables
   for ii=1:n

   fprintf('%s\n','');
   fprintf(fid,'%s\n','');

   endoinfo=['Endogenous: ' endo{ii,1}];
   fprintf('%s\n',endoinfo);
   fprintf(fid,'%s\n',endoinfo);

   temp='fprintf(''%12s';
      for jj=1:Fcperiods-1
      temp=[temp ' %10s'];
      end
   temp=[temp ' %10s\n'','''''];
      for jj=1:Fcperiods
      temp=[temp ',''' stringdates3{jj,1} ''''];
      end
   temp=[temp ');'];
   eval(temp);
   temp='fprintf(fid,''%12s';
      for jj=1:Fcperiods-1
      temp=[temp ' %10s'];
      end
   temp=[temp ' %10s\n'','''''];
      for jj=1:Fcperiods
      temp=[temp ',''' stringdates3{jj,1} ''''];
      end
   temp=[temp ');'];
   eval(temp);

   label='RMSE:       ';
   values=Forecasteval.RMSE(1:Fcperiods,ii)';
   temp='fprintf(''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);
   temp='fprintf(fid,''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);

   label='MAE:        ';
   values=Forecasteval.MAE(1:Fcperiods,ii)';
   temp='fprintf(''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);
   temp='fprintf(fid,''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);

   label='MAPE:       ';
   values=Forecasteval.MAPE(1:Fcperiods,ii)';
   temp='fprintf(''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);
   temp='fprintf(fid,''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);

   label='Theil''s U:  ';
   values=Forecasteval.Ustat(1:Fcperiods,ii)';
   temp='fprintf(''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);
   temp='fprintf(fid,''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);

   label='CRPS:       ';
   values=Forecasteval.CRPS_estimates(1:Fcperiods,ii)';
   temp='fprintf(''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);
   temp='fprintf(fid,''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);

   label='Log score 1:';
   values=Forecasteval.S1_estimates(1:Fcperiods,ii)';
   temp='fprintf(''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);
   temp='fprintf(fid,''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);

   label='Log score 2:';
   values=Forecasteval.S2_estimates(1:Fcperiods,ii)';
   temp='fprintf(''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);
   temp='fprintf(fid,''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);

    end

end

fclose(fid);

author: european-central-bank

repo: BEAR-toolbox


text:
function p00_start_test ( problem_num )

%*****************************************************************************80
%
%% P00_START_TEST prints the norm of the starting point and its function value.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM_NUM, the number of problems.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'P00_START_TEST\n' );
  fprintf ( 1, '  Get norms of starting point X0 and F(X0)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   Problem    Option        ||X0||          ||F(X0)||\n' );
  fprintf ( 1, '\n' );

  for problem = 1 : problem_num

    option_num = p00_option_num ( problem );

    fprintf ( 1, '\n' );

    for option = 1 : option_num

      nvar = p00_nvar ( problem, option );

      x0 = p00_start ( problem, option, nvar );

      x0_norm = norm ( x0(1:nvar) );

      fx0 = p00_fun ( problem, option, nvar, x0 );

      fx0_norm = norm ( fx0(1:nvar-1) );

      fprintf ( 1, '  %8d  %8d  %14f  %14f\n', ...
        problem, option, x0_norm, fx0_norm )

    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% [INPUT]
% dep = A vector of floats of length t representing the dependent variable.
% indep_s = A float t-by-cs matrix (-Inf,Inf) representing the switching independent variables, without intercept because it is internally handled by the model (optional, default=[]).
% indep_ns = A float t-by-cns matrix (-Inf,Inf) representing the non-switching independent variables (optional, default=[]).
% k = An integer [2,4] representing the number of states of the model (optional, default=2).
% vs = A boolean that indicates whether the variance is allowed to switch (optional, default=true).
% finit = A function handle representing a hook for generating a customized minimization problem (optional, by default a standard minimization problem is generated).
%   The function handle must accept the following input arguments (in the same order):
%     - x0: A column vector of floats of length cx representing the initial parameters.
%     - dep, indep_s, indep_s, k, vs, tmm: See the description above.
%     - p0 = A float k-by-k matrix [0,1] representing the initial stochastic row-wise transition matrix.
%     - options = A structure representing the optimization options structure.
%   The function handle must return the following output arguments (in the same order):
%     - x0 = A column vector of floats of length cx = (vs ? k : 1) + (cs * k) + cns representing the initial parameters for variance and independent variables (mandatory).
%     - ai = A float ci-by-cx matrix (-Inf,Inf) representing the "A" element of linear inequality constraints (optional, defaultable to []).
%     - bi = A column vector of floats of length ci representing the "b" element of linear inequality constraints (optional, defaultable to []).
%     - ae = A float ce-by-cx matrix (-Inf,Inf) representing the "Aeq" element of linear equality constraints (optional, defaultable to []).
%     - ce = A column vector of floats of length ce representing the "beq" element of linear equality constraints (optional, defaultable to []).
%     - lb = A column vector of floats of length cx representing the lower bounds (mandatory).
%     - ub = A column vector of floats of length cx representing the upper bounds (mandatory).
% tmm = A float k-by-k matrix [0,1] representing the mask of the stochastic row-wise transition matrix, in which NaNs represent the elements to be computed (optional, by default all the elements are computed).
% fnlcon = A function handle representing a hook for applying non-linear constraints (optional, by default no non-linear constraints are applied).
%   The function handle must accept the following input arguments (in the same order):
%     - x: A column vector of floats representing the current parameters.
%     - dep, indep_s, indep_s, k, vs, tmm: See the description above.
%     - options = A structure representing the optimization options structure.
%   The function handle must return the following output arguments (in the same order):
%     - ci = A column vector of floats representing the "C" element of non-linear constraints (optional, defaultable to []).
%     - ce = A column vector of floats representing the "Ceq" element of non-linear constraints (optional, defaultable to []).
%
% [OUTPUT]
% indep_s_params = A cs-by-1 cell array of row vectors of floats containing the parameters of switching independent variables.
% indep_ns_params = A row vector of floats of length cns representing the parameters of non-switching independent variables.
% s2_params = A row vector of floats of length k representing the variance parameters.
% p = A float k-by-k matrix [0,1] representing the stochastic row-wise transition matrix.
% sprob = A float t-by-k matrix [0,1] representing the smoothed probabilities of each state.
% dur = A row vector of floats of length 4 representing the duration of each state.
% cmu = A column vector of floats of length t representing the conditional means.
% cs2 = A column vector of floats of length t representing the conditional variances.
% e = A column vector of floats of length t representing the standardized residuals.

function [indep_s_params,indep_ns_params,s2_params,p,sprob,dur,cmu,cs2,e] = regime_switching(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('dep',@(x)validateattributes(x,{'double'},{'real' 'finite' 'column' 'nonempty'}));
        ip.addOptional('indep_s',[],@(x)validateattributes(x,{'double'},{'real' 'finite' '2d'}));
        ip.addOptional('indep_ns',[],@(x)validateattributes(x,{'double'},{'real' 'finite' '2d'}));
        ip.addOptional('k',2,@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer' '>=' 2 '<=' 4 'scalar'}));
        ip.addOptional('vs',true,@(x)validateattributes(x,{'logical'},{'scalar'}));
        ip.addOptional('finit',[],@(x)validateattributes(x,{'double' 'function_handle'},{}));
        ip.addOptional('tmm',[],@(x)validateattributes(x,{'double'},{'real' '2d'}));
        ip.addOptional('fnlcon',[],@(x)validateattributes(x,{'double' 'function_handle'},{}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    vs = ipr.vs;
    [dep,indep_s,indep_ns,k,finit,tmm,p0,fnlcon] = validate_input(ipr.dep,ipr.indep_s,ipr.indep_ns,ipr.k,ipr.finit,ipr.tmm,ipr.fnlcon);

    nargoutchk(6,9);

    [indep_s_params,indep_ns_params,s2_params,p,sprob,dur,cmu,cs2,e] = regime_switching_internal(dep,indep_s,indep_ns,k,vs,finit,tmm,p0,fnlcon);

end

function [indep_s_params,indep_ns_params,s2_params,p,sprob,dur,cmu,cs2,e] = regime_switching_internal(dep,indep_s,indep_ns,k,vs,finit,tmm,p0,fnlcon)

    persistent options;

    if (isempty(options))
        options = optimset(optimset(@fmincon),'Algorithm','sqp','Diagnostics','off','Display','off','LargeScale','off','MaxSQPIter',1000,'TolFun',1e-6);
    end

    t = numel(dep);

    [p_params,p_lb,p_ub,p_ae,p_be] = parametrize_p(k,tmm,p0);
    p_params_count = numel(p_params);

    if (isempty(finit))
        x0 = [];
        [ai,bi] = deal([]);
        [ae,be] = deal([]);
        [lb,ub] = deal([]);

        [s2_params,s2_lb,s2_ub] = parametrize_s2(dep,k,vs,1.5,0.75);
        s2_params_count = numel(s2_params);
        x0 = [x0; s2_params];
        lb = [lb; s2_lb];
        ub = [ub; s2_ub];

        [indep_s_params,indep_s_lb,indep_s_ub] = parametrize_indep_s(dep,indep_s,k);
        indep_s_params_count = numel(indep_s_params);
        x0 = [x0; indep_s_params];
        lb = [lb; indep_s_lb];
        ub = [ub; indep_s_ub];

        [indep_ns_params,indep_ns_lb,indep_ns_ub] = parametrize_indep_ns(dep,indep_ns);
        indep_ns_params_count = numel(indep_ns_params);
        x0 = [x0; indep_ns_params];
        lb = [lb; indep_ns_lb];
        ub = [ub; indep_ns_ub];

        params_count = indep_s_params_count + indep_ns_params_count + s2_params_count;
    else
        [x0,ai,bi,ae,be,lb,ub] = finit(dep,indep_s,indep_ns,k,vs,tmm,p0,options);

        if (vs)
            s2_params_count = k;
        else
            s2_params_count = 1;
        end

        indep_s_params_count = size(indep_s,2) * k;
        indep_ns_params_count = size(indep_ns,2);

        params_count = s2_params_count + indep_s_params_count + indep_ns_params_count;

        if (~isa(x0,'double') || ~ismatrix(x0) || (size(x0,2) ~= 1) || (size(x0,1) ~= params_count))
            error(['The function ''finit'' generated an invalid value for ''x0'': it must be a float column vector of length ' num2str(params_count) '.']);
        end

        if (~isempty(ai) || ~isempty(bi))
            if (~isa(ai,'double') || ~ismatrix(ai))
                error(['The function ''finit'' generated an invalid value for ''ai'': it must be a float 2d matrix with no zero-valued dimensions and the number of columns equal to ' num2str(params_count) '.']);
            end

            [ai_r,ai_c] = size(ai);

            if ((ai_r == 0) || (ai_c ~= params_count))
                error(['The function ''finit'' generated an invalid value for ''ai'': it must be a float 2d matrix with no zero-valued dimensions and the number of columns equal to ' num2str(params_count) '.']);
            end

            if (~isa(bi,'double') || ~ismatrix(bi))
                error(['The function ''finit'' generated an invalid value for ''bi'': it must be a a column vector of length ' num2str(ai_r) '.']);
            end

            [bi_r,bi_c] = size(bi);

            if ((bi_r ~= ai_r) || (bi_c ~= 1))
                error(['The function ''finit'' generated an invalid value for ''bi'': it must be a a column vector of length ' num2str(ai_r) '.']);
            end

            ai = [ai zeros(ai_r,p_params_count)];
        end

        if (~isempty(ae) || ~isempty(be))
            if (~isa(ae,'double') || ~ismatrix(ae))
                error(['The function ''finit'' generated an invalid value for ''ae'': it must be a float 2d matrix with no zero-valued dimensions and the number of columns equal to ' num2str(params_count) '.']);
            end

            [ae_r,ae_c] = size(ae);

            if ((ae_r == 0) || (ae_c ~= params_count))
                error(['The function ''finit'' generated an invalid value for ''ae'': it must be a float 2d matrix with no zero-valued dimensions and the number of columns equal to ' num2str(params_count) '.']);
            end

            if (~isa(be,'double') || ~ismatrix(be))
                error(['The function ''finit'' generated an invalid value for ''be'': it must be a a column vector of length ' num2str(ae_r) '.']);
            end

            [be_r,be_c] = size(be);

            if ((be_r ~= ae_r) || (be_c ~= 1))
                error(['The function ''finit'' generated an invalid value for ''be'': it must be a a column vector of length ' num2str(ae_r) '.']);
            end

            ae = [ae zeros(ae_r,p_params_count)];
        end

        if (~isa(lb,'double') || ~ismatrix(lb) || (size(lb,2) ~= 1) || (size(lb,1) ~= params_count))
            error(['The function ''finit'' generated an invalid value for ''lb'': it must be a float column vector of length ' num2str(params_count) '.']);
        end

        if (~isa(ub,'double') || ~ismatrix(ub) || (size(ub,2) ~= 1) || (size(ub,1) ~= params_count))
            error(['The function ''finit'' generated an invalid value for ''ub'': it must be a float column vector of length ' num2str(params_count) '.']);
        end
    end

    x0 = [x0; p_params];
    ae = [ae; zeros(k,params_count) p_ae];
    be = [be; p_be];
    lb = [lb; p_lb];
    ub = [ub; p_ub];

    if (isempty(fnlcon))
        params = fmincon(@(x)likelihood(x,dep,indep_s,indep_ns,k,vs,tmm),x0,ai,bi,ae,be,lb,ub,[],options);
    else
        params = fmincon(@(x)likelihood(x,dep,indep_s,indep_ns,k,vs,tmm),x0,ai,bi,ae,be,lb,ub,@(x)fnlcon(x,dep,indep_s,indep_ns,k,vs,tmm,options),options);
    end

    [~,mu,g] = likelihood(params,dep,indep_s,indep_ns,k,vs,tmm);

    if (vs)
        s2_params = params(1:k).';
        o = k + 1;
    else
        s2_params = ones(1,k) .* params(1);
        o = 2;
    end

    indep_s_count = size(indep_s,2);
    indep_s_params = num2cell(reshape(params(o:o+indep_s_params_count-1),k,indep_s_count).',2);
    o = o + indep_s_params_count;

    if (indep_ns_params_count > 0)
        indep_ns_params = params(o:o+indep_ns_params_count-1).';
        o = o + indep_ns_params_count;
    else
        indep_ns_params = [];
    end

    p = tmm;
    p(isnan(tmm)) = params(o:end);

    pt = p.';

    prob = [ones(1,k) .* (1 / k); zeros(t - 1,k)];

    for i = 2:t
        prob(i,:) = pt * g(i-1,:).';
    end

    sprob = [zeros(t - 1,k); g(t,:)];

    for i = t-1:-1:1
        for j = 1:k
            sprob(i,j) = sum((sprob(i+1,:) .* g(i,j) .* pt(:,j).') ./ prob(i+1,:),'omitnan');
        end
    end

    dur = round(1 ./ (1 - diag(p).'),0);

    cmu = sum(mu .* prob,2);
    cs2 = sum(repmat(sqrt(s2_params),t,1) .* prob,2);

    e = dep - cmu;
    e = (e - mean(e)) ./ std(e);

end

function [ll,mu,g] = likelihood(x,dep,indep_s,indep_ns,k,vs,tmm)

    t = numel(dep);

    if (vs)
        c = x(1:k);
        o = k + 1;
    else
        c = ones(k,1) .* x(1);
        o = 2;
    end

    indep_s_count = size(indep_s,2);
    indep_s_params_count = indep_s_count * k;
    indep_s_params = reshape(x(o:o+indep_s_params_count-1),k,indep_s_count).';
    o = o + indep_s_params_count;

    indep_ns_count = size(indep_ns,2);

    if (indep_ns_count == 0)
        indep_ns = zeros(t,1);
        indep_ns_params = 0;
    else
        indep_ns_params = x(o:o+indep_ns_count-1);
        o = o + indep_ns_count;
    end

    p = tmm;
    p(isnan(tmm)) = x(o:end);

    pt = p.';

    mu = zeros(t,k);
    z = zeros(t,k);
    nc = (2 * pi())^0.5;

    for i = 1:k
        c_i = c(i);

        mu_i = (indep_s * indep_s_params(:,i)) + (indep_ns * indep_ns_params);
        e_i = dep - mu_i;

        mu(:,i) = mu_i;
        z(:,i) = (1 / (nc * sqrt(c_i))) .* exp(-0.5 .* sum((e_i / c_i) .* e_i,2));
    end

    w1 = (pt * (ones(k,1) .* (1 / k))) .* z(1,:).';
    f1 = ones(1,k) * w1;

    f = [f1; zeros(t - 1,1)];
    g = [(w1 ./ f1).'; zeros(t - 1,k)];

    for i = 2:t
        wi = (pt * g(i-1,:).') .* z(i,:).';
        fi = ones(1,k) * wi;

        f(i,1) = fi;
        g(i,:) = wi ./ fi;
    end

    ll_v = log(f(2:end));

    if (any(~isfinite(ll_v)))
        ll = Inf;
    else
        ll = -sum(ll_v);
    end

end

function [params,lb,ub] = parametrize_indep_ns(dep,indep_ns)

    indep_ns_count = size(indep_ns,2);

    if (indep_ns_count > 0)
        params = regress(dep,indep_ns);
        lb = -Inf(indep_ns_count,1);
        ub = Inf(indep_ns_count,1);
    else
        params = [];
        lb = [];
        ub = [];
    end

end

function [params,lb,ub] = parametrize_indep_s(dep,indep_s,k)

    indep_s_count = size(indep_s,2);

    b = regress(dep,indep_s);
    b_factor = 1;

    indep_s_params = zeros(indep_s_count,k);

    for i = 1:k
        indep_s_params(:,i) = b * b_factor;
        b_factor = b_factor * -1;
    end

    params = reshape(indep_s_params.',indep_s_count * k,1);
    lb = -Inf(numel(params),1);
    ub = Inf(numel(params),1);

end

function [params,lb,ub,ae,be] = parametrize_p(k,tmm,p0)

    tmm_nans = isnan(tmm);

    if (all(all(tmm_nans)))
        params = p0(:);

        lb = zeros(k^2,1);
        ub = ones(k^2,1) - 1e-4;

        ae = repmat(eye(k),1,k);
        be = ones(k,1);
    else
        filter = ~tmm_nans(:);

        params = p0(:);
        params(filter) = [];

        lb = zeros(numel(params),1);
        ub = ones(numel(params),1);

        ae = repmat(eye(k),1,k);
        ae(:,filter) = [];

        be = ones(k,1) -  sum(p0 .* ~tmm_nans,2);
    end

end

function [params,lb,ub] = parametrize_s2(dep,k,vs,factor,multiplier)

    mm = (dep - mean(dep)).^2;
    s2 = var(dep);

    if (vs)
        params = zeros(k,1);

        for i = 1:k
            params(i) = s2 * multiplier;
            multiplier = multiplier * factor;
        end

        lb = ones(k,1) .* min(mm);
        ub = ones(k,1) .* max(mm);
    else
        params = s2 * multiplier;

        lb = min(mm);
        ub = max(mm);
    end

end

function [dep,indep_s,indep_ns,k,finit,tmm,p0,fnlcon] = validate_input(dep,indep_s,indep_ns,k,finit,tmm,fnlcon)

    dep = dep(:);
    t = numel(dep);

    if (~isempty(indep_s))
        if (size(indep_s,1) ~= t)
            error(['The value of ''indep_s'' is invalid. Expected input to have ' num2str(t) ' rows.']);
        end

        if (all(indep_s(:,1) == 1))
            error('The value of ''indep_s'' is invalid. Expected input to exclude the intercept because it is internally handled by the model.');
        end

        if (any(sum(indep_s == 0,1) == t))
            error('The value of ''indep_s'' is invalid. Expected input to contain no zero-valued vectors.');
        end
    end

    indep_s = [ones(t,1) indep_s];

    if (~isempty(indep_ns) && (size(indep_ns,1) ~= t))
        if (size(indep_ns,1) ~= t)
            error(['The value of ''indep_ns'' is invalid. Expected input to have ' num2str(t) ' rows.']);
        end

        if (any(sum(indep_ns,1) == t) || any(sum(indep_ns == 0,1) == t))
            error('The value of ''indep_ns'' is invalid. Expected input to contain no zero-valued or one-valued vectors.');
        end
    end

    if (~isempty(finit))
        if (~isa(finit,'function_handle') || ~isscalar(finit))
            error('The value of ''finit'' is invalid. Expected input to be a single function handle.');
        end

        if (nargin(finit) ~= 8)
            error('The value of ''finit'' is invalid. Expected input to accept 8 input arguments.');
        end

        if (nargout(finit) ~= 7)
            error('The value of ''finit'' is invalid. Expected input to accept 7 output arguments.');
        end
    end

    if (isempty(tmm))
        tmm = NaN(k,k);
        p0 = repmat(0.1,k,k) + (eye(k) * (1 - (k * 0.1)));
    else
        if (any(size(tmm) ~= k))
            error(['The value of ''tmm'' is invalid. Expected input to be a matrix of size ' num2str(k) 'x' num2str(k) '.']);
        end

        tmm_nans = isnan(tmm);
        p0 = zeros(k,k);

        for i = 1:k
            tmm_i = tmm(i,:);
            tmm_nans_i = tmm_nans(i,:);
            tmm_nans_i_count = sum(tmm_nans_i);

            if (tmm_nans_i_count == k)
                p0(i,i) = (1 - (k * 0.1)) + 0.1;

                tmm_nans_i(i) = 0;
                p0(i,tmm_nans_i) = 0.1;
            else
                if (any((tmm_i < 0) | (tmm_i >= 1)))
                    error(['The value of ''tmm'' is invalid. Expected input to have a valid definition for row ' num2str(i) ': all the constrained elements must have a value in range [0,1).']);
                end

                if (tmm_nans_i_count == 0)
                    if (sum(tmm_i) ~= 1)
                        error(['The value of ''tmm'' is invalid. Expected input to have a valid definition for row ' num2str(i) ': the sum of elements must be equal to 1.']);
                    end
                elseif (tmm_nans_i_count == 1)
                    if (k == 2)
                        error(['The value of ''tmm'' is invalid. Expected input to have a valid definition for row ' num2str(i) ': when the number of states is equal to 2, only fully constrained or fully unconstrained rows are accepted.']);
                    else
                        error(['The value of ''tmm'' is invalid. Expected input to have a valid definition for row ' num2str(i) ': the number of unconstrained elements must be greater than 1.']);
                    end
                else
                    if (sum(tmm_i(~tmm_nans_i)) >= 1)
                        error(['The value of ''tmm'' is invalid. Expected input to have a valid definition for row ' num2str(i) ': the number of constrained elements must less than 1.']);
                    end
                end

                if (~tmm_nans_i(i) && (tmm_i(i) == 0))
                    error(['The value of ''tmm'' is invalid. Expected input to have a valid definition for row ' num2str(i) ': diagonal elements cannot be constrained to a value equal to 0.']);
                end

                p0(i,~tmm_nans_i) = tmm_i(~tmm_nans_i);

                k_i = sum(tmm_nans_i);
                delta = 1 - sum(tmm_i(~tmm_nans_i));

                if (tmm_nans_i(i))
                    pd = 0.7 * delta;
                    pnd = (delta - pd) / (k_i - 1);

                    p0(i,i) = pd;

                    tmm_nans_i(i) = 0;
                    p0(i,tmm_nans_i) = pnd;
                else
                    p0(i,tmm_nans_i) = delta / sum(tmm_nans_i);
                end
            end
        end

        we = (k - 1)^2 + 1;
        q = p0^we;
        is_ergodic = ~any(q(:) < (k * eps()));

        if (~is_ergodic)
            error('The value of ''tmm'' is invalid. Expected input to produce an ergodic stochastic row-wise transition matrix.');
        end
    end

    if (~isempty(fnlcon))
        if (~isa(fnlcon,'function_handle') || ~isscalar(fnlcon))
            error('The value of ''fnlcon'' is invalid. Expected input to be a single function handle.');
        end

        if (nargin(fnlcon) ~= 8)
            error('The value of ''fnlcon'' is invalid. Expected input to accept 8 input arguments.');
        end

        if (nargout(fnlcon) ~= 2)
            error('The value of ''fnlcon'' is invalid. Expected input to accept 2 output arguments.');
        end
    end

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
function [is,os,str] = strengths_dir(CIJ)
%STRENGTHS_DIR      In-strength and out-strength
%
%   [is,os,str] = strengths_dir(CIJ);
%
%   Node strength is the sum of weights of links connected to the node. The
%   instrength is the sum of inward link weights and the outstrength is the
%   sum of outward link weights.
%
%   Input:      CIJ,    directed weighted connection matrix
%
%   Output:     is,     node instrength
%               os,     node outstrength
%               str,    node strength (instrength + outstrength)
%
%   Notes:  Inputs are assumed to be on the columns of the CIJ matrix.
%
%
%   Olaf Sporns, Indiana University, 2002/2006/2008

% compute strengths
is = sum(CIJ,1);    % instrength = column sum of CIJ
os = sum(CIJ,2)';   % outstrength = row sum of CIJ
str = is+os;        % strength = instrength+outstrength

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function benchmark_tracker_wrapper(tracker_name, runfile_name, do_cleanup)

if nargin < 3
    do_cleanup = true;
end

% *************************************************************
% VOT: Always call exit command at the end to terminate Matlab!
% *************************************************************
if do_cleanup
    cleanup = onCleanup(@() exit() );
else
    [pathstr, ~, ~] = fileparts(mfilename('fullpath'));
    cd_ind = strfind(pathstr, filesep());
    pathstr = pathstr(1:cd_ind(end)-1);
    cleanup = onCleanup(@() cd(pathstr));
end

try

% *************************************************************
% VOT: Set random seed to a different value every time.
% *************************************************************
RandStream.setGlobalStream(RandStream('mt19937ar', 'Seed', sum(clock)));

seq.format = 'vot';

setup_tracker_paths(tracker_name);

disp([runfile_name '(seq, [], []);']);
res = eval([runfile_name '(seq, [], []);']);

catch err
    [wrapper_pathstr, ~, ~] = fileparts(mfilename('fullpath'));
    cd_ind = strfind(wrapper_pathstr, filesep());
    VOT_path = wrapper_pathstr(1:cd_ind(end));

    error_report_path = [VOT_path 'error_reports/'];
    if ~exist(error_report_path, 'dir')
        mkdir(error_report_path);
    end

    report_file_name = [error_report_path tracker_name '_' runfile_name datestr(now,'_yymmdd_HHMM') '.mat'];

    save(report_file_name, 'err')

    rethrow(err);
end

author: martin-danelljan

repo: ECO


text:
function [data,p]=create_data(p)
if nargin==0
    p=struct();
end

gygprrp=load('gygprrp.dat');

p.zss = exp(mean(gygprrp(:,1)));
p.beta = exp(log(p.zss)-mean(gygprrp(:,3)));

start_date='1959Q1';

% per capital GDP growth
gyt = gygprrp(:,1) - log(p.zss);
% Inflation growth
gpt = gygprrp(:,2);
% Real interest rate
rrpt = gygprrp(:,3) - log(p.zss) + log(p.beta);

data=ts(start_date,[gyt,gpt,rrpt],{'GY_HAT','GPAI_HAT','RRPAI_HAT'});
data=pages2struct(data);

end

author: jmaih


text:
function preproc = spm_cfg_preproc8
% Configuration file for 'Combined Segmentation and Spatial Normalisation'
%__________________________________________________________________________
% Copyright (C) 2008-2016 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_cfg_preproc8.m 7629 2019-06-27 12:35:45Z john $

%--------------------------------------------------------------------------
% vols Volumes
%--------------------------------------------------------------------------
vols         = cfg_files;
vols.tag     = 'vols';
vols.name    = 'Volumes';
vols.help    = {
    'Select scans from this channel for processing.'
    'If multiple channels are used (eg T1 & T2), then the same order of subjects must be specified for each channel and they must be in register (same position, size, voxel dims etc..).'
    }';
vols.filter  = 'image';
vols.ufilter = '.*';
vols.num     = [1 Inf];
vols.preview = @(f) spm_check_registration(char(f));

%--------------------------------------------------------------------------
% biasreg Bias regularisation
%--------------------------------------------------------------------------
biasreg         = cfg_menu;
biasreg.tag     = 'biasreg';
biasreg.name    = 'Bias regularisation';
biasreg.help    = {
    'MR images are usually corrupted by a smooth, spatially varying artifact that modulates the intensity of the image (bias). These artifacts, although not usually a problem for visual inspection, can impede automated processing of the images.'
    ''
    'An important issue relates to the distinction between intensity variations that arise because of bias artifact due to the physics of MR scanning, and those that arise due to different tissue properties.  The objective is to model the latter by different tissue classes, while modelling the former with a bias field. We know a priori that intensity variations due to MR physics tend to be spatially smooth, whereas those due to different tissue types tend to contain more high frequency information. A more accurate estimate of a bias field can be obtained by including prior knowledge about the distribution of the fields likely to be encountered by the correction algorithm. For example, if it is known that there is little or no intensity non-uniformity, then it would be wise to penalise large values for the intensity non-uniformity parameters. This regularisation can be placed within a Bayesian context, whereby the penalty incurred is the negative logarithm of a prior probability for any particular pattern of non-uniformity.'
    'Knowing what works best should be a matter of empirical exploration.  For example, if your data has very little intensity non-uniformity artifact, then the bias regularisation should be increased.  This effectively tells the algorithm that there is very little bias in your data, so it does not try to model it.'
    }';
biasreg.labels = {
                  'no regularisation (0)'
                  'extremely light regularisation (0.00001)'
                  'very light regularisation (0.0001)'
                  'light regularisation (0.001)'
                  'medium regularisation (0.01)'
                  'heavy regularisation (0.1)'
                  'very heavy regularisation (1)'
                  'extremely heavy regularisation (10)'
                  }';
biasreg.values = {
                  0
                  1e-05
                  0.0001
                  0.001
                  0.01
                  0.1
                  1
                  10
                  }';
biasreg.val    = {0.001};

%--------------------------------------------------------------------------
% biasfwhm Bias FWHM
%--------------------------------------------------------------------------
biasfwhm         = cfg_menu;
biasfwhm.tag     = 'biasfwhm';
biasfwhm.name    = 'Bias FWHM';
biasfwhm.help    = {
    'FWHM of Gaussian smoothness of bias.'
    'If your intensity non-uniformity is very smooth, then choose a large FWHM. This will prevent the algorithm from trying to model out intensity variation due to different tissue types. The model for intensity non-uniformity is one of i.i.d. Gaussian noise that has been smoothed by some amount, before taking the exponential. Note also that smoother bias fields need fewer parameters to describe them. This means that the algorithm is faster for smoother intensity non-uniformities.'
    }';
biasfwhm.labels = {
                   '30mm cutoff'
                   '40mm cutoff'
                   '50mm cutoff'
                   '60mm cutoff'
                   '70mm cutoff'
                   '80mm cutoff'
                   '90mm cutoff'
                   '100mm cutoff'
                   '110mm cutoff'
                   '120mm cutoff'
                   '130mm cutoff'
                   '140mm cutoff'
                   '150mm cutoff'
                   'No correction'
                   }';
biasfwhm.values = {
                   30
                   40
                   50
                   60
                   70
                   80
                   90
                   100
                   110
                   120
                   130
                   140
                   150
                   Inf
                   }';
biasfwhm.val    = {60};

%--------------------------------------------------------------------------
% write Save Bias Corrected
%--------------------------------------------------------------------------
write         = cfg_menu;
write.tag     = 'write';
write.name    = 'Save Bias Corrected';
write.help    = {
    'Option to save a bias corrected version of your images from this channel, or/and the estimated bias field.'
    'MR images are usually corrupted by a smooth, spatially varying artifact that modulates the intensity of the image (bias). These artifacts, although not usually a problem for visual inspection, can impede automated processing of the images.  The bias corrected version should have more uniform intensities within the different types of tissues.'
    }';
write.labels = {
                'Save Nothing'
                'Save Bias Corrected'
                'Save Bias Field'
                'Save Field and Corrected'
                }';
write.values = {
                [0 0]
                [0 1]
                [1 0]
                [1 1]
                }';
write.val    = {[0 0]};

%--------------------------------------------------------------------------
% channel Channel
%--------------------------------------------------------------------------
channel         = cfg_branch;
channel.tag     = 'channel';
channel.name    = 'Channel';
channel.val     = {vols biasreg biasfwhm write };
channel.help    = {
    'Specify a channel for processing.'
    'If multiple channels are used (eg PD & T2), then the same order of subjects must be specified for each channel and they must be in register (same position, size, voxel dims etc..). The different channels can be treated differently in terms of inhomogeneity correction etc. You may wish to correct some channels and save the corrected images, whereas you may wish not to do this for other channels.'
    }';

%--------------------------------------------------------------------------
% data Data
%--------------------------------------------------------------------------
data         = cfg_repeat;
data.tag     = 'data';
data.name    = 'Data';
data.val     = {channel };
data.help    = {
    'Specify the number of different channels (for multi-spectral classification).'
    'If you have scans of different contrasts for each of the subjects, then it is possible to combine the information from them in order to improve the segmentation accuracy. Note that only the first channel of data is used for the initial affine registration with the tissue probability maps.'
    }';
data.values  = {channel };
data.num     = [1 Inf];

%--------------------------------------------------------------------------
% tpm Tissue probability map
%--------------------------------------------------------------------------
tpm         = cfg_files;
tpm.tag     = 'tpm';
tpm.name    = 'Tissue probability map';
tpm.help    = {
    'Select the tissue probability image for this class.'
    'These should be maps of eg grey matter, white matter or cerebro-spinal fluid probability. A nonlinear deformation field is estimated that best overlays the tissue probability maps on the individual subjects'' image.'
    ''
    'Rather than assuming stationary prior probabilities based upon mixing proportions, additional information is used, based on other subjects'' brain images.  Priors are usually generated by registering a large number of subjects together, assigning voxels to different tissue types and averaging tissue classes over subjects. Three tissue classes are used: grey matter, white matter and cerebro-spinal fluid. A fourth class is also used, which is simply one minus the sum of the first three. These maps give the prior probability of any voxel in a registered image being of any of the tissue classes - irrespective of its intensity.'
    ''
    'The model is refined further by allowing the tissue probability maps to be deformed according to a set of estimated parameters. This allows spatial normalisation and segmentation to be combined into the same model.'
    }';
tpm.filter  = 'image';
tpm.ufilter = '.*';
tpm.num     = [1 1];
tpm.preview = @(f) spm_image('Display',char(f));

%--------------------------------------------------------------------------
% ngaus Num. Gaussians
%--------------------------------------------------------------------------
ngaus         = cfg_menu;
ngaus.tag     = 'ngaus';
ngaus.name    = 'Num. Gaussians';
ngaus.help    = {
    'The number of Gaussians used to represent the intensity distribution for each tissue class can be greater than one.'
    'In other words, a tissue probability map may be shared by several clusters. The assumption of a single Gaussian distribution for each class does not hold for a number of reasons. In particular, a voxel may not be purely of one tissue type, and instead contain signal from a number of different tissues (partial volume effects). Some partial volume voxels could fall at the interface between different classes, or they may fall in the middle of structures such as the thalamus, which may be considered as being either grey or white matter. Various other image segmentation approaches use additional clusters to model such partial volume effects. These generally assume that a pure tissue class has a Gaussian intensity distribution, whereas intensity distributions for partial volume voxels are broader, falling between the intensities of the pure classes. Unlike these partial volume segmentation approaches, the model adopted here simply assumes that the intensity distribution of each class may not be Gaussian, and assigns belonging probabilities according to these non-Gaussian distributions. Typical numbers of Gaussians could be two for grey matter, two for white matter, two for CSF, three for bone, four for other soft tissues and two for air (background).'
    'Note that if any of the Num. Gaussians is set to non-parametric, then a non-parametric approach will be used to model the tissue intensities. This may work for some images (eg CT), but not others - and it has not been optimised for multi-channel data. Note that it is likely to be especially problematic for images with poorly behaved intensity histograms due to aliasing effects that arise from having discrete values on the images.'
    }';
ngaus.labels = {
                '1'
                '2'
                '3'
                '4'
                '5'
                '6'
                '7'
                '8'
                'Nonparametric'
                }';
ngaus.values = {
                1
                2
                3
                4
                5
                6
                7
                8
                Inf
                }';
ngaus.val    = {Inf};

%--------------------------------------------------------------------------
% native Native Tissue
%--------------------------------------------------------------------------
native         = cfg_menu;
native.tag     = 'native';
native.name    = 'Native Tissue';
native.help    = {'The native space option allows you to produce a tissue class image (c*) that is in alignment with the original/* (see Figure \ref{seg1})*/. It can also be used for ``importing'''' into a form that can be used with the Dartel toolbox (rc*).'};
native.labels = {
    'None'
    'Native Space'
    'Dartel Imported'
    'Native + Dartel Imported'
    }';
native.values = {
                 [0 0]
                 [1 0]
                 [0 1]
                 [1 1]
                 }';
native.val    = {[1 0]};

%--------------------------------------------------------------------------
% warped Warped Tissue
%--------------------------------------------------------------------------
warped         = cfg_menu;
warped.tag     = 'warped';
warped.name    = 'Warped Tissue';
warped.help    = {
    'You can produce spatially normalised versions of the tissue class - both with (mwc*) and without (wc*) modulation (see below). These can be used for voxel-based morphometry. All you need to do is smooth them and do the stats.'
    ''
    '``Modulation'''' is to compensate for the effect of spatial normalisation.  When warping a series of images to match a template, it is inevitable that volumetric differences will be introduced into the warped images.  For example, if one subject''s temporal lobe has half the volume of that of the template, then its volume will be doubled during spatial normalisation. This will also result in a doubling of the voxels labelled grey matter.  In order to remove this confound, the spatially normalised grey matter (or other tissue class) is adjusted by multiplying by its relative volume before and after warping.  If warping results in a region doubling its volume, then the correction will halve the intensity of the tissue label. This whole procedure has the effect of preserving the total amount of grey matter signal in the normalised partitions.  Actually, in this version of SPM the warped data are not scaled by the Jacobian determinants when generating the "modulated" data.  Instead, the original voxels are projected into their new location in the warped images.  This exactly preserves the tissue count, but has the effect of introducing aliasing artifacts - especially if the original data are at a lower resolution than the warped images.  Smoothing should reduce this artifact though.'
    'Note also that the "unmodulated" data are generated slightly differently in this version of SPM. In this version, the projected data are corrected using a kind of smoothing procedure. This is not done exactly as it should be done (to save computational time), but it does a reasonable job. It also has the effect of extrapolating the warped tissue class images beyond the range of the original data.  This extrapolation is not perfect, as it is only an estimate, but it may still be a good thing to do.'
    }';
warped.labels = {
                 'None'
                 'Modulated'
                 'Unmodulated'
                 'Modulated + Unmodulated'
                 }';
warped.values = {
                 [0 0]
                 [0 1]
                 [1 0]
                 [1 1]
                 }';
warped.val    = {[0 0]};

%--------------------------------------------------------------------------
% tissue Tissue
%--------------------------------------------------------------------------
tissue         = cfg_branch;
tissue.tag     = 'tissue';
tissue.name    = 'Tissue';
tissue.val     = {tpm ngaus native warped};
tissue.help    = {'A number of options are available for each of the tissues.  You may wish to save images of some tissues, but not others. If planning to use Dartel, then make sure you generate ``imported'''' tissue class images of grey and white matter (and possibly others).  Different numbers of Gaussians may be needed to model the intensity distributions of the various tissues.'};

%--------------------------------------------------------------------------
% tissues Tissues
%--------------------------------------------------------------------------
tissues         = cfg_repeat;
tissues.tag     = 'tissues';
tissues.name    = 'Tissues';
tissues.values  = {tissue};
tissues.num     = [0 Inf];

tissues.val     = {tissue tissue tissue tissue tissue tissue};
tpm_nam = fullfile(spm('dir'),'tpm','TPM.nii');
ngaus   = [1 1 2 3 4 2];
% Change to: ngaus   = [2 2 2 3 4 2];
nval    = {[1 0],[1 0],[1 0],[1 0],[1 0],[0 0]};
for k=1:numel(ngaus),
    tissue.val{1}.val = {{[tpm_nam ',' num2str(k)]}};
    tissue.val{2}.val = {ngaus(k)};
    tissue.val{3}.val = {nval{k}};
    tissues.val{k}    = tissue;
end

tissues.help    = {'The data for each subject are classified into a number of different tissue types.  The tissue types are defined according to tissue probability maps, which define the prior probability of finding a tissue type at a particular location. Typically, the order of tissues is grey matter, white matter, CSF, bone, soft tissue and air/background (if using tpm/TPM.nii).'};
tissues.preview = @(f) spm_check_registration(char([f.tpm]));

%--------------------------------------------------------------------------
% mrf MRF Parameter
%--------------------------------------------------------------------------
mrf         = cfg_entry;
mrf.tag     = 'mrf';
mrf.name    = 'MRF Parameter';
mrf.help    = {'When tissue class images are written out, a few iterations of a simple Markov Random Field (MRF) cleanup procedure are run.  This parameter controls the strength of the MRF. Setting the value to zero will disable the cleanup.'};
mrf.strtype = 'r';
mrf.num     = [1 1];
mrf.val     = {1};

%--------------------------------------------------------------------------
% cleanup Clean up any partitions
%--------------------------------------------------------------------------
cleanup         = cfg_menu;
cleanup.tag     = 'cleanup';
cleanup.name    = 'Clean Up';
cleanup.help    = {
    'This uses a crude routine for extracting the brain from segmented images.'
    'It begins by taking the white matter, and eroding it a couple of times to get rid of any odd voxels.  The algorithm continues on to do conditional dilations for several iterations, where the condition is based upon gray or white matter being present.This identified region is then used to clean up the grey and white matter partitions.  Note that the fluid class will also be cleaned, such that aqueous and vitreous humour in the eyeballs, as well as other assorted fluid regions (except CSF) will be removed.'
    ''
    'If you find pieces of brain being chopped out in your data, then you may wish to disable or tone down the cleanup procedure. Note that the procedure uses a number of assumptions about what each tissue class refers to.  If a different set of tissue priors are used, then this routine should be disabled.'
    }';
cleanup.labels = {
                  'Dont do cleanup'
                  'Light Clean'
                  'Thorough Clean'
}';
cleanup.values = {0 1 2};
cleanup.val    = {1};

%--------------------------------------------------------------------------
% reg Warping Regularisation
%--------------------------------------------------------------------------
reg         = cfg_entry;
reg.tag     = 'reg';
reg.name    = 'Warping Regularisation';
reg.help    = {...
    'Registration involves simultaneously minimising two terms.  One of these is a measure of similarity between the images (mean-squared difference in the current situation), whereas the other is a measure of the roughness of the deformations.  This measure of roughness involves the sum of the following terms:',...
    '* Absolute displacements need to be penalised by a tiny amount.  The first element encodes the amount of penalty on these.  Ideally, absolute displacements should not be penalised, but it is necessary for technical reasons.',...
    '* The `membrane energy'' of the deformation is penalised (2nd element), usually by a relatively small amount. This penalises the sum of squares of the derivatives of the velocity field (ie the sum of squares of the elements of the Jacobian tensors).',...
    '* The `bending energy'' is penalised (3rd element). This penalises the sum of squares of the 2nd derivatives of the velocity.',...
    '* Linear elasticity regularisation is also included (4th and 5th elements).  The first parameter (mu) is similar to that for linear elasticity, except it penalises the sum of squares of the Jacobian tensors after they have been made symmetric (by averaging with the transpose).  This term essentially penalises length changes, without penalising rotations.',...
    '* The final term also relates to linear elasticity, and is the weight that denotes how much to penalise changes to the divergence of the velocities (lambda).  This divergence is a measure of the rate of volumetric expansion or contraction.',...
    'The amount of regularisation determines the tradeoff between the terms. More regularisation gives smoother deformations, where the smoothness measure is determined by the bending energy of the deformations.'};
reg.strtype = 'r';
reg.num     = [1  5];
reg.val     = {[0 0.001 0.5 0.05 0.2]};

%--------------------------------------------------------------------------
% affreg Affine Regularisation
%--------------------------------------------------------------------------
affreg         = cfg_menu;
affreg.tag     = 'affreg';
affreg.name    = 'Affine Regularisation';
affreg.help    = {
    'The procedure is a local optimisation, so it needs reasonable initial starting estimates. Images should be placed in approximate alignment using the Display function of SPM before beginning. A Mutual Information affine registration with the tissue probability maps (D''Agostino et al, 2004) is used to achieve approximate alignment. Note that this step does not include any model for intensity non-uniformity. This means that if the procedure is to be initialised with the affine registration, then the data should not be too corrupted with this artifact.If there is a lot of intensity non-uniformity, then manually position your image in order to achieve closer starting estimates, and turn off the affine registration.'
    ''
    'Affine registration into a standard space can be made more robust by regularisation (penalising excessive stretching or shrinking).  The best solutions can be obtained by knowing the approximate amount of stretching that is needed (e.g. ICBM templates are slightly bigger than typical brains, so greater zooms are likely to be needed). For example, if registering to an image in ICBM/MNI space, then choose this option.  If registering to a template that is close in size, then select the appropriate option for this.'
    }';
affreg.labels = {
                 'No Affine Registration'
                 'ICBM space template - European brains'
                 'ICBM space template - East Asian brains'
                 'Average sized template'
                 'No regularisation'
                 }';
affreg.values = {
                 ''
                 'mni'
                 'eastern'
                 'subj'
                 'none'
                 }';
affreg.val    = {'mni'};

%--------------------------------------------------------------------------
% fwhm Smoothness
%--------------------------------------------------------------------------
smo         = cfg_entry;
smo.tag     = 'fwhm';
smo.name    = 'Smoothness';
smo.help    = {
    'This is used to derive a fudge factor to account for correlations between neighbouring voxels.'
    'Smoother data have more spatial correlations, rendering the assumptions of the model inaccurate.'
    'For PET or SPECT, set this value to about 5 mm, or more if the images have smoother noise.  For MRI, you can usually use a value of 0 mm.'
    }';
smo.strtype = 'r';
smo.num     = [1  1];
smo.val     = {0};

%--------------------------------------------------------------------------
% samp Sampling distance
%--------------------------------------------------------------------------
samp         = cfg_entry;
samp.tag     = 'samp';
samp.name    = 'Sampling distance';
samp.help    = {
    'This encodes the approximate distance between sampled points when estimating the model parameters.'
    'Smaller values use more of the data, but the procedure is slower and needs more memory. Determining the ``best'''' setting involves a compromise between speed and accuracy.'
    }';
samp.strtype = 'r';
samp.num     = [1  1];
samp.val     = {3};

%--------------------------------------------------------------------------
% write Deformation Fields
%--------------------------------------------------------------------------
write         = cfg_menu;
write.tag     = 'write';
write.name    = 'Deformation Fields';
write.help    = {
    'Deformation fields can be saved to disk, and used by the Deformations Utility.'
    'For spatially normalising images to MNI space, you will need the forward deformation, whereas for spatially normalising (eg) GIFTI surface files, you''ll need the inverse. It is also possible to transform data in MNI space on to the individual subject, which also requires the inverse transform. Deformations are saved as .nii files, which contain three volumes to encode the x, y and z coordinates.'
    };
write.labels = {
                'None'
                'Inverse'
                'Forward'
                'Inverse + Forward'
                }';
write.values = {
                [0 0]
                [1 0]
                [0 1]
                [1 1]
                }';
write.val    = {[0 0]};

%--------------------------------------------------------------------------
% bb Bounding box
%--------------------------------------------------------------------------
bb         = cfg_entry;
bb.tag     = 'bb';
bb.name    = 'Bounding box';
bb.help    = {'The bounding box (in mm) of the volume which is to be written (relative to the anterior commissure).'};
bb.strtype = 'r';
bb.num     = [2 3];
bb.val     = {[NaN NaN NaN; NaN NaN NaN]};
bb.hidden  = true;

%--------------------------------------------------------------------------
% vox Voxel sizes
%--------------------------------------------------------------------------
vox         = cfg_entry;
vox.tag     = 'vox';
vox.name    = 'Voxel sizes';
vox.help    = {'The voxel size (isotropic, in mm) of the written normalised or imported images.'};
vox.strtype = 'r';
vox.num     = [1 1];
vox.val     = {NaN};
vox.hidden  = true;

%--------------------------------------------------------------------------
% warp Warping
%--------------------------------------------------------------------------
warp         = cfg_branch;
warp.tag     = 'warp';
warp.name    = 'Warping & MRF';
warp.val     = {mrf cleanup reg affreg smo samp write vox bb};
warp.help    = {
    'A number of warping options.'
    'The main one that you could consider changing is the one for specifying whether deformation fields or inverse deformation fields should be generated.'
    }';

%--------------------------------------------------------------------------
% preproc Segment
%--------------------------------------------------------------------------
preproc      = cfg_exbranch;
preproc.tag  = 'preproc';
preproc.name = 'Segment';
preproc.val  = {data tissues warp};
preproc.help = {
    'Segmentation, bias correction and spatially normalisation - all in the same model.'
    ''
    'This procedure is an extension of the old unified segmentation algorithm (and was known as "New Segment" in SPM8). The algorithm is essentially the same as that described in the Unified Segmentation paper /* \cite{ashburner05}*/, except for (i) a slightly different treatment of the mixing proportions, (ii) the use of an improved registration model, (iii) the ability to use multi-spectral data, (iv) an extended set of tissue probability maps, which allows a different treatment of voxels outside the brain. Some of the options in the toolbox do not yet work, and it has not yet been seamlessly integrated into the SPM8 software.  Also, the extended tissue probability maps need further refinement. The current versions were crudely generated (by JA) using data that was kindly provided by Cynthia Jongen of the Imaging Sciences Institute at Utrecht, NL.'
    ''
    'This function segments, bias corrects and spatially normalises - all in the same model/* \cite{ashburner05}*/. Many investigators use tools within older versions of SPM for a technique that has become known as "optimised" voxel-based morphometry (VBM). VBM performs region-wise volumetric comparisons among populations of subjects. It requires the images to be spatially normalised, segmented into different tissue classes, and smoothed, prior to performing statistical tests/* \cite{wright_vbm,am_vbmreview,ashburner00b,john_should}*/. The "optimised" pre-processing strategy involved spatially normalising subjects'' brain images to a standard space, by matching grey matter in these images, to a grey matter reference.  The historical motivation behind this approach was to reduce the confounding effects of non-brain (e.g. scalp) structural variability on the registration. Tissue classification in older versions of SPM required the images to be registered with tissue probability maps. After registration, these maps represented the prior probability of different tissue classes being found at each location in an image.  Bayes rule can then be used to combine these priors with tissue type probabilities derived from voxel intensities, to provide the posterior probability.'
    ''
    'This procedure was inherently circular, because the registration required an initial tissue classification, and the tissue classification requires an initial registration.  This circularity is resolved here by combining both components into a single generative model. This model also includes parameters that account for image intensity non-uniformity. Estimating the model parameters (for a maximum a posteriori solution) involves alternating among classification, bias correction and registration steps. This approach provides better results than simple serial applications of each component.'
    }';
preproc.prog = @spm_local_preproc_run;
preproc.vout = @vout;

%==========================================================================
function varargout = spm_local_preproc_run(job)
varargout{1} = spm_preproc_run(job);

%==========================================================================
function dep = vout(job)
% This depends on job contents, which may not be present when virtual
% outputs are calculated.

cdep = cfg_dep;
cdep(end).sname      = 'Seg Params';
cdep(end).src_output = substruct('.','param','()',{':'});
cdep(end).tgt_spec   = cfg_findspec({{'filter','mat','strtype','e'}});

for i=1:numel(job.channel)
    if job.channel(i).write(1)
        cdep(end+1)          = cfg_dep;
        cdep(end).sname      = sprintf('Bias Field (%d)',i);
        cdep(end).src_output = substruct('.','channel','()',{i},'.','biasfield','()',{':'});
        cdep(end).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});
    end
    if job.channel(i).write(2)
        cdep(end+1)          = cfg_dep;
        cdep(end).sname      = sprintf('Bias Corrected (%d)',i);
        cdep(end).src_output = substruct('.','channel','()',{i},'.','biascorr','()',{':'});
        cdep(end).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});
    end
end

for i=1:numel(job.tissue)
    if job.tissue(i).native(1)
        cdep(end+1)          = cfg_dep;
        cdep(end).sname      = sprintf('c%d Images',i);
        cdep(end).src_output = substruct('.','tiss','()',{i},'.','c','()',{':'});
        cdep(end).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});
    end
    if job.tissue(i).native(2)
        cdep(end+1)          = cfg_dep;
        cdep(end).sname      = sprintf('rc%d Images',i);
        cdep(end).src_output = substruct('.','tiss','()',{i},'.','rc','()',{':'});
        cdep(end).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});
    end
    if job.tissue(i).warped(1)
        cdep(end+1)          = cfg_dep;
        cdep(end).sname      = sprintf('wc%d Images',i);
        cdep(end).src_output = substruct('.','tiss','()',{i},'.','wc','()',{':'});
        cdep(end).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});
    end
    if job.tissue(i).warped(2)
        cdep(end+1)          = cfg_dep;
        cdep(end).sname      = sprintf('mwc%d Images',i);
        cdep(end).src_output = substruct('.','tiss','()',{i},'.','mwc','()',{':'});
        cdep(end).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});
    end
end

if job.warp.write(1)
    cdep(end+1)          = cfg_dep;
    cdep(end).sname      = 'Inverse Deformations';
    cdep(end).src_output = substruct('.','invdef','()',{':'});
    cdep(end).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});
end

if job.warp.write(2)
    cdep(end+1)          = cfg_dep;
    cdep(end).sname      = 'Forward Deformations';
    cdep(end).src_output = substruct('.','fordef','()',{':'});
    cdep(end).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});
end

dep = cdep;

author: spm


text:
function [hs,probs] = fernsClfApply( data, ferns, inds )
% Apply learned fern classifier.
%
% USAGE
%  [hs,probs] = fernsClfApply( data, ferns, [inds] )
%
% INPUTS
%  data     - [NxF] N length F binary feature vectors
%  ferns    - learned fern classification model
%  inds     - [NxM] cached inds (from previous call to fernsInds)
%
% OUTPUTS
%  hs       - [Nx1] predicted output labels
%  probs    - [NxH] predicted output label probabilities
%
% EXAMPLE
%
% See also fernsClfTrain, fernsInds
%
% Piotr's Image&Video Toolbox      Version 2.50
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]
if( nargin<3 || isempty(inds) )
  inds = fernsInds(data,ferns.fids,ferns.thrs); end
[N,M]=size(inds); H=ferns.H; probs=zeros(N,H);
for m=1:M, probs = probs + ferns.pFern(inds(:,m),:,m); end
if(ferns.bayes==0), probs=probs/M; end; [~,hs]=max(probs,[],2);
end

author: zouchuhang

repo: LayoutNet


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function pathS = MC_BS(S0,r,d,T,sigma,NTime,NSim,NBatches)
%
% Path discretization in a Black-Scholes model
%
%   S0 spot price
%   K the strike price
%   r the riskless rate
%   d the dividend yield
%   sigma the volatility
%   T the maturity
%   NTime number of timesteps
%   NSim number of simulations
%   NBatches number of batches

Delta = T/NTime;                          % The discretisation step

pathS = zeros(NSim,NTime+1,NBatches);
lnS1 = zeros(NSim,NTime+1);               % init the logspot price path
lnS1(:,1)=log(S0*exp(-d*T));              % adjust due to dividend yield

omega = (r-d-0.5*sigma^2) * Delta;

for l = 1:NBatches
    dW = randn(NSim,NTime);               % precompute all randoms

    for i=2:NTime+1
        lnS1(:,i) = lnS1(:,i-1) + omega + sigma* sqrt(Delta)*dW(:,i-1);
    end
    pathS(:,:,l) = exp(lnS1);
end

author: Sable

repo: mcbench-benchmarks