text:
function triangulation_order6_contour ( prefix )

%*****************************************************************************80
%
%% TRIANGULATION_ORDER6_CONTOUR contour plots data at the nodes of a six node triangle.
%
%  Discussion:
%
%    This program can read data files defining a set of nodes, the triangulation of
%    those nodes using 6 node triangles, and the value of a scalar at the nodes.
%    It then rearranges the data to set up calls to MATLAB that display the
%    triangulation, and color contour plots of the scalar.
%
%    The MATLAB plotting routines require a 3-node triangulation.  This program
%    shows two ways to turn a 6-node triangulation into a 3-node triangulation.
%    The default MATLAB plot of data on a triangular mesh uses constant color
%    over each triangle.  After displaying these plots, the program shows how to
%    get a much nicer plot using interpolated color.
%
%  Usage:
%
%    triangulation_order6_contour ( 'prefix' )
%
%    where
%
%    * 'prefix'_nodes.txt contains the node coordinates;
%    * 'prefix'_elements.txt contains the element definitions.
%    * 'prefix'_values.txt contains the nodal values
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 October 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string PREFIX, the common file prefix.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGULATION_ORDER6_CONTOUR:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot a scalar defined on a 6-node triangle triangulation.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This program expects to find three files to read:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  * "nodes.txt",     the node file,\n' );
  fprintf ( 1, '  * "elements.txt",  the element file,\n' );
  fprintf ( 1, '  * "solution.txt",  the solution file,\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  It reads the files, and makes two plots of the solution.\n' );
%
%  The command line argument is the common filename prefix.
%
  if ( nargin < 1 )

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

    prefix = input ( ...
      'Please enter the filename prefix:' );

  end
%
%  Create the filenames.
%
  node_filename = strcat ( prefix, '_nodes.txt' );
  element_filename = strcat ( prefix, '_elements.txt' );
  value_filename = strcat ( prefix, '_values.txt' );
%
%  Read the node data.
%
  [ dim_num, node_num ] = r8mat_header_read ( node_filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the header of "%s".', node_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension DIM_NUM = %d\n', dim_num );
  fprintf ( 1, '  Number of points NODE_NUM = %d\n', node_num );

  if ( dim_num ~= 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_ORDER6_CONTOUR - Fatal error!\n' );
    fprintf ( 1, '  Dataset must have spatial dimension 2.\n' );
    error ( 'TRIANGULATION_ORDER6_CONTOUR - Fatal error!' );
  end

  node_xy = r8mat_data_read ( node_filename, dim_num, node_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the data in "%s".\n', node_filename );

  r8mat_transpose_print_some ( dim_num, node_num, node_xy, 1, 1, dim_num, 5, ...
    '  First 5 nodes:' );
%
%  Read the element data.
%
  [ triangle_order, triangle_num ] = i4mat_header_read ( ...
    element_filename );

  if ( triangle_order ~= 3 && triangle_order ~= 6 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_ORDER6_CONTOUR - Fatal error!\n' );
    fprintf ( 1, '  Data is not for a 3-node or 6-node triangulation.\n' );
    error ( 'TRIANGULATION_ORDER6_CONTOUR - Fatal error!' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the header of "%s".\n', ...
    element_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Triangle order = %d\n', triangle_order );
  fprintf ( 1, '  Number of triangles TRIANGLE_NUM  = %d\n', ...
    triangle_num );

  triangle_node = i4mat_data_read ( element_filename, ...
    triangle_order, triangle_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the data in "%s".\n', element_filename );

  i4mat_transpose_print_some ( triangle_order, triangle_num, ...
    triangle_node, 1, 1, triangle_order, 10, ...
    '  First 10 elements:' );
%
%  Detect and correct 0-based indexing.
%
  triangle_node = mesh_base_one ( node_num, triangle_order, triangle_num, ...
    triangle_node );
%
%  Read the values.
%
  [ value_dim, value_num ] = r8mat_header_read ( value_filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the header of "%s".', value_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension = %d\n', value_dim );
  fprintf ( 1, '  Number of values  = %d\n', value_num );

  if ( value_dim ~= 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_ORDER6_CONTOUR - Fatal error!\n' );
    fprintf ( 1, '  VALUE data must be scalar.\n' );
    error ( 'TRIANGULATION_ORDER6_CONTOUR - Fatal error!' );
  end

  value = r8mat_data_read ( value_filename, value_dim, value_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the data in "%s".\n', value_filename );

  r8mat_transpose_print_some ( value_dim, value_num, value, 1, 1, value_dim, 5, ...
    '  First 5 values:' );
%
%  Make a 3-node triangulation by discarding the midside nodes.
%
  t2 = triangle_node(1:3,:)';
%
%  Display the mesh.
%
  trimesh ( t2, node_xy(1,:), node_xy(2,:), 'Color', 'blue' );

  xlabel ( 'X', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16 );

  ylabel ( 'Y', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16, 'Rotation', 0 );

  title ( 'Dropping midside nodes', 'FontName', 'Helvetica', 'FontWeight', ...
    'bold', 'FontSize', 16 );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Press return...\n' );
  pause
%
%  Display the solution on this crude mesh.
%
  trisurf ( t2, node_xy(1,:), node_xy(2,:), value )

  xlabel ( 'X', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16 );

  ylabel ( 'Y', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16, 'Rotation', 0 );

  zlabel ( 'U', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16, 'Rotation', 0 );

  title ( 'Contours of U(X,Y)', 'FontName', 'Helvetica', 'FontWeight', ...
    'bold', 'FontSize', 16 );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Press return...\n' );
  pause
%
%  Make a 3-node triangulation by breaking each 6 node triangle
%  into 4 3-node triangles.
%
  t3 = triangulation_order6_to_order3 ( triangle_num, triangle_node );
  t3 = t3';
%
%  Display the mesh.
%
  trimesh ( t3, node_xy(1,:), node_xy(2,:), 'Color', 'blue' );

  xlabel ( 'X', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16 );

  ylabel ( 'Y', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16, 'Rotation', 0 );

  title ( 'Splitting 6-Node Triangles', 'FontName', 'Helvetica', 'FontWeight', ...
    'bold', 'FontSize', 16 );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Press return...\n' );
  pause
%
%  Display the solution on this finer mesh.
%
  trisurf ( t3, node_xy(1,:), node_xy(2,:), value );
  xlabel ( 'X', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16 );

  ylabel ( 'Y', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16, 'Rotation', 0 );

  zlabel ( 'U', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16, 'Rotation', 0 );

  title ( 'Contours of U(X,Y)', 'FontName', 'Helvetica', 'FontWeight', ...
    'bold', 'FontSize', 16 );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Press return...\n' );
  pause
%
%  Make a nicer plot on the finer mesh by using color interpolation.
%
  trisurf ( t3, node_xy(1,:), node_xy(2,:), value, 'FaceColor', 'interp', ...
    'EdgeColor', 'interp' )

  xlabel ( 'X', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16 );

  ylabel ( 'Y', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16, 'Rotation', 0 );

  zlabel ( 'U', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16, 'Rotation', 0 );

  title ( 'Contours of U(X,Y)', 'FontName', 'Helvetica', 'FontWeight', ...
    'bold', 'FontSize', 16 );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Press return...\n' );
  pause
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGULATION_ORDER6_CONTOUR:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

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

  return
end
function column_num = file_column_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_COLUMN_COUNT counts the columns in the first line of a file.
%
%  Discussion:
%
%    The file is assumed to be a simple text file.
%
%    Most lines of the file are presumed to consist of COLUMN_NUM words,
%    separated by spaces.  There may also be some blank lines, and some
%    comment lines, which have a "#" in column 1.
%
%    The routine tries to find the first non-comment non-blank line and
%    counts the number of words in that line.
%
%    If all lines are blanks or comments, it goes back and tries to analyze
%    a comment line.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the file.
%
%    Output, integer COLUMN_NUM, the number of columns in the file.
%
  FALSE = 0;
  TRUE = 1;
%
%  Open the file.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_COLUMN_COUNT - Error!' );
  end
%
%  Read one line, but skip blank lines and comment lines.
%  Use FGETL so we drop the newline character!
%
  got_one = FALSE;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( s_len_trim ( line ) == 0 )

    elseif ( line(1) == '#' )

    else
      got_one = TRUE;
      break;
    end

  end

  fclose ( input_unit );

  if ( got_one == FALSE )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Warning!\n' );
    fprintf ( 1, '  The file does not seem to contain any data.\n' );
    column_num = -1;
    return;
  end

  column_num = s_word_count ( line );

  return
end
function row_num = file_row_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_ROW_COUNT counts the number of row records in a file.
%
%  Discussion:
%
%    Each input line is a "RECORD".
%
%    The records are divided into three groups:
%
%    * BLANK LINES (nothing but blanks)
%    * COMMENT LINES (begin with a '#')
%    * DATA RECORDS (anything else)
%
%    The value returned by the function is the number of data records.
%
%    By the way, if the MATLAB routine FGETS is used, instead of
%    FGETL, then the variable LINE will include line termination
%    characters, which means that a blank line would not actually
%    have zero characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the input file.
%
%    Output, integer ROW_NUM, the number of rows found.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_ROW_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_ROW_COUNT - Error!' );
  end

  blank_num = 0;
  comment_num = 0;
  row_num = 0;

  record_num = 0;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    record_num = record_num + 1;
    record_length = s_len_trim ( line );

    if ( record_length <= 0 )
      blank_num = blank_num + 1;
    elseif ( line(1) == '#' )
      comment_num = comment_num + 1;
    else
      row_num = row_num + 1;
    end

  end

  fclose ( input_unit );

  return
end
function table = i4mat_data_read ( input_filename, m, n )

%*****************************************************************************80
%
%% I4MAT_DATA_READ reads data from an I4MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns in the data.
%
%    Output, integer TABLE(M,N), the point coordinates.
%
  table = zeros ( m, n );
%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

  for i = 0 : m
    string = strcat ( string, ' %d' );
  end

  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the input file.\n' );
    error ( 'I4MAT_DATA_READ - Error!' );
  end

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'I4MAT_DATA_READ - Error!\n' );
      fprintf ( 1, '  End of input while reading data.\n' );
      error ( 'I4MAT_DATA_READ - Error!' );
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      [ x, count ] = sscanf ( line, string );

      if ( count == m )
        i = i + 1;
        table(1:m,i) = x(1:m);
      end

    end

  end

  fclose ( input_unit );

  return
end
function [ m, n ] = i4mat_header_read ( input_filename )

%*****************************************************************************80
%
%% I4MAT_HEADER_READ reads the header from an I4MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer M, the spatial dimension.
%
%    Output, integer N, the number of points.
%
  m = file_column_count ( input_filename );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data columns in\n' );
    fprintf ( 1, '  the file %s.\n', input_filename );
  end

  n = file_row_count ( input_filename );

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data rows in\n' );
    fprintf ( 1, '  the file %s\n', input_filename );
  end

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

%*****************************************************************************80
%
%% I4MAT_TRANSPOSE_PRINT_SOME prints some of an I4MAT, transposed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    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, an optional title.
%
  incx = 10;

  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d  ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );
    fprintf ( 1, '\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d  ', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%7d  ', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

  return
end
function element_node = mesh_base_one ( node_num, element_order, ...
  element_num, element_node )

%*****************************************************************************80
%
%% MESH_BASE_ONE ensures that the element definition is one-based.
%
%  Discussion:
%
%    The ELEMENT_NODE array contains nodes indices that form elements.
%    The convention for node indexing might start at 0 or at 1.
%    Since a MATLAB program will naturally assume a 1-based indexing, it is
%    necessary to check a given element definition and, if it is actually
%    0-based, to convert it.
%
%    This function attempts to detect 0-based node indexing and correct it.
%
%    Thanks to Feifei Xu for pointing out that I was subtracting 1 when I
%    should have been adding 1!  29 November 2012.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 November 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer ELEMENT_ORDER, the order of the elements.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input/output, integer ELEMENT_NODE(ELEMENT_ORDE,ELEMENT_NUM), the element
%    definitions.
%
  node_min = min ( min ( element_node(1:element_order,1:element_num) ) );
  node_max = max ( max ( element_node(1:element_order,1:element_num) ) );

  if ( node_min == 0 && node_max == node_num - 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MESH_BASE_ONE:\n' );
    fprintf ( 1, '  The element indexing appears to be 0-based!\n' );
    fprintf ( 1, '  This will be converted to 1-based.\n' );
    element_node(1:element_order,1:element_num) = ...
      element_node(1:element_order,1:element_num) + 1;
  elseif ( node_min == 1 && node_max == node_num )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MESH_BASE_ONE:\n' );
    fprintf ( 1, '  The element indexing appears to be 1-based!\n' );
    fprintf ( 1, '  No conversion is necessary.\n' );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MESH_BASE_ONE - Warning!\n' );
    fprintf ( 1, '  The element indexing is not of a recognized type.\n' );
    fprintf ( 1, '  NODE_MIN = %d\n', node_min );
    fprintf ( 1, '  NODE_MAX = %d\n', node_max );
    fprintf ( 1, '  NODE_NUM = %d\n', node_num );
  end

  return
end
function table = r8mat_data_read ( input_filename, m, n )

%*****************************************************************************80
%
%% R8MAT_DATA_READ reads data from an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns of data.
%
%    Output, real TABLE(M,N), the point coordinates.
%
  table = zeros ( m, n );
%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

  for i = 0 : m
    string = strcat ( string, ' %f' );
  end

  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the file.\n' );
    error ( 'R8MAT_DATA_READ - Error!' );
  end

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      [ x, count ] = sscanf ( line, string );

      if ( count == m )
        i = i + 1;
        table(1:m,i) = x(1:m);
      end

    end

  end

  fclose ( input_unit );

  return
end
function [ m, n ] = r8mat_header_read ( input_filename )

%*****************************************************************************80
%
%% R8MAT_HEADER_READ reads the header from an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer M, the spatial dimension.
%
%    Output, integer N, the number of points.
%
  m = file_column_count ( input_filename );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data columns in\n' );
    fprintf ( 1, '  the file %s.\n', input_filename );
  end

  n = file_row_count ( input_filename );

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data rows in\n' );
    fprintf ( 1, '  the file %s\n', input_filename );
  end

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

%*****************************************************************************80
%
%% R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real 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, an optional title.
%
  incx = 5;

  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d       ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d ', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%12f', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

  return
end
function len = s_len_trim ( s )

%*****************************************************************************80
%
%% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LEN, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end

  return
end
function word_num = s_word_count ( s )

%*****************************************************************************80
%
%% S_WORD_COUNT counts the number of "words" in a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be examined.
%
%    Output, integer WORD_NUM, the number of "words" in the string.
%    Words are presumed to be separated by one or more blanks.
%
  FALSE = 0;
  TRUE = 1;

  word_num = 0;
  s_length = length ( s );

  if ( s_length <= 0 )
    return;
  end

  blank = TRUE;

  for i = 1 : s_length

    if ( s(i) == ' ' )
      blank = TRUE;
    elseif ( blank == TRUE )
      word_num = word_num + 1;
      blank = FALSE;
    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 triangle2 = triangulation_order6_to_order3 ( tri_num1, triangle1 )

%*****************************************************************************80
%
%% TRIANGULATION_ORDER6_TO_ORDER3 linearizes a quadratic triangulation.
%
%  Discussion:
%
%    A quadratic triangulation is assumed to consist of 6-node triangles,
%    as in the following:
%
%    11-12-13-14-15
%     |\    |\    |
%     | \   | \   |
%     6  7  8  9 10
%     |   \ |   \ |
%     |    \|    \|
%     1--2--3--4--5
%
%   This routine rearranges information so as to define the 3-node
%   triangulation:
%
%    11-12-13-14-15
%     |\ |\ |\ |\ |
%     | \| \| \| \|
%     6--7--8--9-10
%     |\ |\ |\ |\ |
%     | \| \| \| \|
%     1--2--3--4--5
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer TRI_NUM1, the number of triangles in the quadratic
%    triangulation.
%
%    Input, integer TRIANGLE1(6,TRI_NUM1), the quadratic triangulation.
%
%    Output, integer TRIANGLE2(3,4*TRI_NUM1), the linear triangulation.
%
  tri2 = 0;

  triangle2 = zeros ( 3, 4 * tri_num );

  for tri1 = 1 : tri_num1

    n1 = triangle1(1,tri1);
    n2 = triangle1(2,tri1);
    n3 = triangle1(3,tri1);
    n4 = triangle1(4,tri1);
    n5 = triangle1(5,tri1);
    n6 = triangle1(6,tri1);

    tri2 = tri2 + 1;
    triangle2(1:3,tri2) = [ n1, n4, n6 ]';
    tri2 = tri2 + 1;
    triangle2(1:3,tri2) = [ n2, n5, n4 ]';
    tri2 = tri2 + 1;
    triangle2(1:3,tri2) = [ n3, n6, n5 ]';
    tri2 = tri2 + 1;
    triangle2(1:3,tri2) = [ n4, n5, n6 ]';

  end

  return
end


text:
function reducedModel=removeMets(model,metsToRemove,isNames,removeUnusedRxns,removeUnusedGenes,removeUnusedComps)
% removeMets
%   Deletes a set of metabolites from a model
%
%   model             a model structure
%   metsToRemove      either a cell array of metabolite IDs, a logical vector
%                     with the same number of elements as metabolites in the model,
%                     of a vector of indexes to remove
%   isNames           true if the supplied mets represent metabolite names
%                     (as opposed to IDs). This is a way to delete
%                     metabolites in several compartments at once without
%                     knowing the exact IDs. This only works if metsToRemove
%                     is a cell array (opt, default false)
%   removeUnusedRxns  remove reactions that are no longer in use (opt,
%                     default false)
%   removeUnusedGenes remove genes that are no longer in use (opt,
%                     default false)
%   removeUnusedComps remove compartments that are no longer in use (opt,
%                     default false)
%
%   reducedModel      an updated model structure
%
%   Usage: reducedModel=removeMets(model,metsToRemove,isNames,...
%           removeUnusedRxns,removeUnusedGenes,removeUnusedComps)
if ~islogical(metsToRemove) && ~isnumeric(metsToRemove)
    metsToRemove=convertCharArray(metsToRemove);
end

if nargin<3
    isNames=false;
end

if nargin<4
    removeUnusedRxns=false;
end

if nargin<5
    removeUnusedGenes=false;
end

if nargin<6
    removeUnusedComps=false;
end

%Check that metsToRemove is a cell array
if isNames==true && ~iscell(metsToRemove)
    error('Must supply a cell array of strings if isNames=true');
end

reducedModel=model;

if isNames==false
    indexesToDelete=getIndexes(model,metsToRemove,'mets');
else
    indexesToDelete=[];
    for i=1:numel(metsToRemove)
        indexesToDelete=[indexesToDelete;find(strcmp(metsToRemove(i),model.metNames))];
    end
end

%Remove metabolites
if ~isempty(indexesToDelete)
    reducedModel.mets(indexesToDelete)=[];
    reducedModel.S(indexesToDelete,:)=[];
    if isfield(reducedModel,'b')
        reducedModel.b(indexesToDelete,:)=[];
    end
    if isfield(reducedModel,'metNames')
        reducedModel.metNames(indexesToDelete)=[];
    end
    if isfield(reducedModel,'metComps')
        reducedModel.metComps(indexesToDelete)=[];
    end
    if isfield(reducedModel,'inchis')
        reducedModel.inchis(indexesToDelete)=[];
    end
    if isfield(reducedModel,'metSmiles')
        reducedModel.metSmiles(indexesToDelete)=[];
    end
    if isfield(reducedModel,'metFormulas')
        reducedModel.metFormulas(indexesToDelete)=[];
    end
    if isfield(reducedModel,'metMiriams')
        reducedModel.metMiriams(indexesToDelete)=[];
    end
    if isfield(reducedModel,'unconstrained')
        reducedModel.unconstrained(indexesToDelete)=[];
    end
    if isfield(reducedModel,'metFrom')
        reducedModel.metFrom(indexesToDelete)=[];
    end
    if isfield(reducedModel,'metCharges')
        reducedModel.metCharges(indexesToDelete)=[];
    end
end

%Remove unused reactions
if removeUnusedRxns==true
    %Get unused reactions
    [~, a]=find(reducedModel.S);
    rxnsToRemove=1:numel(reducedModel.rxns);
    rxnsToRemove(a)=[];
    reducedModel=removeReactions(reducedModel,rxnsToRemove,false,removeUnusedGenes);
end

%Remove unused compartments
if removeUnusedComps==true
    oldComps=reducedModel.comps;
    I=ismember(1:numel(oldComps),reducedModel.metComps);
    if ~all(I)
        reducedModel.comps(~I)=[];
        reducedModel.compNames(~I)=[];
        if isfield(reducedModel,'compOutside')
            reducedModel.compOutside(~I)=[];
        end
        if isfield(reducedModel,'compMiriams')
            reducedModel.compMiriams(~I)=[];
        end
        [~, J]=ismember(oldComps(reducedModel.metComps),reducedModel.comps);
        reducedModel.metComps=J;
    end
end
end

author: SysBioChalmers

repo: RAVEN


text:
%% INS惯导解算
clear;
clc;
close all;

%%
Fs = 100;  %采样频率
N = 1000; %采样次数

dt = 1 / Fs;
gyr = [0.01, 0.02, 0.03]; %单位rad
acc = [0, 0, 9.8]; % 单位m/s^(2)

% 捷联惯导解算
p = zeros(3, 1);
v = zeros(3, 1);
q= [1 0 0 0]';

for i=1:N
    [p ,v, q] = ch_nav_equ_local_tan(p, v, q, acc', gyr' , dt, [0, 0, -9.8]');
    h_pos(i,:) = p;
    h_vel(i,:) = v;
    h_eul(i,:) = ch_q2eul(q);
end

figure;
subplot(2,2,1);
ch_plot_pos3d(h_pos);
subplot(2,2,2);
ch_plot_pos2d(h_pos);
subplot(2,2,3);
ch_plot_att(h_eul);

fprintf('纯积分测试: 陀螺(rad/s):%.3f %.3f %.3f\n', gyr(1), gyr(2), gyr(3));
fprintf('纯积分测试: 加计(m/s^(2)):%.3f %.3f %.3f\n', acc(1), acc(2), acc(3));

fprintf('解算:%d次 总时间:%.3fs\n', N, N /Fs);
fprintf('最终误差(m): %.3f %.3f %.3f\n', h_pos(end, 1),  h_pos(end, 2),  h_pos(end, 3));

author: yandld


text:
function bd=dfun(t,b);
bd=zeros(3,1);
bd(1)=-b(1)+b(2)+b(3)+b(1)*((b(1)^2)+(b(2)^2)+(b(3)^2));
bd(2)=-b(1)-b(2)+b(3)+b(2)*((b(1)^2)+(b(2)^2)+(b(3)^2));
bd(3)=-b(1)-b(2)-b(3)+b(3)*((b(1)^2)+(b(2)^2)+(b(3)^2));

author: sfvsfv

repo: Mathematical-modeling


text:
function endpointsAll = gt2endpoints(annolistFilename)

if isstruct(annolistFilename)
    annolist = annolistFilename;
else
    load(annolistFilename);
end
%if (~exist('annolist','var'))
%    annolist = single_person_annolist;
%end
endpointsAll = cell(length(annolist),1);

skel = [0 1; 1 2; 3 4; 4 5; 6 7; 8 9; 10 11; 11 12; 13 14; 14 15];

for imgidx = 1:length(annolist)
    endpoints = nan(10,4);
    points = annolist(imgidx).annorect.annopoints.point;
    for sidx = 1:size(skel,1)
        p1 = util_get_annopoint_by_id(points,skel(sidx,1));
        p2 = util_get_annopoint_by_id(points,skel(sidx,2));
        if (~isempty(p1) && ~isempty(p2))
            endpoints(sidx,:) = [p1.x p1.y p2.x p2.y];
        end
    end
    endpointsAll{imgidx} = endpoints;
end

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
% This file is part of the following project:
% Oliver Parson, Siddhartha Ghosh, Mark Weal, Alex Rogers.
% Non-intrusive Load Monitoring using Prior Models of General Appliance Types.
% In: 26th AAAI Conference on Artificial Intelligence. Toronto, Canada. 2012.
% Code available for download: https://sites.google.com/site/oliparson/phd-work/research-files/aaai-2012-code.zip?attredirects=0
% Copyright: Oliver Parson et al., University of Southhampton, 2012.

% Modified by Romano Cicchetti, ETH Zurich, in the context of the NILM-Eval project

function [result] = parsonAppliance(evaluation_and_training_days, setup, fid)

        %load parameters
        dataset = setup.dataset;
        household = setup.household;
        granularity = setup.granularity;
        filtering = setup.filtering;
        if isfield(setup,'filtering2')
            filtering2 = setup.filtering2;
        else
            filtering2 = 'noFiltering';
        end
        filtWindowLength = setup.filtWindowLength;
        filtLRatio = setup.filtLRatio;
        phases = setup.phases;

        applianceName = setup.applianceName;
        state_means = [setup.meanOff, setup.meanOn] ;
        state_covs = [setup.varOff, setup.varOn];
        transition_matrix = [1-setup.transOn, setup.transOn;
                            setup.transOff, 1-setup.transOff];
        likelihoodThreshold = setup.likThresh;
        numOfWindows = setup.numOfWindows;
        windowLength = setup.windowLength;
        trainingType = setup.trainingType;

        windowLength = windowLength/granularity;

        % get training and test data
        evaluation_days = evaluation_and_training_days{1};
        training_days = evaluation_and_training_days{2};
        applianceID = getApplianceID(applianceName);
        if strcmpi(phases, 'indv')
            phase_of_appliance = getPhase(household, applianceID);
            read_smartmeter_option = strcat('powerl', num2str(phase_of_appliance));
            test_data = read_smartmeter_data(dataset, num2str(household, '%02d'), evaluation_days, granularity, read_smartmeter_option);
            training_data = read_smartmeter_data(dataset, num2str(household, '%02d'), training_days, granularity, read_smartmeter_option);
        elseif strcmpi(phases, 'all')
            test_data = read_smartmeter_data(dataset, num2str(household, '%02d'), evaluation_days, granularity, 'powerallphases');
            training_data = read_smartmeter_data(dataset, num2str(household, '%02d'), training_days, granularity, 'powerallphases');
        else
            error('wrong value for phases parameter');
        end
        plug_data  = read_plug_data(dataset, household, applianceID, training_days, granularity);

        % apply filtering method to consumption data
        if ~strcmpi(filtering, 'noFiltering')
            function_handle = str2func(filtering);
            if strcmpi(filtering, 'totalVar')
                test_data = function_handle(test_data, filtLRatio);
                training_data = function_handle(training_data, filtLRatio);
                plug_data = function_handle(plug_data, filtLRatio);
            elseif strcmpi(filtering, 'edgeEnhance5')
                test_data = function_handle(test_data);
                training_data = function_handle(training_data);
                plug_data = function_handle(plug_data);
            else
                test_data = function_handle(test_data, filtWindowLength);
                plug_data = function_handle(plug_data, filtWindowLength);
            end
        end
        if ~strcmpi(filtering2, 'noFiltering')
            function_handle = str2func(filtering2);
            if strcmpi(filtering2, 'totalVar')
                test_data = function_handle(test_data, filtLRatio);
                training_data = function_handle(training_data, filtLRatio);
                plug_data = function_handle(plug_data, filtLRatio);
            elseif strcmpi(filtering2, 'edgeEnhance5')
                test_data = function_handle(test_data);
                training_data = function_handle(training_data);
                plug_data = function_handle(plug_data);
            else
                test_data = function_handle(test_data, filtWindowLength);
                training_data = function_handle(training_data, filtWindowLength);
                plug_data = function_handle(plug_data, filtWindowLength);
            end
        end

        % calculate diff model params
        init = ones(1, length(state_means)) / length(state_means);
        num_states = length(init);
        idx = repmat(1:num_states,num_states,1);
        emit_mean = state_means(idx) - state_means(idx');
        emit_cov = state_covs(idx) + state_covs(idx');

        % make difference hmm
        prior_dhmm_bnet = make_dhmm(init, ...
                          emit_mean, ...
                          emit_cov, ...
                          state_means, ...
                          state_covs, ...
                          transition_matrix);

        prior_hmm_bnet = make_hmm(init, ...
                         state_means, ...
                         state_covs, ...
                         transition_matrix);

        % data evidence
        diffs = [0 diff(test_data)];
        evidence = {};
        evidence(2,:) = num2cell(diffs);

        %create engine
        engine = smoother_engine(jtree_2TBN_inf_engine(prior_dhmm_bnet));

        %select windows of training data
        if strcmpi(trainingType, 'noTraining')

        elseif strcmpi(trainingType, 'aggregateTraining')
                [training_windows, ~] = find_training_ranges_generic(training_data, ...
                            windowLength, prior_dhmm_bnet, numOfWindows);

        elseif strcmpi(trainingType, 'plugTraining')
            plug_data_with_noise = awgn(plug_data, 25, 'measured');
            [training_windows, ~] = find_training_ranges_generic(plug_data_with_noise, ...
                            windowLength, prior_dhmm_bnet, numOfWindows);
        end

        %train models
        if strcmpi(trainingType, 'noTraining')
            trained_hmm_bnet = prior_hmm_bnet;
            trained_dhmm_bnet = prior_dhmm_bnet;
        elseif strcmpi(trainingType, 'aggregateTraining') || strcmpi(trainingType, 'plugTraining')
            [trained_dhmm_bnet, ~] = learn_params_generic(prior_dhmm_bnet, diff(training_windows'));
                emission3 = struct(trained_dhmm_bnet.CPD{4});
                dHMM_means = permute(emission3.mean, [2,3,1]);
                dHMM_covs = permute(emission3.cov, [3,4,1,2]);
            try
                [trained_hmm_bnet, ~] = learn_params_generic(prior_hmm_bnet, training_windows');
            catch err
                trained_hmm_bnet = prior_hmm_bnet;
            end
            trained_dhmm_bnet.CPD{1} = tabular_CPD(trained_dhmm_bnet, 1, 'CPT', ones(1, length(init)) / length(init));
            trained_dhmm_bnet.CPD{2} = trained_hmm_bnet.CPD{2};
            trained_dhmm_bnet.CPD{3} = prior_dhmm_bnet.CPD{3};
        end
        hmm_emissions = struct(trained_hmm_bnet.CPD{2});

        % infer power (viterbi)
        evidence(3,:) = num2cell(test_data);
        [mpe] = my_viterbi_diff(trained_dhmm_bnet, evidence, 1, likelihoodThreshold);
        means = hmm_emissions.mean;
        result = struct;
        result.consumption = means(cell2mat(mpe(1,:)));
        result.appliance_names = {applianceName};

        % write to summary text file
        fprintf(fid, 'trained HMM means:\n');
        fprintf(fid, [repmat('%f\t', 1, size(hmm_emissions.mean(:),2)) '\n'], hmm_emissions.mean(:));
        fprintf(fid, '\ntrained HMM covs:\n');
        fprintf(fid, [repmat('%f\t', 1, size(hmm_emissions.cov(:),2)) '\n'], hmm_emissions.cov(:));
        if strcmpi(trainingType, 'aggregateTraining') || strcmpi(trainingType, 'plugTraining')
            fprintf(fid, '\ntrained dHMM means:\n');
            fprintf(fid, [repmat('%f\t', 1, size(dHMM_means,2)) '\n'], dHMM_means');
            fprintf(fid, '\ntrained dHMM covs:\n');
            fprintf(fid, [repmat('%f\t', 1, size(dHMM_covs,2)) '\n'], dHMM_covs');
            fprintf(fid, '\n\n');
        else
            fprintf(fid, '\n');
        end
end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
% This file is part of the Matlab Toolbox for Dimensionality Reduction v0.7.2b.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, 2010
% University California, San Diego / Delft University of Technology

function lnst(hs,he)
ls=get(he,'value');

switch ls
    case 1
        set(hs,'LineStyle','-','marker','x','color','b');
    case 2
        set(hs,'LineStyle','-','marker','o','color','r');
    case 3
        set(hs,'LineStyle','none','marker','o','color','b');
    case 4
        set(hs,'LineStyle','-','marker','none','color','g');
    case 5
        set(hs,'LineStyle','--','marker','d','color','b');

end
drawnow;

repo: cluster


text:
function resFileClean = preprocessResult(resFile, seqName, dataDir, force, minvis)
% reads submitted (raw) MOT16 result from .txt
% and removes all boxes that are associated with ambiguous annotations
% such as sitting people, cyclists or mannequins.
% Also removes partially occluded boxes if minvis>0

% resFile='/home/amilan/research/projects/bmtt-dev/code/bae/res/0001/ADL-Rundle-6.txt';
% seqName = 'MOT16-09';
assert(cleanRequired(seqName),'preproccessing should only be done for MOT16/17')

if nargin<4, force=1; end
if nargin<5, minvis=0; end

fprintf('Preprocessing (cleaning) %s...\n',seqName);

% if file does not exist, do nothing
if ~exist(resFile,'file')
    fprintf('Results file does not exist\n');
    resFileClean = [];
    return;
end

[p,f,e]=fileparts(resFile);
cleanDir = [p,filesep,'clean'];
if ~exist(cleanDir, 'dir'), mkdir(cleanDir); end
resFileClean = [cleanDir,filesep,f,e];

% if clean file already exists and no need to redo, skip
if ~force && exist(resFileClean, 'file')
    fprintf('skipping...\n');
    return;
end

% if file empty, just copy it
tf = dir(resFile);
if tf.bytes == 0
    fprintf('Results file empty\n');
    copyfile(resFile,resFileClean);
    return;
end

if nargin<3, dataDir = getDataDir; end

% [seqName, seqFolder, imgFolder, imgExt, F, dirImages] ...
%     = getSeqInfo(seq, dataDir);

[seqName, seqFolder, imgFolder, frameRate, F, imWidth, imHeight, imgExt] ...
    = getSeqInfoFromFile(seqName, dataDir);

% read in result
resRaw = dlmread(resFile);

%%

% read ground truth
gtFolder = [dataDir,seqName,filesep,'gt',filesep];
gtFile = [gtFolder,'gt.txt'];
gtRaw = dlmread(gtFile);

% make sure we have MOT16 ground truth (= 9 columns)
assert(size(gtRaw,2)==9, 'unknown GT format')

% define which classes should be ignored
distractors = {'person_on_vhcl','static_person','distractor','reflection'};

keepBoxes = true(size(resRaw,1),1);

showVis = 0;

td=0.5; % overlap threshold
for t=1:F
    if ~mod(t,100), fprintf('.'); end

    % find all result boxes in this frame
    resInFrame = find(resRaw(:,1)==t); N = length(resInFrame);
    resInFrame = reshape(resInFrame,1,N);

    % find all GT boxes in frame
    GTInFrame = find(gtRaw(:,1)==t); Ngt = length(GTInFrame);
    GTInFrame = reshape(GTInFrame,1,Ngt);

    % compute all overlaps for current frame
    allisects=zeros(Ngt,N);
    g=0;
    for gg=GTInFrame
        g=g+1; r=0;
        bxgt=gtRaw(gg,3); bygt=gtRaw(gg,4); bwgt=gtRaw(gg,5); bhgt=gtRaw(gg,6);
        for rr=resInFrame
            r=r+1;
            bxres=resRaw(rr,3); byres=resRaw(rr,4); bwres=resRaw(rr,5); bhres=resRaw(rr,6);

            if bxgt+bwgt<bxres, continue; end % ignore if no horizontal overlap
            if bxgt>bxres+bwres, continue; end

            if bygt+bhgt<byres, continue; end % ignore if no vertical overlap
            if bygt>byres+bhres, continue; end

            allisects(g,r)=boxiou(bxgt,bygt,bwgt,bhgt,bxres,byres,bwres,bhres);
        end
    end
%     t

    tmpai=allisects;
    tmpai=1-tmpai;
    tmpai(tmpai>td)=Inf;
    [Mtch,Cst]=Hungarian(tmpai);
    [mGT,mRes]=find(Mtch);
%     pause
    nMtch = length(mGT);
    % go through all matches
    for m=1:nMtch
        g=GTInFrame(mGT(m)); % gt box
        gtClassID = gtRaw(g,8);
        gtClassString = classIDToString(gtClassID);

        % if we encounter a distractor, mark to remove box
        if ismember(gtClassString, distractors)
            r = resInFrame(mRes(m)); % result box
            keepBoxes(r) = false;

            if showVis
                bxgt=gtRaw(g,3); bygt=gtRaw(g,4); bwgt=gtRaw(g,5); bhgt=gtRaw(g,6); idgt=gtRaw(g,2);
                bxres=resRaw(r,3); byres=resRaw(r,4); bwres=resRaw(r,5); bhres=resRaw(r,6); idres=resRaw(r,2);

                clf
                im = imread(fullfile(dataDir,seqName,'img1',sprintf('%06d.jpg',t)));
                imshow(im); hold on
                text(50,50,sprintf('%d',t),'color','w')

                % show GT box
%                 text(bxgt,bygt-20,sprintf('%d',idgt),'color','w')
                classString = insertEscapeChars(classIDToString(gtClassID));
                text(bxgt+50,bygt-20,sprintf('%s',classString),'color','w')
                rectangle('Position',[bxgt,bygt,bwgt,bhgt],'EdgeColor','w');

                % show Res box
%                 text(bxres,byres-20,sprintf('%d',idres),'color','y')
                rectangle('Position',[bxres,byres,bwres,bhres],'EdgeColor','y');

                pause(.01)
            end
        end

        % if we encounter a partially occluded box, mark to remove
        if gtRaw(g,9)<minvis
            r = resInFrame(mRes(m)); % result box
            keepBoxes(r) = false;

            if showVis
                bxgt=gtRaw(g,3); bygt=gtRaw(g,4); bwgt=gtRaw(g,5); bhgt=gtRaw(g,6); idgt=gtRaw(g,2);
                bxres=resRaw(r,3); byres=resRaw(r,4); bwres=resRaw(r,5); bhres=resRaw(r,6); idres=resRaw(r,2);

                clf
                im = imread(fullfile(dataDir,seqName,'img1',sprintf('%06d.jpg',t)));
                imshow(im); hold on
                text(50,50,sprintf('%d',t),'color','w')

                % show GT box
                text(bxgt+50,bygt-20,sprintf('vis %.1f',gtRaw(g,9)*100),'color','w')
                rectangle('Position',[bxgt,bygt,bwgt,bhgt],'EdgeColor','w');

                % show Res box
%                 text(bxres,byres-20,sprintf('%d',idres),'color','y')
                rectangle('Position',[bxres,byres,bwres,bhres],'EdgeColor','y');

                pause(.01)
                pause
            end
        end

    end
end

%%
fprintf('\nRemoving %d boxes from solution...\n',numel(find(~keepBoxes)));
resNew = resRaw;
resNew=resNew(keepBoxes,:);

%% write new file into new dir (clean)
dlmwrite(resFileClean, resNew);

% end

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function [y]=tt_rc2(d,n,elem_fun,eps,varargin)
%[Y]=TT_RC2(D,N,ARR,ELEM_FUN,EPS,[OPTIONS])
%The TT-Renormalization-Cross algorithm
%input is: the dimension of the array d, the
%size vector n, the function that computes the prescribed element
%of an array, accuracy parameter eps
%and also additional arguments specified as 'rmax',10,'nswp',10
%'verb',true and so on
% 'x0', x0
%The elem function has syntax elem_fun(ind); all other additional
%information
%has to be passed as "anonymous" arguments

%The algorithm. It needs left & right index sets to be initialized.
%It needs supecores to be precomputed
%It is a good idea to store also the supercores (there will be d-1,
%i think) & update them; also certain indices maybe removed from the
%index set?? we will not do that;
%Two possibilities: current index contains maxvol submatrix;
%If not, increase, increase & increase (maybe even beyond the rank!)
%rank one subtraction can be made very stable

%Default parameters
if ( numel(n) == 1 )
  n=n*ones(1,d);
end
nswp=10;
rmax=1000;
x0=[];
verb=true;
vec=false; %If there is a vectorization in options, i.e. if elem_fun
%supports vectorized computations of many elements at once
for i=1:2:length(varargin)-1
    switch lower(varargin{i})
        case 'nswp'
            nswp=varargin{i+1};
        case 'rmax'
            rmax=lower(varargin{i+1});
        case 'x0'
            x0=varargin{i+1};
        case 'verb'
            verb=varargin{i+1};
        case 'vec'
            vec=varargin{i+1};
        otherwise
            error('Unrecognized option: %s\n',varargin{i});
    end
end

%The initial setup is done.

%We need initial index sets. Say, the vector of all ones (he-he)
%i_left;
%r_left;
%i_right;
%r_right;
i_left=cell(d,1);
i_right=cell(d,1);
r_left=cell(d,1);
r_right=cell(d,1);
ry=ones(d+1,1); %Initial ranks
%Find fiber maximum
ind=2*ones(d,1);
ind=find_fiber_maximum(elem_fun,ind);
for i=1:d
  i_left{i}=ind(i);
  i_right{i}=ind(i);
  r_left{i}=1;
  r_right{i}=1;
end

%Before we get the multiindex sets (helps to compute s-u-p-e-r-c-o-r-e-s)

ind_left=get_multi_left(i_left,r_left,ry);
ind_right=get_multi_right(i_right,r_right,ry);

%Now we have to compute s-u-p-e-r-c-o-r-e-s (not cores!) using ind_left &
%ind_right
super_core=cell(d-1,1);
for i=1:d-1 %There are d-1 s-u-p-e-r-c-o-r-e-s
   %Compute the index set for the i-th one
   index_set=zeros(d,ry(i),n(i),n(i+1),ry(i+2));
   if ( i == 1 )
       ileft=zeros(1,0);
   else
      ileft=ind_left{i-1};
   end
   if ( i == d - 1 )
      iright = zeros(1,0);
   else
      iright=ind_right{i+2};
   end
   for s1=1:ry(i)
      for s2=1:ry(i+2)
        for i1=1:n(i)
           for i2=1:n(i+1)
              index_set(:,s1,i1,i2,s2)=[ileft(s1,:),i1,i2,iright(s2,:)];
           end
        end
      end
   end
   M=ry(i)*n(i)*n(i+1)*ry(i+2);
   index_set=reshape(index_set,[d,M]);
   if ( vec )
      super_core{i}=reshape(elem_fun(index_set),[ry(i),n(i),n(i+1),ry(i+2)]);
   else
       cur_core=zeros(M,1);
      for k=1:M
         cur_core(k)=elem_fun(index_set(:,k));
      end
      super_core{i}=reshape(cur_core,[ry(i),n(i),n(i+1),ry(i+2)]);
   end

end
%Initialization is done.
%Now to the main iteration

%Through all the s-u-p-e-r-c-o-r-e-s
swp=1;
dir='lr';
i=1;
while ( swp <= nswp )
   %The method acts as follows.
   %1) Test if the current supercore is well approximated by the
   %current cross
   %2) If yes, do nothing.
   %3) If no, add the required indices to the cross and perform
   %modifications of the s-u-p-e-r-c-o-r-e-s

   %Nikrena ne soobrajau
   %Bilo: (i1,i2,i3,i4,i5)
   %Pofiksili (i1,i2), uvelichili r2. Rashirilos' mnojestvo (i2,i3,i4,i5)
   %(nikomu ne nujno)
   %V seredine: (i1,i2) (i3,i4,i5)+ - uvelichilos' mojectvo
   %Kogda uvelichili (i3,i4,i5) viros ry(3), uvelichilos 1 superyadro
   cur_core=super_core{i};
   %we have to convert the pair i_left{i} & r_left{i} to the
   %index set for the matrix cur_core

   cur_core=reshape(cur_core,[ry(i)*n(i),n(i+1)*ry(i+2)]);
   %i1=get_full_ind(i_left{i},r_left{i});
   %i2=get_full_ind(i_right{i+1},r_right{i+1});
    %ry(i)*n(i)

    i1=r_left{i}+(i_left{i}-1)*ry(i);
    i2=i_right{i+1}+(r_right{i+1}-1)*n(i+1);
  % if ( numel(i1)==2)
  %     if ( i1(1) == 2 && i1(2) == 2 )
  %   keyboard;
  %     end
  % end
  [i1,i2]=new_cross(cur_core,i1,i2,eps);
       %Increase i_left{i} & i_right{i+1}
       [radd_left,iadd_left]=ind2sub([ry(i);n(i)],i1);
       [iadd_right,radd_right]=ind2sub([n(i+1);ry(i+2)],i2);
       i_left{i}=iadd_left;
       r_left{i}=radd_left;
       %if ( i == 8 )
       %    keyboard
       %end
       i_right{i+1}=iadd_right;
       r_right{i+1}=radd_right;
       ry(i+1)=numel(i1);
       %Zaglushka
       ind_left=get_multi_left(i_left,r_left,ry);

       ind_right=get_multi_right(i_right,r_right,ry);

      %   fprintf('Index set in question: \n');
      %   ind_left{3}
      % if ( i == 3 && swp == 2 )
      %   keyboard;
      % end
       %keyboard;
       if ( i > 1 )
         %Recompute supercore{i-1}
         super_core{i-1}=compute_supercore(i-1,elem_fun,d,n,ry,ind_left,ind_right,vec);
       end
       if ( i < d - 1)
         super_core{i+1}=compute_supercore(i+1,elem_fun,d,n,ry,ind_left,ind_right,vec);
         %Recompute supercore{i+1};
       end

   %Convert iadd1 to i_left and i_right format; compute new
   %ind_left{i+1},ind_right{i+1}
   if ( strcmp(dir,'lr') )
     if ( i == d-1 )
        dir='rl';
     else
        i=i+1;
     end
   else
     if ( i == 1 )
        dir ='lr';
        swp=swp+1;
     else
        i=i-1;
     end
   end
   fprintf('Step %d sweep %d rank=%3.1f \n',i,swp,mean(ry));
end
fprintf('Done! \n');
keyboard;

return
end

function [super_core]=compute_supercore(i,elem_fun,d,n,ry,ind_left,ind_right,vec)
%[super_core]=compute_supercore(i,elem_fun,n,ry,ind_left,ind_right)
   index_set=zeros(d,ry(i),n(i),n(i+1),ry(i+2));
   if ( i == 1 )
       ileft=zeros(1,0);
   else
      ileft=ind_left{i-1};
   end
   if ( i == d - 1 )
      iright = zeros(1,0);
   else
      iright=ind_right{i+2};
   end
   for s1=1:ry(i)
      for s2=1:ry(i+2)
        for i1=1:n(i)
           for i2=1:n(i+1)
              index_set(:,s1,i1,i2,s2)=[ileft(s1,:),i1,i2,iright(s2,:)];
           end
        end
      end
   end
   M=ry(i)*n(i)*n(i+1)*ry(i+2);
   index_set=reshape(index_set,[d,M]);
   if ( vec )
      super_core=reshape(elem_fun(index_set),[ry(i),n(i),n(i+1),ry(i+2)]);
   else
       cur_core=zeros(M,1);
      for k=1:M
         cur_core(k)=elem_fun(index_set(:,k));
      end
      super_core=reshape(cur_core,[ry(i),n(i),n(i+1),ry(i+2)]);
   end

return
end

function [i1,i2]=new_cross(mat,i1,i2,eps)
%[i1,i2]=enlarge_cross(mat,i1,i2,eps)
%Tests whether the current index set gives a reasonable approximation
%to the matrix, and enlarges the basis if necessary.
%The method acts as follows.

%A-C*A11^{-1} R

%sbm=mat(i1,i2);
%[u,s,v]=svd(sbm,'econ');
%nrm=norm(mat);s=diag(s);
%r=numel(find(
%mat_save=mat;
i1=[];
i2=[];
desirata=eps*norm(mat,'fro');
er=norm(mat,'fro');
while ( er > desirata )
   %Find maximal element in mat
   [~,ind]=max(abs(mat(:)));
   ind=tt_ind2sub(size(mat),ind);
   i=ind(1); j=ind(2);
    i1=[i1;i];
    i2=[i2;j];
    u1=mat(:,j); u2=mat(i,:);
    u1=u1/mat(i,j);
    mat=mat-u1*u2;
    er=norm(mat,'fro');
end
% i0=[i1;iadd1]; j0=[i2;iadd2];
% sbm=mat_save(i0,j0);
% ss=svd(sbm);
% fprintf('%10.10e \n',ss)
% sbm=mat_save(i1,i2);
% fprintf('AND PREVIOUS \n');
% ss=svd(sbm);
% fprintf('%10.10e \n',ss)

return
end

function [ind_left]=get_multi_left(i_left,r_left,ry)
%[ind_left]=get_multi_left(i_left,r_left,ry)
%Computes (all) left multiindex sets for a given compact representation
  d=numel(i_left);
  ind_left=cell(d,1);
  ind_left{1}=i_left{1};
  for i=2:d
      %ind_cur=zeros(
      %ind_cur is an array of size ry(i-1) x i
      ind_cur=zeros(ry(i+1),i);
      r_prev=r_left{i};
      ind_prev=ind_left{i-1};
      i_prev=i_left{i};
      for s=1:ry(i+1)
          %ind_prev(p1(s
         ind_cur(s,:)=[ind_prev(r_prev(s),:),i_prev(s)];
      end
      ind_left{i}=ind_cur;
  end
return
end

function [ind_right]=get_multi_right(i_right,r_right,ry)
%[ind_right]=get_multi_right(i_right,r_right,ry)
%Computes (all) right multiindex sets for a given compact representation
  d=numel(i_right);
  ind_right=cell(d,1);
  ind_right{d}=i_right{d};
  for i=d-1:-1:1
      %ind_cur=zeros(
      %ind_cur is an array of size ry(i) x i
      ind_cur=zeros(ry(i),d-i+1);
      r_prev=r_right{i};
      ind_prev=ind_right{i+1};
      i_prev=i_right{i};
      for s=1:ry(i)
         ind_cur(s,:)=[i_prev(s), ind_prev(r_prev(s),:)];
      end
      ind_right{i}=ind_cur;
  end

return
end

function [ind]=find_fiber_maximum(elem_fun,ind)
%[ind]=find_fiber_maximum(elem_fun,ind)
%Simple ALS method to compute (approximate) maximum in a tensor
%In fact, need some non-zero
%fact=2; %If the new one <= fact times larger than the previous, stop

return
end

author: oseledets

repo: TT-Toolbox


text:
function parseOut = parser(lexIn)
%STRCONVPARSER    LL(1) parser for mathematical expressions
%  PARSEOUT = STRCONVPARSER(LEXIN) returns a syntax tree of expressions so that
%  it can be converted to a format Chebfun is able to work with. The input,
%  LEXIN, is the output of the method STRCONVLEXER(), and is a cell array of
%  strings, containaining the tokens of strings.
%
%  This method implements a LL(1) parser, a technique from compiler theory. For
%  more details, see e.g.,
%   [1] Aho, Sethi, Ullman, Compilers: Principles, Techniques, and Tools,
%       Addison-Wesley, 1986.
%
% See also: STRINGPARSER, STRINGPARSER/STR2ANON, STRINGPARSER/LEXER.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Developers note: Usually, a parser relies on having access to pointers, which
% is not possible in MATLAB. We get around this issue using global variables.
%
% TODO: Can we now get around this? See #515.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialize all global variables
global NEXT
global COUNTER
global LEX
global STACK

% Enter the main routine
parseMain(lexIn);

% Return the stored stack
parseOut = STACK;

% Clear all global variables
NEXT = [];
COUNTER = [];
LEX = [];
STACK = [];

end

function parseMain(lexIn)
%PARSEMAIN  The main parsing routine, starts the recursive parsing.

global NEXT
global STACK
global LEX
global COUNTER

COUNTER = 1;
LEX = lexIn;
NEXT = char(LEX(COUNTER, 2));
STACK = [];

% Our expression can only start with certain labels, make sure we are
% starting with one of them.
validTypes = {'NUM', 'VAR', 'INDVAR', 'PDEVAR', 'LAMBDA',...
            'FUNC1', 'FUNC2', 'FUNC3', 'UN-', 'UN+', 'LPAR'};

if ( any(strcmp(NEXT, validTypes)) )

    % Enter the recursion and check for successful termination.
    parseExpA();

    % We've put a $ at the end of the lexer output. Hence, we have only
    % successfully parsed our string if the only remaining token left is a
    % $ sign.
    success = match('$');

    if ( ~success )
        reportError('Parse:end', ...
            'Input expression ended in unexpected manner.');
    end

else

    reportError('Parse:start', 'Input field started with unaccepted symbol.');

end

end

% The following is the allowed grammar of mathematical expressions in the
% chebgui (¬ denotes an empty symbol):
%
%   ExpA    -> ExpA   COMMA  Exp0
%   ExpA    -> Exp0
%   Exp0    -> Exp0   OP=    Exp05
%   Exp0    -> Exp05
%   Exp05   -> Exp05  OPREL  Exp1       where OPREL can be OP>, OP>=, OP<, OP<=
%   Exp05   -> Exp1
%   Exp1    -> Exp1   OP+    Exp2
%   Exp1    -> Exp1   OP-    Exp2
%   Exp1    -> Exp2
%   Exp2    -> Exp2   OP*    Exp3
%   Exp2    -> Exp2   OP/    Exp3
%   Exp2    -> Exp3
%   Exp3    -> Exp3   OP^    Exp4
%   Exp3    -> Exp4
%   Exp4    ->        UN+    Exp4
%   Exp4    ->        UN-    Exp4
%   Exp4    -> Exp5
%   Exp5    -> NUM
%   Exp5    -> INDVAR
%   Exp5    -> PDEVAR
%   Exp5    -> LAMBDA
%   Exp5    -> VAR
%   Exp5    -> VAR ExpDer
%   Exp5    -> VAR LPAR ExpList RPAR
%   Exp5    -> VAR ExpDer LPAR ExpList RPAR
%   Exp5    -> FUNC1 LPAR Exp1 RPAR
%   Exp5    -> FUNC2 LPAR Exp1 COMMA Exp1 RPAR]
%   Exp5    -> FUNC3 LPAR Exp1 COMMA Exp1 COMMA Exp1 RPAR
%   Exp5    -> LPAR Exp1 RPAR
%   ExpDer  -> DER
%   ExpDer  -> ¬
%   ExpList -> Exp1
%   ExpList -> Exp1 COMMA STR

function parseExpA()

parseExp0();
parseExpApr();

end

function parseExp0()

parseExp05();
parseExp0pr();

end

function parseExp05()

parseExp1();
parseExp05pr();

end

function parseExp1()

parseExp2();
parseExp1pr();

end

function parseExp2()

parseExp3();
parseExp2pr();

end

function parseExp3()

parseExp4();
parseExp3pr();

end

function parseExp4()

parseExp5();

end

function parseExp5()

global NEXT
global COUNTER
global LEX

% We begin by checking whether we have hit a terminal case. In that case,
% we push that into the stack. We need to treat variables with _ in the
% names separately, as we only allow certain operators around the time
% derivative.
if ( strcmp(NEXT, 'VAR') )
    % Create a new tree for later use, corresponding to the dependent
    % variable
    tempLeaf = struct('center', {{char(LEX(COUNTER)), char(NEXT)}}, ...
        'pdeflag', 0);

    % Push the variable onto the stack so that we can take the correct
    % actions if we have derivatives or expressions on the form u(...)
    % involved. If not, we'll simply pop it later.
    push(tempLeaf);

    % Begin by advancing as usual.
    advance();

    % If there follows a derivative symbol, ', we parse that
    parseExpDer();

    % If we then have a u(3) or u(3,left) situation (or u'(3) etc.), we
    % parse that as required
    parseExpList();

elseif ( any(strcmp(NEXT, {'NUM', 'INDVAR', 'LAMBDA', 'STR'})) )
    newLeaf = struct('center', {{char(LEX(COUNTER)), char(NEXT)}}, ...
        'pdeflag', 0);
    push(newLeaf);
    advance();

elseif ( strcmp(NEXT, 'PDEVAR') )
    newLeaf = struct('center', {{char(LEX(COUNTER)), char(NEXT)}}, ...
        'pdeflag', 1);
    push(newLeaf);
    advance();

elseif ( strcmp(NEXT, 'FUNC1') ) % Functions which take one argument
    parseFunction1();

elseif ( strcmp(NEXT, 'FUNC2') ) % Functions which take two arguments
    parseFunction2();

elseif ( strcmp(NEXT, 'FUNC3') ) % Functions which take three arguments
    parseFunction3();

elseif ( strcmp(NEXT, 'LPAR') )
    advance();
    parseExp05();

    % Check if NEXT symbol is ')' as it should be. If not, there is a
    % parenthesis imbalance in the expression and we return an error.
    m = match('RPAR');
    if ( ~m )
        reportError('Parse:parenths', 'Parenthesis imbalance in input fields.')
    end

elseif  ( strcmp(NEXT, 'UN-') || strcmp(NEXT, 'UN+') || ...
          strcmp(NEXT, 'OP-') || strcmp(NEXT,'OP+') )
    % If + or - reaches this far, we have an unary operator.
    % ['UN', char(NEXT(3))] determines whether we have UN+ or UN-.
    newCenterNode = {{char(LEX(COUNTER)), ['UN', char(NEXT(3))]}};
    advance();
    parseExp4();

    rightArg = pop();

    pdeflag = rightArg.pdeflag;

    newTree = struct('center', newCenterNode, 'right', rightArg, ...
        'pdeflag',pdeflag);
    push(newTree);
else

    reportError('Parse:terminal', ...
        ['Unrecognized character in input field:', NEXT]);
end

end

function parseFunction1()

global COUNTER
global LEX

functionName = char(LEX(COUNTER));
advance();

if ( match('LPAR') )
    parseExp1();

    if ( match('COMMA') )
        reportError('Parse:func1', ...
            ['Method ''', functionName, ''' only takes one input argument.']);
    elseif ( ~match('RPAR') )
        reportError('Parse:parenths', 'Parenthesis imbalance in input fields.')
    end

    rightArg =  pop();
    if ( rightArg.pdeflag )
        error('CHEBFUN:STRINGPARSER:parser:PDE', ...
            'Cannot use time derivative as function arguments.')
    end
    % Can assume no pde if we reach here
    newTree = struct('center', {{functionName, 'FUNC1'}}, ...
        'right', rightArg, 'pdeflag', 0);
    push(newTree);
else
    reportError('Parse:parenths', ...
        'Need parenthesis when using functions in input fields.')
end

end

function parseFunction2()

global NEXT
global COUNTER
global LEX

% Function which allow one or two arguments
oneArgAllowed = {'diff', 'cumsum', 'airy', 'mean'};
functionName = char(LEX(COUNTER));
advance();

% Here we need ( as the next symbol
if ( strcmp(NEXT, 'LPAR') )
    advance();
    parseExp1();

    firstArg =  pop();
    if ( firstArg.pdeflag )
        error('CHEBFUN:STRINGPARSER:parser:PDE', ...
            'Cannot use time derivative as function arguments.')
    end

    % Check whether we have a comma, if so, continue as normal
    if ( match('COMMA') )
        parseExp1();
        m = match('RPAR');
        if ( ~m )
            reportError('Parse:parenths', ...
                'Parenthesis imbalance in input fields.')
        end

        secondArg =  pop();
        if ( secondArg.pdeflag )
            error('CHEBFUN:STRINGPARSER:parser:PDE', ...
                'Cannot use time derivative as function arguments.')
        end

        % Can assume no pde if we reach here
        newTree = struct('left', firstArg, ...
            'center', {{functionName, 'FUNC2'}}, ...
            'right', secondArg,'pdeflag',0);
        push(newTree);

    elseif ( match('RPAR') )
        if ( any(strcmp(functionName, oneArgAllowed)) )
            % Have hit a function which allows one or two args.

            % If we only had one argument, we convert the function to type
            % FUNC1.  Can assume no PDE if we reach here
            newTree = struct('center', {{functionName, 'FUNC1'}}, ...
                'right', firstArg, 'pdeflag', 0);
            push(newTree);
        else
            % We tried to call a method which requires two args with only one.
            reportError('Parse:func2', ['Method ''', functionName, ...
                ''' requires two input arguments.']);
        end

    else
        reportError('Parse:parenths', 'Parenthesis imbalance in input fields.')
    end

else
    reportError('Parse:parenths', ...
        'Need parenthesis when using functions in input fields.')
end

end

function parseFunction3()

global NEXT
global COUNTER
global LEX

% Function which allow one or two arguments
oneArgAllowed = {'sum', 'integral'};
twoArgAllowed = {'feval', 'fred', 'volt'};
functionName = char(LEX(COUNTER));
advance();

% Here we need ( as the next symbol
if ( strcmp(NEXT, 'LPAR') )
    advance();
    parseExp1();
    firstArg = pop();
    if ( firstArg.pdeflag )
        error('CHEBFUN:STRINGPARSER:parser:PDE', ...
            'Cannot use time derivative as function arguments.')
    end

    % Check whether we have a comma, if so, continue as normal
    if ( match('COMMA') )
        parseExp1();
        secondArg = pop();
        if ( secondArg.pdeflag )
            error('CHEBFUN:STRINGPARSER:parser:PDE', ...
                'Cannot use time derivative as function arguments.')
        end

        if ( match('COMMA') ) % and again
            parseExp1();
            thirdArg = pop();

            % Define the new branch
            newTree = struct('left', firstArg, ...
                'center', {{functionName, 'FUNC3'}}, ...
                'right', secondArg, 'arg', thirdArg, 'pdeflag', 0);

            if ( ~match('RPAR') ) % Check the final parenthesis
                reportError('Parse:parenths', ...
                    'Parenthesis imbalance in input fields.')
            end
        elseif ( match('RPAR') && any(strcmp(functionName, twoArgAllowed)) )
            % This was actually a two argument function in disguise
            newTree = struct('left', firstArg, ...
                'center', {{functionName, 'FUNC2'}},...
                'right', secondArg, 'pdeflag', 0);
        else
            reportError('Parse:parenths', ...
                'Parenthesis imbalance in input fields.')
        end
        push(newTree);

    elseif ( match('RPAR') )
        if ( any(strcmp(functionName,oneArgAllowed)) )
            % Have hit a function which allows one or two args.

            % If we only had one argument, convert the function to type FUNC1.
            % Can assume no PDE if we reach here.
            newTree = struct('center', {{functionName, 'FUNC1'}}, ...
                'right', firstArg, 'pdeflag', 0);
            push(newTree);
        else
            % We tried to call a method which requires two args with only one.
            reportError('Parse:func2', ['Method ''', functionName, ...
                ''' requires two input arguments.']);
        end

    else
        reportError('Parse:parenths', 'Parenthesis imbalance in input fields.')
    end

else
    reportError('Parse:parenths', 'Need parenthesis when using functions in input fields.')
end

end

function parseExp1pr()

global NEXT

if ( strcmp(NEXT, 'OP+') )

    advance();
    leftArg  = pop();
    parseExp2();
    rightArg = pop();

    pdeflag = leftArg.pdeflag || rightArg.pdeflag;

    newTree = struct('left', leftArg, 'center', {{'+', 'OP+'}}, ...
        'right', rightArg, 'pdeflag', pdeflag);
    push(newTree);
    parseExp1pr();

elseif ( strcmp(NEXT, 'OP-') )
    advance();
    leftArg  = pop();
    parseExp2();
    rightArg = pop();

    pdeflag = leftArg.pdeflag || rightArg.pdeflag;
    newTree = struct('left', leftArg, 'center', {{'-', 'OP-'}}, ...
        'right', rightArg, 'pdeflag', pdeflag);
    push(newTree);
    parseExp1pr();
elseif ( strcmp(NEXT, 'RPAR') || strcmp(NEXT, '$') || strcmp(NEXT, 'OP=') )
	% Do nothing
else % If we don't have ) or the end symbol now something has gone wrong.
%     reportError('Parse:end','Syntax error in input fields.')
end

end

function parseExp2pr()

global NEXT

if ( strcmp(NEXT,'OP*') )
    leftArg  = pop();   % Pop from the stack the left argument
    advance();          % Advance in the input
    parseExp3();
    rightArg = pop();  % Pop from the stack the right argument

    % Check whether we have _ variables
    if ( leftArg.pdeflag || rightArg.pdeflag )
        error('CHEBFUN:STRINGPARSER:parser:PDE','Cannot multiply time derivative')
    end

    % Can assume no PDE if we reach here
    newTree = struct('left', leftArg, 'center', {{'.*', 'OP*'}}, ...
        'right', rightArg, 'pdeflag', 0);
    push(newTree);
    parseExp2pr();

elseif ( strcmp(NEXT,'OP/') )
    leftArg  = pop();   % Pop from the stack the left argument
    advance();          % Advance in the input
    parseExp3();
    rightArg = pop();  % Pop from the stack the right argument

    % Check whether we have _ variables
    if ( leftArg.pdeflag || rightArg.pdeflag )
        error('CHEBFUN:STRINGPARSER:parser:PDE','Cannot divide with time derivatives')
    end

    % Can assume no PDE if we reach here
    newTree = struct('left', leftArg, 'center', {{'./', 'OP/'}}, ...
        'right', rightArg, 'pdeflag', 0);
    push(newTree);
    parseExp2pr();
else
    % Do nothing
end

end

function parseExp3pr()

global NEXT

if ( strcmp(NEXT,'OP^') )
    leftArg  = pop();
    advance();
    parseExp4();
    rightArg = pop();

    % Check whether we have _ variables
    if ( leftArg.pdeflag || rightArg.pdeflag )
        error('CHEBFUN:STRINGPARSER:parser:PDE','Cannot take powers with time derivative')
    end

    % Can assume no pde if we reach here
    newTree = struct('left', leftArg, 'center', {{'.^', 'OP^'}}, ...
        'right', rightArg, 'pdeflag', 0);
    push(newTree);
    parseExp3pr();

elseif ( ~isempty(strfind(NEXT, 'DER')) )
    leftArg  = pop();

    % Check whether we have _ variables
    if ( leftArg.pdeflag )
        error('CHEBFUN:STRINGPARSER:parser:PDE','Cannot differentiate time derivative')
    end

    newTree = struct('center', {{'D', NEXT}}, 'right', leftArg, 'pdeflag', 0);
    push(newTree);
    advance();
    parseExp3pr();
else
    % Do nothing
end

end

function parseExpDer()
% parseExpDer deals with ' denoting derivatives on the dependent variables
% in the problem.

global NEXT

if ( ~isempty(strfind(NEXT, 'DER')) )
    leftArg  = pop();

    % Check whether we have _ variables
    if ( leftArg.pdeflag )
        error('CHEBFUN:STRINGPARSER:parser:PDE','Cannot differentiate time derivative')
    end

    newTree = struct('center', {{'D', NEXT}}, 'right', leftArg, 'pdeflag', 0);
    push(newTree);
    advance();
else
    % Do nothing
end

end

function parseExpList()
% parseExpList deals with potential arguments to the dependent variables in
% the problem, e.g. u(3,left).

global NEXT
global COUNTER
global LEX

if ( match('LPAR') ) % We are in a u(3) or u(3,left) situation
    % The first argument expression can only be of type EXP05 or higher
    % (no = or comma-dividers are allowed)
    parseExp1();
    % Check whether we are in in u(3) or u(3,left) situation. If we
    % have a match with ), we must have a feval with two arguments. If
    % we have a match with a comma, we must have a feval with three
    % arguments, but if there is no match, we have parenthesis
    % imbalance.
    if ( match('RPAR') )
        % Here we only had one argument in u(...), so we create a feval of type
        % FUNC2 (since the feval method has two input arguments)
        secondArg = pop();
        firstArg = pop();

        % Can assume no PDE if we reach here
        newTree = struct('left', firstArg, 'center', {{'feval', 'FUNC2'}}, ...
            'right', secondArg,'pdeflag', 0);
        push(newTree);

    elseif ( match('COMMA') )
        % Here we only had two arguments in u(...), so we create a
        % feval of type FUNC3 (since the feval method has three input
        % arguments)

        % Check whether we have any of the allowed option type as the
        % second argument. If not, throw an error.
        if ( strcmp(NEXT,'STR') )
            % We got a match! But we also need to check for parenthesis
            % balance. Store the option as a temporary leaf
            optLeaf = struct('center', {{char(LEX(COUNTER)), char(NEXT)}}, ...
                'pdeflag', 0);
            % Advance and check for parenth. balance
            advance();
            if ( match('RPAR') )
                % Create a feval of type FUNC3 (since the feval method
                % has three input arguments)
                secondArg = pop();
                firstArg = pop();
                thirdArg = optLeaf;

                % Can assume no pde if we reach here
                newTree = struct('left', firstArg, ...
                    'center', {{'feval', 'FUNC3'}}, 'right', secondArg, ...
                    'arg',thirdArg, 'pdeflag', 0);
                push(newTree);
            else
                reportError('Parse:parenths', ...
                    'Parenthesis imbalance in input fields.')
            end
        else
            reportError('Parse:secondArg', ...
                'Invalid second argument to u(0,...) type of expression.')
        end

    else % There must be a parenthesis imbalance.
        reportError('Parse:parenths', 'Parenthesis imbalance in input fields.')
    end
else % We reach here for the default behaviour, e.g. u' = 1, no parentheses involved
    % Do nothing, since we've already pushed the variable onto the
    % stack.
end

end

function parseExp0pr()

global NEXT

if ( strcmp(NEXT, 'OP=') )

    leftArg  = pop();
    advance();
    parseExp05();

    rightArg = pop();

    pdeflag = leftArg.pdeflag || rightArg.pdeflag;

    newTree = struct('left', leftArg, 'center', {{'=', 'OP='}}, ...
        'right', rightArg, 'pdeflag', pdeflag);
    push(newTree);
else
	% Do nothing
end

end

function parseExpApr()

global NEXT

if ( strcmp(NEXT, 'COMMA') )

    advance();
    leftArg  = pop();
    parseExpA();
    rightArg = pop();

    pdeflag = leftArg.pdeflag || rightArg.pdeflag;

    newTree = struct('left', leftArg, 'center', {{',', 'COMMA'}}, ...
        'right', rightArg, 'pdeflag', pdeflag);
    push(newTree);
else
	% Do nothing
end

end

function parseExp05pr()

global NEXT

if ( any(strcmp(NEXT, {'OP>', 'OP>=', 'OP<', 'OP<='})) )
    tempOpType = NEXT;
    tempOpLabel = tempOpType(3:end);

    advance();
    leftArg  = pop();
    parseExp1();
    rightArg = pop();

    pdeflag = leftArg.pdeflag || rightArg.pdeflag;

    newTree = struct('left', leftArg, 'center', {{tempOpLabel, tempOpType}}, ...
        'right', rightArg, 'pdeflag', pdeflag);
    push(newTree);
else
	% Do nothing
end

end

function advance()
%ADVANCE    Move to the next token in the output from the lexer.
global NEXT
global COUNTER
global LEX

COUNTER = COUNTER + 1;
NEXT = char(LEX(COUNTER,2));

end

function m = match(label)

global NEXT
m = strcmp(label, NEXT);

% If we found a match and are not at the end of output, we want to advance
% to the NEXT symbol
if ( m && ~strcmp(label, '$') )
    advance();
end

end

function push(new)
%PUSH   Push a tree to the stack of syntax trees
global STACK

if ( ~stackRows() )
    STACK = new;
else
    % Ensure number of fields matches
    [STACK, new] = mergeFields(STACK, new);
    % Update the Stack
    STACK = [STACK ; new];
end

end

function p = pop()
%POP    Pop a tree from the stack.
global STACK

% Throw a sensible error if we have an empty stack
if ( isempty(STACK) )
    reportError('Parse:end', 'Syntax error in input expression (Empty stack)');
end

p = STACK(end,:);
STACK(end,:) = [];

end

function m = stackRows()

global STACK
[m, n] = size(STACK); %#ok<NASGU>

end

% TODO:  Can we get rid of this function entirely?
function reportError(id, msg)

error(['CHEBFUN:', id], msg);
% ME = MException(id,msg);
% throw(ME);

end

function [a, b] = mergeFields(a, b)

aFields = fieldnames(a);
bFields = fieldnames(b);

for k = 1:numel(aFields)
    if ( ~isfield(b, aFields{k}) )
        b.(aFields{k}) = [];
    end
end

for k = 1:numel(bFields)
    if ( ~isfield(a, bFields{k}) )
        a.(bFields{k}) = [];
    end
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function C=makeSymbolicStiffnessTensor(opt)

C=sym(zeros([3 3 3 3]));
I=eye(3,3);
II1=dyadicProduct(I,I,1);
II3=dyadicProduct(I,I,3);

for i=1:3;
    for j=1:3;
        for k=1:3;
            for l=1:3;
                cvar=['c',num2str(i),num2str(j),num2str(k),num2str(l)];

                C(i,j,k,l)=sym(cvar);

            end;
        end;
    end;
end;

switch opt
    case 'iso'
        L=(II1+II3)==0;
    case 'transiso'

    case'ortho'

    case'full'
        L=false(size(C));
    case 'empty'
        L=true(size(C));
end
C(L)=0;

end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% dallrgl - DAL with logistic loss and grouped L1 regularization
%
% Overview:
%  Solves the optimization problem:
%   [xx,bias] = argmin sum(log(1+exp(-yy.*(A*x+bias)))) + lambda*||x||_G1
%  where
%   ||x||_G1 = sum(sqrt(sum(xx(Ii).^2)))
%   (Ii is the index-set of the i-th group
%
% Syntax:
%  [xx,bias,status]=dallrgl(xx0, bias0, A, yy, lambda, <opt>)
%
% Inputs:
%  xx0    : initial solution ([nn,1] with opt.blks or [ns nc] with
%           ns*nc=nn for nc groups of size ns)
%  bias0  : initial bias (set [] if bias term is unnecessary)
%  A      : the design matrix A ([mm,nn]) or a cell array {fA, fAT, mm, nn}
%           where fA and fAT are function handles to the functions that
%           return A*x and A'*x, respectively, and mm and nn are the
%           numbers of rows and columns of A.
%  yy     : the target label vector (-1 or +1) ([mm,1])
%  lambda : the regularization constant
%  <opt>  : list of 'fieldname1', value1, 'filedname2', value2, ...
%   blks     : vector that contains the size of the groups.
%              sum(opt.blks)=nn. If omitted, opt.blks = [ns,..., ns]
%              and length(opt.blks)=nc, where nc is the number of groups.
%   stopcond : stopping condition, which can be
%              'pdg'  : Use relative primal dual gap (default)
%              'fval' : Use the objective function value
%           (see dal.m for other options)
% Outputs:
%  xx     : the final solution ([nn,1])
%  bias   : the final bias term (scalar)
%  status : various status values
%
% Example:
% m = 1024; n = [64 64]; k = round(0.1*n(1)); A=randn(m,prod(n));
% w0=randsparse(n,k); yy=sign(A*w0(:)+0.01*randn(m,1));
% lambda=0.1*max(sqrt(sum(reshape(A'*yy/2,n).^2)));
% [ww,bias,stat]=dallrgl(zeros(n), 0, A, yy, lambda);
%
% Copyright(c) 2009 Ryota Tomioka
% This software is distributed under the MIT license. See license.txt

function [ww,bias,status]=dallrgl(ww, bias, A, yy, lambda, varargin)

opt=propertylist2struct(varargin{:});
opt=set_defaults(opt,'solver','cg',...
                     'stopcond','pdg',...
                     'blks',[]);

if ~isequal(unique_bc(yy), [-1;1])
  error('yy must be a column vector of -1''s and 1''s');
end

if isempty(opt.blks)
  opt.blks=size(ww,1)*ones(1,size(ww,2));
  ww = ww(:);
end

prob.floss    = struct('p',@loss_lrp,'d',@loss_lrd,'args',{{yy}});
prob.fspec    = @gl_spec;
prob.dnorm    = @gl_dnorm;
prob.obj      = @objdalgl;
prob.softth   = @gl_softth;
prob.stopcond = opt.stopcond;
prob.ll       = min(0,yy);
prob.uu       = max(0,yy);
prob.Ac       =[];
prob.bc       =[];
prob.info     = struct('blks',opt.blks);

if isequal(opt.solver,'cg')
  prob.hessMult = @hessMultdalgl;
end

if isequal(opt.stopcond,'fval')
  opt.feval = 1;
end

if isnumeric(A)
  A = A(:,:);
  [mm,nn]=size(A);
  At=A';
  fA = struct('times',@(x)A*x,...
              'Ttimes',@(x)At*x,...
              'slice', @(I)A(:,I));
  clear At;
elseif iscell(A)
  mm = A{3};
  nn = A{4};
  fAslice = @(I)fA(sparse(I,1:length(I),ones(length(I),1), nn, length(I)));
  fA = struct('times',A{1},...
              'Ttimes',A{2},...
              'slice',fAslice);
else
  error('A must be either numeric or cell {@(x)A*x, @(y)(A''*y), mm, nn}');
end

prob.mm       = mm;
prob.nn       = nn;

if isempty(bias)
  B = [];
else
  B = ones(mm,1);
end

[ww,bias,status]=dal(prob,ww,bias,fA,B,lambda,opt);

if all(opt.blks==opt.blks(1))
  ns=opt.blks(1);
  nc=length(ww)/ns;
  ww=reshape(ww, [ns,nc]);
end

repo: REST


text:
function [ a, b ] = p19_lim ( )

%*****************************************************************************80
%
%% P19_LIM returns the integration limits for problem 19.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the limits of integration.
%
  a = 0.0;
  b = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function p = hex_grid_points ( nodes_per_layer, layers, n, box )

%*****************************************************************************80
%
%% HEX_GRID_POINTS returns coordinate box hex grid points.
%
%  Discussion:
%
%    This routine determines the coordinates of the elements of
%    a hexagonal grid in the unit square.
%
%    A hexagonal grid is defined in the coordinate box [A,B] x [C,D].
%
%    All nodes of the grid lie on one of LAYERS horizontal lines.
%    The first of these lines is from (A,C) to (B,C), and each
%    successive line is HY units higher.
%
%    On all the odd numbered lines, there are NODES_PER_LAYER points,
%    equally spaced from A to B, with a spacing of HX.
%
%    On the even numbered lines, there are NODES_PER_LAYER-1 points,
%    whose values are the midpoints of successive intervals on
%    an odd numbered line.  (The grid is staggered).
%
%    HY = HX * sqrt ( 3 ) / 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODES_PER_LAYER, the number of grid points on the first
%    horizontal layer of points.
%
%    Input, integer LAYERS, the number of horizontal layers.
%
%    Input, integer N, the total number of hex grid points.
%
%    Input, real BOX(2,2), the values of A, B, C and D
%    that define the coordinate box.
%
%    Output, real P(2,N), the coordinates of the
%    mesh points, listed one horizontal layer at a time.
%
  ndim = 2;

  if ( nodes_per_layer < 1 )
    p = [];
    return
  end

  if ( nodes_per_layer == 1 )
    p(1:ndim,1) = ( box(1:ndim,1) + box(1:ndim,2) ) / 2.0;
    return
  end

  [ hx, hy ] = hex_grid_h ( nodes_per_layer, box );

  k = 0;

  for j = 1 : layers

    y = box(2,1) + hy * ( j - 1 );

    jmod = mod ( j, 2 );

    if ( jmod == 1 )

      for i = 1 : nodes_per_layer
        x = box(1,1) + ( box(1,2) - box(1,1) ) * ( i - 1 ) ...
          / ( nodes_per_layer - 1 );
        k = k + 1;
        if ( k <= n )
          p(1,k) = x;
          p(2,k) = y;
        end
      end

    else

      for i = 1 : nodes_per_layer-1
        x = box(1,1) + ( box(1,2) - box(1,1) ) * ( 2 * i - 1 ) ...
          / ( 2 * nodes_per_layer - 2 );
        k = k + 1;
        if ( k <= n )
          p(1,k) = x;
          p(2,k) = y;
        end
      end

    end

  end

  return
end


text:
function img = flowToColor(flow, varargin)

%  flowToColor(flow, maxFlow) flowToColor color codes flow field, normalize
%  based on specified value,
%
%  flowToColor(flow) flowToColor color codes flow field, normalize
%  based on maximum flow present otherwise

%   According to the c++ source code of Daniel Scharstein
%   Contact: schar@middlebury.edu

%   Author: Deqing Sun, Department of Computer Science, Brown University
%   Contact: dqsun@cs.brown.edu
%   $Date: 2007-10-31 18:33:30 (Wed, 31 Oct 2006) $

% Copyright 2007, Deqing Sun.
%
%                         All Rights Reserved
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for any purpose other than its incorporation into a
% commercial product is hereby granted without fee, provided that the
% above copyright notice appear in all copies and that both that
% copyright notice and this permission notice appear in supporting
% documentation, and that the name of the author and Brown University not be used in
% advertising or publicity pertaining to distribution of the software
% without specific, written prior permission.
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
% INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR ANY
% PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR BROWN UNIVERSITY BE LIABLE FOR
% ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
% WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
% ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
% OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

UNKNOWN_FLOW_THRESH = 1e9;
UNKNOWN_FLOW = 1e10;            %

[height widht nBands] = size(flow);

if nBands ~= 2
    error('flowToColor: image must have two bands');
end;

u = flow(:,:,1);
v = flow(:,:,2);

maxu = -999;
maxv = -999;

minu = 999;
minv = 999;
maxrad = -1;

% fix unknown flow
idxUnknown = (abs(u)> UNKNOWN_FLOW_THRESH) | (abs(v)> UNKNOWN_FLOW_THRESH) ;
u(idxUnknown) = 0;
v(idxUnknown) = 0;

maxu = max(maxu, max(u(:)));
minu = min(minu, min(u(:)));

maxv = max(maxv, max(v(:)));
minv = min(minv, min(v(:)));

rad = sqrt(u.^2+v.^2);
maxrad = max(maxrad, max(rad(:)));

fprintf('max flow: %.4f flow range: u = %.3f .. %.3f; v = %.3f .. %.3f\n', maxrad, minu, maxu, minv, maxv);

if isempty(varargin) ==0
    maxFlow = varargin{1};
    if maxFlow > 0
        maxrad = maxFlow;
    end;
end;

u = u/(maxrad+eps);
v = v/(maxrad+eps);

% compute color

img = computeColor(u, v);

% unknown flow
IDX = repmat(idxUnknown, [1 1 3]);
img(IDX) = 0;

author: shuochsu

repo: DeepVideoDeblurring


text:
% Finite phase difference and Kay frequency estimation
%
% Estimates the instantaneous frequency of the input signal using
% General Phase Difference (FFD, CFD, 4th, 6th) and Kay Frequency
% Estimation
%
% Usage:
%
%     ife = pde( signal, order [, window_length]);
%
% Inputs
%
%    signal
%
%	  Input one dimensional signal to be analysed.
%
%    order
%
%	  Order of the finite phase difference estimator. Available
%	  estimator orders are: 1, 2, 4, 6.
%
%    window_length
%
%	  Kay smoothing window length in the case of weighted phase
%	  difference estimator.  NB Set order = 1 for Kay's method.
%
%
%     For General Phase	Difference:
%
%     window_length:0;
%     order:1, 2, 4 or 6;
%
%     e.g. ife = pdeIF( signal, order);
%
%     for Kay smoothing weighted phase difference:
%
%     window_length: desired window length;
%     order=1;
%
%     e.g. ife = pde( signal, 0, window_length );
%
%
% TFSAP 7.0
% Copyright Prof. B. Boashash
% Qatar University, Doha
% email: tfsap.research@gmail.com

author: HeLiangHIT


text:
%demo_SOFAcalculateITD - Load HRTF and plots ITD.

% #Author: Piotr Majdak
% #Author: Michael Mihocic: bugs fixed (10.2021)
% #Author: Michael Mihocic: header documentation updated (28.10.2021)
%
% SOFA Toolbox - demo script
% Copyright (C) Acoustics Research Institute - Austrian Academy of Sciences
% Licensed under the EUPL, Version 1.2 or – as soon they will be approved by the European Commission - subsequent versions of the EUPL (the "License")
% You may not use this work except in compliance with the License.
% You may obtain a copy of the License at: https://joinup.ec.europa.eu/software/page/eupl
% Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing  permissions and limitations under the License.
%

%% Define parameters
% Subject index of the file to convert
subject=3;

%% load SOFA file
SOFAfn=fullfile(SOFAdbPath, 'database', 'cipic', ['subject_' sprintf('%03d',subject) '.sofa']);
Obj=SOFAload(SOFAfn, 'nochecks');

%% Calculate Interaural time delay
[itd_time, ~, ~, Obj_time] = SOFAcalculateITD(Obj, 'time', 'thr', 20);
[itd_samples, ~, ~, Obj_samples] = SOFAcalculateITD(Obj, 'samples', 'thr', 20);

%% Plot results
h = figure('Name',mfilename);
subplot(211)
plot((Obj_time.Data.Delay(:,1) - Obj_time.Data.Delay(:,2))*1e6)
xlabel('Position')
ylabel('Time (\mus)')
title('ITD (time)')
axis tight

subplot(212)
plot((Obj_samples.Data.Delay(:,1) - Obj_samples.Data.Delay(:,2)))
xlabel('Position')
ylabel(['Samples (Fs:' num2str(Obj.Data.SamplingRate), 'Hz)'])
title('ITD (samples)')
axis tight

%% Polar plot (not working in Octave)
if ~exist('OCTAVE_VERSION','builtin')
    figure('Name',mfilename);
    SOFAplotHRTF(Obj, 'itdhorizontal');
    title('ITD (time, horizontal plane)')
end

author: sofacoustics

repo: SOFAtoolbox


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

clear; clc;
%% Parameters
S = 100;                                   % Spot prices
T = 3;      % maturities

r = 0;                      % discount factors
d = 0;                              % dividends

titS = 'MC Bates - Asset';
titV =' MC Bates - Variance';

legend_base = 'Base scenario';

%% Specify the base model
vInst = 0.04;                  % instantanuous variance of base parameter set
vLong = 0.04;                  % long term variance of base parameter set
kappa = 0.2;                   % mean reversion speed of variance of base parameter set
omega = 0.1;                   % volatility of variance of base parameter set
rho = 0;                       % correlation of base parameter set
sigj = 0.25;
muj = 0.2;
lambda = 0.5;

%% Simulation parameters
NTime = 750; NSim = 1; NBatches = 1;
K = ones(NTime+1,1);
K = T*cumsum(K)/(NTime+1);
%rand('seed', 1);
UV = rand(1,2*NTime);              % precompute all randoms
K1 = K(2:end);

rstream = RandStream('mt19937ar','Seed',12345);
rstreamstate = rstream.State;

[PathS, PathV] = MC_QE_j(S,r,d,T,vInst,vLong,...
    kappa,omega,rho,muj,sigj,lambda,NTime,NSim,NBatches);
%MC_Bates_path(S,r,d,T,vInst,vLong,kappa,omega,rho,muj,sigj,lambda,UV);
%% Changing vInst
muj_low = 0.05;
muj_high = 0.5;

    PathS_low  = MC_QE_j(S,r,d,T,vInst,vLong,kappa,omega,rho,muj_low,sigj,lambda,NTime,NSim,NBatches);
    PathS_high = MC_QE_j(S,r,d,T,vInst,vLong,kappa,omega,rho,muj_high,sigj,lambda,NTime,NSim,NBatches);

legend_low = 'Changing \mu_j low';
legend_high = 'Changing \mu_j high';

createfigure_path(K,PathS_low, PathS, PathS_high, titS, legend_low, legend_base, legend_high);
%% Changing vLong
sigj_low = 0.1;
sigj_high = 0.5;

    PathS_low = MC_QE_j(S,r,d,T,vInst,vLong,kappa,omega,rho,muj,sigj_low,lambda,NTime,NSim,NBatches);
    PathS_high = MC_QE_j(S,r,d,T,vInst,vLong,kappa,omega,rho,muj,sigj_high,lambda,NTime,NSim,NBatches);

    legend_low  = 'Changing \sigma_j low';
    legend_high = 'Changing \sigma_j high';

createfigure_path(K,PathS_low, PathS, PathS_high, titS, legend_low, legend_base, legend_high);

%% Changing kappa
lambda_low = 0.2;
lambda_high = 0.8;

    PathS_low = MC_QE_j(S,r,d,T,vInst,vLong,kappa,omega,rho,muj,sigj,lambda_low,NTime,NSim,NBatches);
    PathS_high = MC_QE_j(S,r,d,T,vInst,vLong,kappa,omega,rho,muj,sigj,lambda_high,NTime,NSim,NBatches);

legend_low  = 'Changing \lambda low';
legend_high = 'Changing \lambda high';

createfigure_path(K,PathS_low, PathS, PathS_high, titS, legend_low, legend_base, legend_high);

author: Sable

repo: mcbench-benchmarks


text:
cd 'nrsfm';
mex computeH.c
cd ..

cd external/SIRFS
compile

cd ../..

repo: CategoryShapes


text:
function vw = initHiddenGray(dataType,scan,roi);
%
% function hiddenGray = initHiddenGray([dataType],[scan],[roi])
%
% djh, sometime in 1999
% ras, 04/05, set name to 'hidden'
% ras, 05/05, added 'curScan' field
% ras, 09/05, ensures session is loaded,
% mrGlobals is run; can specify data type,
% scan, and roi at outset.

%TODO: Decide if we want to add a 'loadAnat' to the initHiddenGray

if ~exist('dataType','var') || isempty(dataType)
    dataType = 1;
end

if ~exist('scan','var') || isempty(scan)
    scan = 1;
end

evalin('base','mrGlobals');
evalin('base','HOMEDIR = pwd;');
evalin('base','loadSession');

verbose = prefsVerboseCheck;
if verbose, disp('Initializing HIDDEN Gray view'); end

vw.name     ='hidden';
vw.viewType ='Gray';
vw.subdir   ='Gray';

% Initialize slots for anat co, amp, and ph
vw.anat         = [];
vw.co           = [];
vw.amp          = [];
vw.ph           = [];
vw.map          = [];
vw.mapName      = '';
vw.mapUnits     = '';
vw.tSeriesSlice = [];
vw.loc          = viewGet(vw,'Size') ./ 2;

% Initialize ROIs
vw.ROIs         = [];
vw.selectedROI  = 0;

% Initialize curDataType
if isnumeric(dataType),
    vw.curDataType = dataType;
elseif ischar(dataType),
    vw.curDataType = existDataType(dataType);
end

% Initialize curScan
vw.curScan = scan;

% need to record mmPerVox
vw.mmPerVox = readVolAnatHeader(getVAnatomyPath);

% Compute/load coords
vw = switch2Gray(vw);

% load any ROIs specified
if exist('roi','var') & ~isempty(roi)
    vw = loadROI(vw,roi);
end

% Dummy UI properties, allows loading of contrast maps on meshes w/o
% mrVista GUI open
vw = viewSet(vw, 'displaymode', 'anat');
vw = viewSet(vw, 'initdisplaymodes');

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [Peg, Peh] = getPe(K, Ec, No, N, pulse);

[mew, sigma] = get_mean_var_G(K, N, Ec, pulse.m1);
Ig = sqrt(N*No/2 + mew);
Th = N*sqrt(Ec);       %Non- Coherent
Peg = qfunc(Th/Ig);
[mew, sigma] = get_mean_var_H(K, N, Ec, pulse.m1, pulse.m2, pulse.w1, pulse.w2);
I1 = sqrt(No*N/2 + mew);
I2 = sqrt(No*N/2 + (mew + sqrt(3) * sigma));
I3 = sqrt(No*N/2 + (mew - sqrt(3) * sigma));
Peh = 2/3*qfunc(Th/I1) + 1/6*qfunc(Th/I2) + 1/6*qfunc(Th/I3);

function [mew_MAI, sigma_MAI] = get_mean_var_H(K, N, P, m1, m2, w1, w2)
T =1;
mew_MAI = (N*P * (K-1)*m1);
sigma_MAI = (N*P) * ((K-1)*( (.375*w1 - m1^2) + ((N-1)/(N^2))*((1.5*w2) + (K-2)*(m2^2) )))^(.5);

function [mew_MAI, sigma_MAI] = get_mean_var_G(K, N, P, m1)
T =1;
mew_MAI = (N*P * (K-1)*m1);
sigma_MAI = 0;


text:
function Rob = motion(Rob, Tim)
%MOTION Robot motion.
%   ROB = MOTION(ROB, TIM) performs one EKF-prediction motion step to robot
%   Rob in the global map Map, following the motion model in Rob.motion.
%   Both Rob and Map are updated. The time information Tim is used only if
%   the motion model requires it, but it has to be provided because MOTION
%   is a generic method.
%
%   The following motion models are supported:
%       'odometry'   uses function odo3()
%       'constVel'   uses function constVel()
%   Edit this file to add new motion models.
%
%   See also SIMMOTION, CONSTVEL, ODO3, UPDATEFRAME.

%   Copyright 2009 David Marquez @ LAAS-CNRS.

global Map

% Update rob and sen info from map
Rob = map2rob(Rob);

% robot state range
r = Rob.state.r;

switch Rob.motion

    case  {'constVel'} % constant velocity

        % motion model of the robot: mean and Jacobians
        [Map.x(r), F_x, F_u] = constVel(Map.x(r),Rob.con.u,Tim.dt);

        % update Rob and Map structures - mean only
        Rob = map2rob(Rob);

        % Covariances matrix update
        predictBlockEkf(r, F_x, Rob.con.U, F_u);

    case  {'odometry'}  % 3D odometry

        % motion model of the robot: mean and Jacobians
        [Rob.frame, F_x, F_u]   = odo3(Rob.frame,Rob.con.u);

        % update Rob and Map structures - mean only
        Map.x(Rob.frame.r) = Rob.frame.x;

        % Covariances matrix update
        predictBlockEkf(r, F_x, Rob.con.U, F_u);

    otherwise

        error('??? Unknown motion model ''%s''.',Rob.motion);
end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb