text:
function mesh = mshSphere(N,rad)
%+========================================================================+
%|                                                                        |
%|                 OPENMSH - LIBRARY FOR MESH MANAGEMENT                  |
%|           openMsh is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : mshSphere.m                                   |
%|    #    |   VERSION    : 0.40                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.03.2018                                    |
%| ( === ) |   SYNOPSIS   : Build uniform mesh for a sphere               |
%|  `---'  |                                                              |
%+========================================================================+

% Fibonnacci rules
or     = (1+sqrt(5))/2;
theta  = (mod((2*pi/or) .* (0:N-1),2*pi))';
phi    = asin( -1 + 2/(N-1) * (0:N-1))';

% Carthesian coordinates
[x,y,z] = sph2cart(theta,phi,rad);
X       = [0 0 0 ; x y z];

% Delaunay triangulation
DT        = delaunayTriangulation(X);
[elt,vtx] = freeBoundary(DT);

% Mesh
mesh = msh(vtx,elt);
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function jacobi_rule ( order, alpha, beta, a, b, filename )

%*****************************************************************************80
%
%% JACOBI_RULE generates a Gauss-Jacobi rule.
%
%  Discussion:
%
%    This program computes a standard Gauss-Jacobi quadrature rule
%    and writes it to a file.
%
%    The user specifies:
%    * the ORDER (number of points) in the rule;
%    * ALPHA, the exponent of (1-x);
%    * BETA, the exponent of (1+x);
%    * A, the left endpoint;
%    * B, the right endpoint;
%    * FILENAME, the root name of the output files.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2010
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'JACOBI_RULE\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compute a Gauss-Jacobi rule for approximating\n' );
  fprintf ( 1, '    Integral ( A <= x <= B ) (B-x)^alpha (x-A)^beta f(x) dx\n' );
  fprintf ( 1, '  of order ORDER.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The user specifies ORDER, ALPHA, BETA and FILENAME.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  ORDER is the number of points\n' );
  fprintf ( 1, '  ALPHA is the exponent of (B-x):\n' );
  fprintf ( 1, '  BETA is the exponent of (x-A);\n' );
  fprintf ( 1, '  A is the left endpoint;\n' );
  fprintf ( 1, '  B is the right endpoint;\n' );
  fprintf ( 1, '  FILENAME is used to generate 3 files:\n' );
  fprintf ( 1, '    filename_w.txt - the weight file\n' );
  fprintf ( 1, '    filename_x.txt - the abscissa file.\n' );
  fprintf ( 1, '    filename_r.txt - the region file.\n' );
%
%  Get ORDER.
%
  if ( nargin < 1 )
    order = input ( '  Enter the rule order ORDER:  ' );
  elseif ( ischar ( order ) )
    order = str2num ( order );
  end
%
%  Get ALPHA.
%
  if ( nargin < 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  ALPHA is the exponent of (B-x) in the integral:\n' );
    fprintf ( 1, '  Note that -1.0 < ALPHA is required.\n' );
    alpha = input ( '  Enter the value of ALPHA:  ' );
  elseif ( ischar ( alpha ) )
    alpha = str2num ( alpha );
  end
%
%  Get BETA.
%
  if ( nargin < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  BETA is the exponent of (x-A) in the integral:\n' );
    fprintf ( 1, '  Note that -1.0 < BETA is required.\n' );
    beta = input ( '  Enter the value of BETA:  ' );
  elseif ( ischar ( beta ) )
    beta = str2num ( beta );
  end
%
%  Get A.
%
  if ( nargin < 4 )
    a = input ( '  Enter the left endpoint A:  ' );
  elseif ( ischar ( a ) )
    a = str2num ( a );
  end
%
%  Get B.
%
  if ( nargin < 5 )
    b = input ( '  Enter the right endpoint B:  ' );
  elseif ( ischar ( b ) )
    b = str2num ( b );
  end
%
%  Get FILENAME.
%
  if ( nargin < 6 )
    fprintf ( 1,  '\n' );
    fprintf ( 1,  '  FILENAME is the ''root name'' of the quadrature files).\n' );
    filename = input ( '  Enter the value of FILENAME as a quoted string:  ' );
  end
%
%  Input summary.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  ORDER = %d\n', order );
  fprintf ( 1, '  ALPHA = %f\n', alpha );
  fprintf ( 1, '  BETA = %f\n', beta );
  fprintf ( 1, '  A = %f\n', a );
  fprintf ( 1, '  B = %f\n', b );
  fprintf ( 1, '  FILENAME = "%s".\n', filename );
%
%  Construct the rule.
%
  kind = 4;
  [ x, w ] = cgqf ( order, kind, alpha, beta, a, b );
%
%  Write the rule.
%
  r = [ a, b ]';
  rule_write ( order, filename, x, w, r );
%
%  Terminate.
%
  fprintf ( 1,  '\n' );
  fprintf ( 1,  'JACOBI_RULE:\n' );
  fprintf ( 1,  '  Normal end of execution.\n' );
  fprintf ( 1,  '\n' );
  timestamp ( );

  return
end
function [ t, wts ] = cdgqf ( nt, kind, alpha, beta )

%*****************************************************************************80
%
%% CDGQF computes a Gauss quadrature formula with default A, B and simple knots.
%
%  Discussion:
%
%    This routine computes all the knots and weights of a Gauss quadrature
%    formula with a classical weight function with default values for A and B,
%    and only simple knots.
%
%    There are no moments checks and no printing is done.
%
%    Use routine EIQFS to evaluate a quadrature computed by CGQFS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev,            (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,inf)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-inf,inf)  |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,inf)     (x-a)^alpha*(x+b)^beta
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Output, real T(NT), the knots.
%
%    Output, real WTS(NT), the weights.
%
  parchk ( kind, 2 * nt, alpha, beta );
%
%  Get the Jacobi matrix and zero-th moment.
%
  [ aj, bj, zemu ] = class_matrix ( kind, nt, alpha, beta );
%
%  Compute the knots and weights.
%
  [ t, wts ] = sgqf ( nt, aj, bj, zemu );

  return
end
function [ t, wts ] = cgqf ( nt, kind, alpha, beta, a, b )

%*****************************************************************************80
%
%% CGQF computes knots and weights of a Gauss quadrature formula.
%
%  Discussion:
%
%    The user may specify the interval (A,B).
%
%    Only simple knots are produced.
%
%    The user may request that the routine print the knots and weights,
%    and perform a moment check.
%
%    Use routine EIQFS to evaluate this quadrature formula.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev Type 1,     (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,+oo)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-oo,+oo)   |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,+oo)     (x-a)^alpha*(x+b)^beta
%    9, Chebyshev Type 2,     (a,b)       ((b-x)*(x-a))^(+0.5)
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Input, real A, B, the interval endpoints.
%
%    Output, real T(NT), the knots.
%
%    Output, real WTS(NT), the weights.
%

%
%  Compute the Gauss quadrature formula for default values of A and B.
%
  [ t, wts ] = cdgqf ( nt, kind, alpha, beta );
%
%  All knots have multiplicity = 1.
%
  mlt = zeros(nt,1);
  mlt(1:nt) = 1;
%
%  NDX(I) = I.
%
  ndx = ( 1 : nt );
%
%  Scale the quadrature rule.
%
  [ t, wts ] = scqf ( nt, t, mlt, wts, nt, ndx, kind, alpha, beta, a, b );

  return
end
function [ aj, bj, zemu ] = class_matrix ( kind, m, alpha, beta )

%*****************************************************************************80
%
%% CLASS_MATRIX computes the Jacobi matrix for a quadrature rule.
%
%  Discussion:
%
%    This routine computes the diagonal AJ and subdiagonal BJ
%    elements of the order M tridiagonal symmetric Jacobi matrix
%    associated with the polynomials orthogonal with respect to
%    the weight function specified by KIND.
%
%    For weight functions 1-7, M elements are defined in BJ even
%    though only M-1 are needed.  For weight function 8, BJ(M) is
%    set to zero.
%
%    The zero-th moment of the weight function is returned in ZEMU.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev,            (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,inf)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-inf,inf)  |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,inf)     (x-a)^alpha*(x+b)^beta
%
%    Input, integer M, the order of the Jacobi matrix.
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Output, real AJ(M), BJ(M), the diagonal and subdiagonal
%    of the Jacobi matrix.
%
%    Output, real ZEMU, the zero-th moment.
%
  temp = eps;

  parchk ( kind, 2 * m - 1, alpha, beta );

  temp2 = 0.5;

  if ( 500.0 * temp < abs ( ( gamma ( temp2 ) )^2 - pi ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CLASS - Fatal error!\n' );
    fprintf ( 1, '  Gamma function does not match machine parameters.\n' );
    error ( 'CLASS - Fatal error!' );
  end

  bj = zeros(m,1);
  aj = zeros(m,1);

  if ( kind == 1 )

    ab = 0.0;

    zemu = 2.0 / ( ab + 1.0 );

    aj(1:m) = 0.0;

    for i = 1 : m
      abi = i + ab * mod ( i, 2 );
      abj = 2 * i + ab;
      bj(i) = abi * abi / ( abj * abj - 1.0 );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 2 )

    zemu = pi;

    aj(1:m) = 0.0;

    bj(1) =  sqrt ( 0.5 );
    bj(2:m) = 0.5;

  elseif ( kind == 3 )

    ab = alpha * 2.0;
    zemu = 2.0^( ab + 1.0 ) * gamma ( alpha + 1.0 )^2 ...
      / gamma ( ab + 2.0 );

    aj(1:m) = 0.0;
    bj(1) = 1.0 / ( 2.0 * alpha + 3.0 );
    for i = 2 : m
      bj(i) = i * ( i + ab ) / ( 4.0 * ( i + alpha )^2 - 1.0 );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 4 )

    ab = alpha + beta;
    abi = 2.0 + ab;
    zemu = 2.0^( ab + 1.0 ) * gamma ( alpha + 1.0 ) ...
      * gamma ( beta + 1.0 ) / gamma ( abi );
    aj(1) = ( beta - alpha ) / abi;
    bj(1) = 4.0 * ( 1.0 + alpha ) * ( 1.0 + beta ) ...
      / ( ( abi + 1.0 ) * abi * abi );
    a2b2 = beta * beta - alpha * alpha;

    for i = 2 : m
      abi = 2.0 * i + ab;
      aj(i) = a2b2 / ( ( abi - 2.0 ) * abi );
      abi = abi^2;
      bj(i) = 4.0 * i * ( i + alpha ) * ( i + beta ) * ( i + ab ) ...
        / ( ( abi - 1.0 ) * abi );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 5 )

    zemu = gamma ( alpha + 1.0 );

    for i = 1 : m
      aj(i) = 2.0 * i - 1.0 + alpha;
      bj(i) = i * ( i + alpha );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 6 )

    zemu = gamma ( ( alpha + 1.0 ) / 2.0 );

    aj(1:m) = 0.0;

    for i = 1 : m
      bj(i) = ( i + alpha * mod ( i, 2 ) ) / 2.0;
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 7 )

    ab = alpha;
    zemu = 2.0 / ( ab + 1.0 );

    aj(1:m) = 0.0;

    for i = 1 : m
      abi = i + ab * mod(i,2);
      abj = 2 * i + ab;
      bj(i) = abi * abi / ( abj * abj - 1.0 );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 8 )

    ab = alpha + beta;
    zemu = gamma ( alpha + 1.0 ) * gamma ( - ( ab + 1.0 ) ) ...
      / gamma ( - beta );
    apone = alpha + 1.0;
    aba = ab * apone;
    aj(1) = - apone / ( ab + 2.0 );
    bj(1) = - aj(1) * ( beta + 1.0 ) / ( ab + 2.0 ) / ( ab + 3.0 );
    for i = 2 : m
      abti = ab + 2.0 * i;
      aj(i) = aba + 2.0 * ( ab + i ) * ( i - 1 );
      aj(i) = - aj(i) / abti / ( abti - 2.0 );
    end

    for i = 2 : m - 1
      abti = ab + 2.0 * i;
      bj(i) = i * ( alpha + i ) / ( abti - 1.0 ) * ( beta + i ) ...
        / ( abti^2 ) * ( ab + i ) / ( abti + 1.0 );
    end

    bj(m) = 0.0;
    bj(1:m) =  sqrt ( bj(1:m) );

  end

  return
end
function [ d, z ] = imtqlx ( n, d, e, z )

%*****************************************************************************80
%
%% IMTQLX diagonalizes a symmetric tridiagonal matrix.
%
%  Discussion:
%
%    This routine is a slightly modified version of the EISPACK routine to
%    perform the implicit QL algorithm on a symmetric tridiagonal matrix.
%
%    The authors thank the authors of EISPACK for permission to use this
%    routine.
%
%    It has been modified to produce the product Q' * Z, where Z is an input
%    vector and Q is the orthogonal matrix diagonalizing the input matrix.
%    The changes consist (essentialy) of applying the orthogonal transformations
%    directly to Z as they are generated.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%    Roger Martin, James Wilkinson,
%    The Implicit QL Algorithm,
%    Numerische Mathematik,
%    Volume 12, Number 5, December 1968, pages 377-383.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real D(N), the diagonal entries of the matrix.
%
%    Input, real E(N), the subdiagonal entries of the
%    matrix, in entries E(1) through E(N-1).
%
%    Input, real Z(N), a vector to be operated on.
%
%    Output, real D(N), the diagonal entries of the diagonalized matrix.
%
%    Output, real Z(N), the value of Q' * Z, where Q is the matrix that
%    diagonalizes the input symmetric tridiagonal matrix.
%
  itn = 30;

  prec = eps;

  if ( n == 1 )
    return
  end

  e(n) = 0.0;

  for l = 1 : n

    j = 0;

    while ( 1 )

      for m = l : n

        if ( m == n )
          break
        end

        if ( abs ( e(m) ) <= prec * ( abs ( d(m) ) + abs ( d(m+1) ) ) )
          break
        end

      end

      p = d(l);

      if ( m == l )
        break
      end

      if ( j == itn )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'IMTQLX - Fatal error!\n' );
        fprintf ( 1, '  Iteration limit exceeded.\n' );
        error ( 'IMTQLX - Fatal error!' );
      end

      j = j + 1;
      g = ( d(l+1) - p ) / ( 2.0 * e(l) );
      r =  sqrt ( g * g + 1.0 );
      g = d(m) - p + e(l) / ( g + r8_sign ( g ) * abs ( r ) );
      s = 1.0;
      c = 1.0;
      p = 0.0;
      mml = m - l;

      for ii = 1 : mml

        i = m - ii;
        f = s * e(i);
        b = c * e(i);

        if ( abs ( f ) >= abs ( g ) )
          c = g / f;
          r =  sqrt ( c * c + 1.0 );
          e(i+1) = f * r;
          s = 1.0 / r;
          c = c * s;
        else
          s = f / g;
          r =  sqrt ( s * s + 1.0 );
          e(i+1) = g * r;
          c = 1.0 / r;
          s = s * c;
        end

        g = d(i+1) - p;
        r = ( d(i) - g ) * s + 2.0 * c * b;
        p = s * r;
        d(i+1) = g + p;
        g = c * r - b;
        f = z(i+1);
        z(i+1) = s * z(i) + c * f;
        z(i) = c * z(i) - s * f;

      end

      d(l) = d(l) - p;
      e(l) = g;
      e(m) = 0.0;

    end

  end

  for ii = 2 : n

     i = ii - 1;
     k = i;
     p = d(i);

     for j = ii : n
       if ( d(j) < p )
         k = j;
         p = d(j);
       end
     end

     if ( k ~= i )
       d(k) = d(i);
       d(i) = p;
       p = z(i);
       z(i) = z(k);
       z(k) = p;
     end

  end

  return
end
function parchk ( kind, m, alpha, beta )

%*****************************************************************************80
%
%% PARCHK checks parameters ALPHA and BETA for classical weight functions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev,            (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,inf)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-inf,inf)  |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,inf)     (x-a)^alpha*(x+b)^beta
%
%    Input, integer M, the order of the highest moment to
%    be calculated.  This value is only needed when KIND = 8.
%
%    Input, real ALPHA, BETA, the parameters, if required
%    by the value of KIND.
%
  if ( kind <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PARCHK - Fatal error!\n' );
    fprintf ( 1, '  KIND <= 0.\n' );
    error ( 'PARCHK - Fatal error!' );
  end
%
%  Check ALPHA for Gegenbauer, Jacobi, Laguerre, Hermite, Exponential.
%
  if ( 3 <= kind && alpha <= -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PARCHK - Fatal error!\n' );
    fprintf ( 1, '  3 <= KIND and ALPHA <= -1.\n' );
    error ( 'PARCHK - Fatal error!' );
  end
%
%  Check BETA for Jacobi.
%
  if ( kind == 4 && beta <= -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PARCHK - Fatal error!\n' );
    fprintf ( 1, '  KIND == 4 and BETA <= -1.0.\n' );
    error ( 'PARCHK - Fatal error!' );
  end
%
%  Check ALPHA and BETA for rational.
%
  if ( kind == 8 )
    tmp = alpha + beta + m + 1.0;
    if ( 0.0 <= tmp || tmp <= beta )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PARCHK - Fatal error!\n' );
      fprintf ( 1, '  KIND == 8 but condition on ALPHA and BETA fails.\n' );
      error ( 'PARCHK - Fatal error!' );
    end
  end

  return
end
function value = r8_sign ( x )

%*****************************************************************************80
%
%% R8_SIGN returns the sign of an R8.
%
%  Discussion:
%
%    The value is +1 if the number is positive or zero, and it is -1 otherwise.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the number whose sign is desired.
%
%    Output, real VALUE, the sign of X.
%
  if ( 0 <= x )
    value = +1.0;
  else
    value = -1.0;
  end

  return
end
function r8mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% R8MAT_WRITE writes an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
%  For smaller data files, and less precision, try:
%
%     fprintf ( output_unit, '  %14.6f', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %24.16f', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function rule_write ( order, filename, x, w, r )

%*****************************************************************************80
%
%% RULE_WRITE writes a quadrature rule to a file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ORDER, the order of the rule.
%
%    Input, string FILENAME, specifies the output files.
%    write files 'filename_w.txt', 'filename_x.txt', 'filename_r.txt' defining
%    weights, abscissas, and region.
%
%    Input, real X(ORDER), the abscissas.
%
%    Input, real W(ORDER), the weights.
%
%    Input, real R(2), the region.
%
  filename_x = strcat ( filename, '_x.txt' );
  filename_w = strcat ( filename, '_w.txt' );
  filename_r = strcat ( filename, '_r.txt' );

  fprintf ( 1, '\n' );
  fprintf ( 1,'  Creating quadrature files.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  "Root" file name is   "%s".\n', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Weight file will be   "%s".\n', filename_w );
  fprintf ( 1, '  Abscissa file will be "%s".\n', filename_x );
  fprintf ( 1, '  Region file will be   "%s".\n', filename_r );

  r8mat_write ( filename_w, 1, order, w' );
  r8mat_write ( filename_x, 1, order, x' );
  r8mat_write ( filename_r, 1, 2,     r' );

  return
end
function [ t, wts ] = scqf ( nt, t, mlt, wts, nwts, ndx, kind, alpha, ...
  beta, a, b )

%*****************************************************************************80
%
%% SCQF scales a quadrature formula to a nonstandard interval.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, real T(NT), the original knots.
%
%    Input, integer MLT(NT), the multiplicity of the knots.
%
%    Input, real WTS(NWTS), the weights.
%
%    Input, integer NWTS, the number of weights.
%
%    Input, integer NDX(NT), used to index the array WTS.
%    For more details see the comments in CAWIQ.
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev Type 1,     (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,+oo)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-oo,+oo)   |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,+oo)     (x-a)^alpha*(x+b)^beta
%    9, Chebyshev Type 2,     (a,b)       ((b-x)*(x-a))^(+0.5)
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Input, real A, B, the interval endpoints.
%
%    Output, real T(NT), the scaled knots.
%
%    Output, real WTS(NWTS), the scaled weights.
%
  temp = eps;

  parchk ( kind, 1, alpha, beta )

  if ( kind == 1 )

    al = 0.0;
    be = 0.0;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 2 )

    al = -0.5;
    be = -0.5;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 3 )

    al = alpha;
    be = alpha;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 4 )

    al = alpha;
    be = beta;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 5 )

    if ( b <= 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  B <= 0.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = a;
    slp = 1.0 / b;
    al = alpha;
    be = 0.0;

  elseif ( kind == 6 )

    if ( b <= 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  B <= 0.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = a;
    slp = 1.0 / sqrt ( b );
    al = alpha;
    be = 0.0;

  elseif ( kind == 7 )

    al = alpha;
    be = 0.0;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 8 )

    if ( a + b <= 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  A + B <= 0.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = a;
    slp = a + b;
    al = alpha;
    be = beta;

  elseif ( kind == 9 )

    al = 0.5;
    be = 0.5;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  end

  p = slp^( al + be + 1.0 );

  for k = 1 : nt

    t(k) = shft + slp * t(k);
    l = abs ( ndx(k) );

    if ( l ~= 0 )
      tmp = p;
      for i = l : l + mlt(k) - 1
        wts(i) = wts(i) * tmp;
        tmp = tmp * slp;
      end
    end

  end

  return
end
function [ t, wts ] = sgqf ( nt, aj, bj, zemu )

%*****************************************************************************80
%
%% SGQF computes knots and weights of a Gauss Quadrature formula.
%
%  Discussion:
%
%    This routine computes all the knots and weights of a Gauss quadrature
%    formula with simple knots from the Jacobi matrix and the zero-th
%    moment of the weight function, using the Golub-Welsch technique.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, real AJ(NT), the diagonal of the Jacobi matrix.
%
%    Input, real BJ(NT), the subdiagonal of the Jacobi
%    matrix, in entries 1 through NT-1.  On output, BJ has been overwritten.
%
%    Input, real ZEMU, the zero-th moment of the weight function.
%
%    Output, real T(NT), the knots.
%
%    Output, real WTS(NT), the weights.
%

%
%  Exit if the zero-th moment is not positive.
%
  if ( zemu <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SGQF - Fatal error!\n' );
    fprintf ( 1, '  ZEMU <= 0.\n' );
    error ( 'SGQF - Fatal error!' );
  end
%
%  Set up vectors for IMTQLX.
%
  wts = zeros ( nt, 1 );
  wts(1) = sqrt ( zemu );
  wts(2:nt) = 0.0;
%
%  Diagonalize the Jacobi matrix.
%
  [ t, wts ] = imtqlx ( nt, aj, bj, wts );

  wts(1:nt) = wts(1:nt).^2;

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [wall_d, x_min, f_min] = pano_line_solver_6(lines, w_res, options)

    pk_loc = lines(1:2:end,1);

    % rotate
    line_label = lines(1:2:end,3);
    line_rot = find(line_label == 2);

    % initialize based on convex/concave
    if line_rot == 1

    else if line_rot == 2
            x_ini = [1.5; 0.5; -1; 2; 1; pk_loc; w_res];
        else if line_rot == 3
                x_ini = [0.5; 1.5; 1; 2; 2; pk_loc; w_res];
            else if line_rot == 4
                    x_ini = [0.25; 0.25; 0.5; 0.5; 1; pk_loc; w_res];
                else if line_rot == 5
                        x_ini = [0.75; 0.25; 1; 0.5; 0.5; pk_loc; w_res];
                        else if line_rot == 6
                            end
                    end
                end
            end
        end
    end

    x_min = Inf;
    f_min = Inf;
    % solve
    for t = 1:400
        %while 1
        f_old = Inf;
        for k = 1:5
            [x,f] = minFunc(@sampleEijOpt_6,x_ini,options);
            x_ini = x;
            % iterate until converge
            if abs(f- f_old) < 1e-5
                break
            end
            f_old = f;
        end
        %disp(f)
        if f < f_min
            f_min = f;
            x_min = x;
        end
        if f < 1e-5
            break
        end
        % TODO: fix initialization here
        rng('shuffle');
        if line_rot == 1
        else if line_rot == 2
                xo = rand(3,1);
                x_ini = [xo(1)+1;xo(2);-xo(3);2;1;pk_loc; w_res];
            else if line_rot == 3
                    xo = rand(2,1)*2;
                    x_ini = [xo(1);xo(2);xo(2)/2;2;2;pk_loc; w_res];
                else if line_rot == 4
                        xo = rand(2,1);
                        x_ini = [xo;(xo(2)+1)/2;(1+xo(1))/2;1;pk_loc; w_res];
                    else if line_rot == 5
                            xo = rand(2,1);
                            x_ini = [xo(1);xo(2);1;xo(1)/2;(1+xo(2))/2;pk_loc; w_res];
                        end
                    end
                end
            end
        end
    end
    disp(f_min)
    %keyboard

    % wall d
    cor = [0 0; 1,0; 1, x_min(3); x_min(4), x_min(3);x_min(4), x_min(5); 0, x_min(5)];
    wall_d = (cor - repmat([x_min(1) x_min(2)], 6, 1));
    wall_d = wall_d(:,1).^2 + wall_d(:,2).^2;
    wall_d = sqrt(wall_d);

author: zouchuhang

repo: LayoutNet


text:
function plotValdata(valdata, fignum)

cthresh = (0.001:0.001:0.5);
err = zeros(size(cthresh));
nregions = zeros(size(cthresh));
for f = 1:numel(valdata)
    for k = 1:numel(cthresh)
        [tmp, minind] = min(abs(valdata(f).mergeCost-cthresh(k)));
        err(k) = err(k) + valdata(f).segError(minind)/numel(valdata);
        nregions(k) = nregions(k) + valdata(f).nregions(minind)/numel(valdata);
    end
end
figure(fignum), subplot(2,2,1), hold off, plot(cthresh, err), title('bias vs. err')
figure(fignum), subplot(2,2,2), plot(err, nregions), title('conf: err vs. nregions')

rthresh = (5:5:1000);
err = zeros(size(rthresh));
nregions = zeros(size(rthresh));
for f = 1:numel(valdata)
    for k = 1:numel(rthresh)
        [tmp, minind] = min(abs(valdata(f).nregions-rthresh(k)));
        err(k) = err(k) + valdata(f).segError(minind)/numel(valdata);
        nregions(k) = nregions(k) + valdata(f).nregions(minind)/numel(valdata);
    end
end
figure(fignum), subplot(2,2,3), plot(rthresh, err),  title('minregion vs. err')
figure(fignum), subplot(2,2,4), plot(err, nregions), title('minreg: err vs. nregions')

author: Cloud-CV

repo: object-proposals


text:
function [newphase,usedmask] = comp_filterbankconstphase(abss,tgrad,fgrad,NEIGH,posInfo,fc,mask,usephase,a,M,N,tol,phasetype)

%sets the mask for the heap integration and acts as a wrapper for
%the corresponding C-function

%% The actual phase construction
sMax = max(abss(:));
absthr = max(sMax)*tol;

if isempty(mask)
    usedmask= zeros(size(abss));
else
    usedmask = cell2mat(mask);
end

%phasetype = 0;
% % Build the phase
if isempty(mask)
     % s is real and positive
     newphase = comp_filterbankheapint(abss,tgrad,fgrad,NEIGH,posInfo,fc,a,M,N,tol(1),phasetype);

     % Find all small coefficients and set the mask
     bigenoughidx = abss>absthr(1);
     usedmask(bigenoughidx) = 1;
 else
     newphase=comp_filterbankmaskedheapint(abss,tgrad,fgrad,NEIGH,posInfo,fc,mask,a,M,N,tol(1),...
                                     phasetype,usephase);
     % Find all small coefficients in the unknown phase area
     missingidx = find(usedmask==0);
     bigenoughidx = abss(missingidx)>absthr(1);
     usedmask(missingidx(bigenoughidx)) = 1;
end
%
% % Do further tol
for ii=2:numel(tol)
     newphase=comp_filterbankmaskedheapint(abss,tgrad,fgrad,NEIGH,posInfo,fc,usedmask,a,M,N,tol(ii),...
                                     phasetype,newphase);
     missingidx = find(usedmask==0);
     bigenoughidx = abss(missingidx)>absthr(ii);
     usedmask(missingidx(bigenoughidx)) = 1;
end

% Convert the mask so it can be used directly for indexing
usedmask = logical(usedmask);
% Assign random values to coefficients below tolerance
zerono = numel(find(~usedmask));
newphase(~usedmask) = rand(zerono,1)*2*pi;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% This script performs our evaluation on simulated datasets obtained from
% our ground truth data.

%% Update MATLAB path.
addpath(genpath('../utility'));
addpath(genpath('../algorithms/OpticalFlow'));
addpath(genpath('../algorithms/SRAlgorithms'));

%% Evaluation settings.

% Path to input data and result directory.
inputDir = '../../data';
datasetConfig = 'baselineDatasets';
resultDir = '../../results/simulatedDatasets';

% Parameters for this evaluation.
initEvaluationParametersForSimulatedDatasets;

%% Start evaluation.
processDatasets(inputDir, datasetConfig, resultDir, binningFactors, numberOfFrames, startReferenceFrame, numberOfFrames_val, sliding_val, sr_method, binning_val, scenes_val, compressions(compress_val), 1, true);

author: thomas-koehler

repo: SupER


text:
% hack to convert autoalign to manual align matrices
% do a regular save from mralign2 after this...
% AW 99
mf=MfA;
trans=transA';
rot=rotA;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function plot_imu_sta(omega, freq, text_st)
% plot_imu_sta: plots static IMU data.
%
% INPUT
%   omega, IMU sensor data (rad/s or m/s^2).%
%   freq, Nx1 IMU sensor sampling frequency (Hz).
%   text_st, title for the figure (string).
%
% OUTPUT
%   a figure.
%
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   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 Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% References:
%
%   M.A. Hopcroft. Allan overlap MATLAB function v2.24.
% https://www.mathworks.com/matlabcentral/fileexchange/13246-allan
%
% Version: 001
% Date:    2021/12/06
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

% Colors
blue    = [0, 0.4470, 0.7410];
orange  = [0.8500, 0.3250, 0.0980];
green   = [0.4660, 0.6740, 0.1880];
yellow  = [0.9290, 0.6940, 0.1250];
light_blue = [0.3010, 0.7450, 0.9330];

% Text size
font_title  = 50;
font_tick   = 15;
font_label  = 20;
font_legend = 15;

% Line width
lw = 3;

M = max(size(omega)) - 1;
dt = 1/freq;

time = (0:dt:dt*M)';

s.linear = polyfit( time(1:length(omega)), omega, 1);

% scale to median for plotting
omega_median=median(omega);
medianfreq=omega-omega_median;

% Screen for outliers using 5x Median Absolute Deviation (MAD) criteria
MAD = median(abs(medianfreq)/0.6745);

% adjust time to remove any starting offset
dtime = time - time(1) + mean(diff(time));

% plot the frequency data, centered on median
% this should not be necessary, but dsplot 1.1 is a little bit brittle
if size(dtime,2) > size(dtime,1), dtime=dtime'; end

% dsplot makes a new figure
if(is_octave)
  hd=plot(dtime,medianfreq);
else
  hd=dsplot(dtime,medianfreq);
end

  set(hd,'Marker','.','LineStyle','none'); % equivalent to '.-'
hold on;

fx = xlim;
% plot([fx(1) fx(2)],[omega_median omega_median],'-k');
plot([fx(1) fx(2)],[0 0],':k','LineWidth', lw);

% show 5x Median Absolute deviation (MAD) values
hm = plot([fx(1) fx(2)],[5*MAD 5*MAD],'-r','LineWidth', lw);
plot([fx(1) fx(2)],[-5*MAD -5*MAD],'-r','LineWidth', lw);

% show linear fit line
hf=plot(xlim,polyval(s.linear,xlim)-omega_median,'-g','LineWidth', lw);
t1 = title(text_st);

hs=plot(xlim,polyval(s.linear,xlim)-omega_median-3*MAD,'--m','LineWidth', lw);
plot(xlim,polyval(s.linear,xlim)-omega_median+3*MAD,'--m','LineWidth', lw);

%set(get(gca,'Title'),'Interpreter','none');
x1=xlabel('Time [sec]');
y1=ylabel('Samples');

if MAD ~= 0.0
    l1=legend([hd hm hs hf],{'data (centered on median)','5x MAD outliers', ...
        '3x MAD outliers', ['Linear Fit (' num2str(s.linear(1),'%g') ')']});
else
    l1=legend([hd hf],{'data (centered on median)', ...
        ['Linear Fit (' num2str(s.linear(1),'%g') ')']});
end

% tighten up
xlim([dtime(1) dtime(end)]);

set(t1,'FontSize', font_title);
set(x1,'FontSize', font_label);
set(y1,'FontSize', font_label);
set(l1,'FontSize', font_legend);
set(gca, 'YTickMode', 'auto', 'FontSize', font_tick);

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
function x = p04_start ( n )

%*****************************************************************************80
%
%% P04_START returns a starting point for optimization for problem 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables X.
%
%    Output, real X(N), a starting point for the optimization.
%
  x = [ 0.0, 1.0 ]';

  return
end


text:
function [K,w_jk,mu_jk,sigma_jk,A1] = estimate_num_comp(Y, A, ...
    sizes, shrink_size, max_num_comp, options)
%ESTIMATE_NUM_COMP Summary of this function goes here
%   Detailed explanation goes here
if nargin < 6
    options = [];
end

thresh = 0.99;
ignore_small_prior = 0.1;

[N,M] = size(A);
K = ones(1,M);
w_jk = cell(1,M);
mu_jk = cell(1,M);
sigma_jk = cell(1,M);

A1 = trim_abundance(A, sizes, shrink_size, thresh);

for j = 1:M
    X = Y(A1(:,j), :);
    num = CVIC(X, max_num_comp, options);
    K(j) = num;

    while true
        try % fit_GMM has a small probability to fail
            obj = fit_GMM(X, num, []);
            break;
        catch me
            disp('Failed to perform fit_GMM.');
        end
    end
    w_jk{j} = obj.PComponents;
    mu_jk{j} = obj.mu;
    sigma_jk{j} = obj.Sigma;
end

disp(['The number of components after CVIC are ',num2str(K)]);

%% remove the components that have too smaller priors
for j = 1:M
    X = Y(A1(:,j), :);
    ignore_masks = w_jk{j} < ignore_small_prior;
    w_jk{j}(ignore_masks) = [];
    mu_jk{j}(ignore_masks,:) = [];
    sigma_jk{j}(:,:,ignore_masks) = [];

    K(j) = length(w_jk{j});
    obj = struct('PComponents',w_jk{j},'mu',mu_jk{j},'Sigma',sigma_jk{j});
    obj = fit_GMM(X, K(j), obj);
    w_jk{j} = obj.PComponents;
    mu_jk{j} = obj.mu;
    sigma_jk{j} = obj.Sigma;
end

disp(['The number of components after removing small priors are ',num2str(K)]);
disp(['Total number of combinations is ',num2str(prod(K))]);

function num = CVIC(X, max_num_comp, options)
[N,B] = size(X);
if N/2 < B % number of samples is too small to allow 2 components
    num = 1;
    return;
end

warning('off','stats:gmdistribution:FailedToConverge');

% p = randperm(N);
% X = X(p,:);
likelihoods = zeros(N,1);
total_lh = zeros(max_num_comp,1);

num_fold = 5;
% set_size = round(N / num_fold);

for K = 1:max_num_comp
    if (N*(num_fold-1)/num_fold) / K < B % number of samples is too small
        total_lh(K) = -Inf;
        continue;
    end
    for i = 1:num_fold
        X1 = X;
%         indices = (i-1)*set_size+1 : min(i*set_size,N);
        indices = i:num_fold:N;
        test_X = X1(indices,:);
        X1(indices,:) = [];
        try
            obj = fit_GMM(X1, K, []);
            value = calc_log_gmm(test_X, obj.PComponents, obj.mu, obj.Sigma);
            likelihoods(indices) = value;
        catch me
            likelihoods(indices) = -Inf;
%             disp(['Error in CVIC: ',me.message]);
        end
    end
    total_lh(K) = sum(likelihoods);
end

[max_lh,num] = max(total_lh);
thresh_CVIC = parse_param(options,'thresh_CVIC',[]);
if ~isempty(thresh_CVIC)
    num = find(abs(total_lh - max_lh) <= max_lh * thresh_CVIC, 1);
end

warning('on','stats:gmdistribution:FailedToConverge');

function obj = fit_GMM(X, num, s)
if num > 1
    if isempty(s)
%         s = kmeans(X, num, 'start', 'cluster');
        s = kmeans(X, num);
    end
    options = statset('Display','off');
    obj = gmdistribution.fit(X, num, 'Start', s, 'Options', options);
else % one Gaussian
    [N,B] = size(X);
    obj = [];
    obj.PComponents = 1;
    obj.mu = mean(X,1);
    X1 = X - repmat(obj.mu, N, 1);
    obj.Sigma = (1/N)*(X1'*X1); % same as in EM (not 1/(N-1))
    if N < B % if the number of samples is too small
        obj.Sigma = obj.Sigma + 1e-9*eye(B);
    end
end

function A_new = trim_abundance(A, sizes, shrink_size, thresh)
if shrink_size == 0
    A_new = A > thresh;
    return
end

[N,M] = size(A);
A_new = false(N,M);
for j = 1:M
    A1 = reshape(A(:,j)>thresh, sizes);
    se = strel('disk', shrink_size);
    erodedBW = imerode(A1, se);

    shrink_size_r = shrink_size;
    while all(erodedBW(:) == 0) % reduce size of structure element
        shrink_size_r = shrink_size_r - 1;
        se_r = strel('disk', shrink_size_r);
        erodedBW = imerode(A1, se_r);
    end

    A_new(:,j) = erodedBW(:);
end

function A_new = trim_abundance1(A, sizes, pure_perc, thresh)
if pure_perc == 1
    A_new = A > thresh;
    return
end

[N,M] = size(A);
A_new = false(N,M);
disk_sizes = (1:100);
for j = 1:M
    A1 = reshape(A(:,j)>thresh, sizes);
    for k = 1:length(disk_sizes)
        se = strel('disk', disk_sizes(k));
        erodedBW = imerode(A1, se);
        if length(find(erodedBW)) / length(find(A1)) <= pure_perc
            A_new(:,j) = erodedBW(:);
            break;
        end
    end
end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function colors = prtPlotUtilDarkenColors(colors)
% Internal function,
% xxx Need Help xxx

colors = colors - 0.1;
colors(colors < 0) = 0;

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
% mcmcTrainScript
% 5-fold cross-validation training

LOAD = 0;

ncv = 5;

imdir = '/IUS/vmr7/dhoiem/context/images/all_images';
outdir = '/IUS/vmr7/dhoiem/context/data';

if LOAD
    load([outdir '/rand_indices.mat']);
    load([outdir '/allimsegs2.mat']);
    load([outdir '/mcmcSuperpixelData.mat']);
    load([outdir '/mcmcSuperpixelClassifier.mat']);
    load([outdir '/mcmcEdgeData.mat']);
    load([outdir '/mcmcEdgeClassifier.mat']);
    load([outdir '/mcmcSegmentationClassifierCV.mat']);
    load([outdir '/mcmcGtSegmentData2.mat']);
    load([outdir '/mcmcSegmentClassifierCV2.mat']);
end

nimages = numel(imsegs);

if ~exist('spfeatures')
    spfeatures = mcmcGetAllSuperpixelData(imdir, imsegs);
    save([outdir '/mcmcSuperpixelData.mat'], 'spfeatures');
end

if ~exist('vclassifierSP')
    [vclassifierSP, hclassifierSP] = ...
        mcmcTrainSuperpixelClassifier(spfeatures(cluster_images), imsegs(cluster_images));
    save([outdir '/mcmcSuperpixelClassifier.mat'], 'vclassifierSP', 'hclassifierSP');
end

if ~exist('efeatures')
    [efeatures, adjlist] = mcmcGetAllEdgeData(spfeatures, imsegs);
    save([outdir '/mcmcEdgeData.mat'], 'efeatures', 'adjlist');
end

if ~exist('eclassifier')
    eclassifier = mcmcTrainEdgeClassifier(efeatures(cluster_images), ...
        adjlist(cluster_images), imsegs(cluster_images));
    save([outdir '/mcmcEdgeClassifier.mat'], 'eclassifier');
end

% if ~exist('gtmaps')
%     [gtmaps, gtlabels] = mcmcGetGroundTruthSegmentations(imsegs, adjlist);
%     gtsegfeatures = mcmcGetAllSegmentFeatures(imsegs, imdir, gtmaps, spfeatures, vclassifierSP, hclassifierSP);
%     save([outdir '/mcmcGtSegmentPerfectData.mat'], 'gtmaps', 'gtlabels', 'gtsegfeatures');
% end
%
% if ~exist('vclassifier')
%     for k = 1:ncv
%         trainind = setdiff([1:numel(imsegs)], ...
%             cv_images(floor((k-1)*numel(cv_images)/ncv)+1:floor(k*numel(cv_images)/ncv)));
%         [vclassifier(k), hclassifier(k)] = ...
%             mcmcTrainSegmentClassifier(imsegs(trainind), gtsegfeatures(trainind), gtlabels(trainind));
%     end
%     save([outdir '/mcmcSegmentClassifierPerfectCV.mat'], 'vclassifier', 'hclassifier');
% end

if ~exist('segclassifier')
    [segdata, seglabels, segimind] = mcmcGenerateRandomSegments3(imsegs(cv_images), imdir,  ...
        adjlist(cv_images), spfeatures(cv_images), efeatures(cv_images), vclassifierSP, hclassifierSP, ...
        eclassifier);

    for k = 1:ncv
        trainind = setdiff([1:numel(cv_images)], ...
            [floor((k-1)*numel(cv_images)/ncv)+1:floor(k*numel(cv_images)/ncv)]);
        segclassifier(k) = mcmcTrainSegmentationClassifier(segdata, seglabels, segimind, trainind);
    end
    save([outdir '/mcmcSegmentationClassifierCV.mat'], 'segclassifier', 'segdata', 'seglabels', 'segimind');
end

if ~exist('gtfeatures')
    [gtfeatures, gtlabels, gtlabelprobs, gtweights] = mcmcGetAllGoodSegments(...
        imsegs(cv_images), imdir, vclassifierSP, hclassifierSP, eclassifier, ...
        segclassifier, spfeatures(cv_images), adjlist(cv_images), efeatures(cv_images), ncv);
    save([outdir '/mcmcGtSegmentData2.mat'], 'gtfeatures', 'gtlabels', 'gtlabelprobs', 'gtweights');
end

if ~exist('vclassifier')
    for k = 1:ncv
        trainind = setdiff([1:numel(cv_images)], ...
            cv_images(floor((k-1)*numel(cv_images)/ncv)+1:floor(k*numel(cv_images)/ncv)));
        [vclassifier(k), hclassifier(k)] = ...
            mcmcTrainSegmentClassifier2(gtfeatures(trainind), gtlabels(trainind), gtweights(trainind));
    end
    save([outdir '/mcmcSegmentClassifierCV2.mat'], 'vclassifier', 'hclassifier');
end

% if ~exist('vregressor')
%     for k = 1:ncv
%         trainind = setdiff([1:numel(cv_images)], ...
%             cv_images(floor((k-1)*numel(cv_images)/ncv)+1:floor(k*numel(cv_images)/ncv)));
%         [vregressor(k), hregressor(k)] = ...
%             mcmcTrainSegmentRegressor(gtfeatures(trainind), gtlabelprobs(trainind), gtweights(trainind));
%     end
%     save([outdir '/mcmcSegmentRegressorCV.mat'], 'vregressor', 'hregressor');
% end


text:
src = gtzan_src('/path/to/gtzan');

N = 5*2^17;

filt1_opt.Q = 8;
filt1_opt.J = 80;
filt1_opt.gabor = 1;

filt2_opt.Q = 1;
filt2_opt.J = 13;

filters1 = morlet_filter_bank_1d(N,filt1_opt);
filters2 = morlet_filter_bank_1d(N,filt2_opt);

scatt_fun = @(x)(scatt_1d(x,{@(x)(wavemod_1d(x,filters1)),@(x)(wavemod_1d(x,filters2)),@(x)(wavemod_1d(x,filters2))}));
scatt_fun = @(x)(format_scatt(log_scatt(renorm_scatt(scatt_fun(x))),'table'));

db = prepare_database(src,{scatt_fun});

[train_set,test_set] = create_partition(src);

model = affine_train(db,train_set);
labels = affine_test(db,model,test_set);

err = classif_err(labels,test_set,src);

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function stringOut=upperFirstLowerRest(stringIn)

stringOut=lower(stringIn);
if ~isempty(stringOut)
  stringOut(1)=upper(stringIn(1));
end

end

author: kristinbranson

repo: JAABA


text:
function [Dt, totalDist] = InterfiberZhangDistance(curve1, curve2, Tt)
%Distance between two 3D paths, defined as two 3xN sets of 3d coordinates.
%
% [Dt, totalDist] = InterfiberZhangDistance(curve1, curve2, [Tt=2])
%
% The terms Dt and totalDist need better definition here.
%
% The metric is described in Zhang et al. (2003) IEEE Transactions on
% Vizualization and Computer Graphics 9(4).  They write:
%
% �In order to emphasize important differences between a pair of
% trajectories, we average the distance between the curves only over the
% region where they are at least Tt apart; smaller differences are assumed
% to be insignificant�.
%
% Although arguable, it is natural to set Tt to data voxel size /not sure,
% before or after resampling/ (as Zhang 2003 did). Here the default value
% is set to 1mm despite the fact that out data are commonly 2x2x2.
%
% Calculated on a point by point basis.
% Note: by definition, if a point on a shorter fiber is closer to the other
% fiber than Tt, this point does NOT contribute to the overall
% curve-to-curve distance measure.
%
%Dt is above-the-threshold (that is, minus Tt) average point-to-curve
%distance  across all points on a shortest fiber to another fiber.
%
%  CLARIFY.
% The output variable totalDist is Dt+Tt, the average point-to-curve
% distance between the two fibers. These output arguments are different by
% a constant,  however included since the first one is the "definition" of
% the distance by Zhang et al., and the second one is easier to interpret.
% Also totalDist is reported as at least at Tt (even for two fibers which
% are closer than Tt).
%
%  Example:
%    curve1 = 10*rand(3,10);
%    curve2 = curve1 + 1;
%    Tt = 0;
%    [Dt, totalDist] = InterfiberZhangDistance(curve1, curve2,Tt)
%    [Dt, totalDist] = InterfiberZhangDistance(curve1, curve2,2)
%
%ER 11/2008
%ER 04/2009 added "totalDist" output and introduced a default value for Tt
%to be 0 (all points contribute to the distance measure).
%
% (c) Stanford VISTA Team

%Pairs of points in closer than Tt are not included when summing the
%distance between fibers.
if (~exist('Tt', 'var')), Tt=0; end

% Check curve sizes.
if size(curve1, 1) ~= 3 || size(curve2, 1) ~= 3
   error('Fiber coordinates must be represented as 3xN');
end

% Use  "nearpoints" to find the nearest point in curve 1 to curve 2.
% [indices, bestSqDist] = nearpoints(src, dest)
%-  src is a 3xM array of points
%-  dest is a 3xN array of points
%-  indices is a 1xM vector, whose elements identifies the closest point
%   in dest for each entry in src.
%   "bestSqDist" is the distance.
%

% For efficiency, choose which is the first term
if size(curve1, 2) < size(curve2, 2),
    [indices, bestSqDist] = nearpoints(curve1, curve2);
else
    [indices, bestSqDist] = nearpoints(curve2, curve1);
end

% Have a look at the data in the two curves.
% mrvNewGraphWin;
%  plot3(curve1(1,:),curve1(2,:),curve1(3,:),'bo'); hold on
%  plot3(curve2(1,:),curve2(2,:),curve2(3,:),'r.');
%

% I don't understand.  If there are no entries bigger than Tt, we set Dt to
% zero, I guess.
nPoints = length(find(bestSqDist > Tt));
if nPoints <= 0
    % Expansive but useless. Precisely,
    % Tt=sum(sqrt(bestSqDist))/length(bestSqDist);
    % We keep it as a constant Tt.
    Dt = 0;
else
    % Otherwise, we set Dt to this average.
    Dt = sum(sqrt(bestSqDist(bestSqDist > Tt)) - Tt) / nPoints;
end

% Then we add Tt back in.
totalDist = Dt + Tt;

return;

%%The slow way (gives equivalent result)
if size(curve1, 2)<size(curve2, 2)
    scurve=curve1;  %shorter
    lcurve=curve2; %longer curve
else
    scurve=curve2;
    lcurve=curve1;
end

countpoints=0;  sumdist=0;

for pnt=1:size(scurve,2)
    distVector = sqrt((lcurve(1, :)-scurve(1, pnt)).^2 + (lcurve(2, :)-scurve(2, pnt)).^2+(lcurve(3, :)-scurve(3, pnt)).^2);
    p2curveDist=min(distVector); %Shortest point to curve distance

    if p2curveDist>Tt
        countpoints=countpoints+1;
        sumdist= sumdist+ (p2curveDist-Tt);
    end

end

if countpoints>0
    Dt=sumdist/countpoints;
else
    Dt=0;
end


text:
function varargout = cstructdecode(buf, varargin);

% CSTRUCTDECODE decodes a structure from a uint8 buffer
%
% See READ_NEURALYNX_NEV for an example

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

if ~isa(buf, 'uint8')
  error('incorrect type of input data, should be uint8');
end

nbytes = numel(buf);
nfield = length(varargin);

wordsize = zeros(1,nfield);
for i=1:nfield
  switch varargin{i}
  case 'uint8'
    wordsize(i) = 1;
  case 'int8'
    wordsize(i) = 1;
  case 'uint16'
    wordsize(i) = 2;
  case 'int16'
    wordsize(i) = 2;
  case 'uint32'
    wordsize(i) = 4;
  case 'int32'
    wordsize(i) = 4;
  case 'uint64'
    wordsize(i) = 8;
  case 'int64'
    wordsize(i) = 8;
  case {'float32' 'single'}
    varargin{i} = 'single';
    wordsize(i) = 4;
  case {'float64' 'double'}
    varargin{i} = 'double';
    wordsize(i) = 8;
  otherwise
    if strncmp(varargin{i}, 'char', 4)
      if length(varargin{i})>4
        % assume a string like 'char128' which means 128 characters
        wordsize(i) = str2num(varargin{i}(5:end));
        varargin{i} = 'char';
      else
        wordsize(i) = 1;
      end
    else
      error('incorrect type specification');
    end
  end
end

pklen = sum(wordsize);
pknum = nbytes/sum(wordsize);

buf = reshape(buf, pklen, pknum);

for i=1:nfield
  rowbeg = sum(wordsize(1:(i-1)))+1;
  rowend = sum(wordsize(1:(i-0)))+0;
  sel = buf(rowbeg:rowend,:);
  if strcmp(varargin{i}, 'char')
    varargout{i} = char(sel)';
  else
    varargout{i} = typecast(sel(:), varargin{i});
  end
end

author: PatternRecognition

repo: OpenBMI


text:
function z = isequal( x, y )
z = x.id_ == y.id_;

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function [cfg, artifact] = ft_artifact_eog(cfg, data)

% FT_ARTIFACT_EOG reads the data segments of interest from file and identifies EOG
% artifacts.
%
% Use as
%   [cfg, artifact] = ft_artifact_eog(cfg)
% with the configuration options
%   cfg.dataset     = string with the filename
% or
%   cfg.headerfile  = string with the filename
%   cfg.datafile    = string with the filename
% and optionally
%   cfg.headerformat
%   cfg.dataformat
%
% Alternatively you can use it as
%   [cfg, artifact] = ft_artifact_eog(cfg, data)
% where the input data is a structure as obtained from FT_PREPROCESSING.
%
% In both cases the configuration should also contain
%   cfg.trl        = structure that defines the data segments of interest. See FT_DEFINETRIAL
%   cfg.continuous = 'yes' or 'no' whether the file contains continuous data
%
% The data is preprocessed (again) with the following configuration parameters,
% which are optimal for identifying EOG artifacts.
%   cfg.artfctdef.eog.bpfilter   = 'yes'
%   cfg.artfctdef.eog.bpfilttype = 'but'
%   cfg.artfctdef.eog.bpfreq     = [1 15]
%   cfg.artfctdef.eog.bpfiltord  = 4
%   cfg.artfctdef.eog.hilbert    = 'yes'
%
% Artifacts are identified by means of thresholding the z-transformed value
% of the preprocessed data.
%   cfg.artfctdef.eog.channel      = Nx1 cell-array with selection of channels, see FT_CHANNELSELECTION for details
%   cfg.artfctdef.eog.cutoff       = z-value at which to threshold (default = 4)
%   cfg.artfctdef.eog.trlpadding   = 0.5
%   cfg.artfctdef.eog.fltpadding   = 0.1
%   cfg.artfctdef.eog.artpadding   = 0.1
%
% The output argument "artifact" is a Nx2 matrix comparable to the
% "trl" matrix of FT_DEFINETRIAL. The first column of which specifying the
% beginsamples of an artifact period, the second column contains the
% endsamples of the artifactperiods.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
% If you specify this option the input data will be read from a *.mat
% file on disk. This mat files should contain only a single variable named 'data',
% corresponding to the input structure.
%
% See also FT_REJECTARTIFACT, FT_ARTIFACT_CLIP, FT_ARTIFACT_ECG, FT_ARTIFACT_EOG,
% FT_ARTIFACT_JUMP, FT_ARTIFACT_MUSCLE, FT_ARTIFACT_THRESHOLD, FT_ARTIFACT_ZVALUE

% Undocumented local options
% cfg.method

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

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init

% ft_preamble provenance is not needed because just a call to ft_artifact_zvalue
% ft_preamble loadvar data is not needed because ft_artifact_zvalue will do this

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamed',    {'datatype', 'continuous'});
cfg = ft_checkconfig(cfg, 'renamedval', {'continuous', 'continuous', 'yes'});

% set default rejection parameters
cfg.artfctdef            = ft_getopt(cfg,               'artfctdef', []);
cfg.artfctdef.eog        = ft_getopt(cfg.artfctdef,     'eog',       []);
cfg.artfctdef.eog.method = ft_getopt(cfg.artfctdef.eog, 'method',    'zvalue');

if isfield(cfg.artfctdef.eog, 'artifact')
  fprintf('eog artifact detection has already been done, retaining artifacts\n');
  artifact = cfg.artfctdef.eog.artifact;
  return
end

if ~strcmp(cfg.artfctdef.eog.method, 'zvalue')
  ft_error('EOG artifact detection only works with cfg.method=''zvalue''');
end

% for backward compatibility
cfg.artfctdef.eog = ft_checkconfig(cfg.artfctdef.eog, 'renamed', {'sgn',     'channel'});
cfg.artfctdef.eog = ft_checkconfig(cfg.artfctdef.eog, 'renamed', {'passbnd', 'bpfreq'});
cfg.artfctdef.eog = ft_checkconfig(cfg.artfctdef.eog, 'renamed', {'padding', 'trlpadding'});
artpadding_oldstyle  = max(ft_getopt(cfg.artfctdef.eog, 'pretim', 0), ft_getopt(cfg.artfctdef.eog, 'psttim', 0));

% settings for preprocessing
cfg.artfctdef.eog.bpfilter    = ft_getopt(cfg.artfctdef.eog, 'bpfilter',   'yes');
cfg.artfctdef.eog.bpfreq      = ft_getopt(cfg.artfctdef.eog, 'bpfreq',     [1 15]);
cfg.artfctdef.eog.bpfiltord   = ft_getopt(cfg.artfctdef.eog, 'bpfiltor',   4);
cfg.artfctdef.eog.bpfilttype  = ft_getopt(cfg.artfctdef.eog, 'bpfilttype', 'but');
cfg.artfctdef.eog.hilbert     = ft_getopt(cfg.artfctdef.eog, 'hilbert',    'yes');

% settings for the zvalue subfunction
cfg.artfctdef.eog.channel    = ft_getopt(cfg.artfctdef.eog, 'channel',    'EOG');
cfg.artfctdef.eog.trlpadding = ft_getopt(cfg.artfctdef.eog, 'trlpadding', 0.5);
cfg.artfctdef.eog.fltpadding = ft_getopt(cfg.artfctdef.eog, 'fltpadding', 0.1);
cfg.artfctdef.eog.artpadding = max(ft_getopt(cfg.artfctdef.eog, 'artpadding', 0.1), artpadding_oldstyle);
cfg.artfctdef.eog.cutoff     = ft_getopt(cfg.artfctdef.eog, 'cutoff',     4);

% construct a temporary configuration that can be passed onto artifact_zvalue
tmpcfg                  = cfg;
tmpcfg.artfctdef.zvalue = cfg.artfctdef.eog;
tmpcfg.artfctdef        = rmfield(tmpcfg.artfctdef, 'eog');

% call the zvalue artifact detection function, where the data is either passed
% into the function by the user or read from file with cfg.inputfile
hasdata = exist('data', 'var');
if ~hasdata
  tmpcfg             = ft_checkconfig(tmpcfg, 'dataset2files', 'yes');
  tmpcfg             = ft_checkconfig(tmpcfg, 'required', {'headerfile', 'datafile'});
  [tmpcfg, artifact] = ft_artifact_zvalue(tmpcfg);
else
  tmpcfg.artfctdef.zvalue.trlpadding = 0;
  tmpcfg.artfctdef.zvalue.fltpadding = 0;
  ft_warning('trlpadding and fltpadding are set to zero to avoid filter problems with NaN, see bug3193 for details');
  [tmpcfg, artifact] = ft_artifact_zvalue(tmpcfg, data);
end
cfg.artfctdef.eog  = tmpcfg.artfctdef.zvalue;

author: spm


text:
function []=panel5disp(n,N,m,p,k,T,d1,d2,d3,d4,d5,Ymat,Xdot,Units,endo,exo,const,Xi,theta_gibbs,theta_median,theta_std,theta_lbound,theta_ubound,sigma_gibbs,sigma_median,D_estimates,gamma_estimates,alpha0,delta0,startdate,enddate,forecast_record,forecast_estimates,Fcperiods,stringdates3,Fstartdate,Fcenddate,Feval,Fcomp,data_endo_c,data_endo_c_lags,data_exo_c,It,Bu,IRF,IRFt,pref,names)

% recover a point estimate (the median) of the VAR coefficients
betatilde=Xi*theta_median;
Btilde=reshape(betatilde,k,N*n);

% check whether the model is stationary
[stationary,eigmodulus]=bear.checkstable(betatilde,N*n,p,k);

% estimate the in-sample evaluation criteria

% obtain a point estimate thetatilde of the structural factors, which is the median
thetatilde=theta_median;
% compute fitted values
Ytilde=full(Xdot*kron(speye(T),thetatilde))';
% reshape for convenience
Ytilde=reshape(Ytilde,T,n,N);
Ymat=reshape(Ymat,T,n,N);

% loop over units
for ii=1:N
% estimate the residuals for this unit
EPS(:,:,ii)=Ymat(:,:,ii)-Ytilde(:,:,ii);

% Compute then the sum of squared residuals
% compute first the RSS matrix, defined in (1.9.5)
RSS(:,:,ii)=EPS(:,:,ii)'*EPS(:,:,ii);
% retain only the diagonal elements to get the vector of RSSi values
rss(:,:,ii)=diag(RSS(:,:,ii));

% Go on calculating R2
% generate Mbar
Mbar=eye(T)-ones(T,T)/T;
% then compute the TSS matrix, defined in (1.9.8)
TSS(:,:,ii)=Ymat(:,:,ii)'*Mbar*Ymat(:,:,ii);
% generate the R2 matrix in (1.9.9)
R2(:,:,ii)=eye(n)-RSS(:,:,ii)./TSS(:,:,ii);
% retain only the diagonal elements to get the vector of R2 values
r2(:,:,ii)=diag(R2(:,:,ii));

% then calculate the adjusted R2, using (1.9.11)
R2bar(:,:,ii)=eye(n)-((T-1)/(T-k))*(eye(n)-R2(:,:,ii));
% retain only the diagonal elements to get the vector of R2bar values
r2bar(:,:,ii)=diag(R2bar(:,:,ii));

end

% estimate the forecast evaluation criteria
% first note that forecast evaluation can only be conducted if it is activated, and if there is some observable data after the beginning of the forecast
if Feval==1 && Fcomp==1

% generate the elements required for the evaluation
beta_gibbs=Xi*theta_gibbs;
forecast_record=reshape(forecast_record,N*n,1);
forecast_estimates=reshape(forecast_estimates,N*n,1);
data_endo_c=reshape(data_endo_c,Fcperiods,N*n);
data_endo_c_lags=reshape(data_endo_c_lags,p,N*n);

% compute forecast evaluation
[RMSE,MAE,MAPE,Ustat,CRPS_estimates,S1_estimates,S2_estimates]=bear.panelfeval(N*n,p,k,beta_gibbs,sigma_gibbs,forecast_record,forecast_estimates,Fcperiods,data_endo_c,data_endo_c_lags,data_exo_c,const,It,Bu);

end

% start displaying and saving the general results

% preliminary task: create and open the txt file used to save the results

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

% print toolbox header

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

% print the list of contributors
bear.printcontributors(fid);

% print then estimation results

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

toolboxinfo='BEAR toolbox estimates';
fprintf('%s\n',toolboxinfo);
fprintf(fid,'%s\n',toolboxinfo);

time=clock;
datestring=datestr(time);
dateinfo=['Date: ' datestring(1,1:11) '   Time: ' datestring(1,13:17)];
fprintf('%s\n',dateinfo);
fprintf(fid,'%s\n',dateinfo);

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

VARtypeinfo='Panel VAR: structural factor (static)';
fprintf('%s\n',VARtypeinfo);
fprintf(fid,'%s\n',VARtypeinfo);

if IRFt==1
SVARinfo='structural decomposition: none';
fprintf('%s\n',SVARinfo);
fprintf(fid,'%s\n',SVARinfo);
elseif IRFt==2
SVARinfo='structural decomposition: choleski factorisation';
fprintf('%s\n',SVARinfo);
fprintf(fid,'%s\n',SVARinfo);
elseif IRFt==3
SVARinfo='structural decomposition: triangular factorisation';
fprintf('%s\n',SVARinfo);
fprintf(fid,'%s\n',SVARinfo);
end

temp='units: ';
for ii=1:N
temp=[temp ' ' Units{ii,1} ' '];
end
unitinfo=temp;
fprintf('%s\n',unitinfo);
fprintf(fid,'%s\n',unitinfo);

temp='endogenous variables: ';
for ii=1:n
temp=[temp ' ' endo{ii,1} ' '];
end
endoinfo=temp;
fprintf('%s\n',endoinfo);
fprintf(fid,'%s\n',endoinfo);

temp='exogenous variables: ';
if const==0 && m==0
temp=[temp ' none'];
elseif const==1 && m==1
temp=[temp ' constant '];
elseif const==0 && m>0
   for ii=1:m-1
   temp=[temp ' ' exo{ii,1} ' '];
   end
elseif const==1 && m>1
temp=[temp ' constant '];
   for ii=1:m-1
   temp=[temp ' ' exo{ii,1} ' '];
   end
end
exoinfo=temp;
fprintf('%s\n',exoinfo);
fprintf(fid,'%s\n',exoinfo);

sampledateinfo=['estimation sample: ' startdate '-' enddate];
fprintf('%s\n',sampledateinfo);
fprintf(fid,'%s\n',sampledateinfo);

samplelengthinfo=['sample size (omitting initial conditions): ' num2str(T)];
fprintf('%s\n',samplelengthinfo);
fprintf(fid,'%s\n',samplelengthinfo);

laginfo=['number of lags included in regression: ' num2str(p)];
fprintf('%s\n',laginfo);
fprintf(fid,'%s\n',laginfo);

hyperparam1='hyperparameters:';
fprintf('%s\n',hyperparam1);
fprintf(fid,'%s\n',hyperparam1);

hyperparam2=['IG shape on residual variance (alpha0):       ' num2str(alpha0)];
fprintf('%s\n',hyperparam2);
fprintf(fid,'%s\n',hyperparam2);

hyperparam3=['IG scale on residual variance (delta0):       ' num2str(delta0)];
fprintf('%s\n',hyperparam3);
fprintf(fid,'%s\n',hyperparam3);

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','');

% display factor estimates

factorinfo=['Structural factors:'];
fprintf('%s\n',factorinfo);
fprintf(fid,'%s\n',factorinfo);

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

factorheader=fprintf('%35s %15s %15s %15s %15s\n','','Median','St.dev','Low.bound','Upp.bound');
factorheader=fprintf(fid,'%35s %15s %15s %15s %15s\n','','Median','St.dev','Low.bound','Upp.bound');

% common component

fprintf('%s\n','theta1 (common component)');
fprintf(fid,'%s\n','theta1 (common component)');
values=[theta_median(1,1) theta_std(1,1) theta_lbound(1,1) theta_ubound(1,1)];
fprintf('%-35s %15.3f %15.3f %15.3f %15.3f\n','common component',values);
fprintf(fid,'%-35s %15.3f %15.3f %15.3f %15.3f\n','common component',values);
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% unit component

fprintf('%s\n','theta2 (unit-specific component)');
fprintf(fid,'%s\n','theta2 (unit component)');
for ii=1:d2
values=[theta_median(d1+ii,1) theta_std(d1+ii,1) theta_lbound(d1+ii,1) theta_ubound(d1+ii,1)];
fprintf('%-35s %15.3f %15.3f %15.3f %15.3f\n',['unit ' int2str(ii) ' component'],values);
fprintf(fid,'%-35s %15.3f %15.3f %15.3f %15.3f\n',['unit ' int2str(ii) ' component'],values);
end
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% variable component

fprintf('%s\n','theta3 (variable-specific component)');
fprintf(fid,'%s\n','theta3 (endogenous variable component)');
for ii=1:d3
values=[theta_median(d1+d2+ii,1) theta_std(d1+d2+ii,1) theta_lbound(d1+d2+ii,1) theta_ubound(d1+d2+ii,1)];
fprintf('%-35s %15.3f %15.3f %15.3f %15.3f\n',['variable ' int2str(ii) ' component'],values);
fprintf(fid,'%-35s %15.3f %15.3f %15.3f %15.3f\n',['variable ' int2str(ii) ' component'],values);
end
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% lag component (if applicable)

if d4~=0
fprintf('%s\n','theta4 (lag-specific component)');
fprintf(fid,'%s\n','theta4 (lag component)');
for ii=1:d4
values=[theta_median(d1+d2+d3+ii,1) theta_std(d1+d2+d3+ii,1) theta_lbound(d1+d2+d3+ii,1) theta_ubound(d1+d2+d3+ii,1)];
fprintf('%-35s %15.3f %15.3f %15.3f %15.3f\n',['lag ' int2str(ii) ' component'],values);
fprintf(fid,'%-35s %15.3f %15.3f %15.3f %15.3f\n',['lag ' int2str(ii) ' component'],values);
end
fprintf('%s\n','');
fprintf(fid,'%s\n','');
end

% exogenous component (if applicable)

if d5~=0
fprintf('%s\n','theta5 (exogenous variable component)');
fprintf(fid,'%s\n','theta5 (exogenous component)');
% initiate equation count
eqcount=0;
% initiate exogenous count
exocount=0;
for ii=1:d5
   if exocount==m
   exocount=0;
   end
   if exocount==0
   eqcount=eqcount+1;
   end
exocount=exocount+1;
values=[theta_median(d1+d2+d3+d4+ii,1) theta_std(d1+d2+d3+d4+ii,1) theta_lbound(d1+d2+d3+d4+ii,1) theta_ubound(d1+d2+d3+d4+ii,1)];
fprintf('%-35s %15.3f %15.3f %15.3f %15.3f\n',['equation ' int2str(eqcount) ', exogenous ' int2str(exocount) ' component'],values);
fprintf(fid,'%-35s %15.3f %15.3f %15.3f %15.3f\n',['equation ' int2str(eqcount) ', exogenous ' int2str(exocount) ' component'],values);
end
fprintf('%s\n','');
fprintf(fid,'%s\n','');
end

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

% display VAR stability results
eigmodulus=reshape(eigmodulus,p,N*n);
stabilityinfo1=['Roots of the characteristic polynomial (modulus):'];
fprintf('%s\n',stabilityinfo1);
fprintf(fid,'%s\n',stabilityinfo1);
for jj=1:p
temp=num2str(eigmodulus(jj,1),'%.3f');
   for kk=2:N*n
   temp=[temp,'  ',num2str(eigmodulus(jj,kk),'%.3f')];
   end
fprintf('%s\n',temp);
fprintf(fid,'%s\n',temp);
end
if stationary==1;
stabilityinfo2=['No root lies outside the unit circle.'];
stabilityinfo3=['The estimated VAR model satisfies the stability condition'];
fprintf('%s\n',stabilityinfo2);
fprintf(fid,'%s\n',stabilityinfo2);
fprintf('%s\n',stabilityinfo3);
fprintf(fid,'%s\n',stabilityinfo3);
else
stabilityinfo2=['Warning: at leat one root lies on or outside the unit circle.'];
stabilityinfo3=['The estimated VAR model will not be stable'];
end

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

% display posterior for sigma
% reshape sigma
sigma_median=reshape(sigma_median,N*n,N*n);
% start displaying
sigmainfo=['sigma (residual covariance matrix): posterior estimates'];
fprintf('%s\n',sigmainfo);
fprintf(fid,'%s\n',sigmainfo);
% calculate the (integer) length of the largest number in sigma, for formatting purpose
width=length(sprintf('%d',floor(max(abs(bear.vec(sigma_median))))));
% add a separator, a potential minus sign, and three digits (total=5) to obtain the total space for each entry in the matrix
width=width+5;
for ii=1:N*n
temp=[];
   for jj=1:N*n
   % convert matrix entry into string
   number=num2str(sigma_median(ii,jj),'% .3f');
      % pad potential missing blanks
      while numel(number)<width
      number=[' ' number];
      end
   number=[number '  '];
   temp=[temp number];
   end
fprintf('%s\n',temp);
fprintf(fid,'%s\n',temp);
end

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

% then display the results for D and gamma, if a structural decomposition was selected

if IRF==1 && IRFt~=1
svarinfo1=['D (structural decomposition matrix): posterior estimates'];
fprintf('%s\n',svarinfo1);
fprintf(fid,'%s\n',svarinfo1);

% recover D
D=reshape(D_estimates,N*n,N*n);
% calculate the (integer) length of the largest number in D, for formatting purpose
width=length(sprintf('%d',floor(max(abs(bear.vec(D))))));
% add a separator, a potential minus sign and three digits (total=5) to obtain the total space for each entry in the matrix
width=width+5;
for ii=1:N*n
temp=[];
   for jj=1:N*n
   % convert matrix entry into string
   number=num2str(D(ii,jj),'% .3f');
      % pad potential missing blanks
      while numel(number)<width
      number=[' ' number];
      end
   number=[number '  '];
   temp=[temp number];
   end
fprintf('%s\n',temp);
fprintf(fid,'%s\n',temp);
end

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

svarinfo2=['gamma (structural disturbances covariance matrix): posterior estimates'];
fprintf('%s\n',svarinfo2);
fprintf(fid,'%s\n',svarinfo2);

% recover gamma
gamma=reshape(gamma_estimates,N*n,N*n);
% calculate the (integer) length of the largest number in D, for formatting purpose
width=length(sprintf('%d',floor(max(abs(bear.vec(D))))));
% add a separator, a potential minus sign and three digits (total=5) to obtain the total space for each entry in the matrix
width=width+5;
for ii=1:N*n
temp=[];
   for jj=1:N*n
   % convert matrix entry into string
   number=num2str(gamma(ii,jj),'% .3f');
      % pad potential missing blanks
      while numel(number)<width
      number=[' ' number];
      end
   number=[number '  '];
   temp=[temp number];
   end
fprintf('%s\n',temp);
fprintf(fid,'%s\n',temp);
end

end

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

% finally display in-sample and forecast evaluation measures
% if forecast evaluation is activated and possible, display the results
if Feval==1 && Fcomp==1
bear.panel5fprint(Units,N,n,endo,rss,r2,r2bar,RMSE,MAE,MAPE,Ustat,CRPS_estimates,S1_estimates,S2_estimates,stringdates3,Fstartdate,Fcenddate,Fcperiods,fid);

% if forecast evaluation is activated but not possible, return a message to signal it
elseif Feval==1 && 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 not activated altogether, do not do anything

end

author: european-central-bank

repo: BEAR-toolbox


text:
% This demo exemplifies the post-dictive interval on a simple GLM
% Let us assume that observed data are given by: Y=aX+b+e, where e are iid
% residuals. We want to portrait our uncertainty regarding "postdicted"
% data, i.e. draw something like a confidence interval around the
% regression line.

clear all
close all

% 0- simulate data under GLM
n = 16;
X = randn(n,1);
a = 4;
b = 1;
e = 1e0*randn(n,1);
y = a*X + b + e;

% 1- fit model (here, using vaggue priors)
options.inG.X = [X,ones(n,1)];
options.priors.SigmaPhi = 1e4*eye(2);
dim.n = 0;
dim.n_phi = 2;
dim.n_theta = 0;
[post,out] = VBA_NLStateSpaceModel(y,[],[],@g_GLM,dim,options);

% 2- get Laplace postdictive density
options.priors = post;
options.priors = rmfield(options.priors,'iQy');
sx = std(X);
X0 = VBA_vec(min(X)-3*sx:1e-1:max(X)+3*sx); % extend the postdiction outside domain of fitted data
dim.p = length(X0);
dim.n_t = 1;
options.inG.X = [X0,ones(dim.p,1)];
[muy,Vy] = VBA_getLaplace([],[],@g_GLM,dim,options);

% 3- plot regression line

% Get classical CI from sampling the posterior density
N = 1e4;
sV = VBA_sqrtm (post.SigmaPhi);
phi = repmat(post.muPhi,1,N) + sV*randn(2,N);
ev = post.b_sigma./post.a_sigma;
E = 0;% sqrt(ev)*randn(length(X0),N); % add in predicted residuals?
Y = [X0,ones(dim.p,1)]*phi + E;
%uY = prctile(Y,95,2);
uY = arrayfun(@(j) VBA_quantile(Y(j,:),.95), 1:size(Y,1))' ;
%lY = prctile(Y,5,2);
lY = arrayfun(@(j) VBA_quantile(Y(j,:),.05), 1:size(Y,1))' ;

vY = var(Y,[],2);

hf = figure('color',[1 1 1]);
ha = subplot(2,1,1,'parent',hf,'nextplot','add');
plotUncertainTimeSeries(muy',diag(Vy)',X0',ha)
plot(ha,X,y,'k.')
plot(ha,X0,uY,'r--')
plot(ha,X0,lY,'r--')
legend(ha,{'postdictive mean','posdictive STD','data','5% and 95% CI on the mean'})
xlabel(ha,'X')
ylabel(ha,'Y')
title(ha,'regression line')
ha = subplot(2,1,2,'parent',hf,'nextplot','add');
plot(ha,X0,sqrt(diag(Vy)),'k')
plot(ha,X0,sqrt(vY),'r--')
xlabel(ha,'X')
ylabel(ha,'STD[Y|y,m]')
title(ha,'postdictive standard deviation')
set(ha,'xlim',[min(X0),max(X0)])

VBA_getSubplots ();

author: MBB-team

repo: VBA-toolbox


text:
function [varargout]=icolorbar(varargin)

% function [h]=icolorbar(cLim,h)
%-------------------------------------------------------------------------
%
%
% Change log:
% 2018/05/10 Fixed bug in relation to older matlab versions which do not
% support h.Colormap=c for axis specific color maps. Older versions use
% colormap(h,c) instead.
% 2019/08/14 Added colorbar handle output
%-------------------------------------------------------------------------

%% parse input

switch nargin
    case 0
        h=[];
        cLim=caxis;
    case 1
        cLim=varargin{1};
        h=[];
    case 2
        cLim=varargin{1};
        h=varargin{2};
end

if isempty(cLim)
    cLim=caxis;
end

if isempty(h)
    %Get current figure or open a new one
    if isempty(findobj('type','figure'))
        hf=cFigure; %Open a new figure
    else
        hf=gcf; %Get the current figure
    end
   figure(hf);
   if isempty(hf.CurrentAxes) %isempty(get(hf,'CurrentAxes'))
       h=axes;
   else
       h=gca;
   end
else
    switch class(h)
        case {'matlab.ui.Figure'} %A figure handle
            hf=h;
            if isempty(hf.CurrentAxes) %isempty(get(hf,'CurrentAxes'))
                h=axes;
            else
                h=gca;
            end
        case {'matlab.graphics.axis.Axes'} %An axes handle
            hf=h.Parent; %Get figure handle
        otherwise
            error('Handle does not represent an axis or figure handle');
    end
end

%% Get colormap
if isprop(h,'Colormap') %Available in MATLAB 2018
    c=h.Colormap;
    matlabNew=1;
else %Older MATLAB versions
    c=colormap(h);
    matlabNew=0;
end

%% Set limits and tick labels
caxis([cLim(1)-0.5 cLim(2)+0.5]);
hc=colorbar;
hc.Ticks=cLim(1):1:cLim(2);
cn=resampleColormap(c,numel(hc.Ticks));
if matlabNew==1 %axis colormap
    h.Colormap=cn;
else %Figure colormap
    colormap(h,cn);
end

%% Collect output
switch nargout
    case 1
        varargout{1}=h;
    case 2
        varargout{1}=h;
        varargout{2}=hc;
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) 2019  Kevin Mattheus Moerman
%
% 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: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
%CORRECTMATCHES  Refines coordinates of corresponding points
%
%     [newPoints1, newPoints2] = cv.correctMatches(F, points1, points2)
%
% ## Input
% * __F__ 3x3 fundamental matrix.
% * __points1__ first set of 2D points. A numeric Nx2/Nx1x2/1xNx2 array or a
%   cell array of 2-element vectors `{[x,y], ...}` (floating-point precision).
% * __points2__ second set of 2D points. Same size and type as `points1`.
%
% ## Output
% * __newPoints1__ The optimized `points1`. Similar in shape to `points1`
%   (either Nx2/1xNx2 numeric array or cell array of 2D points).
% * __newPoints2__ The optimized `points2`.
%
% The function implements the Optimal Triangulation Method (see [Hartley2004]
% for details). For each given point correspondence `points1[i] <-> points2[i]`,
% and a fundamental matrix `F`, it computes the corrected correspondences
% `newPoints1[i] <-> newPoints2[i]` that minimize the geometric error:
%
%     d(points1[i], newPoints1[i])^2 + d(points2[i], newPoints2[i])^2
%
% (where `d(a,b)` is the geometric distance between points `a` and `b`)
% subject to the epipolar constraint
%
%     newPoints2' * F * newPoints1 = 0
%
% ## References
% [Hartley2004]:
% > R.I. Hartley and A. Zisserman. "Multiple View Geometry in Computer Vision"
% > Cambridge University Press, 2004.
%
% See also: cv.findFundamentalMat
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function pass = test_feval( pref )

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = 1e4*pref.techPrefs.chebfuneps;

%% Evaluate at spherical coordinates

% Example 1
f = ballfun(@(r,lam,th)r.*cos(lam).*sin(th),'spherical');
F = feval(f,[0.5,0.7], [0,0], [pi/2,pi/2], 'spherical');
exact = [0.5, 0.7];
pass(1) = norm(F(:)-exact(:)) < tol;

% Example 2
f = ballfun(@(r,lam,th)r.*cos(lam).*sin(th),'spherical');
F = feval(f,[1,1], [pi/4,pi/3], [pi/2,pi/2], 'spherical');
exact = [cos(pi/4); cos(pi/3)];
pass(2) = norm(F(:)-exact(:)) < tol;

% Example 3
f = ballfun(@(r,lam,th)r.*cos(lam).*sin(th),'spherical');
F = feval(f,[1,1], [0,0], [pi/5, pi/7], 'spherical');
exact = zeros(1,1,2);
exact(1,1,1) = sin(pi/5);
exact(1,1,2) = sin(pi/7);
pass(3) = norm(F(:)-exact(:)) < tol;

% Example 4
S = [22,23,25];
f = ballfun(@(x,y,z)x.*y);
r = chebpts(S(1));
lam = pi*trigpts(S(2));
th = pi*trigpts(S(3));
exact = ballfun.coeffs2vals(coeffs3(f,S(1),S(2),S(3)));
F = fevalm(f,r,lam,th);
pass(4) = norm(F(:)-exact(:)) < tol;

%% Evaluate at cartesian coordinates

% Example 5
f = ballfun(@(x,y,z)x);
F = f(1,0,0);
pass(5) = abs(1-F) < tol;

% Example 6
f = ballfun(@(x,y,z)y);
F = f(0,0.5,0);
pass(6) = abs(0.5-F) < tol;

% Example 7
f = ballfun(@(x,y,z)z);
F = f(0,0,-0.3);
pass(7) = abs(-0.3-F) < tol;

%% EXTRACT_SPHEREFUN EXAMPLES:

% Example 8
f = ballfun(@(r,lam,th)r.*cos(lam).*sin(th),'spherical');
g = f(1,:,:,'spherical');
h = spherefun(@(lam,th)cos(lam).*sin(th));
pass(8) = norm( g - h ) < tol;

% Example 9
f = ballfun(@(r,lam,th)r.*cos(th),'spherical');
g = f(.5,:,:,'spherical');
h = spherefun(@(lam,th)0.5.*cos(th));
pass(9) = norm( g - h ) < tol;

%% COMPLEX EVALUATION
f = ballfun(@(x,y,z) cos(x)*1i);
pass(10) = abs(f(0,0,0)-1i);

if (nargout > 0)
    pass = all(pass(:));
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function A = sup(A)
%SUP          Supremum of long (only with error term)
%
%  C = sup(A)
%
%For A being a long scalar or column vector with error term (see long),
%  output C is the long upper bound
%

% written  11/06/99     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 09/20/04     S.M. Rump  E.sign fixed
% modified 11/19/04     S.M. Rump  exponent update for error
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 07/25/05     S.M. Rump  result improved, thanks to
%                                      Nozomu Matsuda and Nobito Yamamoto
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 08/26/12     S.M. Rump  global variables removed
%

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

  INTLAB_LONG_ERROR = getappdata(0,'INTLAB_LONG_ERROR');

  if ~INTLAB_LONG_ERROR
    error('sup called for long without error term')
  end

  n = length(A.sign);
  E = long(A.error.mant);
  E.exponent = E.exponent + A.error.exp;

  A.error.mant = 0;
  A.error.exp = 0;

  A = A + E;
  index = ( A.error.mant~=0 );
  if any(index)
    %VVVV  A(index) = sup(A(index));
    s.type = '()'; s.subs = {index}; A = subsasgn(A,s,sup(subsref(A,s)));
    %AAAA  Matlab bug fix
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function res = imQuantile(img, q, varargin)
%IMQUANTILE Computes value that threshold a given proportion of pixels
%
%   RES = imQuantile(IMG, Q)
%
%   Example
%     % compute the median value in an imag by two different ways
%     img = imread('cameraman.tif');
%     imQuantile(img, .5)
%     ans =
%         144
%     imMedian(img)
%     ans =
%         144
%
%   See also
%     imHistogram, imMedian

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2014-05-26,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2014 INRA - Cepia Software Platform.

if q < 0 || q > 1
    error('Requires a quantile value between 0 and 1');
end

% compute cumulative histogram
[h, x] = imHistogram(img);
cs = cumsum(h);

% iterate for all quantiles to compute
res = zeros(size(q));
for i = 1:length(q)
    % find index of quantile
    qi = sum(h) * q(i);
    ind = find(cs > qi, 1, 'first');

    % coonvert to grey level
    res(i) = x(ind);
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage