text:
function ising_2d_simulation ( m, n, iterations, thresh, seed )

%*****************************************************************************80
%
%% ISING_2D_SIMULATION carries out a 2D Ising simulation.
%
%  Discussion:
%
%    Note that, when all the cells are updated in a single cycle, there is
%    a mathematically stable checkerboard solution, in which all the reds
%    and blacks flip color repeatedly.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 November 2011
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ISING_2D_SIMULATION\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Monte Carlo simulation of a 2D Ising model.\n' );
%
%  Arguments not supplied are set to default values.
%
  if ( nargin < 1 )
    m = 10;
  else
    m = str2num ( m );
  end

  if ( nargin < 2 )
    n = 10;
  else
    n = str2num ( n );
  end

  if ( nargin < 3 )
    iterations = 15;
  else
    iterations = str2num ( iterations );
  end

  if ( nargin < 4 )
    thresh = 0.50;
  else
    thresh = str2num ( thresh );
  end

  if ( nargin < 5 )
    seed = 123456789;
  else
    seed = str2num ( seed );
  end
%
%  Define the probability of flipping if you are in a neighborhood of
%  1, 2, 3, 4, or 5 of the same sign.
%
%  This should really be some kind of exponential dependent on "temperature".
%
  prob = [ 0.98, 0.85, 0.50, 0.15, 0.02 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The row dimension M = %d\n', m );
  fprintf ( 1, '  The column dimension N = %d\n', n );
  fprintf ( 1, '  The number of iterations taken is ITERATIONS = %d\n', iterations );
  fprintf ( 1, '  The threshhold THRESH = %f\n', thresh );
  fprintf ( 1, '  The seed SEED = %d\n', seed );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The transition probability table, based on the number of\n' );
  fprintf ( 1, '  neighbors with the same charge.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '        1           2           3           4           5\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  %10.4f  %10.4f  %10.4f  %10.4f  %10.4f\n', prob(1:5) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  An automatic pause of 2 seconds is imposed between \n;' );
  fprintf ( 1, '  successive displays.\n' );
%
%  Initialize C1.
%
  [ c1, seed ] = ising_2d_initialize ( m, n, thresh, seed );
%
%  Do the simulation.
%
  [ c1, seed ] = transition ( m, n, iterations, prob, seed, c1 );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ISING_2D_SIMULATION\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function i4mat_print ( m, n, a, title )

%*****************************************************************************80
%
%% I4MAT_PRINT prints an I4MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 September 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of the matrix.
%
%    Input, integer A(M,N), an M by N matrix to be printed.
%
%    Input, string TITLE, a title.
%
  i4mat_print_some ( m, n, a, 1, 1, m, n, title );

  return
end
function i4mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% I4MAT_PRINT_SOME prints out a portion of an I4MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of the matrix.
%
%    Input, integer A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, a title.
%
  incx = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );

  if ( m <= 0 || n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (None)\n' );
    return
  end

  for j2lo = max ( jlo, 1 ): incx : min ( jhi, n )

    j2hi = j2lo + incx - 1;
    j2hi = min ( j2hi, n );
    j2hi = min ( j2hi, jhi );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col: ' );

    for j = j2lo : j2hi
      fprintf ( 1, '%7d  ', j );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row\n' );
    fprintf ( 1, '\n' );

    i2lo = max ( ilo, 1 );
    i2hi = min ( ihi, m );

    for i = i2lo : i2hi

      fprintf ( 1, '%5d: ', i );
      for j = j2lo : j2hi
        fprintf ( 1, '%7d  ', a(i,j) );
      end

      fprintf ( 1, '\n' );

    end

  end

  return
end
function [ c, seed ] = i4mat_uniform ( m, n, a, b, seed )

%*****************************************************************************80
%
%% I4MAT_UNIFORM returns a scaled pseudorandom I4MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Pierre L'Ecuyer,
%    Random Number Generation,
%    in Handbook of Simulation,
%    edited by Jerry Banks,
%    Wiley Interscience, page 95, 1998.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input, integer M, N, the row and column dimensions of the matrix.
%
%    Input, integer A, B, the minimum and maximum acceptable values.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, integer C(M,N), the randomly chosen integer vector.
%
%    Output, integer SEED, the updated seed.
%
  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_UNIFORM - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'I4MAT_UNIFORM - Fatal error!' );
  end

  seed = floor ( seed );
  a = round ( a );
  b = round ( b );

  for j = 1 : n

    for i = 1 : m

      seed = mod ( seed, 2147483647 );

      if ( seed < 0 )
        seed = seed + 2147483647;
      end

      k = floor ( seed / 127773 );

      seed = 16807 * ( seed - k * 127773 ) - k * 2836;

      if ( seed < 0 )
        seed = seed + 2147483647;
      end

      r = seed * 4.656612875E-10;
%
%  Scale R to lie between A-0.5 and B+0.5.
%
      r = ( 1.0 - r ) * ( min ( a, b ) - 0.5 ) ...
        +         r   * ( max ( a, b ) + 0.5 );
%
%  Use rounding to convert R to an integer between A and B.
%
      value = round ( r );

      value = max ( value, min ( a, b ) );
      value = min ( value, max ( a, b ) );

      c(i,j) = value;

    end
  end

  return
end
function c5 = ising_2d_agree ( m, n, c1 )

%*****************************************************************************80
%
%% ISING_2D_AGREE returns the number of neighbors agreeing with each cell.
%
%  Discussion:
%
%    The count includes the cell itself, so it is between 1 and 5.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of cells in each spatial dimension.
%
%    Input, integer C1(M,N), an array of 1's and -1's.
%
%    Output, integer C5(M,N), the number of neighbors that agree.
%    1, 2, 3, 4, or 5.
%
  c5 = c1 ...
     + circshift ( c1, [ -1,  0 ] ) ...
     + circshift ( c1, [ +1,  0 ] ) ...
     + circshift ( c1, [  0, -1 ] ) ...
     + circshift ( c1, [  0, +1 ] );

  i = find ( 0 < c1 );
  c5(i) = ( 5 + c5(i) ) / 2;

  i = find ( c1 < 0 );
  c5(i) = ( 5 - c5(i) ) / 2;

  return
end
function ising_2d_display ( step, m, n, c1 )

%*****************************************************************************80
%
%% ISING_2D_DISPLAY displays the current Ising status.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer STEP, the step number.
%
%    Input, integer M, N, the number of rows and columnss.
%
%    Input, integer C1(M,N), the status of each cell:
%    -1, display as a shade of red.
%    +1, display as a shade of blue.
%
  figure ( 1 )
%
%  Clear the graphics frame.
%
  clf
%
%  Determine the plot range.
%
  margin = 0.05;

  x_axes_min = 1.0 - 0.5 - margin;
  x_axes_max = m   + 0.5 + margin;
  y_axes_min = 1.0 - 0.5 - margin;
  y_axes_max = n   + 0.5 + margin;
%
%  Fill in the background with black.
%
  x1 = x_axes_min;
  x2 = x_axes_max;
  y1 = y_axes_min;
  y2 = y_axes_max;

  rgb = [ 0.5, 0.5, 0.5 ];

  fill ( [ x1, x2, x2, x1 ], [ y1, y1, y2, y2 ], rgb );
  hold on
%
%  Draw a square, representing the bed,
%  with most of the length and width, centered at (I,J).
%
  for i = 1 : m
    for j = 1 : n

      x1 = j - 0.47;
      x2 = j + 0.47;
      y1 = ( m + 1 - i ) - 0.47;
      y2 = ( m + 1 - i ) + 0.47;

      if ( c1(i,j) == - 1 )
        rgb = [ 1.0, 0.0, 0.0 ];
      elseif ( c1(i,j) == + 1 )
        rgb = [ 0.0, 0.0, 1.0 ];
      end

      fill ( [ x1, x2, x2, x1 ], [ y1, y1, y2, y2 ], rgb );
      hold on

    end
  end
%
%  Make a title.
%
  title_string = sprintf ( 'Ising charges +1/-1 on step %d', step );

  title ( title_string )
%
%  Choose the aspect ratio and other plot details.
%
  axis ( [ x_axes_min, x_axes_max, y_axes_min, y_axes_max] );
  axis equal
  axis tight
  axis off

  hold off

  pause ( 2 )

  return
end
function [ c1, seed ] = ising_2d_initialize ( m, n, thresh, seed )

%*****************************************************************************80
%
%% ISING_2D_INITIALIZE initializes the Ising array.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real THRESH, the threshhold value, between 0 and 1.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, integer C1(M,N), an array of 1's and -1's.
%
%    Output, integer SEED, a seed for the random number generator.
%
  c1 = ones ( m, n );

  [ r, seed ] = r8mat_uniform_01 ( m, n, seed );

  i = find ( r <= thresh );

  c1(i) = -1;

  return
end
function ising_2d_stats ( step, m, n, c1 )

%*****************************************************************************80
%
%% ISING_2D_STATS prints information about the current step.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer STEP, the step number.
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, integer C1(M,N), the current state of the system.
%
  if ( step == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Step       Positives       Negatives\n' );
    fprintf ( 1, '             #    %%          #    %%\n' );
    fprintf ( 1, '\n' );
  end

  pos_count = sum ( sum ( 0 < c1 ) );
  neg_count = m * n - pos_count;
  pos_percent = ( 100 * pos_count ) / ( m * n );
  neg_percent = ( 100 * neg_count ) / ( m * n );

  fprintf ( 1, '  %4d  %6d  %6.2f  %6d  %6.2f\n', ...
    step, pos_count, pos_percent, neg_count, neg_percent );

  return
end
function neighbor_2d_display ( step, m, n, c1, c5 )

%*****************************************************************************80
%
%% NEIGHBOR_2D_DISPLAY displays the current Ising neighbor status.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer STEP, the step number.
%
%    Input, integer M, N, the number of rows and columnss.
%
%    Input, integer C1(M,N), the status of each cell:
%
%    Input, integer C5(M,N), the number of agreeable neighbors.
%
  figure ( 2 )
%
%  Clear the graphics frame.
%
  clf
%
%  Determine the plot range.
%
  margin = 0.05;

  x_axes_min = 1.0 - 0.5 - margin;
  x_axes_max = m   + 0.5 + margin;
  y_axes_min = 1.0 - 0.5 - margin;
  y_axes_max = n   + 0.5 + margin;
%
%  Fill in the background with black.
%
  x1 = x_axes_min;
  x2 = x_axes_max;
  y1 = y_axes_min;
  y2 = y_axes_max;

  rgb = [ 0.5, 0.5, 0.5 ];

  fill ( [ x1, x2, x2, x1 ], [ y1, y1, y2, y2 ], rgb );
  hold on
%
%  Draw a square, representing the bed,
%  with most of the length and width, centered at (I,J).
%
  for i = 1 : m
    for j = 1 : n

      x1 = j - 0.47;
      x2 = j + 0.47;
      y1 = ( m + 1 - i ) - 0.47;
      y2 = ( m + 1 - i ) + 0.47;

      c = c1(i,j) * c5(i,j);

      if ( c == - 5 )
        rgb = [ 1.0, 0.0, 0.0 ];
      elseif ( c == - 4 )
        rgb = [ 1.0, 0.2, 0.2 ];
      elseif ( c == - 3 )
        rgb = [ 1.0, 0.4, 0.4 ];
      elseif ( c == - 2 )
        rgb = [ 1.0, 0.7, 0.7 ];
      elseif ( c == - 1 )
        rgb = [ 1.0, 0.8, 0.8 ];
      elseif ( c == 0 )
        rgb = [ 1.0, 1.0, 1.0 ];
      elseif ( c == + 1 )
        rgb = [ 0.7, 0.7, 1.0 ];
      elseif ( c == + 2 )
        rgb = [ 0.6, 0.6, 1.0 ];
      elseif ( c == + 3 )
        rgb = [ 0.4, 0.4, 1.0 ];
      elseif ( c == + 4 )
        rgb = [ 0.2, 0.2, 1.0 ];
      elseif ( c == + 5 )
        rgb = [ 0.0, 0.0, 1.0 ];
      end

      fill ( [ x1, x2, x2, x1 ], [ y1, y1, y2, y2 ], rgb );
      hold on

    end
  end
%
%  Make a title.
%
  title_string = sprintf ( 'Ising neighborhoods -5 to +5, step %d to step %d', step-1, step );

  title ( title_string )
%
%  Choose the aspect ratio and other plot details.
%
  axis ( [ x_axes_min, x_axes_max, y_axes_min, y_axes_max] );
  axis equal
  axis tight
  axis off

  hold off

  pause ( 2 )

  return
end
function neighbor_2d_stats ( step, m, n, c1, c5 )

%*****************************************************************************80
%
%% NEIGHBOR_2D_STATS prints neighbor statistics about the current step.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer STEP, the step number.
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, integer C1(M,N), the current state of the system.
%
%    Input, integer C5(M,N), the number of agreeable neighbors.
%
  if ( step == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Step     Neighborhood Charge:\n' );
    fprintf ( 1, ...
      '           -5    -4    -3    -2    -1    +1    +2    +3    +4    +5\n' );
    fprintf ( 1, '\n' );
  end

  fprintf ( 1, '  %4d', step );

  for n = -5 : 5
    if ( n ~= 0 )
      c = sum ( sum ( c1 .* c5 == n ) );
      fprintf ( 1, '  %4d', c );
    end
  end
  fprintf ( 1, '\n' );

  return
end
function [ r, seed ] = r8mat_uniform_01 ( m, n, seed )

%*****************************************************************************80
%
%% R8MAT_UNIFORM_01 returns a unit pseudorandom R8MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns in the array.
%
%    Input, integer SEED, the integer "seed" used to generate
%    the output random number.
%
%    Output, real R(M,N), an array of random values between 0 and 1.
%
%    Output, integer SEED, the updated seed.  This would
%    normally be used as the input seed on the next call.
%
  r = zeros ( m, n );

  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_UNIFORM_01 - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'R8MAT_UNIFORM_01 - Fatal error!' );
  end

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

      seed = floor ( seed );

      seed = mod ( seed, 2147483647 );

      if ( seed < 0 )
        seed = seed + 2147483647;
      end

      k = floor ( seed / 127773 );

      seed = 16807 * ( seed - k * 127773 ) - k * 2836;

      if ( seed < 0 )
        seed = seed + 2147483647;
      end

      r(i,j) = seed * 4.656612875E-10;

    end
  end

  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
function [ c1, seed ] = transition ( m, n, iterations, prob, seed, c1 )

%*****************************************************************************80
%
%% TRANSITION carries out a Monte Carlo simulation of a 2D Ising model.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of cells in each spatial
%    dimension.
%
%    Input, integer ITERATIONS, the number of iterations to carry out.
%
%    Input, real PROB(1:5).  PROB(I) represents the probability
%    that the spin of a given cell will be reversed, given that it has I immediate
%    neighbors (including itself) with spin the same as its own.
%
%    Input/output, integer SEED, a seed for the random number generator.
%
%    Input/output, integer C1(M,N).  On input, the current state of the
%    system.  On output, the state of the system after the iterations.
%
  step = 0;
  ising_2d_stats ( step, m, n, c1 );
  ising_2d_display ( step, m, n, c1 );

  for step = 1 : iterations
%
%  C5 contains 1 through 5, the number of cells that agree with the center cell.
%
    c5 = ising_2d_agree ( m, n, c1 );

    if ( 0 )
      neighbor_2d_stats ( step, m, n, c1, c5 );
    end

    if ( 0 )
      neighbor_2d_display ( step, m, n, c1, c5 );
    end
%
%  Determine the chances of flipping cell (I,J).
%
    threshhold = zeros ( m, n );

    for j = 1 : 5
      i = find ( c5 == j );
      threshhold(i) = prob(j);
    end

    [ r, seed ] = r8mat_uniform_01 ( m, n, seed );

    i = find ( r < threshhold );

    c1(i) = - c1(i);

    ising_2d_stats ( step, m, n, c1 );
    ising_2d_display ( step, m, n, c1 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function ptCloudOut = interpolatePtCloud(camK, camKc, ptCloudIn, imROI, imMask, verbose)
%% Interpolate the input sparse point cloud

%% License
% ACADEMIC OR NON-PROFIT ORGANIZATION NONCOMMERCIAL RESEARCH USE ONLY
% Copyright (c) 2018 Bingyao Huang
% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:

% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.

% If you publish results obtained using this software, please cite our paper.

% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
% SOFTWARE.

%% interpolate point cloud
x = ptCloudIn.Location(:,1);
y = ptCloudIn.Location(:,2);
z = ptCloudIn.Location(:,3);

% # of interpolated points in X Y meshgrid
numPts = 100;

% xRange contains lower and higher bounds
xRange = ptCloudIn.XLimits;
xDiff = xRange(1) - xRange(2);
xStep = -xDiff / numPts;

yRange = ptCloudIn.YLimits;
yDiff = yRange(1) - yRange(2);
yStep = -yDiff / numPts;

%% use griddata to interpolate a mesh from given point cloud
[xq, yq] = meshgrid(xRange(1):xStep:xRange(2), yRange(1):yStep:yRange(2));
zq = griddata(x, y, z, xq, yq, 'v4');

% interpolated denser point cloud
pts3d = [xq(:), yq(:), zq(:)];

%% only keep the interpolated points in mask ROI
% project interpolated mesh points to camera image
pts2d = cv.projectPoints(pts3d, zeros(3,1), zeros(3,1), camK, 'DistCoeffs', camKc);
pts2dColor = pts2d;

% remove points that are out of camera's fov
inlierIdx = (pts2d(:,2) > 0) & (pts2d(:,1) > 0) & (pts2d(:,2) < size(imMask,1)) & (pts2d(:,1) < size(imMask,2));
pts2dColor(~inlierIdx,:) = 1;
pts2d = pts2d(inlierIdx,:);
pts3d = pts3d(inlierIdx,:);

% only keep masked 2d, 3d points
pixVal = interp2(im2single(imMask), pts2d(:, 1), pts2d(:, 2));
pts2dInliers = pts2d(pixVal>0,:);
pts3dInliers = pts3d(pixVal>0,:);

% convert to single for color interpolation
imROI = im2single(imROI);

% assign color to pt3d
r = interp2(imROI(:,:,1), pts2dInliers(:, 1), pts2dInliers(:, 2));
g = interp2(imROI(:,:,2), pts2dInliers(:, 1), pts2dInliers(:, 2));
b = interp2(imROI(:,:,3), pts2dInliers(:, 1), pts2dInliers(:, 2));

% create ptCloud with colors
ptCloudOut = pointCloud(pts3dInliers, 'Color', [r,g,b]);

%% Debug info
if(verbose)
    %% plot point cloud
%     figure;
%     pcshow(ptCloudOut, 'VerticalAxis', 'y', 'VerticalAxisDir', 'down', 'MarkerSize', 60);
%     title('Interpolated point cloud with color')
%     daspect([1 1 1]);
%     % view(3);
%     axis vis3d tight;

    %% plot mesh (edges only)
    figure;
    scatter3(x,y,z, 'bo'); hold on
    mesh(xq, yq, zq);hold off
    title('Interpolated mesh, Nodes are shown as blue circles');
    daspect([1 1 1]);
    % view(3);
    axis vis3d tight;

    %% plot surface
    % interpolate color
    r = reshape(interp2(imROI(:,:,1), pts2dColor(:, 1), pts2dColor(:, 2)), [numPts+1, numPts+1]);
    g = reshape(interp2(imROI(:,:,2), pts2dColor(:, 1), pts2dColor(:, 2)), [numPts+1, numPts+1]);
    b = reshape(interp2(imROI(:,:,3), pts2dColor(:, 1), pts2dColor(:, 2)), [numPts+1, numPts+1]);

    figure;
    h = surf(xq, yq, zq, cat(3, r, g, b), 'FaceColor', 'texturemap', 'EdgeColor', 'texturemap', 'FaceLighting', 'gouraud', 'LineStyle', 'none');
    title('Reconstructed mesh')
    h.XData(pixVal==0) = nan;
    h.YData(pixVal==0) = nan;
    h.ZData(pixVal==0) = nan;
    h.CData(pixVal==0) = nan;
    %     h.EdgeColor = [0.3,0.3,0.3];
    daspect([1 1 1]);
    % view(3);
    axis vis3d tight;
    rotate3d on
end
end

author: BingyaoHuang

repo: single-shot-pro-cam-calib


text:
function H=transferfunction(g,L)
%TRANSFERFUNCTION  The transferfunction of a filter
%   Usage:  H=transferfunction(g,L);
%
%   `transferfunction(g,L)` computes the transferfunction of length *L*
%   of the filter defined by *g*.
%
%   See also: pfilt

complainif_argnonotinrange(nargin,2,2,mfilename);

[g,info] = comp_fourierwindow(g,L,upper(mfilename));

H=comp_transferfunction(g,L);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [k, sk] = biasKernCompute(kern, x, x2)

% BIASKERNCOMPUTE Compute the BIAS kernel given the parameters and X.
% FORMAT
% DESC computes the kernel parameters for the bias
% kernel given inputs associated with rows and columns.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : the input matrix associated with the rows of the kernel.
% ARG x2 : the inpute matrix associated with the columns of the kernel.
% RETURN k : the kernel matrix computed at the given points.
%
% FORMAT
% DESC computes the kernel matrix for the bias
% kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : the kernel matrix computed at the given points.
%
% SEEALSO : biasKernParamInit, kernCompute, kernCreate, biasKernDiagCompute
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% KERN

if nargin< 3
  dims = [size(x, 1), size(x, 1)];
else
  dims = [size(x, 1), size(x2, 1)];
end
k = repmat(kern.variance, dims);
if nargout > 1
  sk = ones(dims);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


% function [ opti_cost ] = get_the_optimal( cutCostsFile, cut_num )
function [ opti_cuts ] = get_the_optimal( cutCostsFile, cut_num )

%GET_THE_OPTIMAL Summary of this function goes here
%   Detailed explanation goes here

    file_t = fopen(cutCostsFile);
    cut_costs = fscanf(file_t, '%d');
    fclose(file_t);

    cut_costs_sorted = sort(cut_costs);

    opti_cost = sum(cut_costs_sorted(1:cut_num-1));

    % get the optimal cuts, sequential~!
    cut_costs_sorted_top = cut_costs_sorted(1:cut_num-1);

    for i=1:(cut_num-1)

        cut_value = cut_costs_sorted_top(i);

        opti_cuts(i,:) = find(cut_costs == cut_value);

    end
end


text:
function d =  testing(a,d)

 d=test(a.child,d);
 d.X=mean(d.X)';

 d=set_name(d,[get_name(d) ' -> ' get_name(a)]);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%% Crystal Orientations
%
%%
%

author: mtex-toolbox

repo: mtex

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


text:
function A = makejcf(n, e, m, X)
%MAKEJCF   A matrix with specified Jordan canonical form.
%          MAKEJCF(N, E, M) is a matrix having the Jordan canonical form
%          whose i'th Jordan block is of dimension M(i) with eigenvalue E(i),
%          and where N = SUM(M).
%          Defaults: E = 1:N, M = ONES(SIZE(E)) with M(1) so that SUM(M) = N.
%          The matrix is constructed by applying a random similarity
%          transformation to the Jordan form.
%          Alternatively, the matrix used in the similarity transformation
%          can be specified as a fifth parameter.
%          In particular, MAKEJCF(N, E, M, EYE(N)) returns the Jordan form
%          itself.
%          NB: The JCF is very sensitive to rounding errors.

if nargin < 2, e = 1:n; end
if nargin < 3, m = ones(size(e)); m(1) = m(1) + n - sum(m); end

if length(e) ~= length(m)
   error('Parameters E and M must be of same dimension.')
end

if sum(m) ~= n, error('Block dimensions must add up to N.'), end

A = zeros(n);
j = 1;
for i=1:max(size(m))
    if m(i) > 1
        Jb = gallery('jordbloc',m(i),e(i));
    else
        Jb = e(i);  % JORDBLOC fails in n = 1 case.
    end
    A(j:j+m(i)-1,j:j+m(i)-1) = Jb;
    j = j + m(i);
end

if nargin < 4
   X = randn(n);
end
A = X\A*X;

author: opencobra

repo: cobratoolbox


text:
function dicomHeaders = MRSI_read_dicom_headers(folderName)
    arguments
        folderName (1, :) char {mustBeFolder}
    end
    dicomFolder = dir(folderName);
    % remove . and .. files from dicomFolder
    dicomFolder = dicomFolder(~ismember({dicomFolder.name}, {'.', '..'}));

    dicomPath = arrayfun(@(listing) fullfile(listing.folder, listing.name), dicomFolder, 'UniformOutput', false);

    dicomHeaders = spm_dicom_headers(dicomPath);
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function upsilon = lfmvvComputeUpsilonDiagVector(gamma, sigma2, t, mode)

% LFMVVCOMPUTEUPSILONDIAGVECTOR Upsilon vector vel. vel. with t1 = t2
% FORMAT
% DESC computes a portion of the LFMVV kernel.
% ARG gamma : Gamma value for system.
% ARG sigma2 : length scale of latent process.
% ARG t : first time input (number of time points x 1).
% ARG mode : operation mode, according to the derivative (mode 0,
% derivative wrt t1, mode 1 derivative wrt t2)
% RETURN upsilon : result of this subcomponent of the kernel for the given values.
%
% COPYRIGHT : Mauricio Alvarez, 2010
%
% SEEALSO : lfmComputeUpsilonMatrix.F, lfmvpComputeUpsilonMatrix.m

% KERN

sigma = sqrt(sigma2);

if mode==0
    upsilon = gamma*lfmvpComputeUpsilonDiagVector(gamma, sigma2, t, mode) ...
        - (2*gamma/(sqrt(pi)*sigma))*exp(-gamma*t).*(exp(-(t.^2)/sigma2));
else
    upsilon = -gamma*lfmvpComputeUpsilonDiagVector(gamma, sigma2, t, mode);
end


text:
function sample_points = tapas_physio_get_sample_points(ons_secs, sqpar, slicenum)
% gets times of slice scan events of a particular slice number in every
% volume acquired
%
% USAGE
%   sample_points = tapas_physio_get_sample_points(ons_secs, sqpar, slicenum)
%
% INPUTS:
%   ons_secs
%   sqpar
%   slicenum    - slice number (1<=slicenum<=Nslices) where signal shall be
%                 sampled; alternative: specify sqpar.onset_slice
%
% OUTPUT:
%   sample_points   - absolute time (in seconds) where the specified slice was
%                     aquired for every volume

% Author: Lars Kasper
%
% Copyright (C) 2013, Institute for Biomedical Engineering, ETH/Uni Zurich.
%
% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% slicenum should be field of onset_slice
if nargin<3
    if isfield(sqpar, 'onset_slice')
        slicenum = sqpar.onset_slice;
    end
end

% default timing: first slice
if isempty(slicenum)
    slicenum = 1;
end

nSampleSlices = length(slicenum);
sample_points = zeros(sqpar.Nscans*nSampleSlices,1);
for n = 1:sqpar.Nscans
    spulse = ons_secs.spulse_per_vol{n + sqpar.Ndummies};
    if length(spulse) < max(slicenum)
        error('scan %d: only %d slice scan events. Cannot resample to slice %d', ...
            n, length(spulse), max(slicenum));
    else
        sample_points((n-1)*nSampleSlices + (1:nSampleSlices)) = spulse(slicenum);
    end
end

end

author: translationalneuromodeling

repo: tapas


text:
function g = mrdivide(f,c)
%/   Right scalar divide for BALLFUN objects.
%   X = B/A or X = mrdivide(B, A) is equivalent to X = B./A.
%
% See also MTIMES.

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

g = f*(1/c);
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function F = stream1d(A,vector)
%% Stream function for any D1Qn Lattice
% There are only 2 basic movements and any other is a combination of the
% basic ones:

% Matlab equivalent solution:
% F = circshift(A,[y,x])

%% Evaluate input data

v = length(vector);

    x = vector(2);
    y = vector(1);

n = length(A);

    b = zeros(1,n);

%% Stream in XY plane

if x == 1       % x+ stream
    b = A(:,n);
    for i = n:-1:2;
        A(:,i) = A(:,i-1);
    end
    A(:,1) = b;
elseif x == -1  % x- stream
    b = A(:,1);
    for i = 1:n-1;
        A(:,i) = A(:,i+1);
    end
    A(:,n) = b;
else            % x=0 stream
    %do nothing
end

F = A;

repo: Aero-matlab


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
function rank_training(params, cache_dir)

n_hiers = length(params.hiers);

%% Load Pareto parameters
n_cands = loadvar(params.files.pareto_point,'n_cands');

% Some checks
assert(n_hiers==size(n_cands,2));
assert(params.n_r_cand==size(n_cands,1));

%% Gather all features and quality of the training set candidates
if exist(params.files.features_file,'file')
    load(params.files.features_file);
    disp(['Loaded: ' params.files.features_file '.'])
else
    disp(['RECOMPUTING: ' params.files.features_file '.'])

    % features = [];
    % jaccards = [];

    % Load which images to consider from the params.database (train, val, etc.)
    im_ids = database_ids(params.database,params.gt_set_ranking);

    num_images = length(im_ids);
    parfor im_id = 1:num_images
        out = cache_mcg_features(params, [], im_ids{im_id}, cache_dir);
        f_lp = out.f_lp; f_ms = full(out.f_ms); feats = out.feats;
        bboxes = out.bboxes; red_cands = out.red_cands;
        b_feats_intersections = full(out.b_feats_intersection);

        if(params.depth_features)
          tt = tic();
          sp = f_lp;
          sp2regC = cands2labels(red_cands, f_ms);
          sp2reg = false(max(sp(:)), length(sp2regC));
          for i = 1:length(sp2regC), sp2reg(sp2regC{i},i) = true; end
          D = getImage(im_ids{im_id}, 'depth');
          D = double(D)./1000;
          C = params.camera_matrix;
          missingMask = getImage(im_ids{im_id}, 'rawdepth') == 0;
          fdepth = depthFeatures(sp, sp2reg, D, missingMask, C);
          feats = cat(2, feats, fdepth');
          fprintf('%s: Time for depth features: %0.3f\n', im_ids{im_id}, toc(tt));
        end

        % Eval candidates
        gt = get_ground_truth(params.database, im_ids{im_id});
        hier = struct('leaves_part', f_lp, 'ms_struct', ms_matrix2struct(f_ms));
        jacc = eval_cands(hier, red_cands, gt);
        max_jacc = max(jacc,[],1)';

        % fprintf('.');
        % Sample candidates
        % Get the optimum candidates for each object
        % to ensure they are on the training set
        [~,ids] = max(jacc,[],2);
        if (length(ids)>params.n_samples)  % More objects than samples asked
            ids = ids(1:params.n_samples);
        end
        ids_rest = setdiff(1:size(jacc,2),ids);
        if (size(jacc,2)-length(ids))>params.n_samples
            ids_rest = ids_rest(randperm(length(ids_rest),params.n_samples-length(ids)));
        end
        if isrow(ids_rest)
            ids_rest = ids_rest';
        end
        sel_ids = [ids; ids_rest];

        % Store
        % features = [features; feats(sel_ids,:)]; %#ok<AGROW>
        % jaccards = [jaccards; max_jacc(sel_ids,:)]; %#ok<AGROW>

        features{im_id} = feats(sel_ids,:);
        jaccards{im_id} = max_jacc(sel_ids,:);
    end
    features = cat(1, features{:});
    jaccards = cat(1, jaccards{:});

    save(params.files.features_file,'features','jaccards');
end

%% Train the random forest
if exist(params.files.trained_classifier,'file')
    disp(['Already trained: ' params.files.trained_classifier '.'])
else
    disp(['TRAINING: ' params.files.trained_classifier '.'])

    % Train and save the result
    rf = regRF_train(features,jaccards,50); %#ok<NASGU>
    disp('Training done')

    % Save the result
    save(params.files.trained_classifier,'rf');
end

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
function test_opt_test ( )

%*****************************************************************************80
%
%% TEST_OPT_TEST tests the TEST_OPT library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
  timestamp (  )
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_OPT_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the TEST_OPT library.\n' );

  test_opt_test01 ( );
  test_opt_test02 ( );
  test_opt_test03 ( );
  test_opt_test04 ( );
  test_opt_test05 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_OPT_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp (  )

  return
end


text:
function [ ml, mu, m ] = bandwidth ( element_order, element_num, element_node )

%*****************************************************************************80
%
%% BANDWIDTH determines the bandwidth associated with the finite element mesh.
%
%  Discussion:
%
%    The quantity computed here is the "geometric" bandwidth determined
%    by the finite element mesh alone.
%
%    If a single finite element variable is associated with each node
%    of the mesh, and if the nodes and variables are numbered in the
%    same way, then the geometric bandwidth is the same as the bandwidth
%    of a typical finite element matrix.
%
%    The bandwidth M is defined in terms of the lower and upper bandwidths:
%
%      M = ML + 1 + MU
%
%    where
%
%      ML = maximum distance from any diagonal entry to a nonzero
%      entry in the same row, but earlier column,
%
%      MU = maximum distance from any diagonal entry to a nonzero
%      entry in the same row, but later column.
%
%    Because the finite element node adjacency relationship is symmetric,
%    we are guaranteed that ML = MU.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ELEMENT_ORDER, the order of the elements.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM);
%    ELEMENT_NODE(I,J) is the global index of local node I in element J.
%
%    Output, integer ML, MU, the lower and upper bandwidths of the matrix.
%
%    Output, integer M, the bandwidth of the matrix.
%
  ml = 0;
  mu = 0;

  for element = 1 : element_num

    for local_i = 1 : element_order
      global_i = element_node(local_i,element);

      for local_j = 1 : element_order
        global_j = element_node(local_j,element);

        mu = max ( mu, global_j - global_i );
        ml = max ( ml, global_i - global_j );

      end
    end
  end

  m = ml + 1 + mu;

  return
end


text:
function B = TTeMPS_op_laplace_to_TT_matrix( A )
    %TTeMPS_to_TT Convert to TT Toolbox matrix format.
    %   TT = TT_to_TTeMPS( A ) takes the TTeMPS Laplace operator A and converts it into
    %   a tt_matrix object using the TT Toolbox 2.x from Oseledets et al.
    %   This toolbox needs to be installed, of course.
    %
    %   See also TTeMPS_to_TT, TTeMPS_op_to_TT, TTeMPS_op_laplace_to_TTeMPS_op.

    %   TTeMPS Toolbox.
    %   Michael Steinlechner, 2013-2016
    %   Questions and contact: michael.steinlechner@epfl.ch
    %   BSD 2-clause license, see LICENSE.txt

    B = TTeMPS_op_to_TT_matrix( TTeMPS_op_laplace_to_TTeMPS_op( A ));

end

author: NicolasBoumal

repo: manopt


text:
function [ r, s, area ] = node_reference_t10 ( )

%*****************************************************************************80
%
%% NODE_REFERENCE_T10 returns the basis nodes for a 10 node triangle.
%
%  Reference Element T10:
%
%    |
%    1  10
%    |  |\
%    |  | \
%    |  8  9
%    |  |   \
%    S  |    \
%    |  5  6  7
%    |  |      \
%    |  |       \
%    0  1--2--3--4
%    |
%    +--0----R---1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real R(10), S(10), the coordinates of the basis nodes.
%
%    Output, real AREA, the area of the element.
%
  r(1) = 0.0;
  s(1) = 0.0;

  r(2) = 1.0 / 3.0;
  s(2) = 0.0;

  r(3) = 2.0 / 3.0;
  s(3) = 0.0;

  r(4) = 1.0;
  s(4) = 0.0;

  r(5) = 0.0;
  s(5) = 1.0 / 3.0;

  r(6) = 1.0 / 3.0;
  s(6) = 1.0 / 3.0;

  r(7) = 2.0 / 3.0;
  s(7) = 1.0 / 3.0;

  r(8) = 0.0;
  s(8) = 2.0 / 3.0;

  r(9) = 1.0 / 3.0;
  s(9) = 2.0 / 3.0;

  r(10) = 0.0;
  s(10) = 1.0;

  area = 0.5;

  return
end


text:
function [col,bol,msz] = spm_MB_col(n)
% FORMAT [col,bol,msz] = spm_MB_col(n)
% Return colours and marker size for number of partitions
% n  - number of partitions
%__________________________________________________________________________
% Copyright (C) 2019-2020 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_MB_col.m 7768 2020-01-07 11:37:19Z spm $

% Marker colour and size
%--------------------------------------------------------------------------
s = rand('twister');
rand('twister',1);
msz   = fix(16 + 64/n);
for k = 1:n
    bol{k} = spm_softmax(log(rand(3,1))*2);
    col{k} = bol{k}*(1 - 1/2) + ones(3,1)/2;
end
rand('twister',s);

author: spm


text:
function [mg] = amu2mg(amu)
% Convert mass from atomic mass units to milligrams.
% Chad Greene 2012
mg = amu*1.6605402e-21;

author: Sable

repo: mcbench-benchmarks


text:
function [pnt, dhk] = mesh_cone(npnt)

% MESH_CONE creates a triangulated cone
%
% Use as
%   [pnt, tri] = mesh_cone(N)
%
% This creates a cone with N-2 vertices on the bottom circle and N vertices in total.
%
% See also MESH_TETRAHEDRON, MESH_OCTAHEDRON, MESH_ICOSAHEDRON, MESH_SPHERE, MESH_CUBE

% Copyright (C) 2002-2019, 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/>.

ncircle = npnt-2;
ndhk = 2*ncircle;

pnt = zeros(npnt, 3);
dhk = zeros(ndhk, 3);

pnt(1:ncircle, 1) = cos(linspace(0, 2*pi, ncircle)');		% bottom plane
pnt(1:ncircle, 2) = sin(linspace(0, 2*pi, ncircle)');		% bottom plane
pnt(npnt-1, :)  = [0 0 0];				% bottom point
pnt(npnt  , :)  = [0 0 1];				% top point

% connect all points on the circle to the bottom point
dhk(1:ncircle, 1) = npnt-1;
dhk(1:ncircle, 2) = (1:ncircle)';
dhk(1:ncircle, 3) = (2:(ncircle+1))';
dhk(  ncircle, 3) = 1;

% connect all points on the circle to the top point
dhk((ncircle+1):ndhk, 1) = npnt;
dhk((ncircle+1):ndhk, 2) = dhk(1:ncircle, 3);
dhk((ncircle+1):ndhk, 3) = dhk(1:ncircle, 2);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function families = compute_families(bnet)
% COMPUTE_FAMILIES
% precomputes the families of nodes in a bnet
%
% The return value is a cell array for now

ss = size(bnet.dag, 1);
families = cell(ss, 1);
for i = 1:ss
  families{i} = family(bnet.dag, i);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [res] = tracking_push_relabel(dres, c_en, c_ex, c_ij, betta, max_it)

dnum = length(dres.x);

dres.c = betta - dres.r; %% cost for each detection window

%%% The mex function works with large numbers.
dres.c  = dres.c  *1e6;
c_en    = c_en    *1e6;
c_ex    = c_ex    *1e6;
c_ij    = c_ij    *1e6;

n_nodes = 2*dnum+2; %% number of nodes in the graph
n_edges = 0;

dat_in = zeros(1e7,3); %% each row represents an edge from node in column 1 to node in column 2 with cost in column 3.
k_dat = 0;
for i = 1:dnum
  k_dat = k_dat+3;
  dat_in(k_dat-2,:) = [1      2*i     c_en      ];
  dat_in(k_dat-1,:) = [2*i    2*i+1   dres.c(i) ];
  dat_in(k_dat,:)   = [2*i+1  n_nodes c_ex      ];
end
for i=1:dnum
  f2 = dres.nei(i).inds;
  for j = 1:length(f2)
    k_dat = k_dat+1;
    dat_in(k_dat,:) = [2*f2(j)+1 2*i c_ij];
  end
end
dat_in = [dat_in repmat([0 1],size(dat_in,1),1)];  %% add two columns: 0 for min capacity in column 4 and 1 for max capacity in column 5 for all edges.

excess_node = [1 n_nodes];  %% push flow in the first node and collect it in the last node.

k = 0;
dat2_old = [0 0 0];

inds_all = [];
tic
lb=1;
ub=1e4;
tr_num_old = 1;

tic

while ub-lb > 1     %% bisection search for the optimum amount of flow. This can be implemented by Golden section search more efficiently.
  tr_num = round((lb+ub)/2);

  [cost_l, dat_l] = cs2_func(dat_in(1:k_dat,:), excess_node, [tr_num -tr_num]);      %% try flow = tr_num
  [cost_u, dat_u] = cs2_func(dat_in(1:k_dat,:), excess_node, [tr_num+1 -tr_num-1]);  %% try flow = tr_num+1
  if cost_u-cost_l > 0
    ub = tr_num;
  else
    lb = tr_num;
  end
  k=k+1;
end

if cost_u < cost_l
  dat1 = dat_u;
else
  cost1 = cost_l;
  dat1 = dat_l;
end

%%%% backtrack tracks to get ids
tmp   = find( dat1(:, 1) == 1);
start = dat1(tmp, 2);       %% starting nodes; is even

tmp   = find( ~mod(dat1(:, 1), 2) .* (dat1(:, 2)-dat1(:, 1) == 11) );
detcs = dat1(tmp, 1);       %% detection nodes; is even

tmp   = find( mod(dat1(:, 1), 2) .* ~mod(dat1(:, 2), 2) .* (dat1(:, 2)-dat1(:, 1) ~= 1) );
links = dat1(tmp, 1:2);     %% links; is [even  odd]

res_inds  = zeros(1, 1e5);
res_ids   = zeros(1, 1e5);

k = 0;
for i = 1:length(start);    %% for each track
  this1 = start(i);
  while this1 ~= n_nodes
    k = k+1;
    res_inds(k) = this1/2;
    res_ids(k) = i;
    this1 = links(find(links(:,1) == this1+1), 2);  %% should have only one value
    if mod(this1, 2) + (length(this1) ~= 1)         %% sanity check
      display('error in the output of solver');
    end
  end
end
res_inds  = res_inds(1:k);    %% only these detection windows are used in tracks.
res_ids   = res_ids(1:k);     %% track id for each detection window

%%% Calculate the cost value to sort tracks
this_cost = zeros(1, max(res_ids));
for i = 1:max(res_ids)  %% for each track
  inds = find(res_ids == i);
  len1= length(inds);
  track_cost(i) = sum(dres.c(res_inds(inds))) + (len1-1) * c_ij + c_en + c_ex;
end
[dummy sort_inds] = sort(track_cost);

for i = 1:length(sort_inds)
  res_ids_sorted(res_ids == sort_inds(i)) = i;
end

res = sub(dres, res_inds);
res.id = res_ids_sorted(:);

% [dummy tmp] = sort(res.id);
% res = sub(res,tmp);
% % res = sub(dres,inds_all_its(1).inds);
% % res.r = res.r/1e6;

author: yuxng


text:
function f = projectOntoBMCIII( f )
% PROJECTONTOBMCI  Projection onto BMC-III symmetry.
%
% f = projectOntoBMCIII(f) is the orthogonal projection of f onto BMC-III
% symmetry, i.e., a function that is
% 1. even in theta for every even wave number in lambda;
% 2. odd in theta for every odd wave number in lambda;
% Additionally, for all but k=0 wavenumber lambda the resulting projection
% enforces the ballfun is zero at the poles.
%
% The projection is orthogonal, i.e., the correction matrix to fix up the
% structure has the smallest possible Frobenius norm.

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

if isempty( f )
   return
end

% Get the tensor of coefficients
F = f.coeffs;

% Permute F
F = permute(F,[1,3,2]);

% Loop over lambda
for l = 1:size(F,3)
    % Update the matrix
    F(:,:,l) = projectOntoRTheta(F(:,:,l));
end

% Permute back
F = permute(F,[1,3,2]);
f = ballfun(F, 'coeffs');
end

function X = projectOntoRTheta( X )
% Project the matrix of Chebyshev--Fourier coefficients onto a BMC-II
% function. The projection is orthogonal, i.e., the correction matrix to
% fix up the structure has the smallest possible Frobenius norm.

% Get the discretization
[m,n] = size(X);

zeroMode = floor(n/2)+1;
evenModes = [fliplr(zeroMode-2:-2:1) zeroMode+2:2:n];  % Not including the zero mode
oddModes = [fliplr(zeroMode-1:-2:1) zeroMode+1:2:n];

% First do the zero-periodic mode in theta.
% Enforce the expansion is even in r.
I = eye( m ); A = I(2:2:end,:);
C = A \ ( A * X(:,zeroMode) );
C = I \ C;

% Update coeff matrix:
X(:,zeroMode) = X(:,zeroMode) - C;

% Second do the even-periodic, non-zero modes in theta.
% Enforce these are zero at the pole and that the expansion is even in
% theta
% Vectors [ T_k(0) ]:
A = real( (-1i).^(0:m-1) );   % A * X should be the zero vector.

% Solution to underdetermined system A*(X + Y) = 0 with smallest Frobenius
% norm:
C = A \ ( A * X(:,evenModes) );
C = I \ C;
X(:,evenModes) = X(:,evenModes) - C;

% Now project onto odd-anti-periodic.  Nothing special has to be done at
% the poles since enforcing the expansion is odd in r guarantees that it
% sums to zero.
A = I(1:2:end,:);
C = A \ ( A * X(:,oddModes) );
C = I \ C;
X(:,oddModes) = X(:,oddModes) - C;
end


text:
function boundary_num = boundary_edge_count_euler_q4_mesh ( node_num, ...
  element_num, hole_num )

%*****************************************************************************80
%
%% BOUNDARY_EDGE_COUNT_EULER_Q4_MESH counts boundary edges using Euler's formula.
%
%  Discussion:
%
%    We assume we are given information about a quadrilateral mesh
%    of a set of nodes in the plane.
%
%    Given the number of nodes, elements and holes, we are going to apply
%    Euler's formula to determine the number of edges that lie on the
%    boundary of the set of nodes.
%
%    The number of faces, including the infinite face and internal holes,
%    is ELEMENT_NUM + HOLE_NUM + 1.
%
%    Let BOUNDARY_NUM denote the number of edges on the boundary.
%    Each of the ELEMENT_NUM quadrilaterals uses four edges.  Every edge
%    occurs in two different elements, so the number of edges must be
%    ( 4 * ELEMENT_NUM + BOUNDARY_NUM ) / 2.
%
%    The number of nodes used in the mesh is NODE_NUM.
%
%    Euler's formula asserts that, for a simple connected figure in the
%    plane with no edge crossings, NODE_NUM nodes, EDGE_NUM edges and
%    FACE_NUM faces:
%
%      NODE_NUM - EDGE_NUM + FACE_NUM = 2
%
%    In our context, this becomes
%
%      NODE_NUM - ( 4 * ELEMENT_NUM + BOUNDARY_NUM ) / 2
%      + ELEMENT_NUM + HOLE_NUM + 1 = 2
%
%    or
%
%      BOUNDARY_NUM = 2 * NODE_NUM + 2 * HOLE_NUM - 2 * ELEMENT_NUM - 2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 February 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Marc de Berg, Marc Krevald, Mark Overmars, Otfried Schwarzkopf,
%    Computational Geometry, Section 9.1,
%    Springer, 2000.
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer HOLE_NUM, the number of internal holes.
%
%    Output, integer BOUNDARY_NUM, the number of edges that
%    lie on the boundary of the mesh.
%
  boundary_num = 2 * node_num + 2 * hole_num - 2 * element_num - 2;

  return
end

author: johannesgerer


text:
function [bic, loglik, penalty] = computebic(pathname, filename, ncomp_vec, dpixc, penalty_type)
% [bic, loglik, penalty] = computebic(pathname, filename, ncomp_vec, dpixc, penalty_type)
% penalty_type = 1 : original NMF
% penalty_type = 2 : original NMF + multidimensinoal data
% penalty_type = 3 : hidden + multidimensinoal data

if ~exist('penalty_type','var')
    penalty_type = 1;
end

for ii=1:length(ncomp_vec)
%     if strcmp(pathname(end), '/')
%         load ([pathname filename num2str(ncomp_vec(ii))])
%     else
%         load ([pathname '/' filename num2str(ncomp_vec(ii))])
%     end
%     load ([pathname filename num2str(ncomp_vec(ii)) '/' filename num2str(ncomp_vec(ii))]);
    load ([pathname '/' filename num2str(ncomp_vec(ii))]);
    loglik(ii) = -ddivergence(reshape(dpixc,peval.numpix, peval.nt),res.w*res.h);
    switch penalty_type
        case 1 % orig
            if ii==1 fprintf ('Penalty: original NMF\n'); end
            penalty(ii) = 0.5*peval.ncomp*(peval.nt+peval.numpix)*log(peval.nt*peval.numpix); %orig
        case 2 % orig + multidiensional data
            if ii==1 fprintf ('Penalty: original NMF + multidimensional data\n'); end
            penalty(ii) = 0.5*peval.ncomp*(peval.nt+peval.numpix)*log(peval.nt);
        case 3 % hidden + multidimensional data
            if ii==1 fprintf ('Penalty: multidimensional + hidden\n'); end
            penalty(ii) = 0.5*peval.ncomp*(peval.numpix)*log(peval.nt);
    end

% penalty(ii) = 0.5*peval.ncomp*(peval.nt+peval.numpix)*log(peval.nt*peval.numpix);
% penalty(ii) = 0.5*peval.ncomp*(peval.nt+peval.numpix)*log(peval.nt*peval.numpix); %orig
% penalty(ii) = 0.5*peval.ncomp*(peval.numpix)*log(peval.nt); %hidden +
% multidimensional data
% penalty(ii) = 0.5*peval.ncomp*(peval.nt+peval.numpix)*log(peval.nt); % orig + multidiensional data

%     bic(ii) = -ddivergence(reshape(dpixc,peval.numpix, peval.nt),res.w*res.h) - 0.5*peval.ncomp*peval.nt*log(peval.nt*peval.numpix);
end
bic =  loglik - penalty;

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function dtiRawReprocess(subDir, dtiDirName, rawBaseName, t1File, rawDirName)
%
% dtiRawReprocessMany(subDir, dtiDirName, dwRawBaseName, t1File, [rawDirName='raw'])
%
% E.g.:
% subDirList = '/biac3/wandell4/data/reading_longitude/dti_adults/ab_090809';
% dtiDirName = 'dti06';
% rawBaseName = 'dti_g13_b800';
% t1File = 'full/path/to/t1.nii.gz';
%
% DY 09/2008: modified from dtiRawReprocessMany to work for just one
% subject at a time in the context of dti_FFA_preprocessScript_fixEddyRT
%
% NOTE: I also set clobber to TRUE just in case (forces recomputation of
% everything). Also, I checked that the eddyCorrect flag is set to default
% (true).

if(~exist('rawDirName','var')||isempty(rawDirName))
    rawDirName = 'raw';
end

rawDir = fullfile(subDir, rawDirName);
d = dir(fullfile(rawDir, [rawBaseName '.*']));
bvec = ''; bval = ''; dwRaw = '';

if(numel(d)>=3)
    tmp = {d.name};
    % now sort bvecs, bvals, nifti file
    try
        bvec = tmp{find(~cellfun('isempty',regexpi(tmp,'\.bvec.$','once')))};
        bval = tmp{find(~cellfun('isempty',regexpi(tmp,'\.bval.$','once')))};
        dwRaw = tmp{find(~cellfun('isempty',regexpi(tmp,'\.nii$|\.nii.gz','once')))};
    catch
    end
end

if(isempty(bvec)||isempty(bval)||isempty(dwRaw))
    disp('FAILED: Can''t find dwRaw/bvec/bval files.');
else
    if(~exist(t1File,'file'))
        disp('FAILED: Can''t find t1 file.');
    else
        outDir = fullfile(subDir, dtiDirName);
        oldDir = fullfile(subDir, [dtiDirName '_OLD' datestr(now,'yymmdd')]);
        % Rename the current data dir and put old raw 'aligned' files
        % in there.
        try
            movefile(outDir, oldDir);
            movefile(fullfile(rawDir, [rawBaseName '_aligned.*']), oldDir);
        catch
        end
        dtiRawPreprocess(fullfile(rawDir,dwRaw), t1File, [], [], true, outDir);
    end
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function julia(arg1,arg2)
%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
%               Julia Fractals Using Matlab
%               Written By Sridharan, Mithun Aiyswaryan
%               Christian Albrechts Universität zu Kiel, Germany
%               Mail Your Comments At: s.mithun@indiatimes.com
%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%               Help :
%               This program takes two arguments and computes the
%               Julia sets using the provided values.
%               For the best viewing results, varhoose the arguments such that:
%
%               Argument 1 > 10     &             Argument 2 > 100
%
%               Note : The generation of fractals is  computationally
%                      very intensive and it may take you a while before
%                      you observe the results on the screen!
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

xaxis=0;
yaxis=0;
l=1.5;
x=linspace(xaxis-l,xaxis+l,arg2);
y=linspace(yaxis-l,yaxis+l,arg2);
[xtrans,ytrans]=meshgrid(x,y);
var= -.745429;
ztrans=xtrans+i*ytrans;
for k=1:arg1;
ztrans=ztrans.^2+var;
t=exp(-abs(ztrans));
end
colormap prism(256)
pcolor(t);
shading flat;
axis('square','equal','off');


text:
function [rhsu] = DispersiveLDGRHS1D(u,time)

% function [rhsu] = DispersiveLDGRHS1D(u,time)
% Purpose  : Evaluate RHS flux in 1D u_xxx using LDG fluxes and periodic BC's

Globals1D;

% Define field differences at faces, incl BC
du   = zeros(Nfp*Nfaces,K); du(:)     = u(vmapM)-u(vmapP);
uin  = u(vmapO);            du (mapI) = u(vmapI) - uin;
uout = u(vmapI);            du (mapO) = u(vmapO) - uout;
fluxu = nx.*(1.0+nx).*du/2.0;

% Compute local variable p, define differences, incl BC
p = rx.*(Dr*u) - LIFT*(Fscale.*fluxu);
dp   = zeros(Nfp*Nfaces,K); dp(:)    = p(vmapM)-p(vmapP);
pin  = p(vmapO);            dp(mapI) = p(vmapI) - pin;
pout = p(vmapI);            dp(mapO) = p(vmapO) - pout;
fluxp = nx.*(1.0-nx).*dp/2.0;

% Compute local variable q, define differences, incl BC
q = rx.*(Dr*p) - LIFT*(Fscale.*fluxp);
dq   = zeros(Nfp*Nfaces,K); dq(:)     = q(vmapM)-q(vmapP);
qin  = q(vmapO);            dq (mapI) = q(vmapI) - qin;
qout = q(vmapI);            dq (mapO) = q(vmapO) - qout;
fluxq = nx.*(1.0-nx).*dq/2.0;

% compute right hand sides of the semi-discrete PDE
rhsu = rx.*(Dr*q) - LIFT*(Fscale.*fluxq);
return


text:
function lattice_rule_test11 ( )

%*****************************************************************************80
%
%% LATTICE_RULE_TEST11 tests MONTE_CARLO;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 November 2008
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LATTICE_RULE_TEST11\n' );
  fprintf ( 1, '  MONTE_CARLO applies a Monte Carlo scheme\n' );
  fprintf ( 1, '  to estimate the integral of a function\n' );
  fprintf ( 1, '  over the unit hypercube.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The spatial dimension DIM_NUM = %d\n', dim_num );

  a(1:dim_num) = 0.0;
  b(1:dim_num) = 1.0;

  exact = e_01_2d ( dim_num, a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '         M    EXACT       ESTIMATE  ERROR\n' );
  fprintf ( 1, '\n' );

  for k = 2 : 5

    m = 10^k;

    quad = monte_carlo ( dim_num, m, @f_01_2d );

    error = abs ( exact - quad );

    fprintf ( 1, '  %8d  %10.6f  %10.6f  %10.6e\n', m, exact, quad, error );

  end

  return
end

repo: jburkardt-m


text:
function dat= proc_movingAverage(dat, ms, varargin)
%PROC_MOVINGAVERAGE - Moving average (low-pass) filter
%
%Usage:
% DAT= proc_movingAverage(DAT, MSEC, <METHOD='causal'>)
%
%Input:
% DAT    - data structure of continuous or epoched data
% MSEC   - length of interval in which the moving average is
%          to be calculated, unit [msec].
% METHOD - 'centered' or 'causal' (default).
%
%Output:
% DAT    - updated data structure

% Author(s): Benjamin Blankertz

misc_checkType(dat, 'STRUCT(x fs)');
misc_checkType(ms, '!DOUBLE[1]');

nSamples = round(ms*dat.fs/1000);
dat.x(:,:)= procutil_movingAverage(dat.x(:,:), nSamples, varargin{:});

author: bbci


text:
function [sys,x0,str,ts]=abcrtoqdor(t,x,u,flag)
%dfsfgsfdgdghdhdfg
switch flag
    case 0
        [sys,x0,str,ts]=mdlInitializeSizes;
    case 3
        sys=mdlOutput(t,x,u);
    case {1,2,4,9}
        sys=[];
    otherwise
        error(['Unhandel flag =',num2str(flag)]);
end;
%========================================================================
%========================================================================
function [sys,x0,str,ts]=mdlInitializeSizes
sizes=simsizes;
sizes.NumContStates= 0;
sizes.NumDiscStates= 0;
sizes.NumOutputs= 1;
sizes.NumInputs= 2;
sizes.DirFeedthrough=1;
sizes.NumSampleTimes=1;
sys=simsizes(sizes);
x0=[];
str=[];
ts=[-1 0];
%=========================================================================
%=========================================================================
function sys=mdlOutput(t,x,u);
%=========================================================================
 vc=u(1);
 vm=u(2);
q=(vc*pi)/(3*vm);
if q>1
    q=1;
end
if q<-1
    q=-1;
end
a=acos(q);
  sys=a*180/pi;


text:
 function display(ob)
%function display(ob)
% "display" method for this class

name = inputname(1);
printf('"%s" is an object of class "%s":', name, class(ob))
ob = struct(ob);
disp(ob)

fnames = fieldnames(ob);
for ii=1:length(fnames)
	fname = fnames{ii};
	if isstruct(ob.(fname))
		printf('%s.%s :', inputname(1), fname)
		disp(ob.(fname))
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% SCRIPT TEST FOR THE KUKA LBR ROBOT KINEMATICS

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
robot = load_robot('KUKA', 'LBR_IIWA_R820_COP');
close all

% test both solutions: based on the transpose an on the Moore-Penrose
q = pi/8*[1,1,1,1,1,1,1]';
q = [0.1 0.2 .3 .4 .5 .6 .7];
drawrobot3d(robot, q)
J = manipulator_jacobian(robot, q);

% rigth
R = J'*inv(J*J');

% left
L=inv(J'*J)*J';

M = pinv(J)

R-L

repo: ARTE


text: %NE Dataset overload

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function variance = r8vec_variance ( n, x )

%*****************************************************************************80
%
%% R8VEC_VARIANCE returns the variance of an R8VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, real X(N), the vector whose variance is desired.
%
%    Output, real VARIANCE, the variance of the vector entries.
%
  mean = r8vec_mean ( n, x );

  variance = sum ( ( x(1:n) - mean ).^2 );

  if ( 1 < n )
    variance = variance / ( n - 1 );
  else
    variance = 0.0;
  end

  return
end


text:
function aboxes_out = post_process_candidate_detections_all_imgs(...
    aboxes_in, nms_iou_thrs, thresh_val, varargin)
% post_process_candidate_detections_all_imgs performs the post-processing
% step of non-maximum-suppression and optionally of box voting on the
% candidate bounding box detections of a set of images.
%
% INPUTS:
% 1) aboxes_in: contains the candidate bounding box detections of each image.
% It can be of two forms:
%   a) (the flag is_per_class must being set to false) a NI x 1 cell array
%   where NI is the number of images. The aboxes_in{i} element contains
%   the candidate bounding box detection of the i-th image in the form of a
%   NB_i x (4 + C) array where C is the number of categories and NB_i
%   is the number of candidate detection boxes of the i-th image. The first
%   4 columns of this array contain the bounding box coordinates in the
%   form of [x0,y0,x1,y1] (where the (x0,y0) and (x1,y1) are the top-left
%   and bottom-right corners) and the rest C columns contain the confidence
%   score of the bounding boxes for each of the C categories.
%   b) (the flag is_per_class must being set to true) a C x 1 cell array
%   where C is the number of categories. The aboxes_in{j} element in this case
%   is NI x 1 cell array, where NI is the number of images, with the candidate
%   bounding box detections of the j-th category for each image. The
%   element aboxes_in{j}{i} is a NB_{i,j} x 5 array, where NB_{i,j} is the
%   number candidate bounding box detections of the i-th image for the j-th
%   category. The first 4 columns of this array are the bounding box
%   coordinates in the form of [x0,y0,x1,y1] (where the (x0,y0) and (x1,y1)
%   are the top-left and bottom-right corners) and the 5-th column contains
%   the confidence score of each bounding box with respect to the j-th category.
% 2) nms_iou_thrs: scalar value with the IoU threshold that will be used
% during the non-max-suppression step.
% 3) thresh_val: is a C x 1 array, where C is the number of categories.
% It must contains the threshold per category that will be used for
% removing candidate boxes with low confidence prior to applying the
% non-max-suppression step.
% OPTIONAL INPUT:
% 4) is_per_class: boolean value that if set to false, then the 1.a) form
% of the aboxes_in input parameter will be expected; otherwise, if set to
% true then the 1.b) form of the aboxes_in input parameter will be expected
% 5) max_per_image: scalar value with the maximum number of detection per
% image and per category. Default value: 200.
% 6) do_bbox_voting: boolean value that if is set to True then the box voting
% step is applied. Default value: false
% 7) box_ave_iou_thresh: scalar value with the minimum IoU threshold that
% is used in order to define the neighbors of bounding box during the box
% voting step. Default value: 0.5
% 8) add_val: scalar value that is added to the confidence score of bounding
% boxes in order to compute the box weight during the box voting step.
% Default value: 1.5
% 9) use_not_static_thr: boolean value that if set to true, then the threshold
% (param thresh_val) that is applied on the candidate bounding box
% detections will be modified (decreased from its input value) such that
% the average number of detections per image and per category to be
% ave_per_image. Default value: true
% 10) ave_per_image: scalar value with the desired average number of
% detections per image and per category. Used only when the use_not_static_thr
% input parameter is set to true. default value: 5
%
%
% OUTPUT:
% 1) aboxes_out: is a C x 1 array where C is the number of
% categories. The aboxes_out{j} element is a NI x 1 cell array, where NI is
% the number of images, with the bounding box detections of the j-th
% category for each image. The element aboxes_out{j}{i} is a ND_{i,j} x 5 array,
% where ND_{i,j} is the number bounding box detections for the i-th image
% and the j-th category. The first 4 columns of this array contain the
% bounding box coordinates in the form of [x0,y0,x1,y1] (where the (x0,y0)
% and (x1,y1) are the top-left and bottom-right corners) and the 5-th
% column contains the confidence score of each bounding box with respect to
% the j-th category.
%
% This file is part of the code that implements the following ICCV2015 accepted paper:
% title: "Object detection via a multi-region & semantic segmentation-aware CNN model"
% authors: Spyros Gidaris, Nikos Komodakis
% institution: Universite Paris Est, Ecole des Ponts ParisTech
% Technical report: http://arxiv.org/abs/1505.01749
% code: https://github.com/gidariss/mrcnn-object-detection
%
% Part of the code in this file comes from the R-CNN code:
% https://github.com/rbgirshick/rcnn
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2015 Spyros Gidaris
%
% "Object detection via a multi-region & semantic segmentation-aware CNN model"
% Technical report: http://arxiv.org/abs/1505.01749
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------
% Copyright (c) 2014, Ross Girshick
%
% This file is part of the R-CNN code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

ip = inputParser;
ip.addParamValue('is_per_class',      false, @islogical);
ip.addParamValue('do_bbox_voting',    false, @islogical);
ip.addParamValue('box_ave_thresh',      0.5, @isnumeric);
ip.addParamValue('use_not_static_thr', true, @islogical);
ip.addParamValue('ave_per_image',         5, @isnumeric);
ip.addParamValue('max_per_image',       200, @isnumeric);
ip.addParamValue('add_val',             1.5, @isnumeric);

ip.parse(varargin{:});
opts = ip.Results;

is_per_class       = opts.is_per_class;
do_bbox_voting     = opts.do_bbox_voting;
box_ave_thresh     = opts.box_ave_thresh;
use_not_static_thr = opts.use_not_static_thr;
max_per_image      = opts.max_per_image;
ave_per_image      = opts.ave_per_image;
add_val            = opts.add_val;

if ~is_per_class
    num_imgs    = length(aboxes_in);
    num_classes = size(aboxes_in{1},2) - 4;
else
    num_classes = length(aboxes_in);
    num_imgs    = length(aboxes_in{1});
end

max_per_set   = ceil(ave_per_image * num_imgs);
aboxes_out    = cell(num_classes, 1);
for i = 1:num_classes, aboxes_out{i} = cell(num_imgs, 1); end

thresh_val_all = ones(num_classes,1,'single') * thresh_val;

for i = 1:num_imgs
    % get the candidate bounding box detection of one image
    if ~is_per_class
        assert(size(aboxes_in{i},2) == (4 + num_classes));
        bbox_cand_dets = aboxes_in{i};
    else
        bbox_cand_dets = cell(num_classes,1);
        for j = 1:num_classes
            assert(size(aboxes_in{j}{i},2) == 5);
            bbox_cand_dets{j} = aboxes_in{j}{i};
        end
    end

    % post-process the candidate bounding box detection of one image
    bbox_detections_per_class = post_process_candidate_detections( ...
        bbox_cand_dets, thresh_val_all, nms_iou_thrs, max_per_image, ...
        do_bbox_voting, box_ave_thresh, add_val );

    for j = 1:num_classes
        aboxes_out{j}{i} = bbox_detections_per_class{j};
    end

    if (mod(i, 1000) == 0) && use_not_static_thr
        for j = 1:num_classes
            [aboxes_out{j}, thresh_val_all(j)] = keep_top_k(aboxes_out{j}, i, max_per_set, thresh_val_all(j));
        end
    end
end

if use_not_static_thr
    disp(thresh_val_all(:)');
    for i = 1:num_classes
        % go back through and prune out detections below the found threshold
        aboxes_out{i} = prune_detections(aboxes_out{i}, thresh_val_all(i));
    end
end

total_num_bboxes = zeros(num_classes, 1);
for j = 1:num_classes
    total_num_bboxes(j) = sum(cellfun(@(x) size(x,1), aboxes_out{j}, 'UniformOutput', true));
end

fprintf('Average number of bounding boxes per class\n');
disp(total_num_bboxes(:)' / num_imgs);
fprintf('Average number of bounding boxes in total\n');
disp(sum(total_num_bboxes) / num_imgs);
end

function [boxes, thresh] = keep_top_k(boxes, end_at, top_k, thresh)
% keep top K detections
X = cat(1, boxes{1:end_at});
if isempty(X), return; end

scores = sort(X(:,end), 'descend');
thresh = scores(min(length(scores), top_k));
for image_index = 1:end_at
    bbox = boxes{image_index};
    keep = find(bbox(:,end) >= thresh);
    boxes{image_index} = bbox(keep,:);
end

end

function bbox_dets = prune_detections(bbox_dets, thresh)
for j = 1:length(bbox_dets)
    if ~isempty(bbox_dets{j})
        bbox_dets{j}(bbox_dets{j}(:,end) < thresh ,:) = [];
    end
end
end

author: gidariss

repo: mrcnn-object-detection


text:
classdef LossNormalized < dagnn.Loss
%   properties
%     loss = 'softmaxlog'
%     ignoreAverage = false
%     opts = {}
%   end
%   properties (Transient)
%     average = 0
%     numAveraged = 0
%   end

  methods
    function outputs = forward(obj, inputs, params)
      outputs{1} = vl_nnloss(inputs{1}, inputs{2}, [], 'loss', obj.loss, obj.opts{:}) ;
      obj.accumulateAverage(inputs, outputs);
      if numel(size(inputs{1}))>3
        bs = size(inputs{1},4) ;
      else
        bs = 1 ;
      end
      outputs{1} = outputs{1} / bs ;
    end

    function accumulateAverage(obj, inputs, outputs)
      if obj.ignoreAverage, return; end;
      n = obj.numAveraged ;
      m = n + size(inputs{1}, 1) *  size(inputs{1}, 2) * size(inputs{1}, 4);
      obj.average = bsxfun(@plus, n * obj.average, gather(outputs{1})) / m ;
      obj.numAveraged = m ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      if numel(size(inputs{1}))>3
        bs = size(inputs{1},4) ;
      else
        bs = 1 ;
      end

      derInputs{1} = vl_nnloss(inputs{1}, inputs{2}, derOutputs{1}, 'loss', obj.loss, obj.opts{:}) / bs;
      derInputs{2} = [] ;
      derParams = {} ;
    end

    function reset(obj)
      obj.average = 0 ;
      obj.numAveraged = 0 ;
    end

    function outputSizes = getOutputSizes(obj, inputSizes, paramSizes)
      outputSizes{1} = [1 1 1 inputSizes{1}(4)] ;
    end

    function rfs = getReceptiveFields(obj)
      % the receptive field depends on the dimension of the variables
      % which is not known until the network is run
      rfs(1,1).size = [NaN NaN] ;
      rfs(1,1).stride = [NaN NaN] ;
      rfs(1,1).offset = [NaN NaN] ;
      rfs(2,1) = rfs(1,1) ;
    end

    function obj = LossNormalized(varargin)
      obj.load(varargin) ;
    end
  end
end

author: hbilen

repo: dynamic-image-nets


text:
function [rar,yclean] = spm_rar (Z,p,m,verbose)
% Bayesian autoregressive modelling with zero-mean Gaussian mixture noise
% function [rar,yclean] = spm_rar (Z,p,m,verbose)
%
% Z          [N x 1] vector of data points
% p          Number of AR coefficients
% m          Number of mixture components (default=2)
% verbose    0/1 to printout inner workings (default=0)
%
% rar        Returned model
% yclean     'Clean' data (ie. with outlier errors removed)
%
% -------------------------------------------------------
% The fields in rar are:
%
% p                The number of AR coefficients
% m                The number of components
% fm               The negative free energy
%
%                  In the field priors:
% lambda_0         Dirichlet parameters for mixing coeffs
% b_0,c_0          Gamma parameters for precisions
%
%                  In the field posts:
% lambda           Dirichlet parameters  for mixing coeffs
% b,c              Gamma parameters for precisions
% a_mean           AR parameters (posterior mean)
% a_cov            AR parameters (posterior cov)
% b_alpha,c_alpha  Gamma parameters for weight precisions
%
%                  Mean posterior values:
% pi               mixing coefficients (lambda/sum(lambda))
% variances        variances (1./(b.*c))
%
% gamma            the responsibilities of each noise component
%
% For details of algorithm see:
%
% S.J. Roberts and W.D. Penny. Variational Bayes for Generalised Autoregressive
% models. IEEE Transactions on Signal Processing, 50(9):2245-2257, 2002
%___________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_rar.m 1276 2008-03-28 18:29:19Z guillaume $

if nargin < 3 | isempty(m)
    m=2;
end

if nargin < 4 | isempty(verbose)
  verbose=0;
end

if (m==1)
  rar=spm_ar(Z,p);
  return;
end

N=length(Z);

% Initialise AR coefficients to maximum likelihood solution
Z=Z(:);
y=Z(p+1:end);
for i=1:p,
    x(:,i)=Z(p-i+1:end-i);
end
x=-x;
y2=y.^2;
xt=x';
N=size(x,1);
wun=ones(N,1);

a_mean = pinv(x)*y;
y_pred = x*a_mean;
err=y-y_pred;
v=mean((y-y_pred).^2);
a_cov = v*inv(x'*x);

% Set mixing priors
lambda_0=5;

% Setting to these values gives updates for mean_alpha
% v. close to evidence framework
bishop_prior=1;
if bishop_prior
  b_alpha_prior=1000;
  c_alpha_prior=0.001;
  b_0=1000;
  c_0=0.001;
end
mean_alpha_prior=b_alpha_prior*c_alpha_prior;

% Cluster on absolute difference from mean
zmix=spm_kmeans1(abs(err-mean(err)),m);
% Posterior for mixers
lambda=100*zmix.pi;
zmean=[zmix.m].^2;
% Posterior for precisions
var_precision=1/std(zmean)^2;
for s=1:m,
      % Set so that b*c=precision and  b^2*c=var_precision
      precision=1/zmean(s);
      b(s)=var_precision/precision;
      c(s)=(precision^2)/var_precision;
end

% Initialise weight precision posterior
E_w=a_mean'*a_mean;
b_alpha=0.5*E_w+0.5*trace(a_cov)+(1/b_alpha_prior);
b_alpha=1/b_alpha;
c_alpha=0.5*p+c_alpha_prior;
mean_alpha=b_alpha*c_alpha;

if verbose
  disp('Init');
  disp('AR Coefficients');
  a_mean(:)'
  for s=1:m,
    disp(sprintf('State %d mix=%1.2f var=%1.2f',s,lambda(s)/sum(lambda),1/(b(s)*c(s))));
  end
end

lik=[];
tol=0.0001;
max_loops=32;
WLOOPS=5;
for loops=1:max_loops,

    % E-step
    lambda_tot=sum(lambda);
    ypred=x*a_mean;
    ypred2=ypred.^2;
    y_err=sum(xt.*(a_cov*xt));
    tv=y2-2*ypred.*y+y_err'+ypred2;
    tv=tv';
    for s=1:m,
        log_tilde_pi(s)=psi(lambda(s))-psi(lambda_tot);
        log_tilde_beta(s)=psi(c(s))+log(b(s));
        tilde_pi(s)=exp(log_tilde_pi(s));
        tilde_beta(s)=exp(log_tilde_beta(s));
        mean_beta(s)=c(s)*b(s);
        tilde_var(s,:)=tv;
        gamma(s,:)=tilde_pi(s)*(tilde_beta(s)^0.5)*exp(-0.5*mean_beta(s)*tv);
    end
    gamma_n=sum(gamma);
    for s=1:m,
        if mean(gamma_n) > eps
            % If component still exists
            gamma(s,:)=gamma(s,:)./gamma_n;
        end
    end

    % M-step
    % Part I
    for s=1:m,
        pi_bar(s)=mean(gamma(s,:));
        N_bar(s)=N*pi_bar(s);
        mean_bar(s)=mean(gamma(s,:)'.*y);
        dg=diag(gamma(s,:));
        x_bar(s,:)=mean(dg*x);
        var_bar_mu(s)=mean(gamma(s,:).*tilde_var(s,:));
    end

    % CALCULATE THE FREE ENERGY
    avg_likelihood=sum(N_bar.*(log_tilde_pi+0.5*log_tilde_beta));
    fit=-0.5*N*sum(mean_beta.*var_bar_mu);
    % Ensure that 0 log 0 = 0
    ent_s=sum(sum(-gamma.*log(gamma+eps)));
    avg_likelihood=avg_likelihood+fit+ent_s;
    lambda_p=lambda_0*ones(1,m);
    kl_dir=spm_kl_dirichlet(lambda,lambda_p,log_tilde_pi);
    kl_gamm=0;
    for s=1:m,
        kl_gamm=kl_gamm+spm_kl_gamma(b(s),c(s),b_0,c_0);
    end
    kl_weights=spm_kl_normal(a_mean,a_cov,zeros(1,p),(1/mean_alpha)*eye(p));
    kl_alpha=spm_kl_gamma(b_alpha,c_alpha,b_alpha_prior,c_alpha_prior);
    fm= avg_likelihood - kl_dir - kl_gamm - kl_weights - kl_alpha;

    % Convergence criterion
    oldlik=lik;
    lik=fm;

    if (mod(loops-1,WLOOPS)==0)
        if (loops>1)
            if abs((lik-oldlik)/lik) < tol
                break;
            end
        end
    end

    % M-Step: Part II
    for s=1:m,
        % Mixers
        lambda(s)=N_bar(s)+lambda_0;
        % Precisions
        b(s)=1/(0.5*N*var_bar_mu(s)+1/b_0);
        c(s)=0.5*N_bar(s)+c_0;
        mean_beta(s)=c(s)*b(s);
    end

    if mod(loops,WLOOPS)==0
        % Weight precisions
        E_w=0.5*a_mean'*a_mean;
        b_alpha=E_w+0.5*trace(a_cov)+(1/b_alpha_prior);
        b_alpha=1/b_alpha;
        c_alpha=0.5*p+c_alpha_prior;
        mean_alpha=b_alpha*c_alpha;

        % AR coefficients
        cc=zeros(p,p);
        cw=zeros(p,1);
        for s=1:m,
            dg=diag(gamma(s,:));
            cc=cc+mean_beta(s)*x'*dg*x;
            cw=cw+mean_beta(s)*x'*dg*y;
        end
        cc=cc+mean_alpha*eye(p);
        a_cov=inv(cc);
        a_mean=a_cov*cw;
    end

    if verbose
        disp(sprintf('It=%d, L_AV =%1.2f, KL Mix=%1.2f, KL Prec=%1.2f, KL-AR=%1.2f, KL-alpha=%1.2f, Fm=%1.2f',loops,avg_likelihood,kl_dir,kl_gamm,kl_weights,kl_alpha,fm));
    end

end

% Put variables into data structure
rar.posts.a_mean=a_mean;
rar.posts.a_cov=a_cov;
rar.m=m;
rar.fm=fm;
rar.priors.lambda_0=lambda_0;
rar.priors.c_0=c_0;
rar.priors.b_0=b_0;

for k=1:m,
  rar.posts.lambda(k)=lambda(k);
end

for k=1:m,
  rar.posts.c(k)=c(k);
  rar.posts.b(k)=b(k);
end

for k=1:m,
  rar.pi(k)=lambda(k)/sum(lambda);
end
rar.variances=1./(b.*c);

% Pre-pad gamma with zeros to get original length time series
rar.gamma=[zeros(m,p),gamma];

% Get 'clean' data
if m > 1
    [tmp,outlier_class]=min(rar.pi);
    e=y-ypred;
    outlier_error=gamma(outlier_class,:)'.*e;
    yclean=ypred+e-outlier_error;
    yclean=[Z(1:p);yclean];  % Pre-padding
else
    yclean=Z;
end


text:
function scimat = scimat_surface2seg(scimat, tri, x)
% scimat_surface2seg  Convert triangular mesh into segmentation.
%
%   scimat_surface2seg sets to "1" those voxels in a segmentation volume
%   that contain or are intersected by the triangles in a triangular mesh.
%   The method works by sampling each triangle in the mesh uniformly, using
%   barycentric coordinates. Samples are close enough to each other so that
%   each voxel will have at least two samples. Then sample coordinates are
%   rounded to the closest voxel centres, and the voxels are set.
%
% SCIMAT2 = scimat_surface2seg(SCIMAT, TRI, X)
%
%   SCIMAT is a struct to hold the segmentation (see "help scimat" for
%   details). SCIMAT.axis must be provided. SCIMAT.axis is important
%   because it tells us where the segmentation volume begings and ends, and
%   the voxel size. SCIMAT.data is not necessary, and if it is not
%   provided, a new one will be created. If it is provided, previously set
%   voxels will not be cleared.
%
%   (TRI, X) describe the triangular mesh.
%
%   TRI is a 3-column matrix. Each row contains the 3 nodes that form one
%   triangular facet in the mesh.
%
%   X is a 2-column matrix. X(i, :) contains the xy-coordinates of the
%   i-th node in the mesh.
%
%   SCIMAT2 is the output segmentation.
%
% See also: surface_tridomain, surface_param, surface_interp.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2014 University of Oxford
% Version: 0.1.1
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
narginchk(3, 3);
nargoutchk(0, 1);

% smallest dimension of the voxel size
lmin = min([scimat.axis.spacing]);

% if the scimat struct has no data field, create one
if (~isfield(scimat, 'data'))
    scimat.data = zeros([scimat.axis.size], 'uint8');
end

% loop triangles
Ntri = size(tri, 1);
for I = 1:Ntri

    % vertices of current triangle
    v = x(tri(I, :), :);

    % sample triangle
    xi = sample_triangle(v, lmin);

    % convert real world coordinates to indices
    idx = round(scimat_world2index(xi, scimat));

    % linearize indices
    idx = sub2ind(size(scimat.data), idx(:, 1), idx(:, 2), idx(:, 3));

    % set voxels that belong to the triangle to 1
    scimat.data(idx) = 1;

end

end

% given a triangle, sample it uniformly using barycentric coordinates
function xi = sample_triangle(x, lmin)

% nomenclature
v1 = x(1, :);
v2 = x(2, :);
v3 = x(3, :);

% compute number of points between triangle vertices. We want to have at
% least 2 samples within each voxel
dmax = max([norm(v1-v2) norm(v2-v3) norm(v1-v3)]);
N = ceil(1 + 2 * dmax / lmin);

% sample the triangle
a1 = linspace(0, 1, N);
a2 = linspace(0, 1, N);

% barycentric coordinate variables
[a1, a2] = meshgrid(a1, a2);
a1 = a1(:);
a2 = a2(:);
idx = a1 + a2 <= 1;
a1 = a1(idx);
a2 = a2(idx);
a3 = 1 - a1 - a2;

% point coordinates
xi = [a1 a1 a1] .* repmat(v1, length(a1), 1) ...
    + [a2 a2 a2] .* repmat(v2, length(a2), 1) ...
    + [a3 a3 a3] .* repmat(v3, length(a3), 1);

% % plot sampled triangle
% hold off
% plot3(xi(:, 1), xi(:, 2), xi(:, 3), '.')
% hold on
% plot3(x([1:end 1], 1), x([1:end 1], 2), x([1:end 1], 3))

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [dat,algo] =  training(algo,dat)

 disp(['training ' get_name(algo) '.... '])
 [algo.input algo.target]= get_xy(dat);
 [n D] = size(algo.input);

 algo.H = zeros(D+2,1);
 [algo.H, fX] = minimize(algo.H, 'gpS00', algo.length, {algo.input, algo.target});

 dat=set_x(dat,get_y(dat));

author: canlab

repo: CanlabCore


text:
function [ vX, paramLambda ] = SolveBp( mA, vB, paramEpsilon )
% ----------------------------------------------------------------------------------------------- %
%[ vX ] = SolveBpAdmm( mA, vB, paramLambda )
% Solve Basis Pursuit (Q1Eps) problem using ADMM.
% Input:
%   - mA                -   Input Matirx.
%                           The model matrix.
%                           Structure: Matrix (m X n).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - vB                -   Input Vector.
%                           The model known data.
%                           Structure: Vector (m X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - paramEpsilon      -   Parameter Epsilon.
%                           Sets the threshold for the Least squares error
%                           of the solution.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: [0, inf).
% Output:
%   - vX                -   Output Vector.
%                           Structure: Vector (n X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - paramLambda       -   Parameter Lambda.
%                           Sets the balance between L1 minimization and
%                           Least Squares minimization.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: [0, inf).
% References
%   1.  A
% Remarks:
%   1.  A
% Known Issues:
%   1.  A
% TODO:
%   1.  B
% Release Notes:
%   -   1.0.000     03/04/2018
%       *   First realease version.
% ----------------------------------------------------------------------------------------------- %

SOLVER_METHOD_ADMM  = 1; %<! ADMM
SOLVER_METHOD_PPGM  = 2; %<! PRoximal Gradient Method
SOLVER_METHOD_CD    = 3; %<! Coordinate Descent

solverMethod = SOLVER_METHOD_CD;

switch(solverMethod)
    case(SOLVER_METHOD_ADMM)
        hOptFun = @(paramLambda) (0.5 * sum(((mA * SolveBpAdmm(mA, vB, paramLambda)) - vB) .^ 2)) - paramEpsilon;
    case(SOLVER_METHOD_PPGM)
        hOptFun = @(paramLambda) (0.5 * sum(((mA * SolveLsL1ProxAccel(mA, vB, paramLambda, 200)) - vB) .^ 2)) - paramEpsilon;
    case(SOLVER_METHOD_CD)
        hOptFun = @(paramLambda) (0.5 * sum(((mA * SolveLsL1Cd(mA, vB, paramLambda, 200)) - vB) .^ 2)) - paramEpsilon;
end

sSolverOptions = optimset('Display', 'off');

paramLambda = fzero(hOptFun, [0.00001, 1000], sSolverOptions);

vX = SolveBpAdmm(mA, vB, paramLambda);

end

author: RoyiAvital

repo: StackExchangeCodes


text:
%% rangeframe
%
%  James Houghton
%  James.P.Houghton@gmail.com
%  December 8 2009
%
% A range frame lets you include information in what would otherwise
% be dumb ink on your plot, such as mean, median, quartiles, or standard deviation.
%
% This code based upon:
%
%  Edward Tufte: Visual Display of Quantitative Information
%  Second Edition
%  Graphics Press
%  pg 132
%
%% ToDo:
%
% # Clean up tick/number positioning
% # Find a better way to name data pulled out of the get return array
%
%% Function Definition
% Inputs:
%
% * *main_axis* - handle to the axis we intend to modify
% * *x_frame*   - a five element vector showing horizontal range frame limits, described below
% * *y_frame*   - a five element vector showing vertical range frame limits, described below
%
% Outputs:
%
% * *main_axis* - handle to the modified axis
%
function [main_axis] = rangeframe(main_axis, x_frame, y_frame)
%% get current axis properties
dims = get(main_axis, 'Position');
    left   = dims(1);           % break these out for clarity
    bottom = dims(2);           % (I wish I knew a better way...)
    width  = dims(3);
    height = dims(4);
x_axis_lims = get(main_axis, 'XLim');
    x_axis_min = x_axis_lims(1);
    x_axis_max = x_axis_lims(2);
y_axis_lims = get(main_axis, 'YLim');
    y_axis_min = y_axis_lims(1);
    y_axis_max = y_axis_lims(2);

%% draw range frame

% set tick length for middle value
l_tick = .01;

% map x_frame array to normalized figure coordinates
x_frame_mapped = (x_frame-x_axis_min)/(x_axis_max-x_axis_min) * width + left;

% map y_frame array to normalized figure coordinates
y_frame_mapped = (y_frame-y_axis_min)/(y_axis_max-y_axis_min) * height + bottom;

% draw horizontal axis range frame
annotation('line', [x_frame_mapped(1),x_frame_mapped(2)], [bottom, bottom]);
annotation('line', [x_frame_mapped(3),x_frame_mapped(3)], [bottom-l_tick, bottom+l_tick], 'Linewidth', 2);
annotation('line', [x_frame_mapped(4),x_frame_mapped(5)], [bottom, bottom]);

% draw vertical axis range frame
annotation('line', [left, left], [y_frame_mapped(1),y_frame_mapped(2)]);
annotation('line', [left-l_tick, left+l_tick], [y_frame_mapped(3),y_frame_mapped(3)], 'Linewidth', 2);
annotation('line', [left, left], [y_frame_mapped(4),y_frame_mapped(5)]);

%% Modify the existing axes to improve the visibility of the range frame

% modify the normal axis
set(main_axis, ...
    'color', 'none', ...    % make sure the background is clear (white would also work)
    'xgrid', 'off', ...     % make sure there is no background grid
    'ygrid','off', ...
    'box', 'off', ...       % get rid of upper and right hand axis frame
    'TickDir', 'out', ...   % set the tick marks to point towards numbers to distringuish from center tick
    'TickLength', [l_tick/2,0]);    % make ths ticks short for the same reason

% cover the axis with a mask to hide the lines
mask_axis = copyobj(main_axis, gcf);
set(mask_axis, ...
    'xcolor', 'w', ...      % maxe the axis lines white
    'ycolor','w', ...
    'Xtick',[], ...         % dont print any ticks or numbers
    'Ytick',[]);


text:
function X = naninterp(X)
% Interpolate over NaNs
% See INTERP1 for more info
X(isnan(X)) = interp1(find(~isnan(X)), X(~isnan(X)), find(isnan(X)),'cubic');
return


text:
function out_backprop()
    global config mem;
    current_layer = config.misc.current_layer;
    mem.deltas{current_layer} = config.DERI_OUT_ACT(config.DERI_COST_FUN(mem.output, mem.GT_output));
    config.EXPAND_DELTA_OUT();
    config.misc.current_layer = config.misc.current_layer - 1;
end

author: jimmy-ren


text:
%% Bivariate Kernel Density Estimation Demonstration
% Several examples show how to use the gkdeb2 function.

%% Distribution with unbounded support
% Distribution with four peaks.
x = [randn(100,1), (randn(100,1)-10)*2;
      randn(100,1)+10, randn(100,1);
      randn(100,1)+10, (randn(100,1)-10)*2;
      randn(100,2)];
gkde2(x);

%% Distribution with upper and lower bounds: uniform distribution
clear
x=rand(10000,2);
% PDF with bounded support
p.xylim=[0 0 1 1];
gkde2(x,p);
% Compare with unbounded PDF estimate
figure;
gkde2(x);

%% Distribution with lower bound only: exponential distribution
clear
x=-log(rand(1000,2));
% PDF with bounded support
p.xylim=[0 0 Inf Inf];
gkde2(x,p);
% Compare with unbounded PDF estimate
figure
gkde2(x);

%% Distribution with lower bound only: log-normal distribution
clear
x=exp(randn(1000,2));
% PDF with bounded support
p.xylim=[0 0 Inf Inf];
gkde2(x,p);
% Compare with unbounded PDF estimate
figure
gkde2(x);

%% Distribution with lower bound only: chi-square distribution
clear
x=randn(1000,2).^2;
% PDF with bounded support
p.xylim=[0 0 Inf Inf];
p.alpha=0.95;
gkde2(x,p);
% Compare with unbounded PDF estimate
figure
p1.alpha=0.95;
gkde2(x,p1);

%%  Distribution with lower bound only: Rayleigh distribution
clear
x=sqrt(randn(2,1000).^2 + randn(2,1000).^2);
% PDF with bounded support
p.xylim=[0 0 Inf Inf];
p.alpha=0.95;
gkde2(x,p);
% Compare with unbounded PDF estimate
figure
p1.alpha=0.95;
gkde2(x,p1);


text:
function prob_test0965 ( )

%*****************************************************************************80
%
%% TEST0965 tests LEVY_CDF, LEVY_CDF_INV, LEVY_PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0965\n' );
  fprintf ( 1, '  For the Levy PDF:\n' );
  fprintf ( 1, '  LEVY_CDF evaluates the CDF;\n' );
  fprintf ( 1, '  LEVY_CDF_INV inverts the CDF.\n' );
  fprintf ( 1, '  LEVY_PDF evaluates the PDF;\n' );

  a = 1.0;
  b = 2.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A =             %14f\n', a );
  fprintf ( 1, '  PDF parameter B =             %14f\n', b );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '       X            PDF           CDF            CDF_INV\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    [ x, seed ] = levy_sample ( a, b, seed );

    pdf = levy_pdf ( x, a, b );

    cdf = levy_cdf ( x, a, b );

    x2 = levy_cdf_inv ( cdf, a, b );

    fprintf ( 1, ' %14f  %14f  %14f  %14f\n', x, pdf, cdf, x2 );

  end

  return
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% DISCRETE VARIANCE SWAP/OPTION PRICER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Descritpion: Script to Price Barrier options in Levy Models
%              using the PROJ method
% Author:      Justin Kirkby
% References:  (1) A General Framework for discretely sampled realized
%              variance derivatives in stocahstic volatility models with
%              jumps, EJOR, 2017
%              (2) Efficient Option Pricing By Frame Duality with The Fast
%              Fourier Transform, SIAM J. Financial Math., 2015
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[folder, name, ext] = fileparts(which( mfilename('fullpath')));
cd(folder);
addpath('../RN_CHF')
addpath('../Helper_Functions')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Step 1) CHOOSE CONTRACT/GENERAL PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
K    = 0.01;  %Strike            %NOTE: no error handling in place for extreme values of W (increase grid if strike falls outside)
r    = .05;  %Interest rate
q    = .00;  %dividend yield
T    = 1;    %Time (in years)
M    = 252;  %number of discrete monitoring points

contract = 1;  %1 for for variance swap, 3 for variance call  (other contracts not yet coded)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Step 2) CHOOSE MODEL PARAMETERS (Levy Models)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

model = 1;   % (CHOOSE) See Models Below (e.g. model 1 is Black Scholes), and choose specific params
params = {};

if model == 1 %BSM (Black Scholes Merton)
    params.sigmaBSM = 0.18;    %CHOOSE

elseif model == 2 %CGMY
    params.C  = 0.02;
    params.G  = 5;
    params.MM = 15;
    params.Y  = 1.2;

elseif model == 3 %NIG
    params.alpha = 15;
    params.beta  = -5;
    params.delta = 0.5;

elseif model == 4 %MJD (Merton Jump Diffusion)
    params.sigma  = 0.12;
    params.lam    = 0.4;
    params.muj    = -0.12;
    params.sigmaj = 0.18;

elseif model == 5 %Kou Double Expo
    params.sigma = 0.15;
    params.lam   = 3;
    params.p_up  = 0.2;
    params.eta1  = 25;
    params.eta2  = 10;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Step 3) CHOOSE PROJ PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
UseCumulant = 1;  %Set to 1 to use the cumulant base rule (Approach 1) to determine gridwidth, else used fixed witdth (Approach 2)

%---------------------
% APPROACH 1: Cumulant Based approach for grid width
% (see "Robust Option Pricing with Characteritics Functions and the BSpline Order of Density Projection")
%---------------------
if UseCumulant ==1  %With cumulant based rule, choose N and Alpha (N = 2^(P+Pbar) based on second approach)
    logN  = 10;   %Uses N = 2^logN  gridpoint
    L1 = 14;  % determines grid witdth (usually set L1 = 8 to 15 for Levy, or 18 for Heston)
%---------------------
% APPROACH 2: Manual GridWidth approach
%---------------------
else %Manually specify resolution and Pbar
    P     = 7;  % resolution is 2^P
    Pbar  = 3;  % Determines density truncation grid with, 2^Pbar
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% PRICE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Note: rnCHF is the risk netural CHF, c1,c2,c4 are the cumulants
modelInput = getModelInput(model, T/M, r, q, params);

if UseCumulant ==1  % Choose density truncation width based on cumulants
    alpha = getTruncationAlpha(T, L1, modelInput, model);
else    % Manually supply density truncation width above
    logN = P + Pbar;
    alpha = 2^Pbar/2;
end
N = 2^logN;

tic
price = PROJ_DiscreteVariance_Swaps_Options( N,alpha,M,r,T,K,modelInput.rnCHF,contract);
toc

fprintf('price: %.8f\n', price)

if contract == 1  % In swap case, the price is known analytically (for Levy process)
    c2 = modelInput.c2; c1 = modelInput.c1;
    analytical = c2 + c1^2*T/M;
    fprintf('analytical SWAP price: %.8f \n', analytical)
    fprintf('----------------------------\n')
    fprintf('Error: %.2e\n', analytical-price)
end


%==============================================================================
% (c) Jan Modersitzki 2011/04/26, see FAIR.2 and FAIRcopyright.m.
% http://www.mic.uni-luebeck.de/people/jan-modersitzki.html
%
% For an extended documentation, see:
% Jan Modersitzki. FAIR: Flexible Algorithms for Image Registration, SIAM, 2009.
% http://www.siam.org/books/fa06/
%
% KERNEL/MATRIXFREE
%
% Contents of FAIR MATRIXFREE
%
%   MGsolver - calls MG solver
%   mfvcyce  - MG solver
%   mfPu.m   - matrix free prolongation operator
%   mfJacobi - MG smoother
%   mfAy     - MG matrix vector operation (M+alpha*B'*B)
%   mfBy     - matrix free B*y for discrete regularizer
%   testMatrixFree - test the files in this folder
%==============================================================================

function debit = contents
if nargout == 0, help(mfilename); return; end;

debit = {
    'contents.m'

    'MGsolver.m'
    'mfvcycle.m'
    'mfAy.m'
    'mfJacobi.m'
    'mfPu.m'

    'testMatrixFree.m'
    };
%==============================================================================