text:
function fem2d_scalar_display_gpl ( prefix )

%*****************************************************************************80
%
%% FEM2D_SCALAR_DISPLAY_GPL creates surface plots of 2D FEM scalar data.
%
%  Discussion:
%
%    This program assumes that you have computed the value of some scalar
%    quantity (such as pressure or temperature) at a set of nodes.
%
%    You may have determined an order 3 triangulation of these nodes,
%    but if you have not, the program will work that out internally.
%
%    This program can read that data, and display a color contour of the
%    solution data.
%
%  Usage:
%
%    fem2d_scalar_display_gpl ( 'prefix' )
%
%    where
%
%    * 'prefix'_nodes.txt contains the node coordinates;
%    * 'prefix'_elements.txt contains the element definitions
%      (this file is optional, and if missing, the elements will be generated
%      by the program);
%    * 'prefix'_values.txt contains the nodal values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 May 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string PREFIX, the common file prefix.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM2D_SCALAR_DISPLAY_GPL:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Surface plot of a scalar U(X,Y) on a triangulated region.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This program expects three input files:\n' );
  fprintf ( 1, '  * a node file,      the node coordinates,\n' );
  fprintf ( 1, '  * an element file,  triples of nodes that form elements,\n' );
  fprintf ( 1, '  * a value file,     solution values.\n' );
%
%  The command line argument is the common filename prefix.
%
  if ( nargin < 1 )

    fprintf ( 1, '\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, 'FEM2D_SCALAR_DISPLAY_GPL - Fatal error!\n' );
    fprintf ( 1, '  Dataset must have spatial dimension 2.\n' );
    error ( 'FEM2D_SCALAR_DISPLAY_GPL - 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 or create the element data.
%
  if ( file_exist ( element_filename ) )

    [ element_order, element_num ] = i4mat_header_read ( ...
      element_filename );

    if ( element_order ~= 3 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'FEM2D_SCALAR_DISPLAY_GPL - Fatal error!\n' );
      fprintf ( 1, '  Data is not for a 3-node triangulation.\n' );
      error ( 'FEM2D_SCALAR_DISPLAY_GPL - Fatal error!' );
    end

%   fprintf ( 1, '\n' );
%   fprintf ( 1, '  Read the header of "%s".\n', ...
%     element_filename );
%   fprintf ( 1, '\n' );
    fprintf ( 1, '  Element order = %d\n', element_order );
    fprintf ( 1, '  Number of elements ELEMENT_NUM  = %d\n', ...
      element_num );

    element_node = i4mat_data_read ( element_filename, ...
      element_order, element_num );

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

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Creating triangulation for data.\n' );
    element_node = delaunayn ( node_xy' );
    element_node = element_node';
    [ element_order, element_num ] = size ( element_node );
    i4mat_write ( element_filename, element_order, element_num, element_node );
    fprintf ( 1, '  Triangulation data written to "%s".\n', element_filename );
  end

% i4mat_transpose_print_some ( element_order, element_num, ...
%   element_node, 1, 1, element_order, 10, ...
%   '  First 10 elements:' );
%
%  Detect and correct 0-based indexing.
%
  element_node = mesh_base_one ( node_num, element_order, element_num, ...
    element_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, 'FEM2D_SCALAR_DISPLAY_GPL - Fatal error!\n' );
    fprintf ( 1, '  VALUE data must be scalar.\n' );
    error ( 'FEM2D_SCALAR_DISPLAY_GPL - 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:' );
%
%  Create the GNUPLOT files.
%
  gnuplot_write ( node_num, node_xy, element_num, element_node, ...
    value, prefix  );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM2D_SCALAR_DISPLAY_GPL:\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 value = file_exist ( file_name )

%*****************************************************************************80
%
%% FILE_EXIST reports whether a file exists.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character FILE_NAME, the name of the file.
%
%    Output, logical FILE_EXIST, is TRUE if the file exists.
%
  fid = fopen ( file_name );

  if ( fid == -1 )
    value = 0;
  else
    fclose ( fid );
    value = 1;
  end

  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 gnuplot_write ( node_num, node_xy, element_num, ...
  element_node, value, prefix )

%*****************************************************************************80
%
%% GNUPLOT_WRITE creates files for a finite element surface plot.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 May 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the coordinates of the nodes.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(3,ELEMENT_NUM), the element
%    definitions.
%
%    Input, integer VALUE(NODE_NUM), the values.
%
%    Input, string PREFIX, the common filename prefix.
%

%
%  Create the data file.
%
  data_filename = strcat ( prefix, '_data.txt' );
  data = fopen ( data_filename, 'wt' );

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

  for element = 1 : element_num
    for vertex = 1 : 3
      node = element_node(vertex,element);
      fprintf ( data, '%f  %f  %f\n', node_xy(1,node), node_xy(2,node), value(node) );
    end
    node = element_node(1,element);
    fprintf ( data, '%f  %f  %f\n', node_xy(1,node), node_xy(2,node), value(node) );
    fprintf ( data, '\n' );
    fprintf ( data, '\n' );
  end

  fclose ( data );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Created the data file "%s"\n', data_filename );
%
%  Create the command file.
%
  command_filename = strcat ( prefix, '_commands.txt' );
  commands = fopen ( command_filename, 'wt' );

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

  fprintf ( commands, 'set term png\n' );
  fprintf ( commands, 'set output "%s.png"\n', prefix );
  fprintf ( commands, 'set grid\n' );
  fprintf ( commands, 'set style data lines\n' );
  fprintf ( commands, 'set timestamp\n' );
  fprintf ( commands, 'unset key\n' );
  fprintf ( commands, 'set xlabel "<---X--->"\n' );
  fprintf ( commands, 'set ylabel "<---Y--->"\n' );
  fprintf ( commands, 'set zlabel "<---U(X,Y)--->"\n' );
  fprintf ( commands, 'set title "%s"\n', prefix );
  fprintf ( commands, 'splot "%s"\n', data_filename );

  fclose ( commands );

  fprintf ( 1, '  Created the command file "%s"\n', command_filename );

  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, a title.
%
  incx = 10;

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

  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.
%
%  Discussion:
%
%    An I4MAT is an array of I4's.
%
%  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 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


text:
classdef TestFastLineDetector
    %TestFastLineDetector

    methods (Static)
        function test_1
            img = cv.imread(fullfile(mexopencv.root(),'test','blox.jpg'), ...
                'Grayscale',true, 'ReduceScale',2);

            obj = cv.FastLineDetector();
            lines = obj.detect(img);
            validateattributes(lines, {'cell'}, {'vector'});
            cellfun(@(v) validateattributes(v, {'numeric'}, ...
                {'vector', 'numel',4, 'real'}), lines);

            img = repmat(img, [1 1 3]);
            out = obj.drawSegments(img, lines);
            validateattributes(out, {class(img)}, {'size',size(img)});
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [x C d] = project_affine(v, A, b, C, d)
% PROJECT_AFFINE    Project a point into an affine set.
%
%   project_affine(v,A,b) is the projection of v onto
%   the affine set { x | Ax = b }.
%
%   You can also call the function as
%
%     [x C d] = project_affine(v,A,b);
%
%   and then call it again with different argument v2 via
%
%     x2 = project_affine(v2,A,b,C,d);
%
%   If calling the function repeatedly with the same A and b,
%   all evaluations after the initial one will be much faster
%   when passing in the cached values C and d.

    if exist('C','var') && ~isempty(C) && exist('d','var') && ~isempty(d)
        x = C*v + d;
    else
        % x = v - pinv(A)*(A*v - b);
        pA = pinv(A);
        C = eye(length(v)) - pA*A;
        d = pA*b;
        x = C*v + d;
    end
end

author: cvxgrp

repo: proximal

save_path: github-repos/MATLAB/cvxgrp-proximal


text:
function newSteps = dtiUnwarpStep(origSteps,param)
% Unwarp (shift and scale) the steps into a new coord system
%
%    newSteps = dtiUnwarpStep(origSteps,param)
%
%    Shifts and scales according to: f(x,p) = ( x/p(1) ) - p(2)
%
%    origSteps - nx1 vector
%    param     - 2x1 vector [scale shift]
%
%    newSteps  - nx1 vector
%
% Example:
%    % original x's
%    x = 1:10;
%    % stretch by 2 times after shifting 5 units to the right
%    nx = dtiWarpStep(x,[2 5]);
%    % go back to the original
%    ox = dtiUnwarpStep(nx,[2 5]);
%
% See also: dtiWarpStep
%
% History:
%    2007/01/17 shc wrote it.
%

if ieNotDefined('origSteps'), error('Require input data!');         end
if ieNotDefined('param'),     error('Require warping parameters!'); end

newSteps = (origSteps/param(1))-param(2);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% A file to demonstrate the power_bounded function
% Author: Muhammad F. Kasim (University of Oxford)

% test == 1 for drawing an ordinary voronoi diagram
% test == 2 for power diagram with random small weights
% test == 3 to see how long does it take for larger number of sites with random small weights
test = 2;

if test < 3
    N = 50; % number of sites
    x = rand(N, 1);
    y = rand(N, 1); % x and y coordinate of the sites
    if test == 1
        wts = zeros(N, 1); % weights of each coordinate (all-zeros weights produce voronoi diagram)
    elseif test == 2
        wts = rand(N, 1) * 0.001;
    end
    crs = [0,0; 0,1; 1,1; 1,0]; % the bounding box in clockwise order

    % get the power diagram
    [V,C] = power_bounded(x,y, wts, crs);

    % draw the resulted power diagram
    plot(x, y, 'r.'); % plot the sites' points
    hold on;
    for i = [1:N]
        xpoly = V(C{i},1);
        ypoly = V(C{i},2);
        xpoly = [xpoly; xpoly(1)];
        ypoly = [ypoly; ypoly(1)];
        plot(xpoly, ypoly, 'b-');
    end
    xlim([min(crs(:,1)), max(crs(:,1))]);
    ylim([min(crs(:,2)), max(crs(:,2))]);
    hold off;

elseif test == 3
    Ns = [100, 1000, 10000, 100000, 1000000];
    times = [];
    for N = Ns
        x = rand(N, 1);
        y = rand(N, 1);
        wts = rand(N, 1) * 0.0001;
        crs = [0,0; 0,1; 1,1; 1,0];

        % time the power diagram
        tic;
        power_bounded(x, y, wts, crs);
        times(end+1) = toc;
    end
    disp(times);

    loglog(Ns, times, 'bo', 'LineWidth', 2);
    hold on;
    loglog(Ns, times, 'b-', 'LineWidth', 2);
    set(gca, 'FontSize', 14)
    xlabel('Number of sites');
    ylabel('Time required (s)');
end

author: optimaltransport


text:
function KC = slcenkernel(K0, K, w)
%SLCENKERNEL Compute the centralized kernel matrix
%
% $ Syntax $
%   - KC = slcenkernel(K0)
%   - KC = slcenkernel(K0, [], w)
%   - KC = slcenkernel(K0, K)
%   - KC = slcenkernel(K0, K, w)
%
% $ Arguments $
%   - K0:       the gram matrix of the referenced samples
%   - K:        the kernel matrix for target samples
%   - w:        the weights for the referenced samples
%   - KC:       the centralized kernel matrix.
%
% $ Description $
%   - KC = slcenkernel(K0) compute the centralized kernel matrix from
%     the original kernel gram matrix K0.
%
%   - KC = slcenkernel(K, [], w) compute the centralized kernel gram
%     matrix from the original kernel gram matrix K0. The mean feature
%     is obtained with the weights for referenced samples, given by w.
%
%   - KC = slcenkernel(K0, K) compute the centralized kernel matrix
%     for target samples, with the original gram matrix for referenced
%     samples K0 and the kernel matrix for target samples w.r.t the
%     referenced samples K given.
%
%   - KC = slcenkernel(K0, K, w) compute the centralied kernel matrix
%     for target samples, with feature mean computed in a weighted
%     manner.
%
% $ Remarks $
%   -# For original kernel matrix K, it is defined as
%      K(i, j) = <phi(i), phi(j)>, given that phi(i) and phi(j) are
%      the feature map of the samples x(i) and x(j) respectively.
%      Then the centralized kernel matrix is defined as
%      KC(i, j) = <phi(i) - mean_phi, phi(j) - mean_phi>, where
%      mean_phi is the mean of all referenced feature maps. If w
%      is specified mean_phi is given by weighted mean.
%      Kernel centralization plays an important role in many
%      kernelized algorithms such as Kernel PCA.
%
%   -# Suppose the mean feature map is defined by
%      mean_phi = sum_i w_i phi(i).
%      Then the centralized kernel gram matrix can be written as
%      KC = K - 1 * (w^T * K) - (K * w) * 1^T + 1 * (w^T * K * w) * 1^T.
%      It can be easily shown that the mean of columns (rows) of KC is
%      a zero vector. Thus, centralize the centralized kernel matrix
%      would keep the input unchanged.
%
%   -# Instead of applying the formula given above, the function
%      implements a more efficient computational routine by reducing
%      the redundant computations.
%
% $ History $
%   - Created by Dahua Lin on May 2nd, 2006
%   - Modified by Dahua Lin on Sep 10, 2006
%       - use sladdrowcols to replace sladd to increase efficiency
%

%% parse and verify input arguments

% for K0
n0 = size(K0, 1);
if ndims(K0) ~= 2 || size(K0, 2) ~= n0;
    error('sltoolbox:invaliddims', ...
        'K0 should be a 2D square matrix');
end

% for K
if nargin < 2 || isempty(K)
    K = K0;
else
    if ndims(K) ~= 2
        error('sltoolbox:invaliddims', ...
            'K should be a 2D matrix');
    end
    if size(K, 1) ~= n0
        error('sltoolbox:sizmismatch', ...
            'Size inconsistency between K0 and K');
    end
end

% for w
if nargin < 3 || isempty(w)
    isweighted = false;
else
    if ~isequal(size(w), [1, n0])
        error('sltoolbox:sizmismatch', ...
            'Size inconsistency between K0 and w');
    end
    isweighted = true;
end

%% compute

% Steps:
% 1. compute v1: mean row vector of K (1 x n)
% 2. compute v2: mean column vector of K0 (n x 1)
% 3. compute s3: mean value of of all elements of K0 (1 x 1)
% 4. KC = K - expand(v1) - expand(v2) + s3

if ~isweighted  % non-weighted case

    v1 = sum(K, 1) * (1 / n0);
    v2 = sum(K0, 2) * (1 / n0);
    s3 = sum(v2) * (1 / n0);

else            % weighted case

    w = w / sum(w);     % normalize the weights

    v1 = w * K;
    v2 = K0 * w';
    s3 = w * v2;

end

KC = sladdrowcols(K, -v1, -v2+s3);

author: lmthang


text:
function test_suite = test_edgePosition
%TEST_EDGEPOSITION  Test case for the file edgePosition
%
%   Test case for the file edgePosition

%   Example
%   test_edgePosition
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@nnates.inra.fr
% Created: 2016-07-19,    using Matlab 8.6.0.267246 (R2015b)
% Copyright 2016 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function test_Simple(testCase) %#ok<*DEFNU>
% Test call of function without argument

edge = [10 10 30 20];
p1 = [10 10];
p2 = [20 15];
p3 = [30 20];

pos1 = edgePosition(p1, edge);
pos2 = edgePosition(p2, edge);
pos3 = edgePosition(p3, edge);

testCase.assertEqual(0, pos1, 'AbsTol', .01);
testCase.assertEqual(.5, pos2, 'AbsTol', .01);
testCase.assertEqual(1, pos3, 'AbsTol', .01);

function test_MultipleEdges(testCase)

edges = [10 10 30 10 ; 10 10 0 0 ; 0 0 20 20];
point = [10 10];

posList = edgePosition(point, edges);

testCase.assertEqual([0 0 .5], posList, 'AbsTol', .01);

function test_MultiplePoints(testCase)

edge = [10 10 30 20];
points = [10 10; 20 15; 30 20; 15 12.5];

posList = edgePosition(points, edge);

testCase.assertEqual([0; .5; 1; .25], posList, 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function mtrSetupDOCCParams(subjID)

% Setting up parameters files for a subject based on ROIs assumes
% 1. You are in the subject directory
% 2. A met_params.txt file is already in the conTrack directory
%subjID = 'ss040804';
localSubjDir = ['/teal/scr1/dti/sisr/' subjID];
dlhSubjDir = ['/radlab_share/home/tony/images/' subjID];
lappySubjDir = ['c:/cygwin/home/sherbond/images/' subjID];
machineListL = {'deep-lumen1','deep-lumen2','deep-lumen4'};
machineListR = {'deep-lumen5','deep-lumen6','deep-lumen8'};
mtrTwoRoiSamplerSISR(localSubjDir, 'CC_FA.mat', 'LV3AB7d_cleaned.mat', 'met_params.txt', 'LDOCC_met_params.txt', 'paths_LDOCC.dat', 'ldoccRoisMask.nii.gz');
mtrTwoRoiSamplerSISR(localSubjDir, 'CC_FA.mat', 'RV3AB7d_cleaned.mat', 'met_params.txt', 'RDOCC_met_params.txt', 'paths_RDOCC.dat', 'rdoccRoisMask.nii.gz');
mtrTwoRoiSamplerSISR(localSubjDir, 'CC_FA.mat', 'LV3AB7d_cleaned.mat', 'met_params.txt', 'LDOCC_met_params_dl.txt', 'paths_LDOCC.dat', 'ldoccRoisMask.nii.gz', 'track_LDOCC.sh', dlhSubjDir, machineListL);
mtrTwoRoiSamplerSISR(localSubjDir, 'CC_FA.mat', 'RV3AB7d_cleaned.mat', 'met_params.txt', 'RDOCC_met_params_dl.txt', 'paths_RDOCC.dat', 'rdoccRoisMask.nii.gz', 'track_RDOCC.sh', dlhSubjDir, machineListR);
mtrTwoRoiSamplerSISR(localSubjDir, 'CC_FA.mat', 'LV3AB7d_cleaned.mat', 'met_params.txt', 'LDOCC_met_params_lap.txt', 'paths_LDOCC.dat', 'ldoccRoisMask.nii.gz', 'temp.sh', lappySubjDir, machineListL);
mtrTwoRoiSamplerSISR(localSubjDir, 'CC_FA.mat', 'RV3AB7d_cleaned.mat', 'met_params.txt', 'RDOCC_met_params_lap.txt', 'paths_RDOCC.dat', 'rdoccRoisMask.nii.gz', 'temp.sh', lappySubjDir, machineListR);

delete(fullfile('conTrack','temp.sh'));


text:
function options = gpsimMapInitParam(y, yvar, options);

% GPSIMMAPINITPARAM Creates a set of options for a GPSIMMAP model as the initial parameters.
% FORMAT
% DESC returns a default initialisation for a GPSIMMAP model.
% ARG y :  the values of each gene at the different time points.
% ARG yvar :  the variances of each gene at the different time points.
% ARG options : the values of each gene at the different time points.
% RETURN options : the options structure.
%
% SEEALSO : gpsimMapCreate
%
% COPYRIGHT : Pei Gao, Neil D. Lawrence and Magnus Rattary, 2008

% GPSIM

numGenes = length(options.S);
options.S = ones(1, numGenes);
initOptions = options;
mu = mean(y, 1);
optionDset = [0.01 0.1 0.4 1];
rand('seed',1);

for j = 1:numGenes
  initOptions.S = options.S(j);
  initOptions.gParam = options.gParam(j);
  %    figure(j);
  %    plot(options.times, y(:,j), 'rx');
  %    hold on;
  for di = 1:length(optionDset)
    initOptions.D = optionDset(di);
    initOptions.B = initOptions.D*mu(j);

    model.gene{j}.init{di} = gpsimMapCreate(1, 1, options.times, y(:,j), ...
                                            yvar(:,j), initOptions);
    if strcmp(options.kern, 'mlp')
      model.gene{j}.init{di}.kern.weightVariance = 30;
      model.gene{j}.init{di}.kern.biasVariance = 20;
      % This forces kernel recompute.
      params = gpsimMapExtractParam(model.gene{j}.init{di});
      model.gene{j}.init{di} = ...
          gpsimMapExpandParam(model.gene{j}.init{di}, params);
    elseif strcmp(options.kern, 'rbf')
      model.gene{j}.init{di}.kern.inverseWidth = 0.001;
      params = gpsimMapExtractParam(model.gene{j}.init{di});
      model.gene{j}.init{di} = ...
          gpsimMapExpandParam(model.gene{j}.init{di}, params);
    end

    updateFoptions = defaultOptions;
    model.gene{j}.init{di} = gpsimMapUpdateF(model.gene{j}.init{di}, ...
                                             updateFoptions);
    %       fmean =  mean(model.gene{j}.init{di}.f);
    %       fsd = sqrt(var(model.gene{j}.init{di}.f));
    %       model.gene{j}.init{di}.f = model.gene{j}.init{di}.f - fmean;
    %       model.gene{j}.init{di}.f = model.gene{j}.init{di}.f/fsd;
    %       model.gene{j}.init{di}.S = model.gene{j}.init{di}.S/exp(fmean);
    %       model.gene{j}.init{di}.gParam = model.gene{j}.init{di}.gParam/exp(fmean);
    ypredMean = mean(model.gene{j}.init{di}.ypred);
    ypredScale = sqrt(var(model.gene{j}.init{di}.ypred(model.gene{j}.init{di}.times_index)));
    model.gene{j}.init{di}.B = (model.gene{j}.init{di}.B - initOptions.D* ...
                                (ypredMean - ypredScale*mean(model.gene{j}.init{di}.y)))/ ...
        ypredScale;
    %       if isfield(model.gene{j}.init{di}, 'alpha')
    %         model.gene{j}.init{di}.alpha = model.gene{j}.init{di}.alpha/ypredScale;
    %       end
    if model.gene{j}.init{di}.B < 0 && isfield(model.gene{j}.init{di}, ...
                                               'alpha')
      model.gene{j}.init{di}.alpha = model.gene{j}.init{di}.alpha/ypredScale;
      model.gene{j}.init{di}.alpha = model.gene{j}.init{di}.alpha + ...
          model.gene{j}.init{di}.B;
      model.gene{j}.init{di}.B = 1e-6;
    end
    model.gene{j}.init{di}.S = model.gene{j}.init{di}.S/ypredScale;
    params = gpsimMapExtractParam(model.gene{j}.init{di});
    model.gene{j}.init{di} = gpsimMapExpandParam(model.gene{j}.init{di}, ...
                                                 params);
    f = gpsimMapFunctionalExtractParam(model.gene{j}.init{di});
    model.gene{j}.init{di} = gpsimMapFunctionalExpandParam(model.gene{j}.init{di}, f);
    model.gene{j}.ll(di) = ...
        gpsimMapLogLikelihood(model.gene{j}.init{di});
    %      plot(model.gene{j}.init{di}.mapt, model.gene{j}.init{di}.ypred);
  end
  [value, index] = max(model.gene{j}.ll);
  %    plot(model.gene{j}.init{index}.mapt, model.gene{j}.init{index}.ypred, 'm');

  options.D(j) = optionDset(index);
  options.S(j) = model.gene{j}.init{index}.S;
  %     if model.gene{j}.init{index}.B < 0
  %       if isfield(model.gene{j}.init{index},'alpha')
  %         options.alpha(j) = model.gene{j}.init{index}.alpha + ...
  %             model.gene{j}.init{di}.B;
  %       end
  %       options.B(j) = 1e-6;
  %     else
  if isfield(model.gene{j}.init{index},'alpha')
    options.alpha(j) = model.gene{j}.init{index}.alpha;
  end
  %     end
  options.B(j) = model.gene{j}.init{index}.B;
  options.gParam(j) = model.gene{j}.init{index}.gParam;
end
options.S = rand(1,numGenes);
options.B = options.D.*mu;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ rel_err, snr ] = get_SNR( X, X0 )

noise = X - X0;
trueNorm = norm(X0);
rel_err = norm( noise ) / trueNorm;
snr = -20*log10( rel_err );

end

author: andrewssobral

repo: lrslibrary


text:
function [Lden,Ld] = dpr1fact(x, d, Lsym, smult, maxu) %#ok
% [Lden,L.d] = dpr1fact(x, d, Lsym, smult, maxu)
%
% DPR1FACT  Factor d[iag] p[lus] r[ank] 1:
%    [Lden,L.d] = dpr1fact(x, d, Lsym, smult, maxu)
%    Computes fi and d such that
%       diag(d_IN) + x*diag(smult)*x' =
%(PI_{i=1}^n L(p_OUT^i,beta_i)) * diag(d_OUT) * (PI_{i=1}^n L(p_OUT^i,beta_i))'
%    where L(p,beta) = eye(n) + tril(p*beta',-1).
%
% Lden.dopiv(k) = 1 if p(:,k) has been reordered, with permutation in
% Lden.pivperm.
% We reorder if otherwise |p(i,k)*beta(j,k)| > maxu.
%
% ******************** INTERNAL FUNCTION OF SEDUMI ********************
%
% See also fwdpr1,bwdpr1,sedumi

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% 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 2 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, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

%Indicate to the user Matlab cannot find the SeDuMi binaries
sedumi_binary_error();

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function stats = rsquare_multiple_regions_multilevel(Y, X, varargin)
% Predict outcome var (y) using data from multiple vars (X, e.g., brain regions)
% Test R-square against permuted data
%
% :Usage:
% ::
%
%     stats = rsquare_multiple_regions_multilevel(Y, X, varargin)
%
% :Inputs: Variable args
%
%   **'colors':**
%        followed by colors cell ({'r' 'g' 'b'}) for each col. of X
%
%   **'nperms:**
%        then num perms
%
% :Examples:
% ::
%
%    % SETUP data
%    cl = [clpos_data clneg_data];
%
%    cl(cat(1, cl.numVox) < min_cluster_size) = [];
%
%    % get brain data cell
%    for i = 1:size(cl(1).all_data, 2)
%        for c = 1:length(cl)
%            data{i}(:,c) = cl(c).all_data(:, i);
%        end
%    end
%
%    % NMDS ANALYSIS
%    OUT = [];
%
%    OUT.ridge = matrix_direct_effects_ridge(data);
%    D = OUT.ridge.mean; D(find(eye(size(D)))) = 1;
%    D = (D' + D) ./ 2;
%    OUT.ridge.D = (1 - D) ./ 2;
%    [OUT.stats_mds.GroupSpace,OUT.stats_mds.obs,OUT.stats_mds.implied_dissim] = shepardplot(OUT.ridge.D,[]);
%
%    OUT.stats_mds = nmdsfig_tools('cluster_solution',OUT.stats_mds, OUT.stats_mds.GroupSpace, 2:max_networks, nperms, []);
%    OUT.stats_mds.colors = {'ro' 'go' 'bo' 'yo' 'co' 'mo' 'ko' 'r^' 'g^' 'b^' 'y^' 'c^' 'm^' 'k^'};
%    create_figure('nmdsfig');
%
%    OUT.stats_mds.names = [];
%    nmdsfig(OUT.stats_mds.GroupSpace,'classes',OUT.stats_mds.ClusterSolution.classes,'names',OUT.stats_mds.names,'sig',OUT.ridge.fdrsig);
%    hh = nmdsfig_fill(OUT.stats_mds);
%    axis image; axis equal
%
%
%    % Multiple regions predict behavior
%
%    % Design matrix with cluster averages
%    classes = OUT.stats_mds.ClusterSolution.classes;
%    clear X
%    for i = 1:length(data)
%        for j = 1:max(classes)
%            X{i}(:, j) = nanmean(data{i}(:, classes == j), 2);
%        end
%    end
%
%    OUT.stats_regression = rsquare_multiple_regions_multilevel(acc_x_dist, X, 'colors', OUT.stats_mds.colors, 'nperms', 100);
%
% ..
%    Tor Wager, June 2008
%
%    var args not all done yet (in development)
% ..

nperms = 100;

for i = 1:length(varargin)
        if ischar(varargin{i})
            switch varargin{i}

                    % functional commands
                case 'colors', colors = varargin{i+1};

                case 'nperms', nperms = varargin{i + 1};

                otherwise, warning(['Unknown input string option:' varargin{i}]);
            end
        end
end

%% initial model fits, all together, then separately to get order
stats = glmfit_multilevel(Y, X, ones(length(X), 1), 'weighted');

nvars = size(X{1}, 2);

for i = 1:nvars
    for s = 1:length(X)
        Xr{s} = [X{s}(:, i)]; % reduced model, with only this regressor.
    end

    stats_ind = glmfit_multilevel(Y, Xr, ones(length(X), 1), 'weighted');
    stats.separate_reg_pvals(i) = stats_ind.p(2);

end

%%% **** Notes: should do best 1, then best combo of 2, best 3...etc ****
%%% permute additional param after best combo from prev. step

%% Predicting using variable # regions
[tmp, myorder] = sort(stats.separate_reg_pvals); % eliminate intercept
stats.myorder = myorder;

for i = 1:length(myorder)

    for s = 1:length(X)

        Xr{s} = [ones(size(X{s}, 1), 1) X{s}(:, myorder(1:i))]; % reduced model, with regressors ordered by predictive accuracy.

        [B,BINT,R,RINT,sstats] = regress(Y{s}, Xr{s});

        r2(s, i) = sstats(1);

    end

end

create_figure('R-squared'); barplot_columns(r2, 'R-squared', [], 'nofig', 'noind');
xlabel('Number of networks in model');
ylabel('Average R^2 value');

if exist('colors', 'var') && ~isempty(colors)

    for i = 1:size(r2, 2)
        hbar(i) = bar(i, nanmean(r2(:, i)), 'FaceColor', colors{myorder(i)}(1));
    end

end

stats.observed_r2 = r2;

drawnow

%% permutation

clear permindx
for i = 1:length(X)
    permindx{i} = permute_setupperms(size(X{i}, 1), nperms);

    %Xrp{s} = [ones(size(X{s}, 1), 1) X{s}(:, myorder(1:i))]; % reduced model, with regressors ordered by predictive accuracy.

end

[r2meanp, b1mean] = deal(zeros(nperms, length(myorder)));

%%
for p = 1:nperms

    fprintf('%3.0f ', p);

    if mod(p, 20) == 0, fprintf(1, '\n'); end

    [r2, bb] = deal(zeros(length(X), length(myorder)));

    for i = 1:length(myorder)

        for s = 1:length(X)

            % Method 1: Permute all variables (permute data)
            Xr{s} = [ones(size(X{s}, 1), 1) X{s}(:, myorder(1:i))]; % reduced model, with regressors ordered by predictive accuracy.

            Yp{s} = Y{s}(permindx{s}(p, :)', :); % permute

% %             % Method 2: Permute only the i:endth regressors
% %             % part to not permute: conditional on this being in model
% %             Xsnp = X{s}(:, myorder(1:i-1));
% %
% %             % part to permute: if everything else were random...
% %             Xsp = X{s}(:, myorder(i:end));
% %             Xsp = Xsp(permindx{s}(p, :)', :); % permute
% %             Xr{s} = [ones(size(X{s}, 1), 1) Xsnp Xsp]; % reduced model, with regressors ordered by predictive accuracy.
% %             Yp{s} = Y{s};

            [B,BINT,R,RINT,sstats] = regress(Yp{s}, Xr{s});

            r2(s, i) = sstats(1);

            bb(s, i) = B(2);

            % %         Xrp_ml{s} = [X{s}(permindx{s}(p, :)', myorder(1:i))]; % reduced model, with regressors ordered by predictive accuracy.

        end

        % %     stats = glmfit_multilevel(Y, Xrp_ml, ones(length(X), 1), 'weighted', 'noverbose');
        % %     multilev_pvals{i}(p, :) = stats.p;

    end

    r2meanp(p, :) = nanmean(r2);

    b1mean(p, :) = nanmean(bb);

end

fprintf(1, '\n');

bar((1:nvars) + .2, mean(r2meanp), 'FaceColor', [.2 .2 .2]);

%% Z-scores and p-values based on perm dist
fp = r2meanp > repmat(nanmean(stats.observed_r2, 1), nperms, 1);

stats.permute.pvals = sum(double(fp), 1) ./ nperms;

yval = nanmean(stats.observed_r2) + .15 * nanmean(stats.observed_r2);

wh = find(stats.permute.pvals <= .001);
if ~isempty(wh), for i = 1:length(wh), text(wh(i), yval(wh(i)), '***', 'FontSize', 18); end, end

wh = find(stats.permute.pvals <= .01 & stats.permute.pvals > .001);
if ~isempty(wh), for i = 1:length(wh), text(wh(i), yval(wh(i)), '**', 'FontSize', 18); end, end

wh = find(stats.permute.pvals <= .05 & stats.permute.pvals > .01);
if ~isempty(wh), for i = 1:length(wh), text(wh(i), yval(wh(i)), '*', 'FontSize', 18); end, end

%%
disp('Added .permute field to stats output strucuture.');
stats.permute.nperms = nperms;
stats.permute.r2meanp = r2meanp;
stats.permute.r2mean_ub = prctile(stats.permute.r2meanp, 95);
stats.permute.permindx = permindx;

stats.permute.b1mean = b1mean;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef UnfittedMesh_Builder_Factory < handle

    methods (Access = public, Static)

        function concreteBuilder = create(type,ndim)
            switch ndim
                case 1
                    concreteBuilder = UnfittedMesh_StraightLine();
                case 2
                    switch type
                        case 'INTERIOR'
                            concreteBuilder = UnfittedMesh_FlatSurface();
                        case 'BOUNDARY'
                            concreteBuilder = UnfittedMesh_FlatCurve();
                    end
                case 3
                    switch type
                        case 'INTERIOR'
                            concreteBuilder = UnfittedMesh_Volumetric();
                        case 'BOUNDARY'
                            concreteBuilder = UnfittedMesh_3DSurface();
                    end
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function err = relaxT1FitFunc(x, m, alpha, tr, te)
%
% function err = relaxT1FitFunc(x, m, alpha, tr, te);
% inputs:
%       x = [r1, pd] where r1 is the inverse of the longitudinal relaxation
%       time (1/T1) and pd is a compsite term representing
%       M_0*G*exp(-TE/T_2star)
%       m: spgr data
%
%       alpha: flip angle (in degrees)
%       tr, te = tr and te
%
%   Output: the error for each measurement
%
%

r1 = x(1);
pd = x(2);

%pd = M_0*G*exp(-te/T_2star);
S_alpha = pd .* sin(alpha) .* (1-exp(-r1.*tr)) ./ (1-cos(alpha) .* exp(-r1.*tr));
err = S_alpha - m;

return;


text:
function C = mtimes(A,B)
%MTIMES Implement A*B (scalar multiply) for ktensor.
%
%   C = mtimes(A,B) computes A * B where A is a Kruskal tensor and B is
%   a scalar (or vice versa). The result C is the same size as A.
%
%   See also KTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

% Note: We can do scalar times a tensor, but anything more complex is
% an error.

if isa(B,'numeric') && isequal(size(B),[1 1])
    C = ktensor(B * A.lambda, A.u);
elseif isa(A,'numeric') && isequal(size(A),[1 1])
    C = ktensor(A * B.lambda, B.u);
else
    error('Use mtimes(full(A),full(B)).');
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function box = minandmax2est(f, N)
%MINANDMAX2EST   Estimates the range of a DISKFUNV.
%   BOX = MINANDMAX2EST(F) returns estimates for the minimum and maximum of each
%   component of the DISKFUNV F over its domain.  BOX is a vector of length
%   four, containing the estimated minimum and maximum of each component.
%
%   BOX = MINANDMAX2EST(F, N) returns estimates for the minimum and maximum of
%   each component of the DISKFUNV F over its domain, based on samples on
%   an N by N grid (N = 33 by default).
%
% See also DISKFUN/MINANDMAX2EST.

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

box = [];

if ( isempty(f) )
    return
end

if ( ( nargin < 2 ) || isempty(N) )
    % Default to N = 33:
    N = 33;
end

for jj = 1:2
    box = [ box, minandmax2est(f.components{jj}, N) ];
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% Synthesis of a stationary time signal given sinusoids parameters
%
% Octave compatible
%
% Inputs
%  sins   : matrix of size(3,N) [freq;amp;phase]
%  fs     : [Hz] Sampling frequency
%  winlen : length of the time signal to synthesize
%  derive : if deriv=1, synthesize the derivative of the signal
%
% Outputs
%  s       : The synthesized signal
%
% Copyright (c) 2012 University of Crete - Computer Science Department
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%

function [s t] = sin2sig(sins, fs, winlen, deriv)
    if nargin<4; deriv=0; end

%      t = (0:winlen-1)/fs;
%      t = (0:winlen-1)/fs;
    % The time reference used for the phase in sins is in the window center
    t = (-(winlen-1)/2:(winlen-1)/2)/fs;

    if deriv==0
        A = cos(2*pi*t'*sins(1,:) + ones(length(t),1)*sins(3,:));
        A = ones(length(t),1)*sins(2,:).*A;
        s = sum(A,2);
    elseif deriv==1
        A = -sin(2*pi*t'*sins(1,:) + ones(length(t),1)*sins(3,:));
        A = ones(length(t),1)*(sins(2,:).*(2*pi*partials(1,:))).*A;
        s = sum(A,2);
    end

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function [a,fps,tps] = aucs(crit,z)
%AUCS Compute area under curve for survival model
%
%  Description
%    A = AUCS(CRIT,Z) Compute are under curve for survival model
%    using criteria vector CRIT (where larger value means larger
%    risk of incidence) and censoring indicator column vector Z
%    (0=event, 1=censored).
%
%    [A,FPS,TPS] = AUCS(CRIT,Z) Return also a vector of false
%    positive rates FPS and a vector of true positive rates TPS.
%
%  Reference
%    L. E. Chambless, C. P. Cummiskey, and G. Cui (2011). Several
%    methods to assess improvement in risk prediction models:
%    Extension to survival analysis. Statistics in Medicine
%    30(1):22-38.

% Copyright (C) 2012 Ernesto Ulloa, Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

ip=inputParser;
ip.addRequired('crit',@(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('z', @(x) isreal(x) && all(isfinite(x(:))))
ip.parse(crit,z)

ye=z;
tps=0;
fps=0;
qs=sort(crit,'descend');

for i2=1:numel(crit)
  tps(i2+1)=mean(crit>=qs(i2) & 1-ye)./mean(1-ye);
  fps(i2+1)=mean(crit>=qs(i2) & ye)./(mean(ye));
end

if (mean(1-ye)==0)||(mean(ye)==0)
  warning('z vector has no different values, function will return 0');
  a=0;
else
  a=sum([diff(fps).*mean([tps(1:end-1);tps(2:end)],1)]);
end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function c = typeadj(a,TYPE)
%TYPEADJ      typecast of  a  to type TYPE
%
%   c = typeadj(a,TYPE)
%
%For details, see intval\typeof and intval\typeadj.
%

% written  12/06/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

% slope\@slope\typeadj:  a  must be slope
  switch TYPE
    case 'slope',          c = a;
  otherwise
    error('invalid type in call of typeadj')
  end

repo: OpenMAS


text:
function [ev,ee,ebound,xyp] = q1p0grid(x,y,xy,mv,bound,mbound);
%q1p0grid   Q1-P0 element grid generator
%   [ev,ee,ebound,xyp] = q1p0grid(x,y,xy,mv,bound,mbound);
%   input
%          x          x coordinate vector
%          y          y coordinate vector
%          xy         nodal coordinate vector
%          mv         Q2 macroelement mapping matrix
%          bound      boundary vertex vector
%          mbound     macroelement boundary vertex vector
%   output
%          ev         element vertex matrix
%          ee         element edge connection matrix
%          ebound     element boundary edge matrix
%          xyp        vertex coordinate vector
%
%   IFISS function: DJS; 28 February 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
xx=xy(:,1); yy=xy(:,2); nvtx=length(xx);
adj=sparse(nvtx,nvtx);
mel=length(mv(:,1)); nel=4*mel;
ev=zeros(nel,4);
%
%% loop over macroelements
k=1:mel;
% first element
ke=4*k-3;
ev(ke,1)=mv(k,1);
ev(ke,2)=mv(k,5);
ev(ke,3)=mv(k,9);
ev(ke,4)=mv(k,8);
% second element
ke=4*k-2;
ev(ke,1)=mv(k,5);
ev(ke,2)=mv(k,2);
ev(ke,3)=mv(k,6);
ev(ke,4)=mv(k,9);
% third element
ke=4*k-1;
ev(ke,1)=mv(k,9);
ev(ke,2)=mv(k,6);
ev(ke,3)=mv(k,3);
ev(ke,4)=mv(k,7);
% fourth element
ke=4*k;
ev(ke,1)=mv(k,8);
ev(ke,2)=mv(k,9);
ev(ke,3)=mv(k,7);
ev(ke,4)=mv(k,4);
%
%% define element edges
ect=1;
% bottom boundary edges
k1=find(mbound(:,2)==1)';
for k=mbound(k1)
   ebound(ect,1)=4*k-3; ebound(ect+1,1)=4*k-2;
   ebound(ect,2)=1    ; ebound(ect+1,2)=1;
   ect=ect+2;
end
% right boundary edges
k2=find(mbound(:,2)==2)';
for k=mbound(k2)
   ebound(ect,1)=4*k-2; ebound(ect+1,1)=4*k-1;
   ebound(ect,2)=2    ; ebound(ect+1,2)=2;
   ect=ect+2;
end
% top boundary edges
k3=find(mbound(:,2)==3)';
for k=mbound(k3)
   ebound(ect,1)=4*k-1; ebound(ect+1,1)=4*k;
   ebound(ect,2)=3    ; ebound(ect+1,2)=3;
   ect=ect+2;
end
% left boundary edges
k4=find(mbound(:,2)==4)';
for k=mbound(k4)
   ebound(ect,1)=4*k; ebound(ect+1,1)=4*k-3;
   ebound(ect,2)=4    ; ebound(ect+1,2)=4;
   ect=ect+2;
end
%%
% centroid coordinates
for ielem=1:nel
xc(ielem)=mean(xx(ev(ielem,:))); yc(ielem)=mean(yy(ev(ielem,:)));
end
xyp=[xc',yc'];
%
%% compute edge to edge connection array ee
      np=nel;
% initialise global matrices
      adj = sparse(nvtx,nvtx);
      ee = zeros(nel,4);
%
% evaluate element number on each edge in turn
% and assemble into adjacency matrix
%% nx= 0, ny=-1
		 adj=adj + sparse(ev(:,1),ev(:,2),1:np,nvtx,nvtx);
%% nx= 1, ny= 0
		 adj=adj + sparse(ev(:,2),ev(:,3),1:np,nvtx,nvtx);
%% nx= 0, ny= 1
		 adj=adj + sparse(ev(:,3),ev(:,4),1:np,nvtx,nvtx);
%% nx=-1, ny= 0
		 adj=adj + sparse(ev(:,4),ev(:,1),1:np,nvtx,nvtx);
%
       for el=1:nel
		   [ii,jj]=find(adj==el);
           ee(el,:)=diag(adj(jj,ii))';
		   end
           ee=ee(:,[2,4,3,1]);
%
% plotting of the grid
%
%if mel <=64,
	adj=sparse(nvtx,nvtx);
    for i=1:nel
	adj(ev(i,1),ev(i,2)) =1;
	adj(ev(i,2),ev(i,3)) =1;
	adj(ev(i,3),ev(i,4)) =1;
	adj(ev(i,4),ev(i,1)) =1;
    end
    figure(30)
    gplot(adj,xy,'b')
    axis('square')
    hold on
    adj=sparse(nvtx,nvtx);
    k1=find(ebound(:,2)==1);
    for k=1:length(k1)
    kk=ebound(k1(k));
    adj(ev(kk,1),ev(kk,2))=1;
    end
    k2=find(ebound(:,2)==2);
    for k=1:length(k2)
    kk=ebound(k2(k));
    adj(ev(kk,2),ev(kk,3))=1;
    end
    k3=find(ebound(:,2)==3);
    for k=1:length(k3)
    kk=ebound(k3(k));
    adj(ev(kk,3),ev(kk,4))=1;
    end
    k4=find(ebound(:,2)==4);
    for k=1:length(k4)
    kk=ebound(k4(k));
    adj(ev(kk,4),ev(kk,1))=1;
    end
%   gplot(adj,xy,'r')
%   axis('off')
plot(xy(:,1),xy(:,2),'ro')
xybd=xy(bound,:);
plot(xybd(:,1),xybd(:,2),'ko')
plot(xyp(:,1),xyp(:,2),'k*')
hold off
drawnow
title('Q1-P0 finite element subdivision')
%end
return

author: johannesgerer

repo: jburkardt-m


text:
function srad = solarradiation(dem,lat,cs,r)
% PUPROSE: Calculate solar radiation for a digital elevation model (DEM)
%          over one year for clear sky conditions in W/m2
% -------------------------------------------------------------------
% USAGE: srad = solarradiation(dem,lat,cs)
% where: dem is the DEM to calculate hillshade for
%        lat is the latitude vector for the DEM - same size as size(dem,2)
%        cs is the cellsize in meters
%        r is the ground reflectance (global value or map, default is 0.2)
%
%       srad is the solar radiation in W/m2 over one year per grid cell
%
% EXAMPLE:
%       srad = solarradiation(peaks(50)*100,54.9:-0.1:50,1000,0.2);
%       - calculates the solar radiation for an example 50x50 peak surface.
%
% See also: GRADIENT, CART2POL
%
% Note: Follows the approach of Kumar et al 1997. Calculates clear sky
%       radiation corrected for the incident angle (selfshading) plus
%       diffuse and reflected radiation. Insolation is depending on time of year (and day),
%       latitude, elevation, slope and aspect.
%       Relief shading is not considered.
%       Script uses simple unweighed gradient of 4 nearest neighbours for slope
%       calculation.
%
% Reference: Kumar, L, Skidmore AK and Knowles E 1997: Modelling topographic variation in solar radiation in
%            a GIS environment. Int.J.Geogr.Info.Sys. 11(5), 475-497
%
%
% Felix Hebeler, Dept. of Geography, University Zurich, May 2008.

%% parameters
%It ;               % total hours of daily sunshine (calculated inline)
%M ;                % air mass ratio parameter (calculated inline)
%r = 0.20;          % ground reflectance coefficient (more sensible to give as input)
%L=lat;             %latitude
n = 1;              % timestep of calculation over sunshine hours: 1=hourly, 0.5=30min, 2=2hours etc
tau_a    = 365;     %length of the year in days
S0 = 1367;          % solar constant W m^-2   default 1367

dr= 0.0174532925;   % degree to radians conversion factor

%%  convert factors
[slop,asp]=get_ders(dem,cs);   % calculate slope and aspect in radians using given cellsize cs
[dummy,L]=meshgrid(1:size(dem,2),lat);   % grid latitude
clear dummy;
L=L*dr;                     % convert to radians
fcirc = 360*dr; % 360 degrees in radians

%% some setup calculations
srad=0;
sinL=sin(L);
cosL=cos(L);
tanL=tan(L);
sinSlop=sin(slop);
cosSlop=cos(slop);
cosSlop2=cosSlop.*cosSlop;
sinSlop2=sinSlop.*sinSlop;
sinAsp=sin(asp);
cosAsp=cos(asp);
term1 = ( sinL.*cosSlop - cosL.*sinSlop.*cosAsp);
term2 = ( cosL.*cosSlop + sinL.*sinSlop.*cosAsp);
term3 = sinSlop.*sinAsp;
%% loop over year
for d = 1:tau_a;
    %display(['Calculating melt for day ',num2str(d)])
    % clear sky solar radiation
    I0 = S0 * (1 + 0.0344*cos(fcirc*d/tau_a)); % extraterr rad per day
    % sun declination dS
    dS = 23.45 * dr* sin(fcirc * ( (284+d)/tau_a ) ); %in radians, correct/verified
    % angle at sunrise/sunset
    % t = 1:It; % sun hour
    hsr = real(acos(-tanL*tan(dS)));  % angle at sunrise
    % this only works for latitudes up to 66.5 deg N! Workaround:
    % hsr(hsr<-1)=acos(-1);
    % hsr(hsr>1)=acos(1);
    It=round(12*(1+mean(hsr(:))/pi)-12*(1-mean(hsr(:))/pi)); % calc daylength
%%  daily loop
    I=0;
    for t=1:n:It % loop over sunshine hours
        % if accounting for shading should be included, calc hillshade here
        % hourangle of sun hs
        hs=hsr-(pi*t/It);               % hs(t)
        %solar angle and azimuth
        %alpha = asin(sinL*sin(dS)+cosL*cos(dS)*cos(hs));% solar altitude angle
        sinAlpha = sinL.*sin(dS)+cosL.*cos(dS).*cos(hs);
        %alpha_s = asin(cos(dS)*sin(hs)/cos(alpha)); % solar azimuth angle
        % correction  using atmospheric transmissivity taub_b
        M=sqrt(1229+((614.*sinAlpha)).^2)-614.*sinAlpha; % Air mass ratio
        tau_b = 0.56 * (exp(-0.65*M) + exp(-0.095*M));
        tau_d = 0.271-0.294*tau_b; % radiation diffusion coefficient for diffuse insolation
        tau_r = 0.271+0.706*tau_b; % reflectance transmitivity
        % correct for local incident angle
        cos_i = (sin(dS).*term1) + (cos(dS).*cos(hs).*term2) + (cos(dS).*term3.*sin(hs));
        Is = I0 * tau_b; % potential incoming shortwave radiation at surface normal (equator)
        % R = potential clear sky solar radiation W m2
        R = Is .* cos_i;
        R(R<0)=0;  % kick out negative values
        Id = I0 .* tau_d .* cosSlop2./ 2.*sinAlpha; %diffuse radiation;
        Ir = I0 .* r .* tau_r .* sinSlop2./ 2.* sinAlpha; % reflectance
        R= R + Id + Ir;
        R(R<0)=0;
        I=I+R;% solar radiation per day (sunshine hours)
     end % end of sun hours in day loop
%%  add up radiation part melt for every day
    srad = srad + I;
end   % end of days in year loop

%%
function [grad,asp] = get_ders(dem,cs)
% calculate slope and aspect (deg) using GRADIENT function
[fx,fy] = gradient(dem,cs,cs); % uses simple, unweighted gradient of immediate neighbours
[asp,grad]=cart2pol(fy,fx); % convert to carthesian coordinates
grad=atan(grad); %steepest slope
asp=asp.*-1+pi; % convert asp 0 facing south

author: Sable

repo: mcbench-benchmarks


text:
function lambda = lambda_search(x, L, delta, xi)
% lambda = phi^{-1}(delta),
% where phi(lambda)=|min{xi/lambda*E, (L/L+lambda)x}|_F and E=ones(m,n)

    treshold = 1e-12;
    norm_x = norm(x);
    if norm_x <= delta
        lambda = 0;
    else
        n = size(x,1);
        x = abs(x);
        x = sort(x, 'descend');

        if norm_x*(1-xi/(x(1)*L))<delta %This also includes the case: x(k)-xi/L<0
            lambda = L*(norm_x/delta - 1);
        else
            partial_sum = norm_x^2 - x(1)^2;
            flag = 0;
            for k=2:n
                xi_lambda_ratio = x(k)-xi/L;
                if xi_lambda_ratio > 0
                    tail_sum = (k-1)*(xi_lambda_ratio)^2;
                    test_sum = sqrt(tail_sum + (1-xi/(x(k)*L))^2*partial_sum);
                    if test_sum<delta
                        %lambda = L*(sqrt(partial_sum/(delta^2-tail_sum))-1);
                        coef = [delta^2, 2*L*delta^2, (delta^2-partial_sum)*L^2-(k-1)*xi^2, -2*L*(k-1)*xi^2, -(k-1)*xi^2*L^2];
                        lambda = roots(coef);
                        ind=(abs(imag(lambda))<treshold).*(lambda>0);
                        lambda = lambda'*ind;
                        flag = 1;
                        break
                    end
                    partial_sum = partial_sum - x(k)^2;
                else
                    coef = [delta^2, 2*L*delta^2, (delta^2-partial_sum)*L^2-(k-1)*xi^2, -2*L*(k-1)*xi^2, -(k-1)*xi^2*L^2];
                    lambda = roots(coef);
                    ind=(abs(imag(lambda))<treshold).*(lambda>0);
                    lambda = lambda'*ind;
                    flag = 1;
                    break
                end
            end
            if flag == 0
                lambda = sqrt(n)*xi/delta;
            end
        end
    end


text:
function check = student_check ( a, b, c )

%*****************************************************************************80
%
%% STUDENT_CHECK checks the parameter of the central Student T CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, shape parameters of the PDF,
%    used to transform the argument X to a shifted and scaled
%    value Y = ( X - A ) / B.  It is required that B be nonzero.
%    For the standard distribution, A = 0 and B = 1.
%
%    Input, real C, is usually called the number of
%    degrees of freedom of the distribution.  C is typically an
%    integer, but that is not essential.  It is required that
%    C be strictly positive.
%
%    Output, logical CHECK, is true if the parameters are legal.
%
  if ( b == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'STUDENT_CHECK - Fatal error!\n' );
    fprintf ( 1, '  B must be nonzero.\n' );
    check = 0;
    return
  end

  if ( c <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'STUDENT_CHECK - Fatal error!\n' );
    fprintf ( 1, '  C must be greater than 0.\n' );
    check = 0;
    return
  end

  check = 1;

  return
end


text:
function [I,oct] = note2interval(N,key)
% MUSIC.NOTE2INTERVAL Converts a note string to a semitone interval and octave.
%   [I,O] = MUSIC.NOTE2INTERVAL(N) returns the interval I and octave O at which
%   note N is found in the key of 'C'. N is a note in scientific pitch notation
%   (e.g. 'A#3') or a simple note string without an octave specification (e.g.
%   'A#'). If the octave is not defined O will be NaN. N may be a cell array of
%   strings.
%
%   [I,O] = MUSIC.NOTE2INTERVAL(N,KEY) returns the interval at which note N is
%   found in the key of KEY. KEY may be a character note (e.g., 'A', 'F#'), or
%   an interval offset from 'C'.
%
%   Example
%      [I,O] = music.note2interval('F#')       % returns 6, NaN
%      [I,O] = music.note2interval('F#5','G')  % returns 11, 5
%
%   See also music.note2tone, music.note2freq, music.interval2note.

%    Author: E. Johnson
%    Copyright 2010 The MathWorks, Inc.

% Resolve character key to offset from C
if nargin < 2
    key = 0;
end
if ischar(key)
    key = music.note2interval(key);
end

% The interval offset from the beginning of the notes octave, can be used to
% index into this table. Which column is used depends on whether or not the key
% is expressed using sharps or flats.
noteTable = { 'C'   'C' ;
              'C#'  'Db';
              'D'   'D' ;
              'D#'  'Eb';
              'E'   'E' ;
              'F'   'F' ;
              'F#'  'Gb';
              'G'   'G' ;
              'G#'  'Ab';
              'A'   'A' ;
              'A#'  'Bb';
              'B'   'B' ; };

N  = cellstr(N);  % ensure N is a cell array
sz = size(N);     % record original size of N
N  = N(:);        % ensure N is a column vector

% Use regular expression to separate octave from note.
C = regexp(N,'([AaBbCcDdEeFfGg][Bb#]?)(\d*$)','tokens');

note = cellfun( @(c)c{1}{1}, C, 'UniformOutput', false );  % simple note
oct  = cellfun( @(c)str2double(c{1}{2}), C );              % octave

ind  = cellfun(@(n)find(strcmpi(noteTable,n),1),note);

[row,col] = ind2sub(size(noteTable),ind); %#ok

I = row - 1 - key;
I = mod(I,12);

I   = reshape(I,  sz);
oct = reshape(oct,sz);


text:
function matrix = genMatCurl3D(pde,mesh,fem)
%% Generate global matrices and load vector of FEM for 3D ellipic eq
%     -div(A grad u)  = f,    x\in \Omega
% INPUTS:
% pde --- given data function from equation, e.g.
%         pde.A --- diffusion coefficient
%         pde.f --- right hand side function
%         pde.gD --- Dirichlet boundary value function
%         pde.one --- constant function 1.
% mesh --- mesh structure.
% fem --- global degree of freedom of FEM
%
% OUTPUTS:
% matrix.S --- stiffness matrix (w/o boundary condition)
% matrix.A --- final FEM matrix (after boundary condition)
% matrix.rhsF --- load vector (w/o boundary condition)
% matrix.f --- final RHS matrix (after boundary condition)

% Last Modified: 08/07/2020 by Xu Zhang

%% 1. Stiffness Matrix
S = globMatrixNed3D(pde.A,1,mesh,fem,fem);
M = globMatrixNed3D(pde.B,0,mesh,fem,fem);

%% 2. Generate the Right Hand Side Vector
rhsF1 = globNedRHS3D(pde.f1, mesh, fem, 0, 1);
rhsF2 = globNedRHS3D(pde.f2, mesh, fem, 0, 2);
rhsF3 = globNedRHS3D(pde.f3, mesh, fem, 0, 3);
rhsF = rhsF1 + rhsF2 + rhsF3;

%% 3. Dirichlet Boundary Conditions
Atotal = S + M;
tu1 = sum(feval(pde.exactu1,fem.gex,fem.gey,fem.gez).*fem.gew,2);
tu2 = sum(feval(pde.exactu2,fem.gex,fem.gey,fem.gez).*fem.gew,2);
tu3 = sum(feval(pde.exactu3,fem.gex,fem.gey,fem.gez).*fem.gew,2);
tgt = mesh.p(mesh.e(:,2),:) - mesh.p(mesh.e(:,1),:);
tgt = tgt./sum(tgt.^2,2).^(1/2);
tu = tu1.*tgt(:,1) + tu2.*tgt(:,2) + tu3.*tgt(:,3);
ub = tu;
ub(fem.mapper) = 0;
rhsB = Atotal*ub;
A = Atotal(fem.mapper,fem.mapper);
f = rhsF(fem.mapper) - rhsB(fem.mapper);

%% Outputs
matrix = struct('A', A, 'f', f, 'S', S, 'rhsF',rhsF,'tu',tu);

repo: ifem


text:
function lambda = kms_eigenvalues ( alpha, n )

%*****************************************************************************80
%
%% KMS_EIGENVALUES returns the eigenvalues of the KMS matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 June 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    William Trench,
%    Spectral decomposition of Kac-Murdock-Szego matrices,
%    Unpublished technical document.
%
%  Parameters:
%
%    Input, real ALPHA, the scalar that defines A.
%    Eigenvalue computations require 0 <= ALPHA <= 1.
%
%    Input, integer N, the order of the matrix.
%
%    Output, real LAMBDA(N,1), the eigenvalues.
%
  theta(1:n,1) = kms_eigenvalues_theta ( alpha, n );

  lambda(1:n,1) = ( 1.0 + alpha ) * ( 1.0 - alpha ) ...
    ./ ( 1.0 - 2.0 * alpha * cos ( theta(1:n,1) ) + alpha * alpha );

  return
end