text:
function stla_display ( input_stla_filename )

%*****************************************************************************80
%
%% STLA_DISPLAY displays the faces of a shape defined by an ASCII STL file.
%
%  Usage:
%
%    stla_display ( 'file.stla' )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 January 2015
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'STLA_DISPLAY\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Reads an object in an ASCII STL file.\n' );
  fprintf ( 1, '  Display it as a MATLAB shape.\n' );
%
%  If at least one command line argument, it's the input file name.
%
  if ( nargin < 1 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'STLA_DISPLAY:\n' );
    input_stla_filename = input ( 'Enter the name of the input file:' );

  end
%
%  Get sizes.
%
  [ solid_num, node_num, face_num, text_num ] = stla_size ( ...
    input_stla_filename );
%
%  Print the sizes.
%
  stla_size_print ( input_stla_filename, solid_num, node_num, face_num, ...
    text_num );
%
%  Get the data.
%
  [ node_xyz, face_node, face_normal, ierror ] = stla_read ( ...
    input_stla_filename, node_num, face_num );

  if ( ierror ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'STLA_DISPLAY - Fatal error!\n' );
    fprintf ( 1, '  STLA_READ returned IERROR = %d\n', ierror );
    return
  end
%
%  Display the shape.
%  The argument pair
%    'EdgeColor', 'None'
%  suppresses the display of the black boundary lines between the patches.
%
  h3 = trisurf ( face_node', node_xyz(1,:), node_xyz(2,:), node_xyz(3,:), ...
      'FaceColor', 'Interp', 'EdgeColor', 'None' );

  axis equal;

  xlabel ( '--X axis--' )
  ylabel ( '--Y axis--' )
  zlabel ( '--Z axis--' )
%
%  The TITLE function will interpret underscores in the title.
%  We need to unescape such escape sequences!
%
  title_string = s_escape_tex ( input_stla_filename );
  title ( title_string )
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'STLA_DISPLAY:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function c2 = ch_cap ( c )

%*****************************************************************************80
%
%% CH_CAP capitalizes a single character.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, the character to capitalize.
%
%    Output, character C2, the capitalized character.
%
  if ( 'a' <= c & c <= 'z' )
    c2 = c + 'A' - 'a';
  else
    c2 = c;
  end

  return
end
function truefalse = ch_eqi ( c1, c2 )

%*****************************************************************************80
%
%% CH_EQI is a case insensitive comparison of two characters for equality.
%
%  Example:
%
%    CH_EQI ( 'A', 'a' ) is TRUE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 July 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C1, C2, the characters to compare.
%
%    Output, logical TRUEFALSE, is TRUE (1) if the characters are equal.
%
  FALSE = 0;
  TRUE = 1;

  if ( ch_cap ( c1 ) == ch_cap ( c2 ) )
    truefalse = TRUE;
  else
    truefalse = FALSE;
  end

  return
end
function truefalse = ch_is_digit ( c )

%*****************************************************************************80
%
% CH_IS_DIGIT returns TRUE if the character C is a digit.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 December 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, a character.
%
%    Output, integer TRUEFALSE, is TRUE (1) if C is a digit, FALSE (0) otherwise.
%
  TRUE = 1;
  FALSE = 0;

  if ( '0' <= c & c <= '9' )
    truefalse = TRUE;
  else
    truefalse = FALSE;
  end

  return
end
function digit = ch_to_digit ( c )

%*****************************************************************************80
%
%% CH_TO_DIGIT returns the integer value of a base 10 digit.
%
%  Example:
%
%     C   DIGIT
%    ---  -----
%    '0'    0
%    '1'    1
%    ...  ...
%    '9'    9
%    ' '    0
%    'X'   -1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, the decimal digit, '0' through '9' or blank
%    are legal.
%
%    Output, integer DIGIT, the corresponding integer value.  If C was
%    'illegal', then DIGIT is -1.
%
  if ( '0' <= c & c <= '9' )

    digit = c - '0';

  elseif ( c == ' ' )

    digit = 0;

  else

    digit = -1;

  end

  return
end
function s3 = s_cat ( s1, s2 )

%*****************************************************************************80
%
%% S_CAT concatenates two strings to make a third string.
%
%  Discussion:
%
%    MATLAB provides the STRCAT function, which you should
%    probably use instead of this function!
%
%    s3 = strcat ( s1, s2 )
%
%    Although STRCAT does not ignore trailing blanks.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 May 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S1, the "prefix" string.
%
%    Input, string S2, the "postfix" string.
%
%    Output, string S3, the string made by
%    concatenating S1 and S2, ignoring any trailing blanks.
%
  l1 = s_len_trim ( s1 );
  l2 = s_len_trim ( s2 );
%
%  The following line essentially "declares" S3 to be a
%  character string.  Omitting this line would cause S3
%  to be regarded as a numeric array.
%
  s3 = '';
%
%  Now copy the strings.
%
  s3(   1:l1   ) = s1(1:l1);
  s3(l1+1:l1+l2) = s2(1:l2);

  return
end
function value = s_eqi ( s1, s2 )

%*****************************************************************************80
%
%% S_EQI is a case insensitive comparison of two strings for equality.
%
%  Example:
%
%    S_EQI ( 'Anjana', 'ANJANA' ) is TRUE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character ( len = * ) S1, S2, the strings to compare.
%
%    Output, logical VALUE, is TRUE if the strings are equal.
%
  FALSE = 0;
  TRUE = 1;

  len1 = length ( s1 );
  len2 = length ( s2 );
  lenc = min ( len1, len2 );

  value = FALSE;

  for i = 1 : lenc

    c1 = ch_cap ( s1(i) );
    c2 = ch_cap ( s2(i) );

    if ( c1 ~= c2 )
      value = FALSE;
      return
    end

  end

  for i = lenc + 1 : len1
    if ( s1(i) ~= ' ' )
      value = FALSE;
      return
    end
  end

  for i = lenc + 1 : len2
    if ( s2(i) ~= ' ' )
      value = FALSE;
      return
    end
  end

  value = TRUE;

  return
end
function s2 = s_escape_tex ( s1 )

%*****************************************************************************80
%
%% S_ESCAPE_TEX de-escapes TeX escape sequences.
%
%  Discussion:
%
%    In particular, every occurrence of the characters '\', '_',
%    '^', '{' and '}' will be replaced by '\\', '\_', '\^',
%    '\{' and '\}'.  A TeX interpreter, on seeing these character
%    strings, is then likely to return the original characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S1, the string to be de-escaped.
%
%    Output, string S2, a copy of the string, modified to avoid TeX escapes.
%
  s1_length = length ( s1 );

  s1_pos = 0;
  s2_pos = 0;
  s2 = [];

  while ( s1_pos < s1_length )

    s1_pos = s1_pos + 1;

    if ( s1(s1_pos) == '\' | ...
         s1(s1_pos) == '_' | ...
         s1(s1_pos) == '^' | ...
         s1(s1_pos) == '{' | ...
         s1(s1_pos) == '}' )
      s2_pos = s2_pos + 1;
      s2 = strcat ( s2, '\' );
    end

    s2_pos = s2_pos + 1;
    s2 = strcat ( s2, s1(s1_pos) );

  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 [ r, lchar, ierror ] = s_to_r8 ( s )

%*****************************************************************************80
%
%% S_TO_R8 reads an R8 from a string.
%
%  Discussion:
%
%    This routine will read as many characters as possible until it reaches
%    the end of the string, or encounters a character which cannot be
%    part of the real number.
%
%    Legal input is:
%
%       1 blanks,
%       2 '+' or '-' sign,
%       2.5 spaces
%       3 integer part,
%       4 decimal point,
%       5 fraction part,
%       6 'E' or 'e' or 'D' or 'd', exponent marker,
%       7 exponent sign,
%       8 exponent integer part,
%       9 exponent decimal point,
%      10 exponent fraction part,
%      11 blanks,
%      12 final comma or semicolon.
%
%    with most quantities optional.
%
%  Example:
%
%    S                 R
%
%    '1'               1.0
%    '     1   '       1.0
%    '1A'              1.0
%    '12,34,56'        12.0
%    '  34 7'          34.0
%    '-1E2ABCD'        -100.0
%    '-1X2ABCD'        -1.0
%    ' 2E-1'           0.2
%    '23.45'           23.45
%    '-4.2E+2'         -420.0
%    '17d2'            1700.0
%    '-14e-2'         -0.14
%    'e2'              100.0
%    '-12.73e-9.23'   -12.73 * 10.0**(-9.23)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string containing the
%    data to be read.  Reading will begin at position 1 and
%    terminate at the end of the string, or when no more
%    characters can be read to form a legal real.  Blanks,
%    commas, or other nonnumeric data will, in particular,
%    cause the conversion to halt.
%
%    Output, real R, the value that was read from the string.
%
%    Output, integer LCHAR, the number of characters of S that were used to form R.
%
%    Output, integer IERROR, is 0 if no error occurred.
%
  s_length = s_len_trim ( s );
  ierror = 0;
  r = 0.0;
  lchar = -1;
  isgn = 1;
  rtop = 0.0;
  rbot = 1.0;
  jsgn = 1;
  jtop = 0;
  jbot = 1;
  ihave = 1;
  iterm = 0;

  while ( 1 )

    lchar = lchar + 1;
    c = s(lchar+1);
%
%  Blank character.
%
    if ( c == ' ' )

      if ( ihave == 2 )

      elseif ( ihave == 6 | ihave == 7 )
        iterm = 1;
      elseif ( 1 < ihave )
        ihave = 11;
      end
%
%  Comma.
%
    elseif ( c == ',' | c == ';' )

      if ( ihave ~= 1 )
        iterm = 1;
        ihave = 12;
        lchar = lchar + 1;
      end
%
%  Minus sign.
%
    elseif ( c == '-' )

      if ( ihave == 1 );
        ihave = 2;
        isgn = -1;
      elseif ( ihave == 6 )
        ihave = 7;
        jsgn = -1;
      else
        iterm = 1;
      end
%
%  Plus sign.
%
    elseif ( c == '+' )

      if ( ihave == 1 )
        ihave = 2;
      elseif ( ihave == 6 )
        ihave = 7;
      else
        iterm = 1;
      end
%
%  Decimal point.
%
    elseif ( c == '.' )

      if ( ihave < 4 )
        ihave = 4;
      elseif ( 6 <= ihave & ihave <= 8 )
        ihave = 9;
      else
        iterm = 1;
      end
%
%  Exponent marker.
%
    elseif ( ch_eqi ( c, 'E' ) | ch_eqi ( c, 'D' ) )

      if ( ihave < 6 )
        ihave = 6;
      else
        iterm = 1;
      end
%
%  Digit.
%
    elseif ( ihave < 11 & ch_is_digit ( c ) )

      if ( ihave <= 2 )
        ihave = 3;
      elseif ( ihave == 4 )
        ihave = 5;
      elseif ( ihave == 6 | ihave == 7 )
        ihave = 8;
      elseif ( ihave == 9 )
        ihave = 10;
      end

      d = ch_to_digit ( c );

      if ( ihave == 3 )
        rtop = 10.0 * rtop + d;
      elseif ( ihave == 5 )
        rtop = 10.0 * rtop + d;
        rbot = 10.0 * rbot;
      elseif ( ihave == 8 )
        jtop = 10 * jtop + d;
      elseif ( ihave == 10 )
        jtop = 10 * jtop + d;
        jbot = 10 * jbot;
      end
%
%  Anything else is regarded as a terminator.
%
    else
      iterm = 1;
    end
%
%  If we haven't seen a terminator, and we haven't examined the
%  entire string, go get the next character.
%
    if ( iterm == 1 | s_length <= lchar + 1 )
      break;
    end

  end
%
%  If we haven't seen a terminator, and we have examined the
%  entire string, then we're done, and LCHAR is equal to S_LENGTH.
%
  if ( iterm ~= 1 & lchar + 1 == s_length )
    lchar = s_length;
  end
%
%  Number seems to have terminated.  Have we got a legal number?
%  Not if we terminated in states 1, 2, 6 or 7!
%
  if ( ihave == 1 | ihave == 2 | ihave == 6 | ihave == 7 )
    ierror = ihave;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'S_TO_R8 - Fatal error!\n' );
    fprintf ( 1, '  IHAVE = %d\n', ihave );
    error ( 'S_TO_R8 - Fatal error!' );
  end
%
%  Number seems OK.  Form it.
%
  if ( jtop == 0 )
    rexp = 1.0;
  else

    if ( jbot == 1 )
      rexp = 10.0^( jsgn * jtop );
    else
      rexp = jsgn * jtop;
      rexp = rexp / jbot;
      rexp = 10.0^rexp;
    end

  end

  r = isgn * rexp * rtop / rbot;

  return
end
function [ node_xyz, face_node, face_normal, ierror ] = stla_read ( ...
  input_file_name, node_num, face_num )

%*****************************************************************************80
%
%% STLA_READ reads graphics information from an ASCII StereoLithography file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    3D Systems, Inc,
%    Stereolithography Interface Specification,
%    October 1989.
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the input file.
%
%    Input, integer NODE_NUM, the number of vertices defined.
%
%    Input, integer FACE_NUM, the number of faces defined.
%
%    Output, real NODE_XYZ(3,NODE_NUM), the coordinates of points.
%
%    Output, integer FACE_NODE(3,FACE_NUM), the nodes that make up each face.
%
%    Output, real FACE_NORMAL(3,FACE_NUM), the normal vector
%    at each face.
%
  ierror = 0;
  state = 0;
  text_num = 0;
  face = 0;
  node = 0;

  node_xyz = zeros ( 3, node_num );
  face_node = zeros ( 3, face_num );
  face_normal = zeros ( 3, face_num );
%
%  Open the file.
%
  iunit = fopen ( input_file_name, 'r' );

  if ( iunit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'STLA_READ - Fatal error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    ierror = 1;
    return
  end
%
%  Read the next line of text.
%
  while ( 1 )

    text = fgetl ( iunit );

    if ( text == -1 )
      if ( state ~= 0 & state ~= 1 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  End-of-file, but model not finished.\n' );
        ierror = 1;
        return
      end
      break;
    end

    text_num = text_num + 1;
    done = 1;
%
%  Read the first word in the line.
%
    [ word1, done ] = word_next_read ( text, done );
%
%  "Doctor" the text, changing a beginning occurrence of:
%
%      END FACET to ENDFACET
%      END LOOP to ENDLOOP
%      END SOLID to ENDSOLID
%      FACET NORMAL to FACETNORMAL
%      OUTER LOOP to OUTERLOOP
%
    if ( s_eqi ( word1, 'END' ) )

      [ word2, done ] = word_next_read ( text, done );

      if ( ~s_eqi ( word2, 'FACET' ) & ...
           ~s_eqi ( word2, 'LOOP' ) & ...
           ~s_eqi ( word2, 'SOLID' ) )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  The tag END was followed by an illegal\n' );
        fprintf ( 1, '  word: "%s", when expecting\n', word );
        fprintf ( 1, '  "FACET", "LOOP", or "SOLID".\n' );
        ierror = 1;
        return
      end

      word1 = s_cat ( word1, word2 );

    elseif ( s_eqi ( word1, 'FACET' ) )

      [ word2, done ] = word_next_read ( text, done );

      if ( ~s_eqi ( word2, 'NORMAL' ) )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  The tag FACET was followed by an illegal\n' );
        fprintf ( 1, '  word: "%s", when expecting\n', word2 );
        fprintf ( 1, '  "NORMAL".\n' );
        ierror = 1;
        return
      end

      word1 = s_cat ( word1, word2 );

    elseif ( s_eqi ( word1, 'OUTER' ) )

      [ word2, done ] = word_next_read ( text, done );

      if ( ~s_eqi ( word2, 'LOOP' ) )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  The tag OUTER was followed by an illegal\n' );
        fprintf ( 1, '  word: "%s", when expecting\n', word2 );
        fprintf ( 1, '  "LOOP".\n' );
        ierror = 1;
        return
      end

      word1 = s_cat ( word1, word2 );

    end
%
%  This first word tells us what to do.
%
%  SOLID - begin a new solid.
%    Valid in state 0, moves to state 1.
%  ENDSOLID - end current solid.
%    Valid in state 1, moves to state 0.
%
%  FACETNORMAL - begin a new facet.
%    Valid in state 0 or 1, moves to state 2.
%  ENDFACET - end current facet.
%    Valid in state 2, moves to state 1.
%
%  OUTERLOOP - begin a list of vertices.
%    Valid in state 2, moves to state 3, sets vertex count to 0.
%  ENDLOOP - end vertex list.
%    Valid in state 3, moves to state 2.
%
%  VERTEX - give coordinates of next vertex.
%    Valid in state 3.
%
%  End of file -
%    Valid in state 0 or 1.
%
    if ( s_eqi ( word1, 'SOLID' ) )

      if ( state ~= 0 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  A new SOLID statement was encountered, but we\n' );
        fprintf ( 1, '  have not finished processing the current solid.\n' );
        ierror = 1
        return
      end

      state = 1;

    elseif ( s_eqi ( word1, 'ENDSOLID' ) )

      if ( state ~= 1 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  An END SOLID statement was encountered, but\n' );
        fprintf ( 1, '  either we have not begun a solid at all, or we\n' );
        fprintf ( 1, '  are not at an appropriate point to finish the\n' );
        fprintf ( 1, '  current solid.\n' );
        ierror = 1;
        return
      end

      state = 0;

    elseif ( s_eqi ( word1, 'FACETNORMAL' ) )

      if ( state ~= 0 & state ~= 1 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  Model not in right state for FACET.\n' );
        ierror = 1;
        return
      end

      state = 2;
      face = face + 1;

      if ( face_num < face )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Warning!\n' );
        fprintf ( 1, '  More faces being read than expected.\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        ierror = 1;
        return
      end

      for i = 1 : 3

        face_normal(i,face) = 0.0;

        [ word2, done ] = word_next_read ( text, done );

        if ( done )
          fprintf ( 1, '\n' );
          fprintf ( 1, 'STLA_READ - Fatal error!\n' );
          fprintf ( 1, '  File line number = %d\n', text_num );
          fprintf ( 1, '  End of information while reading a component\n' );
          fprintf ( 1, '  of the normal vector.\n' );
          ierror = 1;
          return
        end

        [ dval, lchar, ierror ] = s_to_r8 ( word2 );

        if ( ierror ~= 0 )
          fprintf ( 1, '\n' );
          fprintf ( 1, 'STLA_READ - Fatal error!\n' );
          fprintf ( 1, '  File line number = %d\n', text_num );
          fprintf ( 1, '  Reading a component of the normal vector.\n' );
          ierror = 1;
          return
        end

        face_normal(i,face) = dval;

      end

    elseif ( s_eqi ( word1, 'ENDFACET' ) )

      if ( state ~= 2 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  Model not in right state for ENDFACET.\n' );
        ierror = 1;
        return
      end

      state = 1;

    elseif ( s_eqi ( word1, 'OUTERLOOP' ) )

      if ( state ~= 2 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  Model not in right state for OUTERLOOP.\n' );
        ierror = 1;
        return
      end

      state = 3;
      vertex = 0;

    elseif ( s_eqi ( word1, 'ENDLOOP' ) )

      if ( state ~= 3 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  Model not in right state for ENDLOOP.\n' );
        ierror = 1;
        return
      end

      state = 2;

    elseif ( s_eqi ( word1, 'VERTEX' ) )

      if ( state ~= 3 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  Model not in right state for VERTEX.\n' );
        ierror = 1;
        return
      end

      if ( 3 <= vertex )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  More than 3 vertices specified for a face.\n' );
        ierror = 1;
        return
      end

      for i = 1 : 3

        [ word2, done ] = word_next_read ( text, done );

        if ( done )
          fprintf ( 1, '\n' );
          fprintf ( 1, 'STLA_READ - Fatal error!\n' );
          fprintf ( 1, '  File line number = %d\n', text_num );
          fprintf ( 1, '  The value of a vertex coordinate is missing.\n' );
          ierror = 1;
          return
        end

        [ dval, lchar, ierror ] = s_to_r8 ( word2 );

        if ( ierror ~= 0 )
          fprintf ( 1, '\n' );
          fprintf ( 1, 'STLA_READ - Fatal error!\n' );
          fprintf ( 1, '  File line number = %d\n', text_num );
          fprintf ( 1, '  The value of a vertex coordinate makes\n' );
          fprintf ( 1, '  no sense.\n' );
          ierror = 1;
          return
        end

        temp(i) = dval;

      end

      if ( node_num <= node )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Warning!\n' );
        fprintf ( 1, '  More nodes being read than expected.\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        ierror = 1;
        return
      end

      node = node + 1;
      node_xyz(1:3,node) = temp(1:3);

      vertex = vertex + 1;
      face_node(vertex,face) = node;

    else

      fprintf ( 1, '\n' );
      fprintf ( 1, 'STLA_READ - Fatal error!\n' );
      fprintf ( 1, '  File line number = %d\n', text_num );
      fprintf ( 1, '  Unrecognized line in file.\n' );
      ierror = 1;
      return

    end

  end
%
%  Close the file.
%
  fclose ( iunit );

  return
end
function [ solid_num, node_num, face_num, text_num ] = stla_size ( ...
  input_file_name )

%*****************************************************************************80
%
%% STLA_SIZE determines sizes associated with an STLA file.
%
%  Discussion:
%
%    This routine assumes that the file is a legal STLA file.
%
%    To perform checks on the file, call STLA_CHECK first.
%
%    Note that the counts for the number of nodes and edges are
%    overestimates, since presumably, most nodes will be defined several
%    times, once for each face they are part of, and most edges will
%    be defined twice.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    3D Systems, Inc,
%    Stereolithography Interface Specification,
%    October 1989.
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the input file.
%
%    Output, integer SOLID_NUM, the number of solids defined.
%    Presumably, this is 1.
%
%    Output, integer NODE_NUM, the number of vertices defined.
%
%    Output, integer FACE_NUM, the number of faces defined.
%
%    Output, integer TEXT_NUM, the number of lines of text.
%
  ierror = 0;

  state = 0;

  text_num = 0;
  solid_num = 0;
  node_num = 0;
  face_num = 0;
%
%  Open the file.
%
  iunit = fopen ( input_file_name, 'r' );

  if ( iunit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'STLA_SIZE - Fatal error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    return
  end
%
%  Read the next line of text.
%
  while ( 1 )

    text = fgetl ( iunit );

    if ( text == -1 )
      if ( state ~= 0 & state ~= 1 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_SIZE - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  End-of-file, but model not finished.\n' );
        return
      end
      break;
    end

    text_num = text_num + 1;

    done = 1;
%
%  Read the first word in the line.
%
    [ word1, done ] = word_next_read ( text, done );

    if ( done )
      return
    end
%
%  "Doctor" the text, changing a beginning occurrence of:
%
%      END FACET to ENDFACET
%      END LOOP to ENDLOOP
%      END SOLID to ENDSOLID
%      FACET NORMAL to FACETNORMAL
%      OUTER LOOP to OUTERLOOP
%
    if ( s_eqi ( word1, 'END' ) )

      [ word2, done ] = word_next_read ( text, done );

      if ( ~s_eqi ( word2, 'FACET' ) & ...
           ~s_eqi ( word2, 'LOOP' ) & ...
           ~s_eqi ( word2, 'SOLID' ) )
        return
      end

      word1 = s_cat ( word1, word2 );

    elseif ( s_eqi ( word1, 'FACET' ) )

      [word2, done ] = word_next_read ( text, done );

      if ( ~s_eqi ( word2, 'NORMAL' ) )
        return
      end

      word1 = s_cat ( word1, word2 );

    elseif ( s_eqi ( word1, 'OUTER' ) )

      [ word2, done ] = word_next_read ( text, done );

      if ( ~s_eqi ( word2, 'LOOP' ) )
        return
      end

      word1 = s_cat ( word1, word2 );

    end
%
%  This first word tells us what to do.
%
%  SOLID - begin a new solid.
%    Valid in state 0, moves to state 1.
%  ENDSOLID - end current solid.
%    Valid in state 1, moves to state 0.
%
%  FACETNORMAL - begin a new facet.
%    Valid in state 0 or 1, moves to state 2.
%  ENDFACET - end current facet.
%    Valid in state 2, moves to state 1.
%
%  OUTERLOOP - begin a list of vertices.
%    Valid in state 2, moves to state 3.
%  ENDLOOP - end vertex list.
%    Valid in state 3, moves to state 2.
%
%  VERTEX - give coordinates of next vertex.
%    Valid in state 3 if current vertex count is 0, 1 or 2.
%
%  End of file -
%    Valid in state 0 or 1.
%
    if ( s_eqi ( word1, 'SOLID' ) )

      if ( state ~= 0 )
        return
      end

      state = 1;

    elseif ( s_eqi ( word1, 'ENDSOLID' ) )

      if ( state ~= 1 )
        return
      end

      state = 0;
      solid_num = solid_num + 1;

    elseif ( s_eqi ( word1, 'FACETNORMAL' ) )

      if ( state ~= 0 & state ~= 1 )
        return
      end

      state = 2;

      for i = 1 : 3

        [ word2, done ] = word_next_read ( text, done );

        if ( done )
          return
        end

        [ dval, lchar, ierror ] = s_to_r8 ( word2 );

        if ( ierror ~= 0 )
          return
        end

      end

    elseif ( s_eqi ( word1, 'ENDFACET' ) )

      if ( state ~= 2 )
        return
      end

      state = 1;

      face_num = face_num + 1;

    elseif ( s_eqi ( word1, 'OUTERLOOP' ) )

      if ( state ~= 2 )
        return
      end

      state = 3;
      vertex = 0;

    elseif ( s_eqi ( word1, 'ENDLOOP' ) )

      if ( state ~= 3 )
        return
      end

      state = 2;

    elseif ( s_eqi ( word1, 'VERTEX' ) )

      if ( state ~= 3 )
        return
      end

      if ( 3 <= vertex )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_SIZE - Fatal error!\n' );
        fprintf ( 1, '  Too many vertices for a face.\n' );
        ierror = 1;
        return
      end

      for i = 1 : 3

        [ word2, done ] = word_next_read ( text, done );

        if ( done )
          return
        end

        [ dval, lchar, ierror ] = s_to_r8 ( word2 );

        if ( ierror ~= 0 )
          return
        end

      end

      vertex = vertex + 1;
      node_num = node_num + 1;

    else

      return

    end

  end
%
%  Close the file.
%
  fclose ( iunit );

  return
end
function stla_size_print ( input_file_name, solid_num, node_num, face_num, ...
  text_num )

%*****************************************************************************80
%
%% STLA_SIZE_PRINT prints sizes associated with an STLA file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    3D Systems, Inc,
%    Stereolithography Interface Specification,
%    October 1989.
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the input file.
%
%    Input, integer SOLID_NUM, the number of solids defined.
%
%    Input, integer NODE_NUM, the number of vertices defined.
%
%    Input, integer FACE_NUM, the number of faces defined.
%
%    Input, integer TEXT_NUM, the number of lines of text.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Object sizes for STLA file "%s":\n', input_file_name );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solids =                   %d\n', solid_num );
  fprintf ( 1, '  Nodes (may be repeated) =  %d\n', node_num );
  fprintf ( 1, '  Faces (triangular only) =  %d\n', face_num );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of lines of text =  %d\n', text_num );

  return
end
function timestamp ( )

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

  return
end
function [ word, done ] = word_next_read ( s, done )

%*****************************************************************************80
%
%% WORD_NEXT_READ "reads" words from a string, one at a time.
%
%  Special cases:
%
%    The following characters are considered to be a single word,
%    whether surrounded by spaces or not:
%
%      " ( ) { } [ ]
%
%    Also, if there is a trailing comma on the word, it is stripped off.
%    This is to facilitate the reading of lists.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 September 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, a string, presumably containing words
%    separated by spaces.
%
%    Input, logical DONE.
%    TRUE, if we are beginning a new string;
%    FALSE, if we are continuing to process the current string.
%
%    Output, string WORD.
%    If DONE is FALSE, then WORD contains the "next" word read.
%    If DONE is TRUE, then WORD is blank, because there was no more to read.
%
%    Output, logical DONE.
%      FALSE if another word was read,
%      TRUE if no more words could be read.
%
  persistent lenc;
  persistent next;

  tab = char ( 9 );
%
%  We "remember" LENC and NEXT from the previous call.
%
%  An input value of DONE = TRUE signals a new line of text to examine.
%
  if ( done )

    next = 1;
    done = 0;
    lenc = s_len_trim ( s );

    if ( lenc <= 0 )
      done = 1;
      word = ' ';
      return
    end

  end
%
%  Beginning at index NEXT, search the string for the next nonblank,
%  which signals the beginning of a word.
%
  ilo = next;
%
%  ...S(NEXT:) is blank.  Return with WORD = ' ' and DONE = TRUE.
%
  while ( 1 )

    if ( lenc < ilo )
      word = ' ';
      done = 1;
      next = lenc + 1;
      return
    end
%
%  If the current character is blank, skip to the next one.
%
    if ( s(ilo) ~= ' ' & s(ilo) ~= tab )
      break
    end

    ilo = ilo + 1;

  end
%
%  ILO is the index of the next nonblank character in the string.
%
%  If this initial nonblank is a special character,
%  then that's the whole word as far as we're concerned,
%  so return immediately.
%
  if ( s(ilo) == '"' | ...
       s(ilo) == '(' | ...
       s(ilo) == ')' | ...
       s(ilo) == '{' | ...
       s(ilo) == '}' | ...
       s(ilo) == '[' | ...
       s(ilo) == ']' )

    word = s(ilo);
    next = ilo + 1;
    return

  end
%
%  Now search for the last contiguous character that is not a
%  blank, TAB, or special character.
%
  next = ilo + 1;

  while ( next <= lenc )

    if ( s(next) == ' ' )
      break;
    elseif ( s(next) == tab )
      break;
    elseif ( s(next) == '"' )
      break;
    elseif ( s(next) == '(' )
      break;
    elseif ( s(next) == ')' )
      break;
    elseif ( s(next) == '{' )
      break;
    elseif ( s(next) == '}' )
      break;
    elseif ( s(next) == '[' )
      break;
    elseif ( s(next) == ']' )
      break;
    end

    next = next + 1;

  end

  if ( s(next-1) == ',' )
    word = s(ilo:next-2);
  else
    word = s(ilo:next-1);
  end

  return
end


text:
1     % GMRES/Bicgstab(2)
1e-6  % tolerance
100   % maximum number of iterations
3     % lst sqrs commutator preconditioner
1     % ideal
5     % Print to (new) Figure 4
3     % Color

% Data file for it_solve for NS2

author: johannesgerer

repo: jburkardt-m


text:
function eloo = gpla_looe(w, gp, x, y, varargin)
%GPLA_LOOE  Evaluate the mean negative log leave-one-out predictive
%           density with Laplace approximation.
%
%  Description
%    LOOE = GPLA_LOOE(W, GP, X, Y, OPTIONS) takes a parameter
%    vector W, Gaussian process structure GP, a matrix X of input
%    vectors and a matrix Y of targets, and evaluates the mean
%    negative log leave-one-out predictive density
%      LOOE = - 1/n sum log p(Y_i | X, Y_{\i}, th)
%    where th represents the parameters (lengthScale,
%    magnSigma2...), X is inputs and Y is observations.
%
%    Laplace leave-one-out is approximated by using analytic
%    leave-one-out equations for the approximated latent
%    distribution.
%
%    OPTIONS is optional parameter-value pair
%      z - optional observed quantity in triplet (x_i,y_i,z_i)
%          Some likelihoods may use this. For example, in case of
%          Poisson likelihood we have z_i=E_i, that is, expected
%          value for ith case.
%
%  See also
%    GP_LOOE, GPLA_LOOPRED, GPLA_E
%
% Copyright (c) 2011 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.FunctionName = 'GPLA_LOOE';
ip.addRequired('w', @(x) isempty(x) || ...
               isvector(x) && isreal(x) && all(isfinite(x)));
ip.addRequired('gp', @(x) isstruct(x) || iscell(x));
ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.parse(w, gp, x, y, varargin{:});
z=ip.Results.z;

gp=gp_unpak(gp, w);
[tmp, tmp, lpyt] = gpla_loopred(gp, x, y, 'z', z);
eloo=-sum(lpyt);

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function too_young = under_age(age, limit)
    if nargin == 1
        limit = 21;
    end
    too_young = age < limit;
end

author: anishLearnsToCode

repo: introduction-to-programming-with-matlab


text:
function test_tutorial_connectivity2

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY

% contains the code that is in the second part of the ocnnectivity
% tutorial. Purpose: simulate linearly mixed data, and show that dependent
% on the parameter settings the connectivity can go wild.

% create some instantaneously mixed data

% define some variables locally
nTrials  = 100;
nSamples = 1000;
fsample  = 1000;

% mixing matrix
mixing   = [0.8 0.2 0;
              0 0.2 0.8];

data       = [];
data.trial = cell(1,nTrials);
data.time  = cell(1,nTrials);
for k = 1:nTrials
  dat = randn(3, nSamples);
  dat(2,:) = ft_preproc_bandpassfilter(dat(2,:), 1000, [15 25]);
  dat = 0.2.*(dat-repmat(mean(dat,2),[1 nSamples]))./repmat(std(dat,[],2),[1 nSamples]);
  data.trial{k} = mixing * dat;
  data.time{k}  = (0:nSamples-1)./fsample;
end
data.label = {'chan1' 'chan2'}';

figure;plot(dat'+repmat([0 1 2],[nSamples 1]));
title('original ''sources''');

figure;plot((mixing*dat)'+repmat([0 1],[nSamples 1]));
axis([0 1000 -1 2]);
set(findobj(gcf,'color',[0 0.5 0]), 'color', [1 0 0]);
title('mixed ''sources''');

%% do spectral analysis
cfg = [];
cfg.method = 'mtmfft';
cfg.output = 'fourier';
cfg.foilim = [0 200];
cfg.tapsmofrq = 5;
freq = ft_freqanalysis(cfg, data);
fd   = ft_freqdescriptives([], freq);

figure;plot(fd.freq, fd.powspctrm);
set(findobj(gcf,'color',[0 0.5 0]), 'color', [1 0 0]);
title('powerpectrum');

%% compute connectivity
cfg = [];
cfg.method = 'granger';
g = ft_connectivityanalysis(cfg, freq);
cfg.method = 'coh';
c = ft_connectivityanalysis(cfg, freq);

%% visualize
cfg = [];
cfg.parameter = 'grangerspctrm';
figure;ft_connectivityplot(cfg, g);
title('granger causality');
cfg.parameter = 'cohspctrm';
figure;ft_connectivityplot(cfg, c);
title('coherence');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef TestThinning
    %TestThinning

    methods (Static)
        function test_1
            bw = cv.imread(fullfile(mexopencv.root(),'test','bw.png'), ...
                'Grayscale',true, 'ReduceScale',2);
            out = cv.thinning(bw, 'ThinningType','ZhangSuen');
            validateattributes(out, {class(bw)}, {'size',size(bw)});
        end

        function test_error_argnum
            try
                cv.thinning();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [feat, x, y, wid, hgt] = extract_feature(feature, img, c)
%
% Copyright Aditya Khosla http://mit.edu/khosla
%
% Please cite this paper if you use this code in your publication:
%   A. Khosla, J. Xiao, A. Torralba, A. Oliva
%   Memorability of Image Regions
%   Advances in Neural Information Processing Systems (NIPS) 2012
%

if(~exist('c', 'var'))
  c = conf();
end

[feat, x, y, wid, hgt] = feval(['extract_' feature], img, c);

author: adikhosla

repo: feature-extraction


text:
% DEMOIL1 Oil data with fully independent training conditional.

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'oil';
experimentNo = 7;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = fgplvmOptions('fitc');
options.optimiser = 'scg';
latentDim = 2;
d = size(Y, 2);
options.kern = {'gibbs', 'bias', 'white'};
model = fgplvmCreate(latentDim, d, Y, options);

% Optimise the model.
iters = 1000;
display = 1;

model = fgplvmOptimise(model, display, iters);

% Save the results.
capName = dataSetName;;
capName(1) = upper(capName(1));
save(['dem' capName num2str(experimentNo) '.mat'], 'model');

if exist('printDiagram') & printDiagram
  fgplvmPrintPlot(model, lbls, capName, experimentNo);
end

% Load the results and display dynamically.
fgplvmResultsDynamic(dataSetName, experimentNo, 'vector')

% compute the nearest neighbours errors in latent space.
errors = fgplvmNearestNeighbour(model, lbls);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function cvx_optval = square_pos( x )

%SQUARE_POS   Internal cvx version.

narginchk(1,1);
cvx_optval = pow_pos( x, 2 );

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
% pop_resample() - resample dataset (pop up window).
%
% Usage:
%   >> [OUTEEG] = pop_resample( INEEG ); % pop up interactive window
%   >> [OUTEEG] = pop_resample( INEEG, freq);
%
% Graphical interface:
%   The edit box entitled "New sampling rate" contains the frequency in
%   Hz for resampling the data. Entering a value in this box  is the same
%   as providing it in the 'freq' input from the command line.
%
% Inputs:
%   INEEG      - input dataset
%   freq       - frequency to resample (Hz)
%
% Outputs:
%   OUTEEG     - output dataset
%
% Author: Arnaud Delorme, CNL/Salk Institute, 2001
%
% Note: uses the resample() function from the signal processing toolbox
%       if present. Otherwise use griddata interpolation method (it should be
%       reprogrammed using spline interpolation for speed up).
%
% See also: resample(), eeglab()

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% 01-25-02 reformated help & license -ad
% 03-08-02 remove ica activity resampling (now set to []) -ad
% 03-08-02 debug call to function help -ad
% 04-05-02 recompute event latencies -ad

function [EEG, command] = pop_resample( EEG, freq);

command = '';
if nargin < 1
    help pop_resample;
    return;
end;
if isempty(EEG(1).data)
    disp('Pop_resample error: cannot resample empty dataset'); return;
end;

if nargin < 2

	% popup window parameters
	% -----------------------
	promptstr    = {['New sampling rate']};
	inistr       = { num2str(EEG(1).srate) };
	result       = inputdlg2( promptstr, 'Resample current dataset -- pop_resample()', 1,  inistr, 'pop_resample');
	if length(result) == 0 return; end;
	freq         = eval( result{1} );

end;

% process multiple datasets
% -------------------------
if length(EEG) > 1
    [ EEG command ] = eeg_eval( 'pop_resample', EEG, 'warning', 'on', 'params', { freq } );
    return;
end;

% finding the best ratio
[p,q] = rat(freq/EEG.srate, 0.0001); % not used right now

% set variable
% ------------
EEG.data = reshape(EEG.data, EEG.nbchan, EEG.pnts, EEG.trials);
oldpnts   = EEG.pnts;

% resample for multiple channels
% -------------------------
if isfield(EEG, 'event') & isfield(EEG.event, 'type') & isstr(EEG.event(1).type)
    tmpevent = EEG.event;
    bounds = strmatch('boundary', { tmpevent.type });
    if ~isempty(bounds),
        disp('Data break detected and taken into account for resampling');
        bounds = [ tmpevent(bounds).latency ];
        if bounds(1) < 0, bounds(1) = []; end; % remove initial boundary if any
    end;
    bounds = [1 round(bounds-0.5)+1 size(EEG.data,2)+1];
    bounds(find(bounds(2:end)-bounds(1:end-1)==0))=[]; % remove doublet boundary if any
else
    bounds = [1 size(EEG.data,2)+1]; % [1:size(EEG.data,2):size(EEG.data,2)*size(EEG.data,3)+1];
end;
if exist('resample') == 2
     usesigproc = 1;
else usesigproc = 0;
    disp('Signal Processing Toolbox absent: using custom interpolation instead of resample() function.');
    disp('This method uses cubic spline interpolation after anti-aliasing (see >> help spline)');
end;

fprintf('resampling data %3.4f Hz\n', EEG.srate*p/q);
for index1 = 1:size(EEG.data,1)
    fprintf('%d ', index1);
    sigtmp = reshape(EEG.data(index1,:, :), oldpnts, EEG.trials);

    if index1 == 1
        tmpres = [];
        indices = [1];
        for ind = 1:length(bounds)-1
            tmpres  = [ tmpres; myresample( double( sigtmp(bounds(ind):bounds(ind+1)-1,:)), p, q, usesigproc ) ];
            indices = [ indices size(tmpres,1)+1 ];
        end;
        if size(tmpres,1) == 1, EEG.pnts  = size(tmpres,2);
        else                    EEG.pnts  = size(tmpres,1);
        end;
        tmpeeglab = zeros(EEG.nbchan, EEG.pnts, EEG.trials);
    else
        for ind = 1:length(bounds)-1
            tmpres(indices(ind):indices(ind+1)-1,:) = myresample( double( sigtmp(bounds(ind):bounds(ind+1)-1,:) ), p, q, usesigproc );
        end;
    end;
    tmpeeglab(index1,:, :) = tmpres;
end;
fprintf('\n');
EEG.srate   = EEG.srate*p/q;
EEG.data = tmpeeglab;

% recompute all event latencies
% -----------------------------
if isfield(EEG.event, 'latency')
    fprintf('resampling event latencies...\n');
    for index1 = 1:length(EEG.event)
        EEG.event(index1).latency = EEG.event(index1).latency * EEG.pnts /oldpnts;
    end;
    if isfield(EEG, 'urevent') & isfield(EEG.urevent, 'latency')
        try,
            for index1 = 1:length(EEG.event)
                EEG.urevent(index1).latency = EEG.urevent(index1).latency * EEG.pnts /oldpnts;
            end;
        catch,
            disp('pop_resample warning: ''urevent'' problem, reinitializing urevents');
            EEG = rmfield(EEG, 'urevent');
        end;
    end;
    EEG = eeg_checkset(EEG, 'eventconsistency');
end;

% resample for multiple channels ica
EEG.icaact = [];

% store dataset
fprintf('resampling finished\n');

EEG.setname = [EEG.setname ' resampled'];
EEG.pnts    = size(EEG.data,2);
EEG.xmax = EEG.xmin + (EEG.pnts-1)/EEG.srate; % cko: recompute xmax, since we may have removed a few of the trailing samples

command = sprintf('EEG = pop_resample( %s, %d);', inputname(1), freq);
return;

% resample if resample is not present
% -----------------------------------
function tmpeeglab = myresample(data, pnts, new_pnts, usesigproc);

    if usesigproc
        tmpeeglab = resample(data, pnts, new_pnts);
        return;
    end;

    % anti-alias filter
    % -----------------
    data         = eegfiltfft(data', 256, 0, 128*pnts/new_pnts);

    % spline interpolation
    % --------------------
    X            = [1:length(data)];
    nbnewpoints  = length(data)*pnts/new_pnts;
    nbnewpoints2 = ceil(nbnewpoints);
    lastpointval = length(data)/nbnewpoints*nbnewpoints2;
    XX = linspace( 1, lastpointval, nbnewpoints2);

    cs = spline( X, data);
    tmpeeglab = ppval(cs, XX)';

author: PatternRecognition

repo: OpenBMI


text:
function [nautmi] = ft2nautmi(ft)
% Convert length from feet to nautical miles.
% Chad Greene 2012
nautmi = ft*0.0001645788336933;

author: Sable

repo: mcbench-benchmarks


text:
function plot_cloud(M, points, style)
    if isempty(points)
        return
    end
    if nargin==2
        style = 'r*';
    end
    if isempty(M)
        if length(style)==3
            plot3(points(:,1),points(:,2),points(:,3),'*','Color',style)
        else
            plot3(points(:,1),points(:,2),points(:,3),style)
        end
    else
        if length(style)==3
            plot3(M.VERT(points,1),M.VERT(points,2),M.VERT(points,3),'*','Color',style)
        else
            plot3(M.VERT(points,1),M.VERT(points,2),M.VERT(points,3),style)
        end
    end
    axis equal
    rotate3d on
end

author: OshriHalimi


text:
function y = sigmoid(y, fac)
% SIGMOID nonlinear funcion for cochlear model
%	y = sigmoid(y, fac);
%	fac: non-linear factor
%	 -- fac > 0, transister-like function
%	 -- fac = 0, hard-limiter
%	 -- fac = -1, half-wave rectifier
%	 -- else, no operation, i.e., linear
%
%	SIGMOID is a monotonic increasing function which simulates
%	hair cell nonlinearity.
%	See also: WAV2AUD, AUD2WAV

% Auther: Powen Ru (powen@isr.umd.edu), NSL, UMD
% v1.00: 01-Jun-97

if fac > 0,
	%y = exp(y/fac); y = 1/(1+.1)-1./(1+.1*y);
	y = exp(-y/fac); y = 1./(1+y);
elseif fac == 0,
	y = (y > 0);	% hard-limiter
elseif fac == -1,
	y = max(y, 0);	% half-wave rectifier
elseif fac == -3,
	y = halfregu(y);
end;

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function Q = deconvolutionDAMAS(S, e, maxIterations)
%deconvolutionDAMAS - deconvolves the intensity plot with the DAMAS algorithm
%as implemented in "A deconvolution approach for the mapping of acoustic sources
%(DAMAS) determined from phased microphone arrays", Brooks and Humphreys, 2005
%
%Q = deconvolutionDAMAS(S, e, maxIterations)
%
%IN
%S - MxN matrix of delay-and-sum steered response power
%e - MxNxP steering vector/matrix for a certain frequency
%
%OUT
%Q - MxN devonvolved intensity plot
%
%Created by J?rgen Grythe
%Last updated 2017-02-27

if ~exist('maxIterations', 'var')
    maxIterations = 100;
end

Y = real(S);
deps = 0.1;

%M # of y-points, N # of x-points, P number of mics
[M, N, P] = size(e);

%Make the A-matrix square size NxM x NxM x P
ee = reshape(e, M*N, P);
A = (abs(ee*ee').^2)./P^2;

%Initialise final source powers Q
Q = zeros(size(Y));
Q0 = Y;

%Solve the system Y = AQ for Q by Gauss-Seidel iteration where Y is the
%original delay-and-sum plot we want to deconvolve, and Q are the true
%source powers
for i=1:maxIterations;

    %Gauss-Seidel iteration. If the solution is negative set it to zero (to
    %ensure that we only have positive and not negative power)
    for n=1:M*N
        Q(n) = max(0, Y(n) - A(n, 1:n-1)*Q(1:n-1)' ...
            - A(n, n+1:end)*Q0(n+1:end)');
    end

    %Break criterion for convergence
    dX = (Q - Q0);
    maxd = max(abs(dX(:)))/mean(Q0(:));

    if  maxd < deps
        break;
    end

    Q0 = Q;
end

if i == maxIterations
    disp(['Stopped after maximum iterations (' num2str(maxIterations) ')'])
else
    disp(['Converged after ' num2str(i) ' iterations'])
end

author: jorgengrythe

repo: beamforming


text:
function out = isequal(f, g)
%ISEQUAL   Equality test for CHEBFUN3.
%   BOOL = ISEQUAL(F, G) returns logical 1 if F and G are same up to relative
%   machine precision. It returns logical 0 otherwise.
%
% The structure of this code is similar to `dematricize.m` from the HTUCKER
% toolbox of Tobler and Kressner.

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

if ( isempty(f) )
    if ( isempty(g) )
        out = true;
    else
        out = false;
    end
    return
end

% Get low rank representation of f:
[fCore, fCols, fRows, fTubes] = tucker(f);

% Get low rank representation of g:
[gCore, gCols, gRows, gTubes] = tucker(g);

% Test every part:
out = ( isequal(fCore, gCore) & isequal(fCols, gCols) & ...
    isequal(fRows, gRows) & isequal(fTubes, gTubes) );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ cum_ret, cumprod_ret, daily_ret, daily_portfolio] = template( fid, data, varargins, opts )
% This is a template for writing a portfolio selection algorithm
%
% [ cum_ret, cumprod_ret, daily_ret, daily_portfolio] = template( fid, data, varargins, opts )
%
% Please put the description of your algorithm here
% Name of Strategy:
% Author:
% Description:
% The sections labelled as "Static" of the file need not be changed

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors:
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %% Make changes to this section to construct your algorithm

    %% Read Parameters
    p1 = varargins{1};

    %% Initialize variables
    [r c] = size(data);
    b = ones(c,1)/c;
    returns = zeros(r,1);
    portfolio = ones(r,c)/c;

    %% Static
    progress = waitbar(0,'Executing Algorithm...');

    %% The looping over the entire dataset for backtesting
    % The algorithm looping over r time periods
    for t = 1:1:r

        %% Static
        % compute x ofthe optimization problem i.e. todaysRelative
        todaysRelative = data(t,:)';
        % Compute the returns of algorithm
        portfolio(t, :) = b;
        returns(t) = b'*(todaysRelative-1);

        %% Change this section to describe your strategy portfolio selection method

        % Use solver/algorithm to find new portfolio vector at end of time
        % period t

        %%  Static
        % Update Progress
        if mod(t, 50) == 0
            waitbar((t/r));
        end

    end

    %% Static

    % Compute additional statistics (for quick individual run without GUI)
    % Can be deleted
    stats.finalValue = prod(returns+1);
    Y =r/252;
    stats.sharpe = ((stats.finalValue)^(1/Y) - 1.04)  / (std(returns)*sqrt(252));
    stats.averageInTopStock = mean(max(portfolio') );
    stats.averageInTop2Stocks = flipud(sort(portfolio'));
    stats.averageInTop2Stocks = mean(sum(stats.averageInTop2Stocks(1:2, :)));
    stats.variance = (std(returns)*sqrt(252));

    % Conversion of results to algorithm format
    % cum_ret, cumprod_ret, daily_ret, daily_portfolio
    daily_portfolio = portfolio;
    daily_ret = returns+ 1;
    cumprod_ret = cumprod(daily_ret);
    cum_ret = cumprod_ret(end);

    close(progress);
end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


%  Script file: lsqfit.m
%
%  Purpose:
%    To perform a least-squares fit of an input data set
%    to a straight line, and print out the resulting slope
%    and intercept values.  The input data for this fit
%    comes from a user-specified input data file.
%
%  Record of revisions:
%      Date       Programmer          Description of change
%      ====       ==========          =====================
%    03/24/07    S. J. Chapman        Original code
%
% Define variables:
%   count     -- number of values read
%   filename  -- Input file name
%   fid       -- File id
%   msg       -- Open error message
%   n         -- Number of input data pairs (x,y)
%   slope     -- Slope of the line
%   sum_x     -- Sum of all input X values
%   sum_x2    -- Sum of all input X values squared
%   sum_xy    -- Sum of all input X*Y values
%   sum_y     -- Sum of all input Y values
%   x         -- An input X value
%   x_bar     -- Average X value
%   y         -- An input Y value
%   y_bar     -- Average Y value
%   y_int     -- Y-axis intercept of the line

% Initialize sums
n = 0; sum_x = 0; sum_y = 0; sum_x2 = 0; sum_xy = 0;
x1 = []; y1 = [];

% Prompt user and get the name of the input file.
disp('This program performs a least-squares fit of an');
disp('input data set to a straight line. Enter the name');
disp('of the file containing the input (x,y) pairs:  ' );
filename = input(' ','s');

% Open the input file
[fid,msg] = fopen(filename,'rt');

% Check to see if the open failed.
if fid < 0

   % There was an error--tell user.
   disp(msg);

else

   % File opened successfully. Read the (x,y) pairs from
   % the input file.  Get first (x,y) pair before the
   % loop starts.
   [in,count] = fscanf(fid,'%g %g',2);

   while ~feof(fid)
      x = in(1);
      y = in(2);
      n      = n + 1;                  %
      sum_x  = sum_x + x;              % Calculate
      sum_y  = sum_y + y;              %   statistics
      sum_x2 = sum_x2 + x.^2;          %
      sum_xy = sum_xy + x * y;         %
      x1(n) = x;
      y1(n) = y;

      % Get next (x,y) pair
      [in,count] = fscanf(fid,'%f',[1 2]);

   end

   % Close the file
   fclose(fid);

   % Now calculate the slope and intercept.
   x_bar = sum_x / n;
   y_bar = sum_y / n;
   slope = (sum_xy - sum_x*y_bar) / (sum_x2 - sum_x*x_bar);
   y_int = y_bar - slope * x_bar;

   % Tell user.
   fprintf('Regression coefficients for the least-squares line:\n');
   fprintf('   Slope (m)     = %12.3f\n',slope);
   fprintf('   Intercept (b) = %12.3f\n',y_int);
   fprintf('   No of points  = %12d\n',n);

   % Now plot this data set
   figure(1);
   plot(x1,y1,'bo','LineWidth',2);
   hold on;

   % Now plot the fitted line
   x2 = x1;
   y2 = slope .* x2 + y_int;
   plot(x2,y2,'r-','LineWidth',2);
   hold off;

   % Title and labels
   title ('\bfLeast-Squares Fit');
   xlabel ('\bf\itx');
   ylabel ('\bf\ity');
   grid on;

end


text:
function pass = test_guide( pref )
% Test Chebfun2 guide commands. This is not exclusive by checks the main
% commands.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1e3*pref.cheb2Prefs.chebfun2eps;

% commands from guide1:
f = chebfun2(@(x,y) cos(x.*y));

pass(1) = abs( sum2(f) - 3.784332281468732 ) < tol;

x = chebfun2(@(x,y) x, [-2 3 -4 4]);
y = chebfun2(@(x,y) y, [-2 3 -4 4]);

op = @(x,y) 1./( 2 + cos(.25 + x.^2.*y + y.^2) );
g = ( 2 + cos(.25 + x.^2.*y + y.^2) );
f = 1./g;

x = linspace(-2,3);
y = linspace(-4,4);
[xx,yy] = meshgrid(x,y);
pass(2) = ( norm( f(xx,yy) - op(xx,yy), inf) ) < 200*tol;

f = chebfun2(@(z) sin(z));
pass(3) = ( abs( f(1+1i) -  sin(1+1i) ) ) < tol;

f = chebfun2(@(z) sin(z)-sinh(z),2*pi*[-1 1 -1 1]);
x = linspace(-2*pi,2*pi);
[xx,yy] = meshgrid(x);
pass(4) = ( norm( f(xx+1i*yy) - f(xx,yy), inf)) < tol;
pass(5) = ( norm( f(xx+1i*yy) - sin(xx+1i*yy)+sinh(xx+1i*yy), inf)) < 1e3*tol;

% commands from guide2:
f = chebfun2(@(x,y) sin(10*x.*y), [0 pi/4 0 3]);
pass(6) = ( norm( sum(f) - chebfun(@(x) sin(15*x).^2./x/5,[0 pi/4])' )) < tol;
pass(7) = ( norm( sum(f,2) - chebfun(@(y) sin(5*pi*y/4).^2./y/5,[0 3]) )) < tol;
pass(8) = ( abs( sum2(f) - sum(sum(f)) )) < tol;

F = @(x,y) exp(-(x.^2 + y.^2 + cos(4*x.*y)));
I1 = quad2d(F,-1,1,-1,1,'AbsTol',tol);
I2 = sum(sum(chebfun2(F)));
I3 = sum2(chebfun2(F));
pass(9) = ( abs(I1 - I2)) < tol;
pass(10) = ( abs(I2 - I3)) < tol;
pass(11) = ( abs(I1 - I3)) < tol;

f = chebfun2( 'exp(-(x.^2 + y.^2 +4*x.*y))' );
pass(12) = ( abs( norm(f) - sqrt(sum2(f.^2)) ) < tol );

f = chebfun2( @(x,y) exp(-1./( sin(x.*y) + x ).^2) );
g = chebfun2( @(x,y) cos(exp(-1./( sin(x.*y) + x ).^2)) );
h = chebfun2( @(x,y) (exp(-1./( sin(x.*y) + x ).^2)).^5 );
pass(13) = ( norm( cos(f) - g ) ) < tol;
pass(14) = ( norm( f.^5 - h )) < tol;

runge = chebfun2( @(x,y) 1./( .01 + x.^2 + y.^2 )) ;
pass(15) = ( abs(mean2(runge) - 3.796119578934828)) < tol;
pass(16) = ( abs(mean(mean(runge)) - 3.796119578934828)) < tol;

f = chebfun2(@(x,y) exp(-(x.^2 + 3*x.*y+y.^2) ));
pass(17) = ( norm( cumsum(cumsum(f),2) - cumsum2(f) )) < tol;

f = chebfun2(@(x,y) cos(10*x.*y.^2) + exp(-x.^2)); % chebfun2 object
C = chebfun(@(t) t.*exp(10i*t),[0 1]);             % spiral curve

pass(18) = ( abs(sum(f(C)) - 1.613596461872283)) < tol;

% commands from guide4:
f = chebfun2(@(x,y) sin(x+1i*y));   % a holomorphic function
u = real(f); v = imag(f);           % real and imaginary parts

pass(19) = ( norm(diff(u) - (-diff(v,1,2)))) < tol;
pass(20) = ( norm(diff(u,1,2) - diff(v)) ) < tol;        % Do the Cauchy-Riemann eqns hold?

d = [0 1 0 2];
F = chebfun2v(@(x,y) sin(x.*y), @(x,y) cos(y),d);  % calling the constructor
f = chebfun2(@(x,y) sin(x.*y),d); g = chebfun2(@(x,y) cos(y),d);
G = [f;g];
plaw = abs((2*norm(F)^2 + 2*norm(G)^2) - (norm(F+G)^2 + norm(F-G)^2));
pass(21) = ( plaw < tol );

f = chebfun2(@(x,y) cos(10*x.*y.^2) + exp(-x.^2));  % chebfun2
F = gradient(f);                                    % gradient (chebfun2v)
C = chebfun(@(t) t.*exp(10i*t),[0 1]);              % spiral curve
v = integral(F,C);ends = f(cos(10),sin(10))-f(0,0); % line integral
pass(22) = ( abs(v-ends) ) < tol;

% commands from guide6:
r1 = 1; r2 = 1/3;   % inner and outer radius
d = [0 2*pi 0 2*pi];
u = chebfun2(@(u,v) u,d);
v = chebfun2(@(u,v) v,d);
F = [-(r1+r2*cos(v)).*sin(u); (r1+r2*cos(v)).*cos(u); r2*sin(v)];  % torus
G = F./3;  % full 3D divergence of G is 1 because F = [x;y;z].
pass(23) = ( abs( integral2(dot(G,normal(F))) - 2*pi^2*r1*r2.^2 ) < tol );


text:
function errs = Fsampson(F,u);
% FSAMPSON ... first order geometrical error (Sampson Distance)
% errs = Fsampson(F,u);
% F ... 3x3 Fundamental matrix
% u ... 6xN point pairs homogenous
%
% errs ... 1xN error for each point pair
%
% $Id: Fsampson.m,v 1.1 2005/05/23 16:15:59 svoboda Exp $

N = size(u,2);

u1 = u(1:3,:);
u2 = u(4:6,:);

errs = zeros(1,N);
for i=1:N
  Fu1 = F*u1(:,i);
  Fu2 = F'*u1(:,i);
  errs(i) = (u2(:,i)'*F*u1(:,i))^2 / (sum([Fu1(1:2)'.^2,Fu2(1:2)'.^2]));
end

author: strawlab

repo: MultiCamSelfCal


text:
%  Figure 3.32      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
%script to generate Fig. 3.32
clf;
k=7.5;
numCL=[k k];
denCL=[1 5 k-6 k];
t=0:.05:12;
y1=step(numCL,denCL,t);
k=13;
numCL=[k k];
denCL=[1 5 k-6 k];
y2=step(numCL,denCL,t);
k=25;
numCL=[k k];
denCL=[1 5 k-6 k];
y3=step(numCL,denCL,t);
axis([0 12 -1 3]);
plot(t,y1,t,y2,t,y3)
grid;
xlabel('Time (sec)');
ylabel('y(t)');
title('Fig. 3.32 transient responses')
axis('normal');
text(3,2,'K=7.5');
text(2.1,1.5,'K=13');
text(1.4,1.2,'K=25');


text:
function p = predictOneVsAll(all_theta, X)
%PREDICT Predict the label for a trained one-vs-all classifier. The labels
%are in the range 1..K, where K = size(all_theta, 1).
%  p = PREDICTONEVSALL(all_theta, X) will return a vector of predictions
%  for each example in the matrix X. Note that X contains the examples in
%  rows. all_theta is a matrix where the i-th row is a trained logistic
%  regression theta vector for the i-th class. You should set p to a vector
%  of values from 1..K (e.g., p = [1; 3; 1; 2] predicts classes 1, 3, 1, 2
%  for 4 examples)

m = size(X, 1);
num_labels = size(all_theta, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters (one-vs-all).
%               You should set p to a vector of predictions (from 1 to
%               num_labels).
%
% Hint: This code can be done all vectorized using the max function.
%       In particular, the max function can also return the index of the
%       max element, for more information see 'help max'. If your examples
%       are in rows, then, you can use max(A, [], 2) to obtain the max
%       for each row.
%

sig = sigmoid(X * all_theta');
[maxSig, maxSig_2] = max(sig');
p = maxSig_2';

% =========================================================================

end

repo: MachineLearning


text:
function pv = constrain_toverlap(pv,t,numframes)
% pv = constrain_toverlap(pv,t,numframes)
% constrains paramvec to be at least t elements apart

pv2 = cat(2,pv{:});

tmp=[];
for i=1:length(pv), tmp=[tmp i*ones(size(pv{i}))];,end

pv2 = [pv2; tmp];
pv2 = sortrows(pv2',1);

alltons = pv2(:,1);
tmp = [Inf; diff(pv2(:,1))];

wh = find(tmp < t);

while ~isempty(wh) & alltons(wh(1)) < numframes
    alltons(wh(1):end) = alltons(wh(1):end) + (t - tmp(wh(1)));
    tmp = [Inf; diff(alltons)];
    wh = find(tmp < t);
end

pv2 = [alltons pv2]; pv2(pv2(:,1) > numframes,:) = [];

for i=1:max(pv2(:,3))
    pv{i} = pv2(pv2(:,3)==i,1)';
end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
 function [mass_atten, kev, mtype, file] = ...
		xray_read_atten(mtype, kev_in, varargin)
%function [mass_atten, kev, mtype, file] = ...
%|		xray_read_atten(mtype, kev_in, [options])
%|
%| Read mass attenuation coefficients for a given material type.
%| Optionally interpolate onto desired energies.
%|
%| in
%|	mtype			'aluminum', 'copper', 2, '2', '02-helium', ...
%|				See xray_material_file_name.m
%|		(Optionally mtype can be a cell array of several materials.
%|		If so, kev_in is mandatory and the output will be an array.)
%|	kev_in		[N 1]	optional desired energies [in keV]
%|
%| options
%|	'units'		cm | mm	default: cm
%|	'interp'	{}	interpolator type.  default {'pchip', 'extrap'}
%|	shortfile	0|1	return short file name instead of full path
%|
%| out
%|	mass_atten	[L 1]	mass attenuation coefficients [cm^2/g],
%|				If "energy" input is provided, then [N L].
%|	kev		[N 1]	at these corresponding energies.
%|	mtype		{L}	material type
%|	file		{L}	file name(s) for each material
%|
%| Copyright 2004-05-1, Jeff Fessler, University of Michigan

% default is to show example
if nargin < 1, ir_usage, end
if nargin == 1 && streq(mtype, 'test'), xray_read_atten_test, return, end
if ~isvar('kev_in'), kev_in = []; end

if isnumeric(mtype) && length(mtype) > 1
	mtype = num2cell(mtype);
end

if iscell(mtype)
	if isempty(kev_in), error 'kev_in required', end
	mass_atten = zeros(length(kev_in), length(mtype)); % array
	for ll=1:length(mtype)
		[mass_atten(:,ll) dum mtype{ll} file{ll}] = ...
			xray_read_atten(mtype{ll}, kev_in, varargin{:});
	end
	kev = kev_in;
return
end

% defaults
arg.units = 'cm';
%arg.interp = {'linear', 'extrap'};
%arg.interp = {'spline', 'extrap'};
arg.interp = {'pchip', 'extrap'};
arg.shortfile = false;

arg = vararg_pair(arg, varargin);

is_mm = 0;
if streq(arg.units, 'mm')
	is_mm = 1;
elseif ~streq(arg.units, 'cm')
	error 'bad units'
end

file = xray_material_file_name(mtype);
tmp = load_ascii_skip_header(file); % read uncommented lines
kev = tmp(:,1) * 1000;	% keV
mass_atten = tmp(:,2);	% mass atten. coeff. [cm^2/g]

if arg.shortfile
	file = regexprep(file, '.*\/', '');
end

if 0 % look at k-edges
	jump = find(diff(kev) == 0);
	if ~isempty(jump)
		file = regexprep(file, '.*\/', '');
		jump = num2str(kev(jump)', '%5.1f');
		printf('%s %s', file, jump)
	end
end

% if input energies are specified, then interpolate onto those.
% trick: allow for the k-edge jumps!
if isvar('kev_in') && ~isempty(kev_in)
	mass_atten = log(mass_atten); % interpolate on a log scale
	mass_atten = interp1_jump(kev, mass_atten, kev_in, ...
		'monodown', arg.interp{:});
	mass_atten = exp(mass_atten);
	kev = kev_in;
end

if is_mm
	mass_atten = mass_atten * 100;
end

% xray_read_atten_test()
% example usage
function xray_read_atten_test
mtypes = {'lead', 'aluminum', 'water', 'lexan'};
mtypes = {'water', 'bone'};
mtypes = {'iodine', 'cesium', 'csi'};
arg1 = {};
kev = [10:510]';
for ll=1:length(mtypes)
	mtype = mtypes{ll};
	[mac1 kev1] = xray_read_atten(mtype);
	ie = min(kev) <= kev1 & kev1 <= max(kev);
	arg1 = {arg1{:}, kev1(ie), mac1(ie), 'o'};
end
mac2 = xray_read_atten(mtypes, kev); % test multiple types (cell)
if im
	clf, semilogy(arg1{:})
	legend(mtypes{:})
	hold on, semilogy(kev, mac2, '-'), hold off
	xlabel 'KeV', ylabel 'mass atten. coef.', axis tight
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function h = twVisualizeRFs(view, rois, dt, scans);
%
% h = twVisualizeRFs(view, [rois=all], [dt='Averages'], [scans=1:2];
%
% Visualize the location of RF esimates for each voxel in the
% specified ROIs, from a traveling wave analysis.
%
% scans: 2x1 matrix of [polar angle, eccentricity] scans, from which to
% take estimates. Retinotopic mapping parameters should be set for each of
% these (retinoSetParams; menu ColorMap | Set Retinotopy Parameters...).
%
% ras, 03/2007.
if notDefined('view'),	view = getCurView;			end
if notDefined('rois'),	rois = 1:length(view.ROIs);	end
if notDefined('dt'),	dt = 'Averages';			end
if notDefined('scans'),	scans = [1 2];				end

rois = tc_roiStruct(view, rois);

anal = twEstimateRFs(view, rois, dt, scans);

figure('Color', 'w', 'Name', ' Traveling Wave pRFs');
nrows = ceil(sqrt(length(rois)));
ncols = ceil(length(rois)/nrows);
for r = 1:length(rois)
	subplot(nrows, ncols, r);
	retinoPlot([], []);  % puts up grid
	plot(anal.x0{r}, anal.y0{r}, 'k.', 'MarkerSize', 1);
	grid on, axis([-10 10 -14 14]);
	axis square; axis equal;
	title(rois(r).name, 'FontSize', 14, 'FontName', 'Helvetica');
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% mri_superres.m
% NOT DONE!
%
% Explore super-resolution in MRI with rotated propeller sequences
% Inspired by Benoit Desjardins' work.
% Copyright 2003-7-23, Jeff Fessler, The University of Michigan

%
% create k-space samples, in units of 1/mm
%
if ~isvar('kspaces'), printm 'setup kspace'
	fov = 256;	% [mm] typical brain FOV
	N0 = 32;	% nominal image size
	kmax = N0/2*(1/fov);	% for display axes

	k1 = ([-N0/2:N0/2-1]+0.0)/fov;
	[kk1, kk2] = ndgrid(k1, k1);
	kcart = [kk1(:), kk2(:)];

	Nprops = [1 5];	% # of propellers

	for il=1:length(Nprops)
		Nprop = Nprops(il);

		kspace = zeros(Nprop*N0^2, 2);
		for ip=0:Nprop-1
			ang = ip/Nprop * pi;
			rot = [cos(ang) sin(ang); -sin(ang) cos(ang)];
			kspace(ip*N0^2+[1:N0^2],:) = kcart * rot;
			if ip > 0, % avoid repeated DC sampling
				kspace((N0/2+1)^2-1,:) = [];
			end
		end

		kspaces{il} = kspace;
	end, clear ang rot ip il k1 kk1 kk2 kcart

	im plc 2 3
	im('subplot', 1)
	if im
		plot(kspace(:,1), kspace(:,2), '.')
		axis(1.1*[-1 1 -1 1]*kmax), axis square
		xlabel 'k_1 [mm^{-1}]', ylabel 'k_2 [mm^{-1}]'
		title(sprintf('%d k-space samples (%d)', size(kspace,1), Nprop))
	end

	clear kmax kspace
prompt
end

%
% true object
%
if 0 || ~isvar('xtrue'), printm 'setup object'

	% display images with many pixels...
	Ndisp = 256;
	x1d = [-Ndisp/2:Ndisp/2-1] / Ndisp * fov;
	[x1dd x2dd] = ndgrid(x1d, x1d);

	% parameter units all in [mm]
	obj = mri_objects('case1');
	xtrue = obj.image(x1dd, x2dd);

	clim = [0 2];
	im(2, x1d, x1d, xtrue, 'x true', clim), cbar

	clear x1dd x2dd xpar
prompt
end

%
% analytical k-space data
%
if 0 || ~isvar('yd'), printm 'data'

	% noiseless data
	rng(0)
	for ik=1:length(Nprops)
		kspace = kspaces{ik};
		ytrue{ik} = obj.kspace(kspace(:,1), kspace(:,2));
		yd{ik} = ytrue{ik} + 0 * randn(size(ytrue{ik}));
	end
	clear x1dd x2dd xpar kspace ytrue

	[xhatg, yd_g, xg, kg] = mri_grid_linear(kspaces{1}, yd{1}, N0, fov);
	im(3, kg{1}, kg{2}, abs(yd_g), '|y_d|'), cbar

	clear kg yd_g xhatg ik kspace
prompt
end

return

%
% iterative recon
%
if ~isvar('xpcg')

Nlist = [2.^[6:6]];
niter = 20;
beta = 2^-10 * size(kspace,1);
xpcg = {};

for in=1:length(Nlist)
	N = Nlist(in);

	% gridding estimate to initialize iterations
	[xhatg, yhatg, xg] = mri_grid_linear(kspace, yd, N, fov);
	tmp = sprintf('|x| gridding %d', N);
	im(5, xg{1}, xg{2}, abs(xhatg), tmp, clim), cbar

	%
	% create Gnufft class object
	%
	if 1 || ~isvar('G'), printm 'setup Gnufft object'
		omega = 2*pi*kspace*fov/N;
%		minmax(omega(:,1))
		G = Gnufft({omega, [N N], [6 6], 2*[N N], [N/2 N/2]});
	end

	%
	% reconstruct by PCG
	%
	if 1 || ~isvar('xpcg'), printm 'PCG with quadratic penalty'
		mask = true(N);
		R = Robject(mask, 'beta', beta);
		ytmp = yd(:) * (1/fov)^2 * N^2;	% scaling!
		xiter = qpwls_pcg(xhatg(:), G, 1, ytmp, 0, R.C, 1, niter);
		xpcg{in} = embed(xiter(:,end), mask);
		im(6, xg{1}, xg{2}, abs(xpcg{in}), '|x| pcg quad'), cbar
	end
end % fof
end % if

return

if 0
	im plc 2 3
	im(1, x1d, x1d, xtrue, clim, 'x true'), cbar
	xtick([-fov/2 0 fov/2-1])
	ytick([-fov/2 0 fov/2-1])
	for in=1:length(Nlist)
		x1g = [-N/2:N/2-1]'/N * fov;
		tmp = sprintf('N=%d', Nlist(in));
		im(in+1, x1g, x1g, abs(xpcg{in}), tmp, clim), cbar
		xtick([-fov/2 0 fov/2-1])
		ytick([-fov/2 0 fov/2-1])
	end
%	ir_savefig fig_mri_pixel_size_image
return
end

if 1
	clf
	h(1) = plot(x1d, xtrue(:,Ndisp/2+1), 'c-');
	hold on
	llist = {'r:', 'y--', 'g-.', 'r.', 'm-'};
	arg = {'true'};
	for in=1:length(Nlist)
		N = Nlist(in);
		x1g = [-N/2:N/2-1]'/N * fov;
		tmp = xpcg{in};
		h(in+1) = plot(x1g, abs(tmp(:,N/2+1)), llist{in});
		arg = {arg{:}, sprintf('N=%d', N)};
	end
	hold off
	axis([-fov/2 fov/2 0 2.1])
	xlabel 'horizontal position [mm]'
	ylabel '|f(x,0)|'
	set(0, 'DefaultTextFontSize', 12)
%	set(0, 'DefaultAxesFontSize', 14)
	legend(h, arg{:})
%	ir_savefig fig_mri_pixel_size_profile
end


text:
function bessel_k0_int_values_test ( )

%*****************************************************************************80
%
%% BESSEL_KO_INT_VALUES_TEST demonstrates the use of BESSEL_KO_INT_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BESSEL_KO_INT_VALUES_TEST:\n' );
  fprintf ( 1, '  BESSEL_K0_INT_VALUES stores values of \n' );
  fprintf ( 1, '  the integral of the Bessel function K0.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X           FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = bessel_k0_int_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function Pvalue=myfisher23(x)
%P=MYFISHER23(X)- Fisher's Exact Probability Test on 2x3 matrix.
%Fisher's exact test of 2x3 contingency tables permits calculation of
%precise probabilities in situation where, as a consequence of small cell
%frequencies, the much more rapid normal approximation and chi-square
%calculations are liable to be inaccurate. The Fisher's exact test involves
%the computations of several factorials to obtain the probability of the
%observed and each of the more extreme tables. Factorials growth quickly,
%so it's necessary use logarithms of factorials. This computations is very
%easy in Matlab because x!=gamma(x+1) and log(x!)=gammaln(x+1). This
%function is now fully vectorized to speed up the computation.
%
% Syntax: 	myfisher23(x)
%
%     Inputs:
%           X - 2x3 data matrix
%     Outputs:
%           - Three p-values
%
%   Example:
%
%                A   B   C
%           -------------------
%      X         0   3   2
%           -------------------
%      Y         6   5   1
%           -------------------
%
%
%   x=[0 3 2; 6 5 1];
%
%   Calling on Matlab the function:
%             myfisher23(x)
%
%   Answer is:
%
% 2x3 matrix Fisher's exact test: 18 tables were evaluated
% -----------------------------------------------------------------
% 		 p-value (2-tails): 0.0882352941
% -----------------------------------------------------------------
%
%           Created by Giuseppe Cardillo
%           giuseppe.cardillo-edta@poste.it
%
% To cite this file, this would be an appropriate format:
% Cardillo G. (2007) MyFisher23: a very compact routine for Fisher's exact
% test on 2x3 matrix
% http://www.mathworks.com/matlabcentral/fileexchange/15399

%Input Error handling
if ~isequal(size(x),[2 3])
    if isequal(size(x),[3 2])
        x=x';
    else
        error('Input matrix must be a 2x3 matrix')
    end
end
if ~all(isfinite(x(:))) || ~all(isnumeric(x(:)))
    error('Warning: all X values must be numeric and finite')
end
if ~isequal(x(:),round(x(:)))
    error('Warning: X data matrix values must be whole numbers')
end

Rs=sum(x,2); %rows sum
Cs=sum(x); %columns sum
N=sum(Rs); %Total observations

%If necessed, rearrange matrix
if ~issorted(Cs)
    [Cs,ind]=sort(Cs);
    x=x(:,ind);
    clear ind
end
if ~issorted(Rs)
    [Rs,ind]=sort(Rs);
    x=x(ind,:);
    clear ind
end

%recall that Fisher's P=[Prod(Rs!)*Prod(Cs!)]/[N!*prod(X(i,j)!)]
%Log(A*B)=Log(A)+Log(B) and Log(A/B)=Log(A)-Log(B)
%Costruct all possible tables
%A 2x3 matrix has 2 degrees of freedom...
A=0:1:min(Rs(1),Cs(1)); %all possible values of X(1,1)
B=min(Cs(2),Rs(1)-A); %max value of X(1,2) given X(1,1)
et=sum(B+ones(size(B))); %tables to evaluate
Tables=zeros(et,6); %Matrix preallocation
%compute the index
stop=cumsum(B+1);
start=[1 stop(1:end-1)+1];
%In the first round of the for cycle, Column 1 assignment should be skipped
%because it is already zero. So, modify the cycle...
Tables(start(1):stop(1),2)=0:1:B(1); %Put in the Column2 all the possible values of X(1,2) given X(1,1)
for I=2:length(A)
    Tables(start(I):stop(I),1)=A(I); %replicate the A(I) value for B(I)+1 times
    %Put in the Column2 all the possible values of X(1,2) given X(1,1)
    Tables(start(I):stop(I),2)=0:1:B(I);
end
clear A B start stop
%The degrees of freedom are finished, so complete the table...
%...Put all the possible values of X(1,3) given X(1,1) and X(1,2)
Tables(:,3)=Rs(1)-sum(Tables(:,1:2),2);
%Complete the second row given the first row
Tables(:,4:6)=repmat(Cs,et,1)-Tables(:,1:3);

%Compute log(x!) using the gammaln function
zf=gammaln(Tables+1); %compute log(x!)
K=sum(gammaln([Rs' Cs]+1))-gammaln(N+1); %The costant factor K=log(prod(Rs!)*prod(Cs!)/N!)
np=exp(K-sum(zf,2)); %compute the p-value of each possible matrix
[tf,obt]=ismember(x(1,:),Tables(:,1:3),'rows'); %Find the observed table
clear zf K tf
%Finally compute the probability for 2-tailed test
P=sum(np(np<=np(obt)));

%display results
tr=repmat('-',1,65); %Set up the divisor
disp(' ')
fprintf('2x3 matrix Fisher''s exact test: %0.0f tables were evaluated\n',et)
disp(tr)
fprintf('\t\t p-value (2-tails): %0.10f\n',P);
disp(tr)
fprintf('Mid-p correction: %0.10f\n',0.5*np(obt)+sum(np(np<np(obt))));
disp(tr)

if nargout
    Pvalue=P;
end


text:
function [ a, b, c ] = sphere_triangle_vertices_to_angles ( r, v1, v2, v3 )

%*****************************************************************************80
%
%% SPHERE_TRIANGLE_VERTICES_TO_ANGLES computes the angles of a spherical triangle.
%
%  Discussion:
%
%    A sphere centered at 0 in 3D satisfies the equation:
%
%      X * X + Y * Y + Z * Z = R * R
%
%    A spherical triangle is specified by three points on the surface
%    of the sphere.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 August 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the sphere.
%
%    Input, real V1(3), V2(3), V3(3), the vertices of the triangle.
%
%    Output, real A, B, C, the angles of the spherical triangle.
%

%
%  Compute the lengths of the sides of the spherical triangle.
%
  [ as, bs, cs ] = sphere_triangle_vertices_to_sides ( r, v1, v2, v3 );
%
%  Get the spherical angles.
%
  [ a, b, c ] = sphere_triangle_sides_to_angles ( r, as, bs, cs );

  return
end


text:
function [PV,DV] = VariableClustering(Problem,Population,nSel,nPer)
% Detect the kind of each decision variable

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    [N,D] = size(Population.decs);
    ND    = NDSort(Population.objs,1) == 1;
    fmin  = min(Population(ND).objs,[],1);
    fmax  = max(Population(ND).objs,[],1);
    if any(fmax==fmin)
        fmax = ones(size(fmax));
        fmin = zeros(size(fmin));
    end

    %% Calculate the proper values of each decision variable
    Angle  = zeros(D,nSel);
    RMSE   = zeros(D,nSel);
    Sample = randi(N,1,nSel);
    for i = 1 : D
        drawnow('limitrate');
        % Generate several random solutions by perturbing the i-th dimension
        Decs      = repmat(Population(Sample).decs,nPer,1);
        Decs(:,i) = unifrnd(Problem.lower(i),Problem.upper(i),size(Decs,1),1);
        newPopu   = Problem.Evaluation(Decs);
        for j = 1 : nSel
            % Normalize the objective values of the current perturbed solutions
            Points = newPopu(j:nSel:end).objs;
            Points = (Points-repmat(fmin,size(Points,1),1))./repmat(fmax-fmin,size(Points,1),1);
            Points = Points - repmat(mean(Points,1),nPer,1);
            % Calculate the direction vector of the determining line
            [~,~,V] = svd(Points);
            Vector  = V(:,1)'./norm(V(:,1)');
            % Calculate the root mean square error
            error = zeros(1,nPer);
            for k = 1 : nPer
                error(k) = norm(Points(k,:)-sum(Points(k,:).*Vector)*Vector);
            end
            RMSE(i,j) = sqrt(sum(error.^2));
            % Calculate the angle between the line and the hyperplane
            normal     = ones(1,size(Vector,2));
            sine       = abs(sum(Vector.*normal,2))./norm(Vector)./norm(normal);
            Angle(i,j) = real(asin(sine)/pi*180);
        end
    end

    %% Detect the kind of each decision variable
    VariableKind = (mean(RMSE,2)<1e-2)';
    result       = kmeans(Angle,2)';
    if any(result(VariableKind)==1) && any(result(VariableKind)==2)
        if mean(mean(Angle(result==1&VariableKind,:))) > mean(mean(Angle(result==2&VariableKind,:)))
            VariableKind = VariableKind & result==1;
        else
            VariableKind = VariableKind & result==2;
        end
    end
    PV = find(~VariableKind);
    DV = find(VariableKind);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% FDLA and FMMC solutions for an 8-node, 13-edge graph
% S. Boyd, et. al., "Convex Optimization of Graph Laplacian Eigenvalues"
% ICM'06 talk examples (www.stanford.edu/~boyd/cvx_opt_graph_lapl_eigs.html)
% Written for CVX by Almir Mutapcic 08/29/06
% (figures are generated)
%
% In this example we consider a graph described by the incidence matrix A.
% Each edge has a weight W_i, and we optimize various functions of the
% edge weights as described in the referenced paper; in particular,
%
% - the fastest distributed linear averaging (FDLA) problem (fdla.m)
% - the fastest mixing Markov chain (FMMC) problem (fmmc.m)
%
% Then we compare these solutions to the heuristics listed below:
%
% - maximum-degree heuristic (max_deg.m)
% - constant weights that yield fastest averaging (best_const.m)
% - Metropolis-Hastings heuristic (mh.m)

% small example (incidence matrix A)
A = [ 1  0  0  1  0  0  0  0  0  0  0  0  0;
     -1  1  0  0  1  1  0  0  0  0  0  0  1;
      0 -1  1  0  0  0  0  0 -1  0  0  0  0;
      0  0 -1  0  0 -1  0  0  0 -1  0  0  0;
      0  0  0 -1  0  0 -1  1  0  0  0  0  0;
      0  0  0  0  0  0  1  0  0  0  1  0  0;
      0  0  0  0  0  0  0 -1  1  0 -1  1 -1;
      0  0  0  0 -1  0  0  0  0  1  0 -1  0];

% x and y locations of the graph nodes
xy = [ 1 2   3 3 1 1 2   3 ; ...
       3 2.5 3 2 2 1 1.5 1 ]';

% Compute edge weights: some optimal, some based on heuristics
[n,m] = size(A);

[ w_fdla, rho_fdla ] = fdla(A);
[ w_fmmc, rho_fmmc ] = fmmc(A);
[ w_md,   rho_md   ] = max_deg(A);
[ w_bc,   rho_bc   ] = best_const(A);
[ w_mh,   rho_mh   ] = mh(A);

tau_fdla = 1/log(1/rho_fdla);
tau_fmmc = 1/log(1/rho_fmmc);
tau_md   = 1/log(1/rho_md);
tau_bc   = 1/log(1/rho_bc);
tau_mh   = 1/log(1/rho_mh);

fprintf(1,'\nResults:\n');
fprintf(1,'FDLA weights:\t\t rho = %5.4f \t tau = %5.4f\n',rho_fdla,tau_fdla);
fprintf(1,'FMMC weights:\t\t rho = %5.4f \t tau = %5.4f\n',rho_fmmc,tau_fmmc);
fprintf(1,'M-H weights:\t\t rho = %5.4f \t tau = %5.4f\n',rho_mh,tau_mh);
fprintf(1,'MAX_DEG weights:\t rho = %5.4f \t tau = %5.4f\n',rho_md,tau_md);
fprintf(1,'BEST_CONST weights:\t rho = %5.4f \t tau = %5.4f\n',rho_bc,tau_bc);

% Plot results
figure(1), clf
plotgraph(A,xy,w_fdla);
text(0.55,1.05,'FDLA optimal weights')

figure(2), clf
plotgraph(A,xy,w_fmmc);
text(0.55,1.05,'FMMC optimal weights')

figure(3), clf
plotgraph(A,xy,w_md);
text(0.5,1.05,'Max degree optimal weights')

figure(4), clf
plotgraph(A,xy,w_bc);
text(0.5,1.05,'Best constant optimal weights')

figure(5), clf
plotgraph(A,xy,w_mh);
text(0.46,1.05,'Metropolis-Hastings optimal weights')


text:
% InitializationforGUI_step1
% Processing nested data: one-time initialization for GUI

clear all; close all

behav_fluo_shift = -2; % shift correction between fictive and fluo, manually chosen

%% Load data
data_masterdir = GetNestedDataDir();
range_fish = GetFishRange();

for i_fish = range_fish,

    disp(['fish ' num2str(i_fish) ' loading...']);
    tic
    data_dir = fullfile(data_masterdir,['subject_' num2str(i_fish)]);

    TimeSeries = h5read(fullfile(data_dir,'TimeSeries.h5'),'/TimeSeries');

    filename = fullfile(data_dir,'CoreInfo.mat');
    load(filename);% {'periods','timelists','timelists_names','stimuluskey_raw','CellXYZ','anat_stack','fpsec'};

    filename = fullfile(data_dir,'OptionalInfo.mat');
    load(filename); % {'Behavior_raw','numcell_full','CellXYZ_norm','IX_inval'};

    filename = fullfile(data_dir,'AdditionalInfo.mat');
    load(filename); % {'frame_keys','IX_inval_anat'} and if exist, 'stimset';
    toc

    %% generate absolute cell index
    % (this indexing is stable; any later cell selection would be a subset of this)
    absIX_full = (1:numcell_full)';
    absIX = setdiff(absIX_full,IX_inval_anat);

    %% Anatomy stack projections
    % x-y view
    im = max(anat_stack,[],3);
    out=imNormalize99(im);
    anat_yx = repmat(out,[1 1 3]);

    % y-z view
    im = squeeze(max(anat_stack,[],2));
    out=imNormalize99(im);
    anat_yz = repmat(out,[1 1 3]);

    % x-z view
    im = squeeze(max(anat_stack,[],1));
    out=imNormalize99(im);
    out = flipud(out'); %%%% empirically necessary...
    anat_zx = repmat(out,[1 1 3]);

%     dimv_yx = size(anat_yx);
%     dimv_yz = size(anat_yz);
%     dimv_zx = size(anat_zx);

    %% get full-length time-index from timelists
    if length(timelists_names)==1, % fish with single stimulus set
        IX_all_raw = timelists_raw;
        % get relative time-index (relative to full-length list), ~ frame number
        timelists = {1:length(IX_all_raw)};

        % stim & behavior
        stim_full = stimuluskey_raw(IX_all_raw);
        Behavior_full = Behavior_raw(:,IX_all_raw);

        % trimming & shift
        CellResp_full = TimeSeries(:,IX_all_raw);
        CellResp_full = circshift(CellResp_full,behav_fluo_shift,2);

        % averages
        period = periods;
        numcell = size(CellResp_full,1);
        CellRespAvr_full = mean(reshape(CellResp_full,numcell,period,[]),3);
        stimAvr = stim_full(1:period);
        BehaviorAvr = mean(reshape(Behavior_full,size(Behavior_full,1),period,[]),3);

    else % fish with multiple stimulus sets

        temp = cat(2, timelists_raw{:});
        IX_all_raw = sort(temp);

        % get relative time-index (relative to full-length list), ~ frame number
        nTypes = length(timelists_names);
        timelists = cell(1,nTypes);
        for i_type = 1:nTypes,
            [~,IX] = ismember(timelists_raw{i_type},IX_all_raw);
            timelists{i_type} = IX(find(IX));
        end

        % stim & behavior
        stim_full = stimuluskey_raw(IX_all_raw);
        Behavior_full = Behavior_raw(:,IX_all_raw);

        % trimming & shift
        CellResp_full = TimeSeries(:,IX_all_raw);
        CellResp_full = circshift(CellResp_full,behav_fluo_shift,2); % using circshift as a shortcut for padding on one end

        % averages
        numcell = size(CellResp_full,1);
        CellRespAvr_full = [];
        stimAvr = [];
        BehaviorAvr = [];
        for i = 1:nTypes,
            M = CellResp_full(:,timelists{i});
            if mod(numel(M),numcell*periods(i))==0,
                avr = mean(reshape(M,numcell,periods(i),[]),3);
            else % for patterns with dummy periods where the stimulus is constant, like 'spont'
                nrep = floor(numel(M)/(numcell*periods(i)));
                M2 = M(:,1:nrep*periods(i));
                avr = mean(reshape(M2,numcell,periods(i),[]),3);
            end
            CellRespAvr_full = horzcat(CellRespAvr_full,avr); %#ok<AGROW>

            m = stim_full(timelists{i});
            stimAvr = horzcat(stimAvr,m(1:periods(i))); %#ok<AGROW>

            M_behav = Behavior_full(:,timelists{i});
            if mod(numel(M_behav),size(M_behav,1)*periods(i))==0,
                avr = mean(reshape(M_behav,size(M_behav,1),periods(i),[]),3);
            else % for patterns with dummy periods where the stimulus is constant, like 'spont'
                nrep = floor(numel(M_behav)/(size(M_behav,1)*periods(i)));
                M2_behav = M_behav(:,1:nrep*periods(i));
                avr = mean(reshape(M2_behav,size(M_behav,1),periods(i),[]),3);
            end
            BehaviorAvr = horzcat(BehaviorAvr,avr); %#ok<AGROW>
        end
    end

    %% Behavior normalizations
    for i = 1:3,
        m = BehaviorAvr(i,:);
        BehaviorAvr(i,:) = (m-min(m))/(max(m)-min(m));
        m = Behavior_full(i,:);
        Behavior_full(i,:) = (m-min(m))/(max(m)-min(m));
    end
    m = BehaviorAvr(4:5,:);
    BehaviorAvr(4:5,:) =  (m-min(min(m)))/(max(max(m))-min(min(m)));
    m = Behavior_full(4:5,:);
    Behavior_full(4:5,:) =  (m-min(min(m)))/(max(max(m))-min(min(m)));

    %% compute z-score
    disp('compute z-score...')
    tic
    CellRespZ_full = zscore(CellResp_full')';
    CellRespAvrZ_full = zscore(CellRespAvr_full')';
    toc

    %% save
    disp('saving to files...')
    tic
    % directory to save data formatted for distribution:
    save_masterdir = GetCurrentDataDir();
    save_dir = fullfile(save_masterdir,['subject_' num2str(i_fish)]);
    if ~exist(save_dir, 'dir')
        mkdir(save_dir);
    end

    % save time-series
    filename = fullfile(save_dir,'TimeSeries.h5');
    h5create(filename,'/CellResp',size(CellResp_full(absIX,:)),'Datatype','single','ChunkSize',[1000 100]);
    h5write(filename,'/CellResp',CellResp_full(absIX,:));

    h5create(filename,'/CellRespZ',size(CellRespZ_full(absIX,:)),'Datatype','single','ChunkSize',[1000 100]);
    h5write(filename,'/CellRespZ',CellRespZ_full(absIX,:));

    h5create(filename,'/CellRespAvr',size(CellRespAvr_full(absIX,:)),'Datatype','single','ChunkSize',[1000 100]);
    h5write(filename,'/CellRespAvr',CellRespAvr_full(absIX,:));

    h5create(filename,'/CellRespAvrZ',size(CellRespAvrZ_full(absIX,:)),'Datatype','single','ChunkSize',[1000 100]);
    h5write(filename,'/CellRespAvrZ',CellRespAvrZ_full(absIX,:));

    h5create(filename,'/absIX',size(absIX),'Datatype','single');
    h5write(filename,'/absIX',absIX);
    toc
    %%
    data = [];
    names = {'periods','timelists_names','stimuluskey_raw','CellXYZ','anat_stack','fpsec',...
        'Behavior_raw','numcell_full','CellXYZ_norm','IX_inval_anat',...
        'anat_yx','anat_yz','anat_zx',...
        'timelists','stim_full','stimAvr','Behavior_full','BehaviorAvr'};
%         'absIX'}; % absIX now stored in hdf5
    if length(timelists_names)>1, % M_stimset(i_fish) > 1,
        names = [names,{'stimset'}]; %#ok<AGROW>
    end

    for i = 1:length(names), % use loop to save variables into fields of 'data'
        eval(['data.',names{i},'=', names{i},';']);
    end

    save(fullfile(save_dir,'data_full.mat'),'data');

    toc

end

author: xiuyechen

repo: FishExplorer


text:
function dtiInitReorientBvecs(dwParams,dwDir,doResamp,doBvecs,bvecs,bvals)
%
%  doRawResample = dtiInitResample(dwParams,dwDir)
%
% Based on user selected params and the outputs of dtiInitAlign and
% dtiInitEddyCC we reorient and save bvecs and bvals files
%
% WEB resources:
%   http://white.stanford.edu/newlm/index.php/DTI_Preprocessing
%   mrvBrowseSVN('dtiInitReorientBvecs');
%
% (C) Stanford VISTA, 8/2011 [lmp]
%

%%
%
% If bvecs file was computed earlier doBvecs will = true. If clobber = true
% or the aligned BVs files do not exist we reorient and save.
if doBvecs || (doResamp && dwParams.clobber == -1) ...
           || dwParams.clobber == 1 ...
           || ~exist(dwDir.alignedBvecsFile,'file') ...
           || ~exist(dwDir.alignedBvalsFile,'file')
    % Reorient and save the b-vectors and bvals
    dtiRawReorientBvecs(bvecs, dwDir.ecFile, dwDir.acpcFile, dwDir.alignedBvecsFile);
    dlmwrite(dwDir.alignedBvalsFile,bvals,' ');
else
    % If clobber = 'ask' prompt the user to overwrite
    if dwParams.clobber == 0
        resp = questdlg([dwDir.alignedBvecsFile ' exists- would you like to overwrite it?'],...
            'Clobber Bvecs/bvals', 'Overwrite','Use Existing File','Abort','Use Existing File');
        if(strcmpi(resp,'Abort')), error('User aborted.'); end
        if(strcmpi(resp,'Overwrite'))
            % If the user selects to overwrite we reorient and save the
            % b-vectors and bvals
            dtiRawReorientBvecs(dwDir.bvecsFile, dwDir.ecFile, dwDir.acpcFile, dwDir.alignedBvecsFile);
            dlmwrite(dwDir.alignedBvalsFile,bvals,' ');
        end
    end
end


text:
function err = getL2error3ND2(node,elem,exactE,Eh,markedElem)
%% GETL2ERROR3ND2 L2 norm of the quadratic (1st type) Nedelec edge element.
%
%  error = getL23errorND2(node,elem,exactE,Eh,markedElem);
%
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Sort elem to ascend ordering
elem = sort(elem,2);

%% Construct Data Structure
[Dlambda,volume] = gradbasis3(node,elem);
% elem2dof
[elem2edge,edge] = dof3edge(elem);
[elem2face,face] = dof3face(elem);
NT = size(elem,1);  NE = size(edge,1); NF = size(face,1);
elem2dof = [elem2edge elem2edge+NE elem2face+2*NE elem2face+2*NE+NF];
% local indices
locBasesIdx = [1 2 0; 1 3 0; 1 4 0; 2 3 0; 2 4 0; 3 4 0; ... % phi
               1 2 0; 1 3 0; 1 4 0; 2 3 0; 2 4 0; 3 4 0; ... % psi
               3 2 4; 3 1 4; 2 1 4; 2 1 3; ...
               4 2 3; 4 1 3; 4 1 2; 3 1 2]; % chi
%% Compute square of the L2 error
err = zeros(NT,1);
[lambda,w] = quadpts3(4); % quadrature order is 4
nQuad = size(lambda,1);
for p = 1:nQuad
    % quadrature points in the x-y-z coordinate
    pxy = lambda(p,1)*node(elem(:,1),:) ...
        + lambda(p,2)*node(elem(:,2),:) ...
        + lambda(p,3)*node(elem(:,3),:) ...
        + lambda(p,4)*node(elem(:,4),:);
    Ep = exactE(pxy);
    Ehp = zeros(NT,3);
    % compute Ehp at quadrature points
    for k = 1:20
        k1 = locBasesIdx(k,1);
        k2 = locBasesIdx(k,2);
        k3 = locBasesIdx(k,3);
        % evaluate basis at quadrature points
        if k<=6
        % phi_k = lambda_{k1}Dlambda_{k2} - lambda_{k2}Dlambda_{k1};
            basis_k = (lambda(p,k1)*Dlambda(:,:,k2) ...
                      -lambda(p,k2)*Dlambda(:,:,k1));
        elseif k<=12
        % phi_k = lambda_{k1}Dlambda_{k2} + lambda_{k2}Dlambda_{k1};
            basis_k = (lambda(p,k1)*Dlambda(:,:,k2) ...
                      +lambda(p,k2)*Dlambda(:,:,k1));
        else
        % chi_k = lambda_{k1}phi_{k2,k3};
            basis_k = lambda(p,k1)*(lambda(p,k2)*Dlambda(:,:,k3) ...
                                   -lambda(p,k3)*Dlambda(:,:,k2));
        end
        Ehp = Ehp + repmat(Eh(elem2dof(:,k)),1,3).*basis_k;
    end
    err = err + w(p)*volume.*sum((Ep - Ehp).^2,2);
end
% modify the error
err(isnan(err)) = 0; % remove the singular part
if (nargin == 5) && ~isempty(markedElem)
    err = err(markedElem); % L2 error on some marked region
end
err = sqrt(sum(err));
%% TODO add example in M-lint

repo: ifem


text:
classdef TestAGAST
    %TestAGAST

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','tsukuba_l.png');
        kfields = {'pt', 'size', 'angle', 'response', 'octave', 'class_id'};
    end

    methods (Static)
        function test_detect_img
            img = imread(TestAGAST.im);
            kpts = cv.AGAST(img, ...
                'Threshold',10, 'NonmaxSuppression',true);
            validateattributes(kpts, {'struct'}, {'vector'});
            assert(all(ismember(TestAGAST.kfields, fieldnames(kpts))));
        end

        function test_detect_types
            types = {'AGAST_5_8', 'AGAST_7_12d', 'AGAST_7_12s', 'OAST_9_16'};
            img = imread(TestAGAST.im);
            for i=1:numel(types)
                kpts = cv.AGAST(img, 'Type',types{i});
                validateattributes(kpts, {'struct'}, {'vector'});
            end
        end

        function test_error_argnum
            try
                cv.AGAST();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv


text:
  function [fw, hr, hl] = fwhm1(psfs, varargin)
%|function [fw, hr, hl] = fwhm1(psfs, options)
%|
%| compute fwhm of point-spread function centered at pixel imid
%| and half-right width and half-left width (fw = hr + hl)
%|
%| in
%|	psfs	[np nc]	column psf(s)
%|
%| option
%|	imid	[1]	which pixel is the middle (default: peak)
%|	dx	[1]	pixel size (default: 1)
%|	min0	0|1	if 1, then negative psf value set to 0 (default: 1)
%|	chat	0|1	if 1, then plot
%| out
%|	fw	[nc]	fwhm of each psf column
%|
%| 2012-09-13, added 'min0' option to avoid fwhm < 1

if ~nargin, ir_usage, end
if streq(psfs, 'test'), fwhm1_test, return, end

arg.chat = 0;
arg.dx = 1;
arg.imid = [];
arg.min0 = true;
arg = vararg_pair(arg, varargin);

psfs = squeeze(psfs); % in case [1 1 np]
[np nc] = size(psfs);
if (np == 1) % single row
	psfs = psfs';
	[np nc] = size(psfs);
end

warned = false;
for ic = 1:nc
	psf = psfs(:,ic);
	if isempty(arg.imid)
		imid = imax(psf);
	else
		imid = arg.imid;
	end

	% normalize
	psf = psf / psf(imid);
	if ~warned && (1 ~= max(psf))
		warn 'peak not at center'
		warned = true;
	end

	% right
	ir = sum(cumprod(double(psf((imid+1):np) >= 0.5)));
	if (imid + ir == np)
		hr(ic,1) = ir;
	else
		high = psf(imid + ir);
		ilow = imid + ir + 1;
		low = psf(ilow);
		if arg.min0 && low < 0
			warn('psf(%d) = %g < 0, so using 0', ilow, low)
			low = 0;
		end
		hr(ic,1) = ir + (high - 1/2) / (high-low);
	end

	% left
	il = sum(cumprod(double(psf((imid-1):-1:1) >= 0.5)));
	if (il == imid-1)
		hl(ic,1) = il;
	else
		high = psf(imid - il);
		ilow = imid - il - 1;
		low = psf(ilow);
		if arg.min0 && low < 0
			warn('psf(%d) = %g < 0, so using 0', ilow, low)
			low = 0;
		end
		hl(ic,1) = il + (high - 1/2) / (high-low);
	end
end

hr = hr * arg.dx;
hl = hl * arg.dx;
fw = hr + hl;

if arg.chat && im
	plot(([1:np]-imid)*arg.dx, psf, '-o', ...
	[-hl hr], [0.5 0.5], '-')
	xlabel 'x', ylabel 'psf(x)'
	title(sprintf('fwhm=%g', fw))
end

% fwhm1_test
function fwhm1_test
dx = 3;
nx = 100;
xx = [-nx/2:nx/2-1]' * dx;
fx = 30;
sx = fx / sqrt(log(256));
psf = exp(-((xx/sx).^2)/2);
fw = fwhm1(psf, 'dx', dx, 'chat', 1);

psf = zeros(7,1); psf(5) = 2; psf([4 6]) = -0.2;
fw = fwhm1(psf, 'dx', dx, 'chat', 1);