text:
function triangulation_mask ( prefix, triangle_mask )

%*****************************************************************************80
%
%% MAIN is the main program for TRIANGULATION_MASK.
%
%  Discussion:
%
%    TRIANGULATION_MASK allows a user to remove triangles from a triangulation.
%
%    The user supplies a node file and a triangle file, containing
%    the coordinates of the nodes, and the indices of the nodes that
%    make up each triangle.  Either 3-node or 6-node triangles may
%    be used.
%
%    The program reads the data.  Then, for each triangle, it calls
%    a user routine
%
%      function mask = triangle_mask ( dim_num, triangle_order, nodes, coord )
%
%    The input to this routine is
%
%      Input, integer DIM_NUM, the spatial dimension (should be 2),
%
%      Input, integer TRIANGLE_ORDER, the triangle order (3 or 6),
%
%      Input, integer NODES(TRIANGLE_ORDER), the indices of the nodes,
%
%      Input, real COORD(DIM_NUM,TRIANGLE_ORDER), the
%      coordinates of the nodes.
%
%      Output, logical MASK, is set by the user to be
%      TRUE if this triangle should be "masked", that is, omitted;
%      FALSE if this triangle should be retained.
%
%    Once the program has determined which triangles are to be retained,
%    it then examines the node data.  Any node which only occurs in
%    deleted triangles is also marked for deletion.
%
%    An output copy of the remaining nodes is made.  An output copy
%    of the remaining triangles is made, with the node indices renumbered
%    to access the updated node file.
%
%  Usage:
%
%    triangulation_mask ( 'prefix', @triangle_mask )
%
%    where 'prefix' is the common filename prefix:
%
%    * prefix_nodes.txt contains the node coordinates,
%    * prefix_elements.txt contains the element definitions.
%    * prefix_mask_nodes.txt will contain the masked node coordinates;
%    * prefix_mask_elements.txt will contain the masked element definitions.
%
%    and
%
%    * @triangle_mask is a pointer to the user triangle mask function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 October 2009
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGULATION_MASK\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Read a node file of NODE_NUM points in 2 dimensions.\n' );
  fprintf ( 1, '  Read an associated triangulation file of TRIANGLE_NUM\n' );
  fprintf ( 1, '  triangles using 3 or 6 nodes.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  For each triangle, call a user masking routine\n' );
  fprintf ( 1, '  to see if the triangle should be MASKED (not shown).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Write new node and triangulation files corresponding\n' );
  fprintf ( 1, '  to the unmasked data only.\n' );
%
%  The command line argument is the common filename prefix.
%
  if ( nargin < 1 )

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

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

  end
%
%  Create the filenames.
%
  node_filename = strcat ( prefix, '_nodes.txt' );
  element_filename = strcat ( prefix, '_elements.txt' );
  node_mask_filename = strcat ( prefix, '_mask_nodes.txt' );
  element_mask_filename = strcat ( prefix, '_mask_elements.txt' );
%
%  Read the data.
%
  [ dim_num, input_node_num ] = r8mat_header_read (  node_filename );

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

  input_node_xy = r8mat_data_read ( node_filename, ...
    dim_num, input_node_num );

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

  r8mat_transpose_print_some ( dim_num, input_node_num, input_node_xy, ...
    1, 1, dim_num, 5, '  First 5 nodes:' );
%
%  Read the element file.
%
  [ triangle_order, input_triangle_num ] = ...
    i4mat_header_read ( element_filename );

  if ( triangle_order ~= 3 && triangle_order ~= 6)
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_MASK - Fatal error!\n' );
    fprintf ( 1, '  Data is not for a 3-node or 6-node triangulation.\n' );
    error ( 'TRIANGULATION_MASK - 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', input_triangle_num );

  input_triangle_node = i4mat_data_read ( ...
    element_filename, triangle_order, input_triangle_num );

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

  i4mat_transpose_print_some ( triangle_order, input_triangle_num, ...
    input_triangle_node, 1, 1, triangle_order, 5, '  First 5 triangles:' );
%
%  Detect and correct 0-based indexing.
%
  input_triangle_node = mesh_base_one ( input_node_num, triangle_order, ...
    input_triangle_num, input_triangle_node );
%
%  Mask the triangles.
%
  output_triangle_num = 0;

  coord = zeros ( dim_num, triangle_order );
  input_triangle_mask = zeros ( input_triangle_num, 1 );
  for triangle = 1 : input_triangle_num

    nodes(1:triangle_order) = input_triangle_node(1:triangle_order,triangle);
    for i = 1 : dim_num
      coord(i,1:triangle_order) = input_node_xy(i,nodes(1:triangle_order));
    end

    mask = feval ( triangle_mask, dim_num, triangle_order, nodes, coord );

    input_triangle_mask(triangle) = mask;

    if ( ~mask )
      output_triangle_num = output_triangle_num + 1;
    end

  end

  if ( 0 )
    lvec_print ( input_triangle_num, input_triangle_mask, ...
      '  INPUT_TRIANGLE_MASK' );
  end
%
%  Determine which nodes are being used.
%
  input_node_mask(1:input_node_num) = 1;

  for triangle = 1 : input_triangle_num
    if ( ~input_triangle_mask(triangle) )
      for order = 1 : triangle_order
        node = input_triangle_node(order,triangle);
        input_node_mask(node) = 0;
      end
    end
  end

  if ( 0 )
    lvec_print ( input_node_num, input_node_mask, '  INPUT_NODE_MASK' );
  end

  input_to_output_node(1:input_node_num) = -1;

  output_node_num = 0;
  for node = 1 : input_node_num
    if ( ~input_node_mask(node) )
      output_node_num = output_node_num + 1;
      input_to_output_node(node) = output_node_num;
    end
  end
%
%  Write the unmasked triangles.
%
  output_triangle_num = 0;
  for triangle = 1 : input_triangle_num

    if ( ~input_triangle_mask(triangle) )
      output_triangle_num = output_triangle_num + 1;
      for order = 1 : triangle_order
        output_triangle_node(order,output_triangle_num) = ...
          input_to_output_node( input_triangle_node(order,triangle) );
      end
    end

  end

  i4mat_write ( element_mask_filename, triangle_order, ...
    output_triangle_num, output_triangle_node );

  fprintf ( 1, '\n');
  fprintf ( 1, '  The masked triangle data was written to "%s".\n', ...
    triangle_mask_filename );
%
%  Write the unmasked nodes.
%
  for input_node = 1 : input_node_num

    if ( ~input_node_mask(input_node) )

      output_node = input_to_output_node(input_node);

      output_node_xy(1:dim_num,output_node) = ...
        input_node_xy(1:dim_num,input_node);

    end

  end

  r8mat_write ( node_mask_filename, dim_num, ...
    output_node_num, output_node_xy );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The masked node data was written to "%s".\n', ...
    node_mask_filename );
%
%  Summary report.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of input triangles =   %d\n', input_triangle_num );
  fprintf ( 1, '  Number of output triangles =  %d\n', output_triangle_num );
  fprintf ( 1, '  Number of deleted triangles = %d\n', ...
    input_triangle_num - output_triangle_num );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of input nodes =   %d\n', input_node_num );
  fprintf ( 1, '  Number of output nodes =  %d\n', output_node_num );
  fprintf ( 1, '  Number of deleted nodes = %d\n', ...
    input_node_num - output_node_num );

  i4mat_transpose_print_some ( triangle_order, output_triangle_num, ...
    output_triangle_node, 1, 1, triangle_order, 5, ...
    '  First 5 output triangles:' );

  r8mat_transpose_print_some ( dim_num, output_node_num, ...
    output_node_xy, 1, 1, dim_num, 5, ...
    '  First 5 output nodes:' );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGULATION_MASK\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 i4mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% I4MAT_WRITE writes an I4MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer TABLE(M,N), the points.
%
%    Input, logical HEADER, is TRUE if the header is to be included.
%

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

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'I4MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %12d', round ( table(i,j) ) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function lvec_print ( n, a, title )

%*****************************************************************************80
%
%% LVEC_PRINT prints a logical vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the vector.
%
%    Input, logical A(N), the vector to be printed.
%
%    Input, string TITLE, a title to be printed first.
%    TITLE may be blank.
%
  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  fprintf ( 1, '\n' );
  for i = 1 : n
    value = ( a(i) ~= 0 );
    fprintf ( 1, '%6d  %1d\n', i, value );
  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 r8mat_write ( output_filename, m, n, table )

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

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

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

  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

author: johannesgerer

repo: jburkardt-m


text:
function p = equivalentPerimeter(grains)
% returns the equivalent perimeter of grain-polygon
%
% Description
% The equivalent perimeter of grain-polygon is defined as
%
% $$ p = 2 \pi ER $$,
%
% where $ER$ is the <grain2d.equivalentRadius.html equivalent radius> of a
% grain
%
% Input
%  grains - @grain2d
%
% Output
%  p - equivalent perimeter in measurement units
%
% See also
% grain2d/deltaarea grain2d/paris grain2d/equivalentRadius

p = 2*pi*equivalentRadius(grains);

author: mtex-toolbox

repo: mtex

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


text:
function [f,relres,iter]=isgramreal(s,g,a,M,varargin)
%ISGRAMREAL  Spectrogram inversion (real signal)
%   Usage:  f=isgramreal(s,g,a,M);
%           f=isgramreal(s,g,a,M,Ls);
%           [f,relres,iter]=isgramreal(...);
%
%   Input parameters:
%         c       : Array of coefficients.
%         g       : Window function.
%         a       : Length of time shift.
%         M       : Number of channels.
%         Ls      : length of signal.
%   Output parameters:
%         f       : Signal.
%         relres  : Vector of residuals.
%         iter    : Number of iterations done.
%
%   `isgramreal(s,g,a,M)` attempts to invert a spectrogram computed by ::
%
%     s = abs(dgtreal(f,g,a,M)).^2;
%
%   by an iterative method.
%
%   `isgramreal(s,g,a,M,Ls)` does as above but cuts or extends *f* to length *Ls*.
%
%   If the phase of the spectrogram is known, it is much better to use
%   |dgtreal|
%
%   `f,relres,iter]=isgramreal(...)` additionally returns the residuals in a
%   vector *relres* and the number of iteration steps done.
%
%   Generally, if the spectrogram has not been modified, the iterative
%   algorithm will converge slowly to the correct result. If the
%   spectrogram has been modified, the algorithm is not guaranteed to
%   converge at all.
%
%   `isgramreal` takes the following parameters at the end of the line of
%   input arguments:
%
%     'lt',lt      Specify the lattice type. See the help on
%                  |matrix2latticetype|. Only the rectangular or quinqux
%                  lattices can be specified.
%
%     'zero'       Choose a starting phase of zero. This is the default
%
%     'rand'       Choose a random starting phase.
%
%     'int'        Construct a starting phase by integration. Only works
%                  for Gaussian windows.
%
%     'griflim'    Use the Griffin-Lim iterative method, this is the
%                  default.
%
%     'bfgs'       Use the limited-memory Broyden Fletcher Goldfarb
%                  Shanno (BFGS) method.
%
%     'tol',t      Stop if relative residual error is less than the specified tolerance.
%
%     'maxit',n    Do at most n iterations.
%
%     'print'      Display the progress.
%
%     'quiet'      Don't print anything, this is the default.
%
%     'printstep',p
%                  If 'print' is specified, then print every p'th
%                  iteration. Default value is p=10.
%
%   The BFGS method makes use of the minFunc software. To use the BFGS method,
%   please install the minFunc software from:
%   `<http://www.cs.ubc.ca/~schmidtm/Software/minFunc.html>`_.
%
%   See also:  dgtreal, idgtreal
%
%   References: griffin1984sem decorsiere2011 liu1989limited

%   AUTHOR : Remi Decorsiere and Peter L. Søndergaard.
%   REFERENCE: OK

% Check input paramameters.

  if nargin<3
    error('%s: Too few input parameters.',upper(mfilename));
  end;

  if numel(g)==1
    error('g must be a vector (you probably forgot to supply the window function as input parameter.)');
  end;

  definput.keyvals.Ls=[];
  definput.keyvals.lt=[0 1];
  definput.keyvals.tol=1e-6;
  definput.keyvals.maxit=100;
  definput.keyvals.printstep=10;
  definput.flags.method={'griflim','bfgs'};
  definput.flags.print={'print','quiet'};
  definput.flags.startphase={'zero','rand','int'};

  [flags,kv,Ls]=ltfatarghelper({'Ls','tol','maxit'},definput,varargin);

  N=size(s,2);
  W=size(s,3);

  % Make a dummy call to test the input parameters
  Lsmallest=dgtlength(1,a,M,kv.lt);

  M2=floor(M/2)+1;

  if M2~=size(s,1)
      error('Mismatch between the specified number of channels and the size of the input coefficients.');
  end;

  L=N*a;

  if rem(L,Lsmallest)>0
      error('%s: Invalid size of coefficient array.',upper(mfilename));
  end;

  %% ----- step 3 : Determine the window

  [g,info]=gabwin(g,a,M,L,kv.lt,'callfun',upper(mfilename));

  if L<info.gl
      error('%s: Window is too long.',upper(mfilename));
  end;

  if ~isreal(g)
      error('%s: Window must be real-valued.',upper(mfilename));
  end;

  %% Actual computation

  sqrt_s=sqrt(s);

  if flags.do_zero
    % Start with a phase of zero.
    c=sqrt(s);
  end;

  if flags.do_rand
    c=sqrt_s.*exp(2*pi*1i*rand(size(s)));
  end;

  if flags.do_int
      if kv.lt(2)>1
          error(['%s: The integration initilization is not implemented for ' ...
                 'non-sep lattices.'],upper(mfilename));
      end;

    s2=zeros(M,N);
    s2(1:M2,:)=s;
    if rem(M,2)==0
      s2(M2+1:M,:)=flipud(s(2:end-1,:));
    else
      s2(M2+1:M,:)=flipud(s(2:end));
    end;
    c=constructphase(s2,g,a);
    c=c(1:M2,:);
  end;

  gd = gabdual(g,a,M);

  % For normalization purposes
  norm_s=norm(s,'fro');

  relres=zeros(kv.maxit,1);
  if flags.do_griflim
    for iter=1:kv.maxit
      f=comp_idgtreal(c,gd,a,M,kv.lt,0);
      c=comp_dgtreal(f,g,a,M,kv.lt,0);

      relres(iter)=norm(abs(c).^2-s,'fro')/norm_s;

      c=sqrt_s.*exp(1i*angle(c));

      if flags.do_print
        if mod(iter,kv.printstep)==0
          fprintf('ISGRAMREAL: Iteration %i, residual = %f.\n',iter,relres(iter));
        end;
      end;

      if relres(iter)<kv.tol
        relres=relres(1:iter);
        break;
      end;

    end;

  end;

  if flags.do_bfgs
    if exist('minFunc')~=2
      error(['To use the BFGS method in ISGRAMREAL, please install the minFunc ' ...
             'software from http://www.cs.ubc.ca/~schmidtm/Software/minFunc.html.']);
    end;

    % Setting up the options for minFunc
    opts = struct;
    opts.display = kv.printstep;
    opts.maxiter = kv.maxit;
    opts.usemex = 0;

    % Don't limit the number of function evaluations, just the number of
    % time-steps.
    opts.MaxFunEvals = 1e9;

    f0 = comp_idgtreal(c,gd,a,M,kv.lt,0);
    [f,fval,exitflag,output]=minFunc(@objfun,f0,opts,g,a,M,s,kv.lt);
    % First entry of output.trace.fval is the objective function
    % evaluated on the initial input. Skip it to be consistent.
    relres = sqrt(output.trace.fval(2:end))/norm_s;
    iter = output.iterations;
  end;

  % Cut or extend f to the correct length, if desired.
  if ~isempty(Ls)
    f=postpad(f,Ls);
  else
    Ls=L;
  end;

  f=comp_sigreshape_post(f,Ls,0,[0; W]);

%  Subfunction to compute the objective function for the BFGS method.
function [f,df]=objfun(x,g,a,M,s,lt);
  c=comp_dgtreal(x,g,a,M,lt,0);

  inner=abs(c).^2-s;
  f=norm(inner,'fro')^2;

  df=4*real(conj(comp_idgtreal(inner.*c,g,a,M,lt,0)));

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%SPATIALGRADIENT  Calculates the first order image derivative in both x and y using a Sobel operator
%
%     [dx, dy] = cv.spatialGradient(src)
%     [...] = cv.spatialGradient(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __src__ input image, 8-bit single-channel.
%
% ## Output
% * __dx__ output `int16` image with first-order derivative in x.
% * __dy__ output `int16` image with first-order derivative in y.
%
% ## Options
% * __KSize__ size of Sobel kernel. It must be 3 in the current implementation.
%   default 3
% * __BorderType__ Pixel extrapolation method, see cv.copyMakeBorder. Only
%   'Default', 'Reflect101', and 'Replicate' are supported. default 'Default'
%
% Equivalent to calling:
%
%     dx = cv.Sobel(src, 'DDepth','int16', 'XOrder',1, 'YOrder',0, 'KSize',3);
%     dy = cv.Sobel(src, 'DDepth','int16', 'XOrder',0, 'YOrder',1, 'KSize',3);
%
% See also: cv.Sobel, imgradientxy, imgradient
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function SaveAsNIfTI(data, target, output)

% function SaveAsNIfTI(data, nifti)
%
% Input:
%
% data: the data array to be saved to disk
%
% target: the NIfTI object specifying the target volume specification
%
% output: the filename for the output NIfTI file
%

% following the example in
% http://niftilib.sourceforge.net/mat_api_html/README.txt
%
% author: Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%

dat = file_array;
dat.fname = output;
dat.dim = target.dim;
dat.dtype = 'FLOAT64-LE';
dat.offset = ceil(348/8)*8;

N = nifti;
N.dat = dat;
N.mat = target.mat;
N.mat_intent = target.mat_intent;
N.mat0 = target.mat0;
N.mat0_intent = target.mat0_intent;

create(N);

N.dat(:,:,:,:) = data;

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function Out = metastockread(fullpath)

% METASTOCKREAD Read metastock files (symbols index files: master, emaster, xmaster; data files: .dat and .mwd)
%
%   METASTOCKREAD Select with UIGETFILE the metastock file to read.
%                 maetastockread() remembers the last selected file.
%
%   METASTOCKREAD(FULLPATH) Read metastock file FULLPATH to output struct OUT.
%                           The FULLPATH input is a string enclosed in single quotes.
%
%   OUT = metastockread(...)
%       OUT is a "m by 1" NON-scalar structure, where "m" is the # of .dat/.mwd files in the same directory
%       of the symbol index files.
%
%       The (sub)scalar structure has the following fields:
%           - datNum : # of the .dat/.mwd file with the data                           UINT16
%           - symbol : security symbol                                                 CHAR
%           - name   : security name                                                   CHAR
%           - inDate : initial date of the data in datestr format 29 ('yyyy-mm-dd')    CHAR
%           - fiDate : final   date of the data in datestr format 29 ('yyyy-mm-dd')    CHAR
%           - freq   : time freq. of the data I (intraday), D (daily), W, M, Q or Y    CHAR
%           - idFreq : intraday time frequency in minutes                              UINT8
%           - data   : time series data with variable number of fields (columns)       DOUBLE
%
%       The data field columns follow the schema:
%
%       | Date/Time* | Open | High | Low | Close | Volume | OpenInterest |
%       * The date is represented in datenum format. The precision depends on the field idFreq.
%
%       If data has only 5 columns then Volume and OpenInterest columns are absent.
%
% Examples:
%
%     % Read metastock sample data supplied with the submission (unzip the MSdatasample.zip)
%     Out1 = metastockread;                                   % Chose any symbols index file with UIGETFILE
%
%     Out2 = metastockread('MSdatasample directory\Emaster'); % Read in data associated with the Emaster
%     Out3 = metastockread('MSdatasample directory\Xmaster'); % Read in data associated with the Xmaster
%     Out  = [Out2; Out3];                                    % Concatenate results
%
%     whos Out2
%     Name        Size              Bytes  Class     Attributes
%     Out1       255x1            1911591  struct
%
%     Out1(1,1)
%     ans =
%         datNum: 1
%         symbol: 'ASR'
%           name: 'A.S. ROMA'
%         inDate: '2010-01-04'
%         fiDate: '2010-07-09'
%           freq: 'D'
%         idFreq: 0
%           data: [127x7 double]
%
% Additional features:
% - FEX metastockread page
%
% See also: FREAD, TYPECAST, BITGET, DATESTR, DATENUM

% Author: Oleg Komarov (oleg.komarov@hotmail.it)
% Tested on R14SP3 (7.1) and on R2009b. In-between compatibility is assumed.
% Platform tested: WIN. If the fcn doesn't work on other platforms please report.
% 06 sep 2010 - Created
% 09 sep 2010 - Added link to FEX page and edit example.
% 03 aug 2011 - Added fix for buggy *.dat files as pointed out by Davide Dalmasso

% Use a global variable to remember last opened file
global alreadyOpenedFile

% IF no inputs
if nargin == 0
    % Select a symbols index file with uigetfile
    [filename, pathname] = uigetfile({'Emaster;Master;Xmaster','Metastock files'},'Load metastock file', alreadyOpenedFile);
    fullpath = [pathname filename];
    alreadyOpenedFile = fullpath;
    % No file has been selected
    if filename == 0
        disp([char(13) 'The user selected CANCEL.' char(13)])
        alreadyOpenedFile = [];
        return
    end
% IF one input
elseif nargin == 1
    % Check
    if ~ischar(fullpath)
        error('metastockread:strFullpath','FULLPATH should be a string.')
    end
% ELSE error
else
    error('metastockread:ninput','Too many input arguments.')
end

% Open fid
fid = fopen(fullpath,'r','b');
if fid == -1
    error('metastockread:invalidFullpath','Invalid FULLPATH. File not found.')
end

% Read symbols index file: MASTER, EMASTER or XMASTER
if     regexpi(fullpath,'\w*\\master')
    Out = readMaster(fid);
elseif regexpi(fullpath,'\w*\\emaster')
    Out = readEmaster(fid);
elseif regexpi(fullpath,'\w*\\xmaster')
    Out = readXmaster(fid);
else
    error('metastockread:namefile','The name of the file should be either MASTER, EMASTER or XMASTER.')
end

% Close fid
fclose(fid);

% Retrieve data files that reside on the same path as the symbols index file selected
pathname = fullpath(1:find(fullpath == '\',1,'last'));
pathname = strrep(pathname, '\','\\');

% Create list of filenames
if Out(1).datNum < 256; fmt = 'F%d.DAT\n';
else                    fmt = 'F%d.MWD\n'; end
names = dataread('string',sprintf([pathname fmt],cat(1,Out.datNum)),'%s','delimiter','\n');

% Files to process
numOut = numel(Out);
for n = 1:numOut

    % Open file
    fid = fopen(names{n},'r','b');

    % If file not found leave NaN and skip to next security
    if fid == -1
        Out(n).data = NaN;
        continue
    end
    % Retrieve the number of records
    fseek(fid,2,'bof');
    numRecs = double(typecast(fread(fid,2,'*uint8'),'uint16'));
    % Retrieve the number of fields
    fseek(fid,0,'eof');
    nfields = ftell(fid)/(numRecs*4);
    % Fix nfields in case .dat file is buggy: second last timestamp-entry
    % is repeated as last
    nfields = fix(nfields);

    % Skip header
    fseek(fid,nfields*4,'bof');

    % Retrieve all data at once; the FULL info set would be organized as follows:
    % | Date / Time | Open | High | Low | Close | Volume | OpenInterest |
    Out(n).data = reshape(mbf2ieee(fread(fid,(numRecs-1)*nfields*4,'*uint8')), nfields, numRecs-1).';

    % Format the date and time into one datenum column (the first)
    [d,m,y] = partNumDate(Out(n).data(:,1));
    if nfields == 8
        h = fix(Out(n).data(:,2)/1e4);
        m = fix(mod(Out(n).data(:,2)/1e2,1e2));
        Out(n).data = Out(n).data(:,[1 3:8]);
        Out(n).data(:,1) = datenummx(y,m,d,h,m,0);
    else
        Out(n).data(:,1) = datenummx(y,m,d);
    end

    % Close fid
    fclose(fid);
end
end

% ------------------------------------------------------------------------------------------------------------
% READMASTER  - Read master file
% ------------------------------------------------------------------------------------------------------------
function Out = readMaster(fid)
% Skip header (each record is 53 bytes long)
fseek(fid,53,'bof');
% Read in the whole file
tmp = fread(fid,[53,inf],'*uint8');
% Select and convert useful info
Out = struct('datNum', num2cell(uint16(tmp(1,:).'))         ,... % # of the .dat file with the data
             'symbol', cellstr(char(tmp(37:52,:).'))        ,... % security symbol
             'name'  , cellstr(char(tmp(8:23,:).'))         ,... % security name
             'inDate', fmtStrDate(mbf2ieee(tmp(26:29,:)))   ,... % initial date as 'yyyy-mm-dd' (datestr format 29)
             'fiDate', fmtStrDate(mbf2ieee(tmp(30:33,:)))   ,... % final   date as 'yyyy-mm-dd' (datestr format 29)
             'freq'  , cellstr(char(tmp(34,:).'))           ,... % time frequency of the data: I, D, M, W, Q, Y
             'idFreq', num2cell(tmp(35,:).'));                   % intraday time frequency in minutes
end

% ------------------------------------------------------------------------------------------------------------
% READEMASTER - Read emaster file
% ------------------------------------------------------------------------------------------------------------
function Out = readEmaster(fid)
% Skip header (each record is 150 bytes long)
fseek(fid,192,'bof');
% Read in the whole file
tmp = fread(fid,[192,inf],'*uint8');
% Select and convert useful info
Out = struct('datNum', num2cell(uint16(tmp(3,:).'))                                 ,...
             'symbol', cellstr(char(tmp(12:25,:).'))                                ,...
             'name'  , cellstr(char(tmp(33:48,:).'))                                ,...
             'inDate', fmtStrDate(typecast(reshape(tmp(65:68,:),[],1),'single'))    ,...
             'fiDate', fmtStrDate(typecast(reshape(tmp(73:76,:),[],1),'single'))    ,...
             'freq'  , cellstr(char(tmp(61,:).'))                                   ,...
             'idFreq', num2cell(tmp(63,:).'));

end

% ------------------------------------------------------------------------------------------------------------
% READXMASTER - Read xmaster file
% ------------------------------------------------------------------------------------------------------------
function Out = readXmaster(fid)
% Skip header (each record is 192 bytes long)
fseek(fid,150,'bof');
% Read in the whole file
tmp = fread(fid,[150,inf],'*uint8');
% Names (trim unuseful chars; observed behavior only in Xmaster)
names = cellfun(@(x) char(x(1:find(x == 0,1,'first'))), num2cell(tmp(17:39,:).',2),'un',0) ;
idx   = cellfun('isempty',names);
names(idx) = cellstr(char(tmp(17:39,idx).'));
% Select and convert useful info
Out = struct('datNum', num2cell(typecast(reshape(tmp(66:67,:),[],1),'uint16'))              ,...
             'symbol', cellstr(char(tmp(2:15,:).'))                                         ,...
             'name'  , names                                                                ,...
             'inDate', fmtStrDate(single(typecast(reshape(tmp(109:112,:),[],1),'uint32')))  ,...
             'fiDate', fmtStrDate(single(typecast(reshape(tmp(117:120,:),[],1),'uint32')))  ,...
             'freq'  , cellstr(char(tmp(63,:).'))                                           ,...
             'idFreq', num2cell(tmp(65,:).'));
end

% ------------------------------------------------------------------------------------------------------------
% MBF2IEEE - Converts MBF to IEEE
% ------------------------------------------------------------------------------------------------------------
function ieee = mbf2ieee(mbf)
% Convert a 32 bit Microsoft Basic Float into a IEEE double of the form (C)YYMMDD (ex: 990123 or 1010123)

% Layout of a MBF:
%
% MBF __________|__Byte 4 __|__Byte 3 __|__Byte 2 __|__Byte 1 __|
% Position bit__| 32 ... 24 | 23 ... 16 | 15 .. . 8 | 7 . . . 0 |
% Value in bit__| EEEE EEEE | XMMM MMMM | MMMM MMMM | MMMM MMMM |
%
% Components:
% X = sign bit
% E = 8 bit exponent
% M = 23 bit mantissa

% Cast into uint32 combining the 4 bytes and convert to double
mbf           = double(typecast(mbf(:),'uint32'));
mantissa_mask =   16777215;                                    % (2^24-1) bottom 24 bits holds the mantissa
exponent_mask = 4278190080;                                    % bitxor(2^32-1,mask1) top 8 bits holds the exponent
mantissa      = bitset  (bitand(mbf,mantissa_mask), 24)    ;   % restore hidden bit
exponent      = bitshift(bitand(mbf,exponent_mask),-24)-152;   % scale exponent
theSign       = double(bitget(mbf,24) == 0)*2 -1 ;             % +1 for zero or positive, -1 for negative
zeroVal       = double(mbf ~= 0) ;                             % 0 for zero values else +1
ieee          = zeroVal .* theSign .* 2.^exponent .* mantissa; % Convert to ieee
end

% ------------------------------------------------------------------------------------------------------------
% FORMATSTRDATE - Converts (C)YYMMDD or YYYYMMDD to datestring format 29
% ------------------------------------------------------------------------------------------------------------
function strdate = fmtStrDate(n)
% Add 1900 if is YYMMDD (ex: 990123 --> 1999-01-23) or CYYMMDD (ex:1010123 --> 2001-01-23) or leave YYYYMMDD (ex: 19990123)
strdate = sprintf('%d-%02d-%02d',[1900*(~fix(n/1e7)) + fix(n/1e4),fix(mod(n/1e2,1e2)),mod(n,1e2)].');
strdate = cellstr(reshape(strdate(:),10,[]).');
% Place 'none' if 0
none = n == 0;
strdate(none) = {'none'};
end

% ------------------------------------------------------------------------------------------------------------
% PARTNUMDATE - Parts (C)YYMMDD or YYYYMMDD to year YYYY, month MM and day DD
% ------------------------------------------------------------------------------------------------------------
function [d,m,y] = partNumDate(n)
% Day, month and year
d = mod(n,1e2);
m = fix(mod(n/1e2,1e2));
y = fix(n/1e4);
% If the number is supplied as YYMMDD (ex: 990123) or CYYMMDD (ex:1010123 --> 2001-01-23) instead of YYYYMMDD (ex: 19990123)
idx = ~fix(n/1e7);
y(idx) = y(idx) +1900;
end

author: Sable

repo: mcbench-benchmarks


text:
function avg = MovAvgFilter2(x)
%
%
persistent n xbuf
persistent firstRun

if isempty(firstRun)
  n = 10;
  xbuf = x*ones(n, 1);

  firstRun = 1;
end

for m=1:n-1
  xbuf(m) = xbuf(m+1);
end
xbuf(n) = x;

avg = sum(xbuf) / n;

author: philbooks

repo: Kalman-Filter-for-Beginners


text:
function t = iseuclidean(D)
%ISEUCLIDEAN Is a distance matrix Euclidean?
%   T = ISEUCLIDEAN(D) returns a logical indicating whether or not the
%   dissimilarity matrix D is a Euclidean distance matrix, i.e., whether
%   there exist n points in p-dimensional space (for some p < n) such that
%   their Euclidean distances are given by D.  D may be specified as either
%   a full (square, symmetric) dissimilarity matrix, or as the lower
%   triangle (e.g., output by PDIST).
%
%   This algorithm is essentially classical multidimensional scaling.
%
%   See also CMDSCALE, PDIST, LINKAGE.
%
%   References:
%     [1] Seber, G.A.F., Multivariate Observations, Wiley, 1984

%   Copyright 1993-2002 The MathWorks, Inc.
%   $Revision: 1.4 $ $Date: 2002/02/04 19:25:44 $

[n m] = size(D);
del = 10*eps;

% lower triangle form for D
if n == 1
    % make sure it's a valid dissimilarity matrix
    n = (1+sqrt(1+8*m))/2;
    if n == fix(n) & all(D >= 0)
        D = squareform(D);
    else
        warning('Not a valid dissimilarity or distance matrix.')
        t = logical(0);
        return
    end

% full matrix form, make sure it's valid dissimilarity matrix
elseif n ~= m | any(any(D < 0 | abs(D - D') > del*max(max(D)))) | any(diag(D) > del)
    warning('Not a valid dissimilarity or distance matrix.')
    t = logical(0);
    return
end

P = eye(n) - repmat(1/n,n,n);
B = P * (-.5 * D .* D) * P;
g = eig((B+B')./2); % guard against spurious complex e-vals from roundoff
t = all(-eps^(3/4)*max(abs(g)) <= g); % all non-negative eigenvals (within roundoff)?

author: zouchuhang

repo: LayoutNet


text:
function [fx,dfdx,dfdp] = f_HRF3(Xt,P,ut,in)
% Balloon (HRF) model evolution function in log-space
% function [fx,dfdx,dfdp] = f_HRF2(Xt,P,ut,in)
% This function evaluates the evolution function derived from the balloon
% model for the hemodynamic response function. It can be called in two
% ways: (i) as a "stand-alone" evolution function, whereby the system's
% states are hemodynamic states of the balloon model, or (ii) as a
% generalized observation function, where the real system's states are the
% neuronal states of a DCM for fMRI model. Note that the hemodynamic states
% are in log-space, for positivity constrints.

deltat = in.deltat;
n = size(Xt,1);

% Get parameters
[E0,V0,tau0,kaf,kas,epsilon,alpha] = BOLD_parameters;
nreg = numel(in.n5); %n./5;
ind1 = in.ind1;
ind2 = in.ind2;
ind3 = in.ind3;
ind4 = in.ind4;
n1 = in.n1;
n2 = in.n2;
n3 = in.n3;
n4 = in.n4;
n5 = in.n5;
epsilon = 1;
alpha = alpha.*exp(P(in.ind5));

E0 = 1./(1+exp(-(P(ind1)-0.6633)));
dsdp = E0.*(1-E0);
tau0 = tau0.*exp(P(ind2));
kaf = kaf.*exp(P(ind3));
kas = kas.*exp(P(ind4));

% hemodynamic states ...
if in.linearized
    x1 = zeros(nreg,1);
    x2 = ones(nreg,1);
    x3 = ones(nreg,1);
    x4 = ones(nreg,1);
else
    % vasodilatory signal s(t)
    x1 = Xt(n1,:);
    % blood inflow f(t)
    if ~in.logx2
        x2 = Xt(n2,:) + 1; % deviation to steady-state!
    else
        x2 = exp(Xt(n2,:)) + in.xshift;
    end
    % blood volume v(t)
    x3 = exp(Xt(n3,:)) + 0.*in.xshift;
    % dHb content q(t)
    x4 = exp(Xt(n4,:)) + 0.*in.xshift;
end
% blood outflow
fv = x3.^(1./alpha);
% d[blood flow]/dXt(3)
dfvdx = (1./alpha).*fv;
% oxygen extraction
ff = (1-(1-E0).^(1./x2))./E0;
% d[O2 extraction]/dXt(2)
if ~in.logx2
    dffdx = log(1-E0).*(1-E0).^(1./x2)./(E0.*x2.^2);
else
    dffdx = log(1-E0).*(1-E0).^(1./x2)./(E0.*x2);
end
% ... and flow field, derivatives, etc...
f = zeros(n,1);
J = zeros(n,n);
dfdp = zeros(size(P,1),n);

% Evaluate flow field
f(n1) = (epsilon.*ut - kas.*x1 - kaf.*(x2 - 1));
if ~in.logx2
    f(n2) = x1;
else
    f(n2) = x1./x2;
end
f(n3) = (x2 - fv)./(tau0.*x3);
f(n4) = (x2.*ff./x4 - fv./x3)./tau0;

% Evaluate jacobian and gradients wrt parameters
for i=1:nreg

    if ~in.logx2
        J(n1(i),n1(i):n1(i)+3) = [ -kas(i) , 1, 0, 0 ];
        J(n2(i),n1(i):n1(i)+3) = [ -kaf(i), 0, ...
            1./(tau0(i).*x3(i)), ...
            (ff(i)+x2(i).*dffdx(i))./(x4(i).*tau0(i))];
    else
        J(n1(i),n1(i):n1(i)+3) = [ -kas(i) , 1./x2(i), 0, 0 ];
        J(n2(i),n1(i):n1(i)+3) = [ -kaf(i).*x2(i), -x1(i)./x2(i), ...
            x2(i)./(tau0(i).*x3(i)), ...
            x2(i).*(ff(i)+dffdx(i))./(tau0(i).*x4(i))];
    end
    J(n3(i),n1(i):n1(i)+3) = [ 0, 0,...
        -f(n3(i)) - dfvdx(i)./(tau0(i).*x3(i)) , ...
        (fv(i)-dfvdx(i))./(tau0(i).*x3(i))];
    J(n4(i),n1(i):n1(i)+3) = [ 0, 0, 0, ...
        -(x2(i).*ff(i))./(tau0(i).*x4(i)) ];
    J(n5(i),n1(i)) = epsilon;

    if in.linearized

        tmp = log(1-E0(i))./E0(i);

        dfdp(ind1(i),n1(i):n1(i)+3) = [0,0,0,...
            -dsdp(i).*(1+tmp).*Xt(n2(i))./(tau0(i).*E0(i))];
        dfdp(ind2(i),n1(i):n1(i)+3) = [0,0,...
            (-Xt(n2(i))+Xt(n3(i))./alpha(i))./tau0(i),...
            -((ff(i)+dffdx(i)).*Xt(n2(i)) + (fv(i)-dfvdx(i)).*Xt(n3(i)) ...
            - ff(i).*Xt(n4(i)))./tau0(i)];
        dfdp(ind3(i),n1(i):n1(i)+3) = kaf(i).*[-Xt(n2(i)),0,0,0];
        dfdp(ind4(i),n1(i):n1(i)+3) = kas(i).*[-Xt(n1(i)),0,0,0];
        dfdp(in.ind5(i),n1(i):n1(i)+3) = [0,0,...
            Xt(n3(i))./(tau0(i).*alpha(i)),...
            Xt(n3(i))./(tau0(i).*alpha(i))];

    else

        % gradient wrt parameters
        dfdp(ind1(i),n1(i):n1(i)+3) = [0,0,0,...
            (((1-E0(i)).^(-1+1./x2(i)))-x2(i).*ff(i))...
            .*dsdp(i)./(tau0(i).*x4(i).*E0(i))];
        dfdp(ind2(i),n1(i):n1(i)+3) = [0,0,...
            -(x2(i) - fv(i))./(tau0(i).*x3(i)),...
            -(x2(i).*ff(i)./x4(i) - fv(i)./x3(i))./tau0(i)];
        dfdp(ind3(i),n1(i):n1(i)+3) = kaf(i).*[-x2(i)+1,0,0,0];
        dfdp(ind4(i),n1(i):n1(i)+3) = kas(i).*[-x1(i),0,0,0];
        dfdp(in.ind5(i),n1(i):n1(i)+3) = [0,0,...
            log(x3(i)).*fv(i)./(tau0(i).*x3(i).*alpha(i)),...
            log(x3(i)).*fv(i)./(tau0(i).*x3(i).*alpha(i))];

    end

end

% Apply Euler discretization
dfdp = deltat.*dfdp;
dfdx = eye(n) + deltat.*J;
if in.linearized
    fx = dfdx'*Xt;
else
    fx = Xt + deltat.*f;
end

author: MBB-team

repo: VBA-toolbox


text:
function [ t_stop, y_stop ] = p09_stop ( neqn )

%*****************************************************************************80
%
%% P09_STOP returns the stopping point for problem p09.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_STOP, Y_STOP(NEQN), the final data.
%
  y_stop = zeros ( neqn, 1 );

  t_stop = 20.0;
  y_stop = [ ...
    9.826950928005993e-01; ...
    2.198447081694832; ...
    9.129452507276399e-01 ];

  return
end


text:
function [y]=gnanmean(varargin)
% function [y]=gnanmean(A,dimDir)
%------------------------------------------------------------------------
% this function is a replacement for nanmean which is part of the
% statistics and machine learning toolbox.
% The use of mean with the omitnan flag is used or if this does not work
% (for old MATLAB versions) a custom implementation is used.
%
% Change log:
% 2020/01/09 Created
%------------------------------------------------------------------------

%%
try
    %Use 'omitnan' flag
    y = mean(varargin{:},'omitnan');
catch
    %Use custom version
    A=varargin{1};
    if nargin==2
        dimDir=varargin{2};
    else
        dimDir=1;
    end
    L=isnan(A);
    Q=sum(~L,dimDir);
    B=A;
    B(L)=0;
    S=sum(B,dimDir);
    y=S./Q;
    y(all(L,dimDir))=NaN;
end

end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [ colorspace ] = get_colorspace( im, fparam, gparam )

    [im_height, im_width, num_im_chan, num_images] = size(im);

    single_im = single(im)/255;

    if strcmp(fparam.colorspace,'gray')
        if num_im_chan == 3
            if num_images == 1
                t_colorspace = rgb2gray(single_im) - 0.5;
            else
                t_colorspace = zeros(im_height, im_width, 1, num_images, 'single');
                for k = 1:num_images
                    t_colorspace(:,:,:,k) = rgb2gray(single_im(:,:,:,k)) - 0.5;
                end
            end
        elseif num_im_chan == 1
            t_colorspace = single_im - 0.5;
        else
            except = MException('get_colorspace','Invalid input data, must have 1 or 3 dimensions');
            throw(except);
        end;
    elseif strcmp(fparam.colorspace,'rgb')
        if num_im_chan == 3
            t_colorspace = single_im - 0.5;
        else
            except = MException('get_colorspace','Invalid input data, must have 3 dimensions for rgb');
            throw(except);
        end;
    end;

    if gparam.use_gpu
        t_colorspace = gpuArray(t_colorspace);
    end

    if gparam.cell_size > 1
        colorspace = average_feature_region(t_colorspace,gparam.cell_size);
    else
        colorspace = t_colorspace;
    end;
end

repo: STRCF


text:
function test_suite = testContainsRegexp
%testContainsRegexp Unit tests for containsRegexp

%   Steven L. Eddins
%   Copyright 2008 The MathWorks, Inc.

initTestSuite;

function testOneStringContains
assertTrue(xunit.utils.containsRegexp('MATLAB is great', '[A-Z]'));

function testOneStringDoesntContain
assertTrue(~ xunit.utils.containsRegexp('no upper-case letters', '[A-Z]'));

function testCellArray
strs = {'MATLAB is great', 'no upper-case letters'};
assertEqual(xunit.utils.containsRegexp(strs, '[A-Z]'), [true false]);
assertEqual(xunit.utils.containsRegexp(strs', '[A-Z]'), [true; false]);

author: ricedsp


text:
function [ x, y, z ] = r8_unswap3 ( x, y, z )

%*****************************************************************************80
%
%% R8_UNSWAP3 unswaps three R8's.
%
%  Example:
%
%    Input:
%
%      X = 2, Y = 3, Z = 1
%
%    Output:
%
%      X = 1, Y = 2, Z = 3
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, Y, Z, three values to be swapped.
%
%    Output, real X, Y, Z, the swapped values.
%
  w = z;
  z = y;
  y = x;
  x = w;

  return
end


text:
function model = rbfOptimise(model, X, Y, display, iters);

% RBFOPTIMISE Optimise RBF for given inputs and outputs.

% MLTOOLS

if nargin < 4
  display = 1;
  if nargin < 5
    iters = 500;
  end
end

options = optOptions;
options(14) = iters;
options(1) = display;
model = netopt(model, options, X, Y, 'scg');

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function class = gmm_classification(xtest, mu_train, sigma_train, c_train)
    % Classify new points in the 'xtest' set using GMM matching
    % Each column in xtest represents a point to be classified.
    % Each column in xtrain contains a training sample.
    % ytrain - column vector.

%     global NUM_OF_GAUSSIANS;
%     global NUM_OF_ITERATIONS;
%     [mu_train, sigma_train, c_train] = ...
%         gmm_training(xtrain, ytrain, NUM_OF_GAUSSIANS, NUM_OF_ITERATIONS);

    NUM_OF_TEST_SAMPLES = size(xtest, 2);
    class = zeros(NUM_OF_TEST_SAMPLES, 1);
    for k = 1:NUM_OF_TEST_SAMPLES
        [class(k), score] = GMM.gmm_classify_sample(mu_train, sigma_train, ...
            c_train, xtest(:, k));
    end;
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function FMatrix=kannumfcc(num,s,Fs)
%Author:        Olutope Foluso Omogbenigun
%Email:         olutopeomogbenigun at hotmail.com
%University:    London Metropolitan University
%Date:          11/09/07
%Syntax:        M=mfccf(num,s, Fs);
%Computes and returns the mfcc coefficients for a speech signal s
%where num is the required number of MFCC coefficients. It utilises the
%function 'melbankm' from the toolbox Voicebox by Mike Brooks copyright(c)
%1997 (GNU General Public License), freely available on the internet,
%to implement the triangular mel filter bank

n=512;              %Number of FFT points
Tf=0.025;           %Frame duration in seconds
N=Fs*Tf;            %Number of samples per frame
fn=24;              %Number of mel filters
l=length(s);        %total number of samples in speech
Ts=0.01;            %Frame step in seconds
FrameStep=Fs*Ts;    %Frame step in samples
a=1;
b=[1, -0.97];       %a and b are high pass filter coefficients

noFrames=floor(l/FrameStep);    %Maximum no of frames in speech sample
FMatrix=zeros(noFrames-2, num); %Matrix to hold cepstral coefficients
lifter=1:num;                   %Lifter vector index
lifter=1+floor((num)/2)*(sin(lifter*pi/num));%raised sine lifter version

if mean(abs(s)) > 0.01
    s=s/max(s);                     %Normalises to compensate for mic vol differences
end

%Segment the signal into overlapping frames and compute MFCC coefficients
for i=1:noFrames-2
    frame=s((i-1)*FrameStep+1:(i-1)*FrameStep+N);  %Holds individual frames
    Ce1=sum(frame.^2);          %Frame energy
    Ce2=max(Ce1,2e-22);         %floors to 2 X 10 raised to power -22
    Ce=log(Ce2);
    framef=filter(b,a,frame);   %High pass pre-emphasis filter
    F=framef.*hamming(N);       %multiplies each frame with hamming window
    FFTo=fft(F,N);              %computes the fft
    melf=melbankm(fn,n,Fs);     %creates 24 filter, mel filter bank
    halfn=1+floor(n/2);
    spectr1=log10(melf*abs(FFTo(1:halfn)).^2);%result is mel-scale filtered
    spectr=max(spectr1(:),1e-22);
    c=dct(spectr);              %obtains DCT, changes to cepstral domain
    c(1)=Ce;                    %replaces first coefficient
    coeffs=c(1:num);            %retains first num coefficients
    ncoeffs=coeffs.*lifter';    %Multiplies coefficients by lifter value
    FMatrix(i, :)=ncoeffs';     %assigns mfcc coeffs to succesive rows i
end


text:
function tf=isNonFailure(obj)
% return whether this test was a non-failure
%
% tf=isSuccess(obj)
%
% Input:
%   obj                     MoxUnitFailedTestOutcome object
%
% Output:
%   tf                      set to false, because a failed test is a
%                           failure
%
% See also: @MoxUnitTestOutcome

    tf=false;


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

% Cleaning
clear all
close all
clc

% Gypsilab path
run('../../addpathGypsilab.m')

% Parameters
N   = 1e2;
gss = 3;

% Spherical mesh
sphere = mshSphere(N,1);

% Square mesh
square = mshSquare(2*N,[3 3]);

% Graphical representation
figure
plot(sphere)
hold on
plot(square)
axis equal

% Domain
sigma = dom(sphere,gss);

% Finite elements
u = fem(sphere,'P0');
v = fem(sphere,'P1');
w = fem(sphere,'RWG');

% Gren kernel
Gxy   = @(X,Y) femGreenKernel(X,Y,'[1/r]',[]);
dyGxy = @(X,Y) femGreenKernel(X,Y,'grady[1/r]1',[]);

% Single radiation P0
M  = integral(square.vtx,sigma,Gxy,u);
Mr = regularize(square.vtx,sigma,'[1/r]',u);
norm(M+Mr,'inf')

% Single radiation P1
M  = integral(square.vtx,sigma,Gxy,v);
Mr = regularize(square.vtx,sigma,'[1/r]',v);
norm(M+Mr,'inf')

% Single radiation RWG
M  = integral(square.vtx,sigma,Gxy,div(w));
Mr = regularize(square.vtx,sigma,'[1/r]',div(w));
norm(M+Mr,'inf')

% Single layer P0
M  = integral(sigma,sigma,u,Gxy,u);
Mr = regularize(sigma,sigma,u,'[1/r]',u);
norm(M+Mr,'inf')

% Single layer P1
M  = integral(sigma,sigma,v,Gxy,v);
Mr = regularize(sigma,sigma,v,'[1/r]',v);
norm(M+Mr,'inf')

% Single layer RWG
M  = integral(sigma,sigma,div(w),Gxy,div(w));
Mr = regularize(sigma,sigma,div(w),'[1/r]',div(w));
norm(M+Mr,'inf')

% Double radiation P0
M  = integral(square.vtx,sigma,dyGxy,u);
Mr = regularize(square.vtx,sigma,'grady[1/r]1',u);
norm(M+Mr,'inf')

% Double layer P0
M   = integral(sigma,sigma,u,dyGxy,u);
Mr  = regularize(sigma,sigma,u,'grady[1/r]1',u);
norm(M+Mr,'inf')

% Stokeslet P0
for i = 1:3
    for j = 1:3
        name  = ['[ij/r+rirj/r^3]',num2str(i),num2str(j)];
        green = @(X,Y) femGreenKernel(X,Y,name,[]);
        M     = integral(sigma,sigma,u,green,u);
        Mr    = regularize(sigma,sigma,u,name,u);
        norm(M+Mr,'inf')
    end
end

disp('~~> Michto gypsilab !')

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% Color Mixture Non-local Pixel Affinities
% This function implements the color-mixture information flow in
% Yagiz Aksoy, Tunc Ozan Aydin, Marc Pollefeys, "Designing Effective
% Inter-Pixel Information Flow for Natural Image Matting", CVPR, 2017
% when the input parameter 'useXYinLLEcomp' is false (default), and
% the affinity definition used in
% Xiaowu Chen, Dongqing Zou, Qinping Zhao, Ping Tan, "Manifold
% preserving edit propagation", ACM TOG, 2012
% when 'useXYinLLEcomp' is true.
% All parameters other than image are optional. The output is a sparse
% matrix which has non-zero element for the non-local neighbors of
% the pixels given by binary map inMap.
% - K defines the number of neighbors from which LLE weights are
%   computed.
% - outMap is a binary map that defines where the nearest neighbor
%   search is done.
% - xyWeight determines how much importance is given to the spatial
%   coordinates in the nearest neighbor selection.

function Wcm = colorMixtureAffinities(image, K, inMap, outMap, xyWeight, useXYinLLEcomp)

    [h, w, ~] = size(image);
    N = h * w;

    if ~exist('K', 'var') || isempty(K)
        K = 20;
    end
    if ~exist('inMap', 'var') || isempty(inMap)
        inMap = true(h, w);
    end
    if ~exist('outMap', 'var') || isempty(outMap)
        outMap = true(h, w);
    end
    if ~exist('xyWeight', 'var') || isempty(xyWeight)
        xyWeight = 1;
    end
    if ~exist('useXYinLLEcomp', 'var') || isempty(useXYinLLEcomp)
        useXYinLLEcomp = false;
    end

    [inInd, neighInd, features] = findNonlocalNeighbors(image, K, xyWeight, inMap, outMap);

    if ~useXYinLLEcomp
        features = features(:, 1 : end - 2);
    end
    flows = zeros(size(inInd, 1), size(neighInd, 2));

    for i = 1 : size(inInd, 1)
        flows(i, :) = localLinearEmbedding(features(inInd(i), :)', features(neighInd(i, :), :)', 1e-10);
    end
    flows = flows ./ repmat(sum(flows, 2), [1, K]);

    inInd = repmat(inInd, [1, K]);
    Wcm = sparse(inInd(:), neighInd(:), flows, N, N);
end

author: yaksoy

repo: AffinityBasedMattingToolbox


text:
function arvar=arloop(data_endo,const,p,n)

% function arvar=arloop(data_endo,const,p,n)
% computes individual OLS estimations of autoregressive models and record their residual variances, as stated p16 of technical guide
% inputs:  - matrix 'data_endo': matrix of endogenous data used for model estimation
%          - integer 'const': 0-1 value to determine if a constant is included in the model
%          - integer 'p': number of lags included in the model (defined p 7 of technical guide)
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
% outputs: - vector 'arvar': residual variance of individual AR models estimated for each endogenous variable

arvar = zeros(n, 1);
%loop over the columns of data_endo
for ii=1:n
% estimate an AR model with p lags for each series
% and record residual variance in vector ARvar
[~,~,arvar(ii,1),~,~,~,~,~,~,~,~,~,~,~,~]=bear.olsvar(data_endo(:,ii),[],const,p);
end

author: european-central-bank

repo: BEAR-toolbox


text:
function info(DEM)

%INFO detailed information on GRIDobj instance
%
% Syntax
%
%     info(DEM)
%
% Description
%
%     info displays detailed information about an instance of an GRIDobj in
%     the command window.
%
% Input arguments
%
%     DEM    instance of GRIDobj
%
% See also: DISP
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 8. January, 2013

[ul] = [1 1 1] * DEM.refmat;
[lr] = [DEM.size 1] * DEM.refmat;

disp(' ')
disp('TopoToolbox GRIDobj')
disp(['  name:                  ' DEM.name])
disp(['  data type:             ' class(DEM.Z)])
disp(['  number of rows:        ' num2str(DEM.size(1))])
disp(['  number of columns:     ' num2str(DEM.size(2))])
disp(['  cellsize:              ' num2str(DEM.cellsize)])
disp(['  extent in x-direction: ' num2str(ul(1)) ' -- ' num2str(lr(1))])
disp(['  extent in y-direction: ' num2str(ul(2)) ' -- ' num2str(lr(2))])
disp(['  maximum z-value:       ' num2str(max(DEM.Z(:)))])
disp(['  minimum z-value:       ' num2str(min(DEM.Z(:)))])
disp(['  z-unit:                ' DEM.zunit])

if ~isempty(DEM.georef);
    try
        disp(['  coordinate system:     ' DEM.georef.GeoKeyDirectoryTag.GTCitationGeoKey])
    catch
        try
            disp(['  coordinate system:     ' DEM.georef.GeoKeyDirectoryTag.GeogCitationGeoKey]);
        catch
        end
    end
else
    disp(['  coordinate system:     ' 'undefined'])
end

disp(' ')

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
function diri = sb_find_elec(vol,sens)

% SB_FIND_ELEC
%
% $Id$

diri = zeros(size(sens.elecpos,1),1);
dist = zeros(size(sens.elecpos,1),1);
for i=1:size(sens.elecpos,1)
    [dist(i), diri(i)] = min(sum(bsxfun(@minus,vol.pos,sens.elecpos(i,:)).^2,2));
end
if ~all(dist < 1e-8)
    error('Electrode positions are not located on mesh nodes! This should not happen, please contact support.');
end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


% Estimates the integer frequency offset using the time domain samples of
% the first ZC sequence (root index 600)
%
% @param received_zc_td Time domain of the received ZC sequence with no
%                       cyclic prefix (should be row vector)
% @param sample_rate Sample rate of the samples in Hz
% @param max_carrier_offset How much integer frequency offset to search for
%                           (must be positive non-zero value)
% @return integer_offset_hz Best estimate for the integer frequency offset
%                           in Hz
function [integer_offset_hz] = est_integer_freq_offset(received_zc_td, sample_rate, max_carrier_offset)
    fft_size = get_fft_size(sample_rate);

    assert(isrow(received_zc_td), "ZC time domain samples must be row vector");
    assert(length(received_zc_td) == fft_size, "ZC time domain samples must be FFT samples wide");
    assert(isnumeric(sample_rate) && sample_rate > 0, "Sample rate must be positive and non-zero");
    assert(isnumeric(max_carrier_offset) && max_carrier_offset > 0, "Max carrier offset must be positive and non-zero");

    % Calculate all of the cyclic shifts that will be done to search for
    % the integer frequency offset value
    shift_range = -max_carrier_offset:1:max_carrier_offset;

    % Allocate storage for the correlation scores from each shift
    scores = zeros(1, length(shift_range));

    % Move the ZC sequence to the frequency domain
    zc_fd = fft(received_zc_td);

    for idx=1:length(shift_range)
        % Circularly shift the FFT
        zc_fd_shifted = circshift(zc_fd, shift_range(idx));

        % Take the upper and lower `(fft_size / 2) - 1` samples
        left = zc_fd_shifted(1:(fft_size / 2) - 1);
        right = zc_fd_shifted((fft_size / 2) + 2:end);

        % Correlate the lower samples with the mirror of the upper samples
        scores(idx) = xcorr(left, fliplr(right), 0);
    end

    % Find the index of the max score
    [~, index] = max(abs(scores));

    % Calculate the number of Hz between each FFT bin
    carrier_spacing = sample_rate / fft_size;

    % Get how many carriers the best score was shifted by and multiply that
    % value by the carrier spacing to get the integer frequency offset
    integer_offset_hz = shift_range(index) * carrier_spacing;
end


text:
function value = he_double_product_integral ( i, j )

%*****************************************************************************80
%
%% HE_DOUBLE_PRODUCT_INTEGRAL: integral of He(i,x)*He(j,x)*e^(-x^2/2).
%
%  Discussion:
%
%    VALUE = integral ( -oo < x < +oo ) He(i,x)*He(j,x) exp(-x^2/2) dx
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, J, the polynomial indices.
%
%    Output, real VALUE, the value of the integral.
%
  if ( i ~= j )
    value = 0.0;
  else
    value = r8_factorial ( i );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [y]=dmrg_cross(d,n,fun,eps,varargin)
%DMRG-cross method for the approximation of TT-tensors
%   [A]=DMRG_CROSS(D,N,FUN,EPS,OPTIONS) Computes the approximation of a
%   given tensor via the adaptive DMRG-cross procedure. The input is a pair
%   (D,N) which determines the size of the tensor (N can be either a
%   number, or array of mode sizes). FUN is the function to compute
%   a prescribed element of a tensor (FUN(IND)), or it can be vectorized to
%   compute series of elements of a tensor (see OPTIONS) To pass parameters
%   to FUN please use anonymous function handles. EPS is the accuracy
%   of the approximation.Options are provided in form
%   'PropertyName1',PropertyValue1,'PropertyName2',PropertyValue2 and so
%   on. The parameters are set to default (in brackets in the following)
%   The list of option names and default values are:
%       o nswp - number of DMRG sweeps [10]
%       o vec  - Fun is vectorized [ true | {false} ]
%       o verb - output debug information [ {true} | false ]
%       o y0   - initial approximation [random rank-2]
%       o radd - minimal rank change [0]
%       o rmin - minimal rank that is allows [1]
%       o kickrank - stabilization parameter [2]
%
%   Example:
%       d=10; n=2; fun = @(ind) sum(ind);
%       tt=dmrg_cross(d,n,fun,1e-7);
%
%  This code implements the algorithm from
%  D. Savostyanov, I. Oseledets, Fast adaptive interpolation of multi-dimensional arrays in tensor train format,
%  http://dx.doi.org/10.1109/nDS.2011.6076873
%  Please cite this paper if your research benefits from the use of this code.
%
%  Vectorized version contributed by Prof. Le Song (http://www.cc.gatech.edu/~lsong/)
%
%
% TT-Toolbox 2.2, 2009-2013
%
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------
%Default parameters
rmin=1;
verb=true;
radd=0;
kickrank=2;
nswp=10;
y=[];
vectorized=false;
maxr = 20;
for i=1:2:length(varargin)-1
    switch lower(varargin{i})
        case 'nswp'
            nswp=varargin{i+1};
        case 'y0'
            y=varargin{i+1};
        case 'verb'
            verb=varargin{i+1};
        case 'rmin'
            rmin=varargin{i+1};
        case 'radd'
            radd=varargin{i+1};
        case 'vec'
            vectorized=varargin{i+1};
        case 'kickrank'
            kickrank=varargin{i+1};
      case 'maxr'
        maxr = varargin{i+1};

        otherwise
            error('Unrecognized option: %s\n',varargin{i});
    end
end

if ( numel(n) == 1 )
   n=n*ones(d,1);
end

sz=n;
if (isempty(y) )
    y=tt_rand(sz,d,2);
end
if ( ~vectorized )
    elem=@(ind) my_vec_fun(ind,fun);
else
    elem = fun;
end
y=round(y,0); %To avoid overranks
ry=y.r;
[y,rm]=qr(y,'rl');
y=rm*y;
%Warmup procedure: orthogonalization from right to left of the initial
%approximation & computation of the index sets & computation of the
%right-to-left R matrix
swp=1;
rmat=cell(d+1,1);
rmat{d+1}=1;
rmat{1}=1; %These are R-matrices from the QR-decomposition.
index_array{d+1}=zeros(0,ry(d+1));
index_array{1}=zeros(ry(1),0);
r1=1;
for i=d:-1:2
    cr=y{i}; cr=reshape(cr,[ry(i)*n(i),ry(i+1)]);
    cr = cr*r1; cr=reshape(cr,[ry(i),n(i)*ry(i+1)]); cr=cr.';
    [cr,rm]=qr(cr,0);
    [ind]=maxvol2(cr);
    ind_old=index_array{i+1};
    rnew=min(n(i)*ry(i+1),ry(i));
    ind_new=zeros(d-i+1,rnew);
    for s=1:rnew
       f_in=ind(s);
       w1=tt_ind2sub([ry(i+1),n(i)],f_in);
       rs=w1(1); js=w1(2);
       ind_new(:,s)=[js,ind_old(:,rs)'];
    end
    index_array{i}=ind_new;
    r1=cr(ind,:);
    cr=cr/r1;
    r1=r1*rm;
    r1=r1.';

    cr=cr.';
    y{i}=reshape(cr,[ry(i),n(i),ry(i+1)]);
    cr=reshape(cr,[ry(i)*n(i),ry(i+1)]);
    cr=cr*rmat{i+1}; cr=reshape(cr,[ry(i),n(i)*ry(i+1)]);
    cr=cr.';
    [~,rm]=qr(cr,0);
    rmat{i}=rm; %The R-matrix
end
%Forgot to put r1 onto the last core
cr=y{1}; cr=reshape(cr,[ry(1)*n(1),ry(2)]);
y{1}=reshape(cr*r1,[ry(1),n(1),ry(2)]);
not_converged = true;
dir = 1; %The direction of the sweep
i=1; %Current position
er_max=0;
while ( swp < nswp && not_converged )
    % A sweep through the cores
    %Compute the current index set, compute the current supercore
    %(right now without any 2D cross inside, but it is trivial to
    %implement). The supercore is (i,i+1) now.
    %Left index set is index_array{i}, right index set is index_array{i+2}
    %We will modify ry(i+1) at this step and use rmat{i} and rmat{i+2}
    %as "weighting" matrices for the low-rank approximation. The initial
    %approximation is simply rmax{i}*u{i}*u{i+1}*rmat{i+2} (hey!)
    %We also have to store the submatrix in the current factors
    %Then the algorithm would be as follows: Computex sets, compute
    %supercore. Compute rmax{i}*Phi*rmax{i+2} = U*V by SVD, then split
    rm1=rmat{i}; rm2=rmat{i+2};
    cr1=y{i}; cr2=y{i+1};
    ind1=index_array{i};
    ind2=index_array{i+2};

%     big_index=zeros(ry(i),n(i),n(i+1),ry(i+2),d);
%     for i1=1:n(i)
%         for i2=1:n(i+1)
%             for s1=1:ry(i)
%                 for s2=1:ry(i+2)
%                     ind=[ind1(s1,:),i1,i2,ind2(:,s2)'];
%                     big_index(s1,i1,i2,s2,:)=ind;
%                 end
%             end
%         end
%     end
%     big_index=reshape(big_index,[numel(big_index)/d,d]);

    big_index = [ ...
      ind1(repmat((1:ry(i))', n(i)*n(i+1)*ry(i+2), 1),:), ...
      kron(repmat((1:n(i))', n(i+1)*ry(i+2), 1), ones(ry(i),1)), ...
      kron(repmat((1:n(i+1))', ry(i+2), 1), ones(ry(i)*n(i),1)), ...
      ind2(:, kron((1:ry(i+2))', ones(ry(i)*n(i)*n(i+1),1)))' ...
      ];

    score=elem(big_index);

    %Now plug in the rmax matrices
    score=reshape(score,[ry(i),n(i)*n(i+1)*ry(i+2)]);
    score=rmat{i}*score;
    ry(i)=size(score,1);
    score=reshape(score,[ry(i)*n(i)*n(i+1),ry(i+2)]);
    score=score*rmat{i+2};
    ry(i+2)=size(score,2);

    %Do the SVD splitting (later on we can replace it by cross for large
    %mode sizes)
    score=reshape(score,[ry(i)*n(i),n(i+1)*ry(i+2)]);

%     [u,s,v]=svd(score,'econ');
    [u,s,v]=svds(score, maxr);
    s=diag(s);
    r=my_chop2(s,norm(s)*eps/sqrt(d-1)); %Truncation
    u=u(:,1:r); v=v(:,1:r); s=diag(s(1:r));

    %Kick rank
    if ( dir == 1 )
%         v=v*diag(s);
        v = v * s';

        ur=randn(size(u,1),kickrank);
        u=reort(u,ur);
        radd=size(u,2)-r;
        if ( radd > 0 )
            vr=zeros(size(v,1),radd);
            v=[v,vr];
        end
        r=r+radd;
    else
%          u=u*diag(s);
         u = u * s;

         vr=randn(size(v,1),kickrank);
         v=reort(v,vr);
         radd=size(v,2)-r;
         if ( radd > 0 )
             ur=zeros(size(u,1),radd);
             u=[u,ur];
         end
         r=r+radd;
    end

    v=v';

%     size(v)

    %Compute the previous approximation
    appr=reshape(cr1,[numel(cr1)/ry(i+1),ry(i+1)])*reshape(cr2,[ry(i+1),numel(cr2)/ry(i+1)]);
    appr=reshape(appr,[ry(i),n(i)*n(i+1)*ry(i+2)]);
    appr=rmat{i}*appr;
    appr=reshape(appr,[ry(i)*n(i)*n(i+1),ry(i+2)]);
    appr=appr*rmat{i+2};
    er_loc=norm(score(:)-appr(:))/norm(score(:));
    er_max=max(er_max,er_loc);
    if ( verb )
        fprintf('swp=%d block=%d new_rank=%d local_er=%3.1e\n',swp,i,r,er_loc);
    end
    ry(i+1)=r;

    u = reshape(u,[ry(i),n(i)*r]);
    u = rmat{i}\u; %Hope it is stable blin
    v=reshape(v,[r*n(i+1),ry(i+2)]);
    u=reshape(u,[ry(i)*n(i),ry(i+1)]);
    v=v/rmat{i+2}; v=reshape(v,[r,n(i+1)*ry(i+2)]);
    if ( dir == 1 )
        [u,rm]=qr(u,0);
        ind=maxvol2(u);
        r1=u(ind,:);
        u=u/r1; y{i}=reshape(u,[ry(i),n(i),ry(i+1)]);
        r1=r1*rm;
        v=r1*v; y{i+1}=reshape(v,[ry(i+1),n(i+1),ry(i+2)]);
        %Recalculate rmat
        u1=reshape(u,[ry(i),n(i)*ry(i+1)]);
        u1=rmat{i}*u1;
        u1=reshape(u1,[ry(i)*n(i),ry(i+1)]);
        [~,rm]=qr(u1,0);
        rmat{i+1}=rm;
        %Recalculate index array
        ind_old=index_array{i};
        ind_new=zeros(ry(i+1),i);
        for s=1:ry(i+1)
            f_in=ind(s);
            w1=tt_ind2sub([ry(i),n(i)],f_in);
            rs=w1(1); js=w1(2);
            ind_new(s,:)=[ind_old(rs,:),js];
        end
        index_array{i+1}=ind_new;
        if ( i == d - 1 )
            dir = -dir;
        else
            i=i+1;
        end
    else %Reverse direction
         v=v.'; %v is standing
        [v,rm]=qr(v,0);
        ind=maxvol2(v);
        r1=v(ind,:);
        v=v/r1; v2=reshape(v,[n(i+1),ry(i+2),ry(i+1)]); y{i+1}=permute(v2,[3,1,2]);
        r1=r1*rm; r1=r1.';
        u=u*r1; y{i}=reshape(u,[ry(i),n(i),ry(i+1)]);
        %Recalculate rmat
        v=v.';
        v=reshape(v,[ry(i+1)*n(i+1),ry(i+2)]);
        v=v*rmat{i+2};
        v=reshape(v,[ry(i+1),n(i+1)*ry(i+2)]); v=v.';
        [~,rm]=qr(v,0);
        rmat{i+1}=rm;
        %Recalculate index array
        ind_old=index_array{i+2};
        ind_new=zeros(d-i,ry(i+1));
        for s=1:ry(i+1);
            f_in=ind(s);
            w1=tt_ind2sub([n(i+1),ry(i+2)],f_in);
            rs=w1(2); js=w1(1);
            ind_new(:,s)=[js,ind_old(:,rs)'];
        end
        index_array{i+1}=ind_new;
        if ( i == 1 )
            dir=-dir;
            swp = swp + 1;
            if ( er_max < eps )
                not_converged=false;
            else
                er_max=0;
            end
        else
            i=i-1;
        end
    end
end
return
end
function val = my_vec_fun(ind, fun)
%Trivial vectorized computation of the elements of a tensor
%   [VAL]=MY_VEC_FUN(IND,FUN) Given a function handle FUN, compute all
%   elements of a tensor given in the index array IND. IND is a M x d
%   array, where M is the number of indices to be computed.
M = size(ind, 1);
val = zeros(1, M);
for i = 1:M
    ind_loc = ind(i,:); ind_loc=ind_loc(:);
%   ind_loc = ind(:,i);
   val(i) = fun(ind_loc);
end
return
end

author: oseledets

repo: TT-Toolbox


text:
clear

if(isunix)
    features_exe = '"../../../build/bin/FeatureExtraction"';
else
    features_exe = '"../../../x64/Release/FeatureExtraction.exe"';
end

find_FERA2011;

output_dir = 'E:\datasets\face_datasets_processed\fera2011';
FERA2011_dir = [FERA2011_dir, '/../'];

% Go two levels deep
fera_dirs = dir([FERA2011_dir, '*']);
fera_dirs = fera_dirs(3:end);
parfor f1=1:numel(fera_dirs)

    fera_dirs_level_2 = dir([FERA2011_dir, fera_dirs(f1).name]);
    fera_dirs_level_2 = fera_dirs_level_2(3:end);

    for f2=1:numel(fera_dirs_level_2)

        vid_files = dir([FERA2011_dir, fera_dirs(f1).name, '/', fera_dirs_level_2(f2).name, '/*.avi']);

        for v=1:numel(vid_files)

            input_file = [FERA2011_dir, fera_dirs(f1).name, '/', fera_dirs_level_2(f2).name, '/', vid_files(v).name];
            [~, name, ~] = fileparts(vid_files(v).name);
            out_name = [fera_dirs(f1).name '_' name];
            command = sprintf('%s -f "%s" -out_dir "%s" -hogalign -pdmparams -of %s', features_exe, input_file, output_dir, out_name );

            dos(command);

        end

    end

end

author: TadasBaltrusaitis

repo: OpenFace