text:
function d = commandlines(d, commandname, varargin)

%tstoolbox/@description/addcommandlines
%   adds new commandline to list of commands that have been applied to
%   that signal
%   example 1
%          addcommandlines(s, 's = spec2(s', 512, 'Hanning' )) will add 's
%          = spec2(s, 512, 'Hanning');' to the list of applied commands
%   example 2
%          len = 512; text = 'Hanning'; addcommandlines(s, 's = spec2(s',
%          512, 'Hanning' )) will add 's = spec2(s, 512, 'Hanning');' to
%          the list of applied commands
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

commandline = commandname;

for i = 1:nargin-2
	commandline = [commandline ', ' tostring(varargin{i})];
end

commandline = [commandline ');'];
d.commandlines = TSTL_append(d.commandlines, commandline);

function out = tostring(in)
switch class(in)
	case 'char'
		out = ['''' in ''''];
	case 'double'
		out = num2str(in);
	otherwise
		error('Wrong type of argument(s) given');
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
%SF_LINE_H3 Third order 1D C1 Hermite shape functions for lines.
%
%   [ VBASE, NLDOF, XLDOF, SFUN ] = SF_LINE_H3( I_EVAL, N_SDIM, N_VERT, I_DOF, XI, AINVJAC, VBASE )
%   Evaluates C1 Hermite shape functions on 1D line elements with value and
%   first derivatives defined in the nodes. XI are Barycentric coordinates.
%
%       Input       Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       i_eval      scalar:  1             Evaluate function values
%                           >1             Evaluate values of derivatives
%       n_sdim      scalar: 1              Number of space dimensions
%       n_vert      scalar: 2              Number of vertices per cell
%       i_dof       scalar: 1-4            Local basis function to evaluate
%       xi          array  [2,1]           Local coordinates of evaluation point
%       aInvJac     [n,3]                  Inverse of transformation Jacobian
%       vBase       [n]                    Preallocated output vector
%                                                                                .
%       Output      Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       vBase       [n]                    Evaluated function values
%       nLDof       [2,4]                  Number of local degrees of freedom on
%                                          vertices, edges, faces, cell interiors,
%                                          and vertices without boundary conditions
%       xLDof       [2,n_ldof]             Local coordinates of local dofs
%       sfun        string                 Function name of called shape function
%
%   See also SF_LINE_P3

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
classdef ElementSetEnum < matlab.mixin.SetGet
    %ElementSetEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        CartesianElements('Cartesian Elements', ...
                         {'Rx','Ry','Rz','Vx','Vy','Vz'}, ...
                         {'km','km','km','km/s','km/s','km/s'});
        KeplerianElements('Keplerian Orbit Elements', ...
                         {'SMA','Eccentricity','Inclination','RAAN','Arg Peri','True Aomaly'}, ...
                         {'km','','deg','deg','deg','deg'});
        GeographicElements('Geographical Elements', ...
                          {'Latitude','Longitude','Altitude','Velocity Az','Velocity El','Velocity Mag'}, ...
                          {'degN','degE','km','deg','deg','km/s'});
        UniversalElements('Universal Orbit Elements', ...
                         {'C3','Rp','Inclination','RAAN','Arg Peri','Time Past Peri.'}, ...
                         {'km^2/s^2','km','deg','deg','deg','sec'});
    end

    properties
        name(1,:) char
        elemNames(6,1) cell
        unitNames(6,1) cell
    end

    methods
        function obj = ElementSetEnum(name, elemNames, unitNames)
            obj.name = name;
            obj.elemNames = elemNames;
            obj.unitNames = unitNames;
        end
    end

    methods(Static)
        function listBoxStr = getListBoxStr()
            m = enumeration('ElementSetEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('ElementSetEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('ElementSetEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% Fig. 6.48   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%

clear all;
%close all;
clf

num=[20 1];
den=[1 0 0];
w=logspace(-3,1,100);
[m,p]=bode(num,den,w);
loglog(w,m);
xlabel('\omega (rad/sec)');
ylabel('Magnitude');
title('Fig. 6.48 Compensated open-loop transfer function.');
bodegrid;

author: Sable

repo: mcbench-benchmarks


text:
function stla_to_tri_surface ( stla_file_name, node_file_name, triangle_file_name )

%*****************************************************************************80
%
%% STLA_TO_TRI_SURFACE converts an ASCII STL file to a TRI_SURFACE file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string STLA_FILE_NAME, the name of the ASCII STL file to be read.
%
%    Input, string NODE_FILE_NAME, the name of the TRI_SURFACE node file
%    to be written.
%
%    Input, string TRIANGLE_FILE_NAME, the name of the TRI_SURFACE triangle
%    file to be written.
%

%
%  Get the size of the STLA data.
%
  [ solid_num, node_num, triangle_num, text_num ] = stla_size ( stla_file_name );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  STLA file "%s".\n', stla_file_name );
  fprintf ( 1, '  Number of solids   = %d\n', solid_num );
  fprintf ( 1, '  Number of nodes    = %d\n', node_num );
  fprintf ( 1, '  Number of faces    = %d\n', triangle_num );
  fprintf ( 1, '  Number of textures = %d\n', text_num );
%
%  Read the STLA data.
%
  [ node_xyz, triangle_node, triangle_normal ] = stla_read ( stla_file_name, ...
    node_num, triangle_num );
%
%  Write the data to the TRI_SURFACE files.
%
  dim_num = size ( node_xyz, 1 );
  order_num = 3;

  tri_surface_write ( node_file_name, triangle_file_name, dim_num, node_num, ...
    order_num, triangle_num, node_xyz, triangle_node );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  TRI_SURFACE data written to:\n' );
  fprintf ( 1, '  node file     "%s".\n', node_file_name );
  fprintf ( 1, '  triangle_file "%s".\n', triangle_file_name );

  return
end
function c2 = ch_cap ( c )

%*****************************************************************************80
%
%% CH_CAP capitalizes a single character.
%
%  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.
%
%  Examples:
%
%    CH_EQI ( 'A', 'a' ) is TRUE.
%
%  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.
%
%  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
%
%  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 i4mat_write0 ( output_filename, m, n, table )

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

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

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

  return
end
function r8mat_write0 ( output_filename, m, n, table )

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

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

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

  return
end
function 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.
%
%  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.
%
%  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 len = s_len_trim ( s )

%*****************************************************************************80
%
% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LENGTH, 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.
%
%  Examples:
%
%    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)
%
%  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 ] = stla_read ( ...
  input_file_name, node_num, face_num )

%*****************************************************************************80
%
%% STLA_READ reads graphics information from an ASCII StereoLithography file.
%
%  Modified:
%
%    17 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 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.
%
  node_xyz = [];
  face_node = [];
  face_normal = [];
%
%  If no file input, try to get one from the user.
%
  if ( nargin < 1 )
    input_file_name = input ( 'Enter the name of the ASCII STL file.' );
    if ( isempty ( input_file_name ) )
      return
    end
  end
%
%  If no size input, scan the file.
%
  if ( nargin < 3 )
    [ solid_num, node_num, face_num, text_num ] = stla_size ( ...
      input_file_name );
  end
%
%  Open the file.
%
  input_file_unit = fopen ( input_file_name, 'r' );

  if ( input_file_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'STLA_READ - Fatal error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'STLA_READ - Fatal error!' )
    return
  end
%
%  Allocate space.
%
  state = 0;
  text_num = 0;

  solid = 0;
  face = 0;
  node = 0;

  node_xyz = zeros ( 3, node_num );
  face_node = zeros ( 3, face_num );
  face_normal = zeros ( 3, face_num );
%
%  Read the next line of text.
%
  while ( 1 )

    text = fgetl ( input_file_unit );
%
%  On end of input, FGETL returns a -1 instead of a character array.
%
    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' );
        error ( 'STLA_READ - Fatal error!' );
        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 )
      break
    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' ) )
        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' );
        error ( 'STLA_READ - Fatal error!' );
        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' );
        error ( 'STLA_READ - Fatal error!' );
        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' );
        error ( 'STLA_READ - Fatal error!' );
        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' );
        error ( 'STLA_READ - Fatal error!' );
        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' );
        error ( 'STLA_READ - Fatal error!' );
        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' );
        error ( 'STLA_READ - Fatal error!' );
        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 );
        error ( 'STLA_READ - Fatal error!' );
        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' );
          error ( 'STLA_READ - Fatal error!' );
          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' );
          error ( 'STLA_READ - Fatal error!' );
          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' );
        error ( 'STLA_READ - Fatal error!' );
        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' );
        error ( 'STLA_READ - Fatal error!' );
        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' );
        error ( 'STLA_READ - Fatal error!' );
        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' );
        error ( 'STLA_READ - Fatal error!' );
        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' );
        error ( 'STLA_READ - Fatal error!' );
        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' );
          error ( 'STLA_READ - Fatal error!' );
          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' );
          error ( 'STLA_READ - Fatal error!' );
          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 );
        error ( 'STLA_READ - Fatal error!' );
        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' );
      error ( 'STLA_READ - Fatal error!' );
      return

    end

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

  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.
%
%  Modified:
%
%    16 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;

  solid_num = -1;
  node_num = -1;
  face_num = -1;
  text_num = -1;
%
%  If no file input, try to get one from the user.
%
  if ( nargin < 1 )
    input_file_name = input ( 'Enter the name of the ASCII STL file.' );
    if ( isempty ( input_file_name ) )
      return
    end
  end
%
%  Open the file.
%
  input_file_unit = fopen ( input_file_name, 'r' );

  if ( input_file_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'STLA_SIZE - Fatal error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'STLA_SIZE - Fatal error!' );
    return
  end

  solid_num = 0;
  node_num = 0;
  face_num = 0;
  text_num = 0;
%
%  Read the text, one line at a time.
%
  while ( 1 )

    text = fgetl ( input_file_unit );

    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 )
      break
    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 ( input_file_unit );

  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 tri_surface_write ( node_file_name, triangle_file_name, dim_num, ...
  node_num, order_num, triangle_num, node_xyz, triangle_node )

%*****************************************************************************80
%
%% TRI_SURFACE_WRITE writes graphics information to a pair of TRI_SURFACE files.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string NODE_FILE_NAME, the name of the node file.
%
%    Input, string TRIANGLE_FILE_NAME, the name of the triangle file.
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer NODE_NUM, the number of points.
%
%    Input, integer ORDER_NUM, the order of the triangles.
%
%    Input, integer TRIANGLE_NUM, the number of triangles.
%
%    Input, real NODE_XYZ(DIM_NUM,NODE_NUM), the node coordinates.
%
%    Input, integer TRIANGLE_NODE(ORDER_NUM,TRIANGLE_NUM),
%    the nodes that form the triangles.
%
  r8mat_write0 ( node_file_name, dim_num, node_num, node_xyz );

  i4mat_write0 ( triangle_file_name, order_num, triangle_num, ...
    triangle_node );

  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.
%
%  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

author: johannesgerer

repo: jburkardt-m


text:
% This script compiles Matlab's MEX interfaces for LIBQP solvers.
%

mex libqp_splx_mex.c -largeArrayDims ../lib/libqp_splx.c
mex libqp_gsmo_mex.c -largeArrayDims ../lib/libqp_gsmo.c

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
% SUMMARY:  Log Forward Backward, to solve overflow problem
% AUTHOR:   QIUQIANG KONG
% Created:  18-11-2015
% Modified: 19-11-2015 modify max to sum (exact solution)
%           20-11-2015 add a ceiling for logbeta
%           21-11-2015 use log instead of all to avoid overfit
%           25-11-2015 Add annotation
% -----------------------------------------------------------
% input:
%   p_xn_given_zn  p(xn|zn), size: N*Q
%   p_start        p(z1), size: Q
%   A              p(zn|zn-1), size: Q*Q
% output:
%   loggamma       ln p(zn|X), size: N*Q
%   logksi         ln p(zn,zn-1|X), size: N*Q*Q
%   loglik         ln p(X), to monitor convergence
% ===========================================================
function [loggamma, logksi, loglik] = LogForwardBackward(logp_xn_given_zn, p_start, A)
    % reserve space
    [N,Q] = size(logp_xn_given_zn);
    logalpha = zeros(N,Q);
    logbeta = zeros(N,Q);
    logc = zeros(N,1);
    loggamma = zeros(N,Q);
    logksi = zeros(N,Q,Q);

    % init log alpha(z1), log beta(zN), c(1)
    tmp = logp_xn_given_zn(1,:) + log(p_start);
    logc(1) = log( sum( exp( tmp - max(tmp) ) ) ) + max(tmp);
    logalpha(1,:) = -logc(1) + logp_xn_given_zn(1,:) + log(p_start);
    logbeta(N,:) = 0;

    % calculate logalpha, c
    for n = 2:N
        tmp = bsxfun(@plus, bsxfun(@plus, log(A), logalpha(n-1,:)'), logp_xn_given_zn(n,:));
        logc(n) = log ( sum( sum ( exp ( tmp - max(tmp(:)) ) ) ) ) + max(tmp(:));
        for q = 1:Q
            tmp2 = logalpha(n-1,:) + log(A(:,q)');
            if (isinf(max(tmp2)))
                logalpha(n,q) = -inf;
            else
                logalpha(n,q) = -logc(n) + logp_xn_given_zn(n,q) + log( sum( exp( tmp2 - max(tmp2) ) ) ) + max(tmp2);
            end
        end
    end

    % calculate logbeta
    for n = N-1:-1:1
        for q = 1:Q
            tmp = logbeta(n+1,:) + logp_xn_given_zn(n+1,:) + log(A(q,:));
            logbeta(n,q) = -logc(n+1) + log( sum( exp( tmp - max(tmp) ) ) ) + max(tmp);
        end
    end

    % calculate loggamma
    loggamma = logalpha + logbeta;

    % calculate logksi
    for n = 2:N
        logksi(n,:,:) = -logc(n) + bsxfun(@plus, bsxfun(@plus, log(A), logalpha(n-1,:)'), logp_xn_given_zn(n,:) + logbeta(n,:));
    end
    logksi(1,:,:) = [];

    % calculate likelihood
    loglik = sum(logc);
end

author: qiuqiangkong

repo: matlab-hmm


text:
function f = tsfilter(ts,flt,mode)

% function f = tsfilter(ts,flt,mode)
%
% <ts> is a b x t matrix with time-series data oriented along the rows
% <flt> is a row vector with the filter (in either the Fourier domain or space domain)
% <mode> (optional) is
%   0 means interpret <flt> as a magnitude filter in the Fourier domain,
%     and do the filtering in the Fourier domain.
%   [1 sz A] means interpret <flt> as a magnitude filter in the Fourier domain,
%     but do the filtering in the space domain using imfilter.m and 'replicate'.
%     in order to convert the Fourier filter to the space domain, we use
%     fouriertospace1D.m and sz and use A as the <mode> input to fouriertospace1D.m.
%     you can omit A, in which case we default A to 1 (which means to ensure that
%     the space filter sums to 1).
%   2 means interpret <flt> as a space-domain filter and do the filtering in the
%     space domain using imfilter.m and 'replicate'.
%   default: 0.
%
% return the filtered time-series data.  we force the output to be real-valued.
% in general, beware of wraparound and edge issues!
%
% history:
% - 2017/10/30 - speed-ups
%
% example:
% flt = zeros(1,100);
% flt(40:60) = 1;
% flt = ifftshift(flt);
% figure; plot(tsfilter(randn(1,100),flt));

% SEE ALSO IMAGEFILTER.M

% constants
num = 1000;  % number to do at a time

% input
if ~exist('mode','var') || isempty(mode)
  mode = 0;
end
if length(mode)==2
  mode = [mode 1];
end

% construct space filter if necessary
if mode(1)==1
  flt = fouriertospace1D(flt,mode(2),[],mode(3));
end

% do it
switch mode(1)
case 0
  f = zeros(size(ts),class(ts));
  for p=1:ceil(size(ts,1)/num)
%    statusdots(p,ceil(size(ts,1)/num));
    mn = (p-1)*num+1;
    mx = min(size(ts,1),(p-1)*num+num);
    f(mn:mx,:) = real(ifft(fft(ts(mn:mx,:),[],2) .* repmat(flt,[mx-mn+1 1]),[],2));
%SLOW:
%    f = cat(1,f,real(ifft(fft(ts(mn:mx,:),[],2) .* repmat(flt,[mx-mn+1 1]),[],2)));
  end
case {1 2}
  f = processmulti1D(@imfilter,ts,flt,'replicate','same','conv');
end

author: cvnlab

repo: GLMsingle

save_path: github-repos/MATLAB/cvnlab-GLMsingle


text:
function varargout = process_evt_combine( varargin )
% PROCESS_EVT_COMBINE: Combine different categories events into one (stimulus / response)
%
% USAGE:  OutputFiles = process_evt_combine('Run', sProcess, sInputs)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2012-2023

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Combine stim/response';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Events';
    sProcess.Index       = 53;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/EventMarkers?highlight=%28Combine+stim%2Fresponse%29#Events';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'raw', 'matrix'};
    sProcess.OutputTypes = {'data', 'raw', 'matrix'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Event name
    sProcess.options.combine.Comment = ['Example: We have one stim (A) and two responses (B and C).
' ...
                                        'We want to create two new pairs of event categories:
' ...
                                        '    A_AB: Event A (followed by B)
' ...
                                        '    A_AC: Event A (followed by C)
' ...
                                        '    B_AB: Event B (preceded by A)
' ...
                                        '    C_AC: Event C (preceded by A)
' ...
                                        'For that, we use the following classification:
' ...
                                        '    A_AB , B_AB , A , B
' ...
                                        '    A_AC , C_AC , A , C
' ...
                                        'To prevent one category to be created, use "ignore". Example:
' ...
                                        '    A_AB , ignore , A , B
' ...
                                        'To create an extended event AB between A and B, use "extend":
' ...
                                        '    AB , extend , A , B
' ...
                                        'You may add as many combinations as needed, one per line.

'];
    sProcess.options.combine.Type    = 'textarea';
    sProcess.options.combine.Value   = ['A_AB, B_AB , A , B' 10 'A_AC, C_AC , A , C'];
    % Maximum duration between simulateous events
    sProcess.options.dt.Comment = 'Maximum delay between grouped events: ';
    sProcess.options.dt.Type    = 'value';
    sProcess.options.dt.Value   = {1, 'ms', 0};
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    % Return all the input files
    OutputFiles = {};

    % ===== GET OPTIONS =====
    dt = sProcess.options.dt.Value{1};
    % Combination string
    combineStr = strtrim(sProcess.options.combine.Value);
    % Split in lines
    combine_lines = str_split(combineStr, [10 13]);
    % Split each line with the ";"
    combineCell = {};
    for i = 1:length(combine_lines)
        % No information on this line: skip
        combine_lines{i} = strtrim(combine_lines{i});
        if isempty(combine_lines{i})
            continue;
        end
        % Split with ";"
        lineCell = str_split(combine_lines{i}, ';,');
        % If there are not 4 elements on this line: ignore
        if (length(lineCell) ~= 4)
            continue;
        end
        % Add combination entry
        iComb = size(combineCell,1) + 1;
        combineCell{iComb,1} = strtrim(lineCell{1});
        combineCell{iComb,2} = strtrim(lineCell{2});
        combineCell{iComb,3} = strtrim(lineCell{3});
        combineCell{iComb,4} = strtrim(lineCell{4});
    end
    % If no combination available
    if isempty(combineCell)
        bst_report('Error', sProcess, [], 'Invalid combinations format.');
        return;
    end

    % For each file
    for iFile = 1:length(sInputs)
        % ===== GET FILE DESCRIPTOR =====
        isRaw = strcmpi(sInputs(iFile).FileType, 'raw');
        % Load the raw file descriptor
        if isRaw
            DataMat = in_bst_data(sInputs(iFile).FileName, 'F');
            sFile = DataMat.F;
        else
            DataMat = in_bst_data(sInputs(iFile).FileName, 'Events');
            sFile.events = DataMat.Events;
            sFile.epochs = [];
        end

        % Process only continuous files
        if ~isempty(sFile.epochs)
            bst_report('Error', sProcess, sInputs(iFile), 'This function can only process continuous recordings (no epochs). Skipping File...');
            continue;
        end
        % If no markers are present in this file
        if isempty(sFile.events)
            bst_report('Warning', sProcess, sInputs(iFile), 'This file does not contain any event. Skipping File...');
            continue;
        end

        % ===== PROCESS EACH COMBINATION =====
        isModified = 0;
        for iComb = 1:size(combineCell,1)
            % === GET EVENTS INVOLVED ===
            % Find all the events in the structure
            evtNewStimName = combineCell{iComb,1};
            evtNewRespName = combineCell{iComb,2};
            iEvtNewStim = find(strcmpi(evtNewStimName, {sFile.events.label}));
            iEvtNewResp = find(strcmpi(evtNewRespName, {sFile.events.label}));
            iEvtStim    = find(strcmpi(combineCell{iComb,3}, {sFile.events.label}));
            iEvtResp    = find(strcmpi(combineCell{iComb,4}, {sFile.events.label}));
            % If any event not available
            if isempty(iEvtStim)
                bst_report('Warning', sProcess, sInputs(iFile), ['Event "' combineCell{iComb,3} '" does not exist in file. Skipping combination...']);
                continue;
            end
            if isempty(iEvtResp)
                bst_report('Warning', sProcess, sInputs(iFile), ['Event "' combineCell{iComb,4} '" does not exist in file. Skipping combination...']);
                continue;
            end
            % If events are extended events: skip
            if (size(sFile.events(iEvtStim).times,1) > 1) || (size(sFile.events(iEvtResp).times,1) > 1)
                bst_report('Error', sProcess, sInputs(iFile), 'Cannot process extended events. Skipping combination...');
                continue;
            end

            % === FIND OTHER RESPONSES ===
            % Get all the other possible responses to this stim
            iOtherResp = find(strcmpi(combineCell(:,3), combineCell{iComb,3}));
            % Do not consider the target response
            iOtherResp = setdiff(iOtherResp, iComb);
            % Loop to get all the times of all the other responses
            tOtherResp = [];
            for iOther = 1:length(iOtherResp)
                % Find event
                evtOtherRespName = combineCell{iOtherResp(iOther),4};
                iEvtOtherResp = find(strcmpi(evtOtherRespName, {sFile.events.label}));
                if isempty(iEvtOtherResp) || isempty(sFile.events(iEvtOtherResp).times) || (size(sFile.events(iEvtOtherResp).times,1) > 1)
                    continue;
                end
                % Get times for this response event
                tOtherResp = [tOtherResp, sFile.events(iEvtOtherResp).times];
            end

            % === DETECT RESPONSE TO STIM ===
            % Get stim and resp times
            tStim = sFile.events(iEvtStim).times;
            tResp = sFile.events(iEvtResp).times;
            iNewStim = [];
            iNewResp = [];
            % For each stim event, look for a response
            for is = 1:length(tStim)
                % Look a TARGET response event immediately following a stim event
                if (is == length(tStim))
                    ir = find(tResp > tStim(is), 1);
                else
                    ir = find((tResp > tStim(is)) & (tResp < tStim(is+1)), 1);
                end
                % Look an OTHER response event immediately following a stim event
                if (is == length(tStim))
                    ior = find(tOtherResp > tStim(is), 1);
                else
                    ior = find((tOtherResp > tStim(is)) & (tOtherResp < tStim(is+1)), 1);
                end
                % If a response was found, and before any other response type
                if ~isempty(ir) && (isempty(ior) || (tResp(ir) <= tOtherResp(ior))) && (tResp(ir) - tStim(is) < dt)
                    iNewStim(end+1) = is;
                    iNewResp(end+1) = ir;
                end
            end

            % === CREATE NEW EVENTS ===
            % If some couples stim/response were detected: create or udpate the correpsonding events
            if ~isempty(iNewStim)
                isModified = 1;

                % Should the events include channels
                if ~isempty(sFile.events(iEvtStim).channels) || ~isempty(sFile.events(iEvtResp).channels)
                    isChannels = 1;
                    if isempty(sFile.events(iEvtStim).channels)
                        sFile.events(iEvtStim).channels = cell(1, size(sFile.events(iEvtStim).times,2));
                    end
                    if isempty(sFile.events(iEvtResp).channels)
                        sFile.events(iEvtResp).channels = cell(1, size(sFile.events(iEvtResp).times,2));
                    end
                else
                    isChannels = 0;
                end

                % Should the events include notes
                if ~isempty(sFile.events(iEvtStim).notes) || ~isempty(sFile.events(iEvtResp).notes)
                    isNotes = 1;
                    if isempty(sFile.events(iEvtStim).notes)
                        sFile.events(iEvtStim).notes = cell(1, size(sFile.events(iEvtStim).times,2));
                    end
                    if isempty(sFile.events(iEvtResp).notes)
                        sFile.events(iEvtResp).notes = cell(1, size(sFile.events(iEvtResp).times,2));
                    end
                else
                    isNotes = 0;
                end

                % === NEW STIM ===
                % Only create new stim event if name is not "ignore"
                if ~strcmpi(evtNewStimName, 'ignore')
                    % If new stim event does not exist: create
                    if isempty(iEvtNewStim)
                        % Initialize new event
                        iEvtNewStim = length(sFile.events) + 1;
                        sEvent = db_template('event');
                        sEvent.label = evtNewStimName;
                        % Get the default color for this new event
                        sEvent.color = panel_record('GetNewEventColor', iEvtNewStim, sFile.events);
                        % Add new event to list
                        sFile.events(iEvtNewStim) = sEvent;
                    end
                    % Add simple events
                    if ~strcmpi(evtNewRespName, 'extend')
                        sFile.events(iEvtNewStim).times    = [sFile.events(iEvtNewStim).times,   sFile.events(iEvtStim).times(iNewStim)];
                        sFile.events(iEvtNewStim).epochs   = [sFile.events(iEvtNewStim).epochs,  sFile.events(iEvtStim).epochs(iNewStim)];
                        % Sort
                        [sFile.events(iEvtNewStim).times, indSort] = unique(sFile.events(iEvtNewStim).times);
                        sFile.events(iEvtNewStim).epochs   = sFile.events(iEvtNewStim).epochs(indSort);
                        % Channels
                        if isChannels
                            sFile.events(iEvtNewStim).channels = [sFile.events(iEvtNewStim).channels,sFile.events(iEvtStim).channels(iNewStim)];
                            sFile.events(iEvtNewStim).channels = sFile.events(iEvtNewStim).channels(indSort);
                        end
                        % Notes
                        if isNotes
                            sFile.events(iEvtNewStim).notes = [sFile.events(iEvtNewStim).notes,   sFile.events(iEvtStim).notes(iNewStim)];
                            sFile.events(iEvtNewStim).notes = sFile.events(iEvtNewStim).notes(indSort);
                        end
                    % Add extended events
                    else
                        sFile.events(iEvtNewStim).times    = [sFile.events(iEvtNewStim).times,   [sFile.events(iEvtStim).times(iNewStim); sFile.events(iEvtResp).times(iNewResp)]];
                        sFile.events(iEvtNewStim).epochs   = [sFile.events(iEvtNewStim).epochs,  sFile.events(iEvtStim).epochs(iNewStim)];
                        if isChannels
                            sFile.events(iEvtNewStim).channels = [sFile.events(iEvtNewStim).channels,sFile.events(iEvtStim).channels(iNewStim)];
                        end
                        if isNotes
                            sFile.events(iEvtNewStim).notes    = [sFile.events(iEvtNewStim).notes,   sFile.events(iEvtStim).notes(iNewStim)];
                        end
                    end
                end

                % === NEW RESP ===
                % Only create new response event if name is not "ignore" or "extend"
                if ~strcmpi(evtNewRespName, 'ignore') && ~strcmpi(evtNewRespName, 'extend')
                    % If new stim event does not exist: create
                    if isempty(iEvtNewResp)
                        % Initialize new event
                        iEvtNewResp = length(sFile.events) + 1;
                        sEvent = db_template('event');
                        sEvent.label = evtNewRespName;
                        % Get the default color for this new event
                        sEvent.color = panel_record('GetNewEventColor', iEvtNewResp, sFile.events);
                        % Add new event to list
                        sFile.events(iEvtNewResp) = sEvent;
                    end
                    % Add occurrences
                    sFile.events(iEvtNewResp).times    = [sFile.events(iEvtNewResp).times,    sFile.events(iEvtResp).times(iNewResp)];
                    sFile.events(iEvtNewResp).epochs   = [sFile.events(iEvtNewResp).epochs,   sFile.events(iEvtResp).epochs(iNewResp)];
                    % Sort
                    [sFile.events(iEvtNewResp).times, indSort] = unique(sFile.events(iEvtNewResp).times);
                    sFile.events(iEvtNewResp).epochs   = sFile.events(iEvtNewResp).epochs(indSort);
                    % Channels
                    if isChannels
                        sFile.events(iEvtNewResp).channels = [sFile.events(iEvtNewResp).channels, sFile.events(iEvtResp).channels(iNewResp)];
                        sFile.events(iEvtNewResp).channels = sFile.events(iEvtNewResp).channels(indSort);
                    end
                    % Notes
                    if isNotes
                        sFile.events(iEvtNewResp).notes = [sFile.events(iEvtNewResp).notes,    sFile.events(iEvtResp).notes(iNewResp)];
                        sFile.events(iEvtNewResp).notes = sFile.events(iEvtNewResp).notes(indSort);
                    end
                end
            end
        end

        % ===== SAVE RESULT =====
        % Report results
        if isRaw
            DataMat.F = sFile;
        else
            DataMat.Events = sFile.events;
        end
        % Only save changes if something was change
        if isModified
            bst_save(file_fullpath(sInputs(iFile).FileName), DataMat, 'v6', 1);
        end
        % Return all the input files
        OutputFiles{end+1} = sInputs(iFile).FileName;
    end
end

author: brainstorm-tools


text:
function [qt] = m32qt(m3)
% Convert volume from cubic meters to US liquid quarts.
% Chad Greene 2012
qt = m3*1056.6882094;


text:
function M = permutation(P)
% PERMUTATION
%
% M = permutation(X)
%
% Creates the model [sum(P,1)==1,sum(P,2)==1,binary(P)]

M = [sum(P,1)==1,sum(P,2)==1,binary(P)];

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [ nb, ni, np, ns, pr_bud, pr_bud_angular, r, seed ] = ...
  parameter_initialize ( )

%*****************************************************************************80
%
%% PARAMETER_INITIALIZE sets the initial values of parameters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 December 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer NB, the number of boundary points.
%
%    Output, integer NI, the number of interior points.
%
%    Output, integer NP, the total number of points.
%
%    Output, integer NS, the number of sample points.
%
%    Output, real PR_BUD, the probability that a boundary point will bud.
%
%    Output, real PR_BUD_ANGULAR, the probability that a boundary point bud
%    will be constrained to the boundary.
%
%    Output, real R, the radius of the circle.
%
%    Output, integer SEED, a seed for the random number generator.
%
  nb = 30;
  ni = 60;
  np = ni + nb;
  ns = 10000;

  pr_bud = 0.05;
  pr_bud_angular = 0.30;

  r = 1.0;
%
%  Call RNG to initialize the random number generator with a random seed,
%  and record that value in SEED.
%
  rng_struct = rng ( 'shuffle' );
  seed = rng_struct.Seed;

  return
end


text:
function MCMC_State = stat_ubsplsm_init(varargin)
% Initialize MCMC estimator for univariate B-Spline smoother
%
% Author: Tim Mullen and Wes Thompson, 2010-13, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

arg_define([0 Inf],varargin, ...
    arg_norep({'Y','TSData'},mandatory,[],sprintf(['Cell array of data to smooth.\n' ...
              'Generally, Y is a cell array where Y{i} is the T x 1 vector of time-varying (or freq-varying) connectivity for the kth channel pair of the sth subject.\n' ...
              'e.g. Y = {s1(1,1) s1(1,2) s1(1,3) ... s1(N,1) s1(N,2) s1(N,3) ... \n' ...
              '          s2(1,1) s2(1,2) s2(1,3) ... } \n'])), ...
    arg({'K','Knots'},5,[],'Positions of spline knots. If K is a scalar, then we automatically assign K evenly spaced knots. A good heuristic is one knot every 5%','shape','row'), ...
    arg({'Q','FPCABasisDim','fpcaBasisDim'},4,[0 Inf],'Number of FPCA basis functions.'), ...
    arg({'order','BSplineOrder'},4,[1 Inf],'B-Spline order. 4 is a good choice'), ...
    arg({'initNoiseVariance'},0.1,[eps Inf],'Initial noise variance'), ...
    arg({'verb','VerbosityLevel'},2,{int32(0) int32(1) int32(2)},'Verbosity level. 0 = no output, 1 = text, 2 = graphical') ...
);
% arg({'T','TimeSeriesLength'},mandatory,[],'Time series length. This is the length of the time-series (vector) to smooth'), ...
% arg({'R','NumTimeSeries'},mandatory,[],'Number of observations. This is the total number of time-series to smooth. For instance, this might be the total number of connectivity edges across all subjects.'), ...

T = size(Y{1},1);
R = length(Y);

% Determine knot locations
% -------------------------------------------------------------------------
if isscalar(K)
    idx   = 1:T;
    K = K+order-1;
    knots = quantile(idx,(0:1:(K-order+1)) / (K-order+1));
    knots = knots(:)';
else
    % knots are already provided
    knots = K(:);
    K = length(knots)+2;  % add endpoints
end

% Construct orthonormal basis functions
% -------------------------------------------------------------------------
phi_t = stat_ubsplsm_mkspl(knots,T,order,verb);

% Initialize FPCA
% -------------------------------------------------------------------------
if verb==2
    multiWaitbar('Initializing FPCA','Reset','Color',hlp_getNextUniqueColor);
end

% Initialize the FPCA components to random, orthonormal vectors
% [!] should be able to replace this whole section with a single line of orth(mvrnd(...)). We only need orthonormal random vectors
THETA(:,1,1) = mvnrnd(zeros(K,1),eye(K));
THETA(:,1,1) = THETA(:,1,1)/norm(THETA(:,1,1));
for q = 2:Q
    if verb==2
        multiWaitbar('Initializing FPCA',q/Q);
    end
    THETA(:,q,1) = mvnrnd(zeros(K,1),eye(K));

    % orthonormalization of random vector
    for r = 1:(q-1)

        THETA(:,q,1) = THETA(:,q,1) ...
                       -(THETA(:,q,1)'*THETA(:,q-r,1)/norm(THETA(:,q-r,1).^2)) ...
                       * THETA(:,q-r,1);
    end
    THETA(:,q,1) = THETA(:,q,1)/norm(THETA(:,q,1));
end

% Construct MCMC state object
% -------------------------------------------------------------------------
MCMC_State.THETA     = THETA;
MCMC_State.ALPHA     = zeros(Q,R);
MCMC_State.ALPHA_BAR = zeros(Q,1);
MCMC_State.SIGMA_EPS = initNoiseVariance;
MCMC_State.phi_t     = phi_t;
MCMC_State.initstate = true;

if verb==2
    multiWaitbar('Initializing FPCA', 'Close');
end

repo: REST


text:
function [text,AIC,BIC]=agarch_display(parameters,ll,vcv,epsilon,p,q,modelType,errorType)
% Display parameters, tstats, pvals, log-likelihood and AIC/BIC
% from estimates of a AGARCH(P,Q) or NAGARCH(P,Q) produced using agarch
%
% USAGE:
%   [TEXT] = agarch_display(PARAMETERS,LL,VCV,DATA,P,Q)
%   [TEXT,AIC,BIC] = agarch_display(PARAMETERS,LL,VCV,DATA,P,Q,MODELTYPE,ERRORTYPE)
%
% INPUTS:
%   PARAMETERS    - A 2+p+q column vector of parameters with
%                   [omega alpha(1) ... alpha(p) gamma beta(1) ... beta(q) [nu lambda]]'
%   LL            - The log likelihood at the optimum
%   VCV           - Variance-covariance matrix (Robust or inverse Hessian)
%   EPSILON       - Column of mean zero data used to fit the model
%   P             - Positive, scalar integer representing the number of
%                   symmetric innovations
%   Q             - Non-negative, scalar integer representing the number
%                   of lags of conditional variance (0 for ARCH)
%   MODEL_TYPE    - [OPTIONAL] The type of variance process, either
%                     'AGARCH'  - Asymmetric GARCH, Engle (1990) [DEFAULT]
%                     'NAGARCH' - Nonlinear Asymmetric GARCH, Engle & Ng (1993)
%   ERRORTYPE     - [OPTIONAL] The error distribution used, valid types are:
%                     'NORMAL'    - Gaussian Innovations [DEFAULT]
%                     'STUDENTST' - T distributed errors
%                     'GED'       - Generalized Error Distribution
%                     'SKEWT'     - Skewed T distribution
%
% OUTPUTS:
%   TEXT          - Character matrix with the formatted parameters of the model
%   AIC           - Aikake Information Criteria computed from the LL
%   BIC           - Schwartz/Bayesian Information Criteria computed from the LL
%
%  See also AGGARCH

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 7/12/2005

%%%%%%%%%%%%%%%%%%%%%%%%%%
% PARAMETER CHECKING
%%%%%%%%%%%%%%%%%%%%%%%%%%
switch nargin
    case 5
        q=0;
        modelType = 'AGARCH';
        errorType = 'NORMAL';
    case 6
        modelType = 'AGARCH';
        errorType = 'NORMAL';
    case 7
        errorType = 'NORMAL';
    case 8
    otherwise
        error('5 to 8 inputs required')
end

% parameters, N by 1, real
if any(~isreal(parameters)) || size(parameters,2)~=1
    error('PARAMETERS must be a column vector.')
end
% LL
if ~isscalar(ll) || ~isreal(ll)
    error('LL must be a scalar.')
end
% VCV
if size(vcv,2)~=size(vcv,1) || any(min(eig(vcv))<=0) || size(vcv,1)~=length(parameters)
    error('VCV must be a square positive definite matrix compatible with PARAMETERS.')
end
% epsilon
if any(~isreal(epsilon)) || size(epsilon,2)~=1
    error('EPSILON must be a T by 1 column vector.')
end
% p
if ~isscalar(p) || p<1 || floor(p)~=p
    error('P must be a postitive scalar integer')
end
% q
if ~isscalar(q) || q<0 || floor(q)~=q
    error('Q must be a non-negative scalar integer')
end
% errorType
if ~ischar(errorType)
    errorType=[];
end
if isempty(errorType)
    errorType='NORMAL';
end
errorType = upper(errorType);
switch errorType
    case 'NORMAL'
        errorType=1;
        extraP=0;
    case 'STUDENTST'
        errorType=2;
        extraP=1;
    case 'GED'
        errorType=3;
        extraP=1;
    case 'SKEWT'
        errorType=4;
        extraP=2;
    otherwise
        error('ERRORTYPE is not one of the supported distributions')
end
% modelType
modelType
if ~ischar(modelType)
    modelType=[];
end
modelType = upper(modelType);
if isempty(modelType)
    modelType='AGARCH';
end
switch modelType
    case 'AGARCH'
        modelType=1;
    case 'NAGARCH'
        modelType=2;
    otherwise
        error('MODELTYPE must be either ''AGARCH'' or ''NAGARCH''.')
end

if length(parameters)~=(2+p+q+extraP)
    error('Size of PARAMETERS is not compatible with input P and Q')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%
% PARAMETER CHECKING
%%%%%%%%%%%%%%%%%%%%%%%%%%

if modelType == 1
    model_name = ['AGARCH(' num2str(p) ',' num2str(q) ')'];
else
    model_name = ['NAGARCH(' num2str(p) ',' num2str(q) ')'];
end

parameters_text=num2str(parameters,'%4.4f');
stderr=sqrt(diag(vcv));
stderr_text=num2str(stderr,'%4.4f');
tstats=parameters./stderr;
tstats_text=num2str(tstats,'%4.4f');
pvals=2-2*normcdf(abs(tstats));
pvals_text=num2str(pvals,'%4.4f');

T=length(epsilon);
AIC = -ll/T+2*length(parameters)/T;
BIC = -ll/T+log(T)*length(parameters)/T;

% Format the output
text=[];
text{1,1}= ' ';
text{2,1}= ' ';
text{3,1}=repmat('-',1,50);
text{4,1}=model_name;
text{5,1}=repmat('-',1,50);
text{6,1}=' ';
text{7,1}=['Loglikelihood: ' sprintf('%1.2f',ll)];
text{8,1}=['AIC: ' sprintf('%1.4f',AIC)];
text{9,1}=['BIC: ' sprintf('%1.4f',BIC)];
text{10,1}=' ';

for i=1:size(text,1);
    disp(text{i,1})
end

% Format the parameter table, need to right align everything
K=size(parameters_text,1);
for i=1:K
    N=size(parameters_text,2);
    if any(parameters_text(i,:)==' ')
        numSpace=sum(parameters_text(i,:)==' ');
        parameters_text(i,:)=[repmat(' ',1,numSpace) parameters_text(i,1:N-numSpace)];
    end
    N=size(stderr_text,2);
    if any(stderr_text(i,:)==' ')
        numSpace=sum(stderr_text(i,:)==' ');
        stderr_text(i,:)=[repmat(' ',1,numSpace) stderr_text(i,1:N-numSpace)];
    end
    N=size(tstats_text,2);
    if any(tstats_text(i,:)==' ')
        numSpace=sum(tstats_text(i,:)==' ');
        tstats_text(i,:)=[repmat(' ',1,numSpace) tstats_text(i,1:N-numSpace)];
    end
    N=size(pvals_text,2);
    if any(pvals_text(i,:)==' ')
        numSpace=sum(pvals_text(i,:)==' ');
        pvals_text(i,:)=[repmat(' ',1,numSpace) pvals_text(i,1:N-numSpace)];
    end
end
% Append column labels
labels={' Parameters','   Std. Err.','     T-stat','      P-val'};
cols = {parameters_text,stderr_text,tstats_text,pvals_text};
for i=1:length(labels);
    text1=labels{i};
    text2=cols{i};
    maxcols=max(size(text1,2),size(text2,2));
    if size(text1,2)<maxcols
        labels{i} = [repmat(' ',1,maxcols-size(text1,2)) text1 ];
    end
    if size(text2,2)<maxcols
        cols{i} = [repmat(' ',size(text2,1),maxcols-size(text2,2)) text2 ];
    end
end

% Finally variable names
variable_names=cell(length(parameters),1);
variable_names{1}='omega';
index=2;
for i=1:p
    variable_names{index}=['alpha(' num2str(i) ')'];
    index=index+1;
end
variable_names{index}=['gamma'];
index=index+1;
for i=1:q
    variable_names{index}=['beta(' num2str(i) ')'];
    index=index+1;
end

maxVarNameLength=max(cellfun('length',variable_names));
for i=1:length(variable_names)
    variable_names{i} = [repmat(' ',1,maxVarNameLength-length(variable_names{i})) variable_names{i}];
end
variable_names=cell2mat(variable_names(:));

% Output the parameters
outmat=strvcat(' ',variable_names);
for i=1:4
    outmat=[outmat repmat('  ',size(outmat,1),1) strvcat(labels{i},cols{i})];
end

text2=[];
for i=1:length(text)
    text2=strvcat(text2,text{i});
end
text=strvcat(text2,outmat);

disp(outmat)

author: bashtage

repo: mfe-toolbox


text:
function Y = multiherm(X)
% Returns the Hermitian parts of the matrices in the 3D matrix X
%
% function Y = multiherm(X)
%
% Y is a 3D matrix the same size as X. Each slice Y(:, :, i) is the
% Hermitian part of the slice X(:, :, i).
%
% See also: multiprod multitransp multihconj multiscale multiskew

% This file is part of Manopt: www.manopt.org.
% Original author: Hiroyuki Sato, April 27, 2015.
% Contributors:
% Change log:

    Y = .5*(X + multihconj(X));

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function [u, s, U_p, U_k, U_d] = pinHole(p,k,d)

% PINHOLE Pin-hole camera model, with optional radial distortion.
%   U = PINHOLE(P) gives the projected pixel U of a point P in a canonical
%   pin-hole camera, that is, with calibration parameters
%     u0 = 0
%     v0 = 0
%     au = 1
%     av = 1
%   It uses reference frames {RDF,RD} (right-down-front for the 3D world
%   points and right-down for the pixel), according to this scheme:
%
%         / z (forward)
%        /
%       +------- x                 +------- u
%       |                          |
%       |      3D : P=[x;y;z]      |     image : U=[u;v]
%       | y                        | v
%
%   U = PINHOLE(P,K) allows the introduction of the camera's calibration
%   parameters:
%     K = [u0 v0 au av]'.
%
%   U = PINHOLE(P,K,D) allows the introduction of the camera's radial
%   distortion parameters:
%     D = [K2 K4 K6 ...]'
%   so that the new pixel is distorted following the distortion equation:
%     U_D = U * (1 + K2*R^2 + K4*R^4 + ...)
%   with R^2 = sum(U.^2), being U the projected point in the image plane
%   for a camera with unit focal length.
%
%   [U,S] = PINHOLE(...) returns the depth S from the camera center.
%
%   If P is a points matrix, PINHOLE(P,...) returns a pixel matrix U and a
%   depths row-vector S. P, U and S are defined as
%     P = [P1 ... Pn];   Pi = [xi;yi;zi]
%     U = [U1 ... Un];   Ui = [ui;vi]
%     S = [S1 ... Sn]
%
%   [U,S,U_p,U_k,U_d] returns the Jacobians of U wrt P, K and D. It only
%   works for single points P=[x;y;z].
%
%   See also PERSP_PROJECT, DISTORT, PIXELLISE, INVPINHOLE, PINHOLEIDP, INVDISTORTION.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

% Point's depth
s = p(3,:);

if nargout <= 2 % only pixel

    switch nargin
        case 1
            u = persp_project(p);
        case 2
            u = pixellise(persp_project(p),k);
        case 3
            if ~isempty(d)
                u = pixellise(distort(persp_project(p),d),k);
            else
                u = pixellise(persp_project(p),k);
            end
    end

else % Jacobians

    if size(p,2) > 1
        error('Jacobians not available for multiple points')
    else

        switch nargin
            case 1
                [u, U_p] = persp_project(p);

            case 2
                [u1, U1_p]     = persp_project(p);
                [u, U_u1, U_k] = pixellise(u1,k);
                U_p            = U_u1*U1_p;

            case 3
                if ~isempty(d)
                    [u1, U1_p]        = persp_project(p);
                    [u2, U2_u1, U2_d] = distort(u1,d);
                    [u, U_u2, U_k]    = pixellise(u2,k);
                    U_d               = U_u2*U2_d;
                    U_p               = U_u2*U2_u1*U1_p;
                else
                    [u1, U1_p]     = persp_project(p);
                    [u, U_u1, U_k] = pixellise(u1,k);
                    U_p            = U_u1*U1_p;
                    U_d            = zeros(2,0);
                end

        end

    end

end

return

%% jacobians
syms x y z u0 v0 au av d2 d4 d6 real
p = [x;y;z];
k = [u0;v0;au;av];
d = [d2;d4;d6];

[u, s, U_p, U_k, U_d] = pinHole(p,k,d);

simplify(U_p - jacobian(u,p))
simplify(U_k - jacobian(u,k))
simplify(U_d - jacobian(u,d))

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

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [coeffs,termMat]=polyFitMultiDim(varVals,funVals,degree)
%%POLYFITMULTIDIM Obtain a least-squares multivariate polynomial fit of a
%                 specified order to a given set of data.
%
%INPUTS: varVals This is a numDimXnumPoints set oif dtapoints where the
%                function was evaluated.
%        funVals This is a numPointsX1 or 1XnumPoints vector of the values
%                of the functions at the specified data points.
%        degree The degree of the desired interpolating polynomial. The
%               degree must be less than the total number of monomial
%               terms, which is
%               (1+degree)*binomial(degree+numDim,numDim-1)/numDim.
%
%OUTPUTS: coeffs A hypermatrix of the coefficients for the multivariate
%               polynomial that can be used in the polyValMultiDim
%               function. These are arranged such that
%               coeffs(a1,a2,a3...an) corresponds to the coefficient of an
%               x1^(a1-1)*x2^(a2-1)*x3^(a3-1)...xn^(an-1) term.  Thus, the
%               number of indices coeffs takes is equal to the
%               dimensionality of x (not counting singleton dimensions at
%               the end of coeffs). Note that this ordering is the reverse
%               that used in the 1D polyval function that is built into
%               Matlab. The number of elements for each index in coeffs is
%               the maximum order of that dimension +1.
%       termMat An (n+1)XnumTerms matrix such that
%               termMat(:,i)=[c,a1,a2,...,an] is a monomial term where c is
%               the value of of the monomial coefficient and the monomial
%               is x1^(a1-1)*x2^(a2-1)*x3^(a3-1)...xn^(an-1).
%
%Multivariate polynomials consist of a sum of monomial terms of the form
%c*x1^(a1)*x2^(a2)*x3^(a3)...xn^(an)
%Here, the variables x are known as are the values z of the polynomial
%evaluated at those points. Thus, the unknowns are the c terms, the
%coefficients of the monomials. This function evaluates the monomial terms
%and build a matrix with all of the monomial term values for all of the
%equations. The solution to the coefficient values is thus the solution to
%a linear system of equations. The
%
%EXAMPLE:
%Here, we have a two-dimensional function evaluated at a number of random
%points and we want to fit a polynomial to it. In this problem, we do not
%add noise to the function.
% points=linspace(-1,1,10);
% [x,y]=meshgrid(points,points);
% z=2*x.^3+2*x.*y-3*y.^3+4*x.^2.*y-5*y.^2.*x.^3;
% varVals=[x(:)';y(:)'];
% funVals=z(:);
% degree=5;
% [~,termMat]=polyFitMultiDim(varVals,funVals,degree);
% %One will get coefficients that are almost exact. Given that we know that
% %the coefficients are integers, we can eliminate all small coefficients by
% %rounding.
% termMat=round(termMat);
% %Get rid of terms that are numerically zero.
% sel=termMat(1,:)~=0;
% termMat=termMat(:,sel);
% terms2String(termMat,{'x','y'})
%
%November 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numDim=size(varVals,1);
numEq=size(varVals,2);

%The total number of monomial terms for the multivariate polynomial. This
%is the sum of binomial(degCur+numDim-1,numDim-1) from degCur=0 to degree,
%which is the number of compositions of the degree degCur into n parts.
totalNumMonomials=(1+degree)*binomial(degree+numDim,numDim-1)/numDim;

if(totalNumMonomials>numEq)
   error('The total number of monomial coefficients to find is greater than the number of data points provided.')
end

termMat=zeros(numDim+1,totalNumMonomials);

V=zeros(numEq,totalNumMonomials);

%x^0,y^0, etc. value.
V(:,1)=1;
curTerm=2;
for curDegree=1:degree
    numMonomials=binomial(curDegree+numDim-1,numDim-1);

    for j=0:(numMonomials-1)
        curMonomial=unrankTComposition(j,numDim,curDegree+numDim,true)-1;
        termMat(2:end,curTerm)=curMonomial;
        for curEq=1:numEq
            V(curEq,curTerm)=prod(varVals(:,curEq).^curMonomial);
        end
        curTerm=curTerm+1;
    end
end

termMat(1,:)=V\funVals(:);
%Get rid of terms that are numerically zero.
sel=termMat(1,:)~=0;
termMat=termMat(:,sel);
coeffs=terms2MultiDimPolyMat(termMat);

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function creds = loadplotlycredentials()

userhome = getuserdir();

plotly_credentials_file = fullfile(userhome,'.plotly','.credentials');

% check if credentials exist
if ~exist(plotly_credentials_file, 'file')
    error('Plotly:CredentialsNotFound',...
        ['It looks like you haven''t set up your plotly '...
        'account credentials yet.\nTo get started, save your '...
        'plotly username and API key by calling:\n'...
        '>>> saveplotlycredentials(username, api_key)\n\n'...
        'For more help, see https://plot.ly/MATLAB or contact '...
        'chris@plot.ly.']);
end

fileIDCred = fopen(plotly_credentials_file, 'r');

if(fileIDCred == -1)
    error('plotly:loadcredentials', ...
        ['There was an error reading your credentials file at '...
        plotly_credentials_file '. Contact chris@plot.ly for support.']);
end

creds_string_array = fread(fileIDCred, '*char');
creds_string = sprintf('%s',creds_string_array);
creds = loadjson(creds_string);

end

author: plotly


text:
function statplot(c)

   % Private method. See ../plot for details.

   % Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   Rmean  = c.stat(:,1);
   Rrmshi = c.stat(:,2);
   Rrmslo = c.stat(:,3);
   T      = c.stat(:,4);
   Trms   = c.stat(:,5);

   % PLOT RESULTS
   figure('Color','w','Position',[50 50 1000 600]);
   set(gcf,'DefaultLineLineWidth',1);
   set(gcf,'DefaultAxesFontSize',14);

   subplot(2,1,1);
   hold on; grid on; box on;
   ylabel('Time (s)');
   title('Best fit delay time and 1\sigma error');
   xfill = [ [1:length(T)]' ; [length(T):-1:1]' ];
   yfill = [T+Trms ; flipud(T-Trms)];
   fill(xfill,yfill,[.9 .9 .9]);
   plot(c.lags','b.','MarkerSize',5);
   plot(T,'ko-','LineWidth',2);
   xlim([1 length(T)]);
   maxval = 1.05*max(max(abs(c.lags)));
   ylim([-1*maxval maxval]);

   subplot(2,1,2);
   hold on; grid on; box on;
   xfill = [ [1:length(T)]' ; [length(T):-1:1]' ];
   yfill = [Rrmshi ; flipud(Rrmslo)];
   fill(xfill,yfill,[.9 .9 .9]);
   plot(Rmean,'ko-','LineWidth',2);
   xlabel('Trace no.');
   ylabel('Mean R');
   title('Mean max correlation and 1\sigma error');
   xlim([1 length(T)]);
   ylim([0 1]);

   % PREP PRINT OUT
   set(gcf, 'paperorientation', 'landscape');
   set(gcf, 'paperposition', [1 3 9 4] );
   print(gcf, '-depsc2', 'FIG_stat.ps');
end

author: geoscience-community-codes

repo: GISMO


text:
function extOut = fileExtension(filename)
%
%
%        extOut = fileExtension(filename)
%
%
%        Description: return the extension of a file
%
%        Input:
%           -filename: the name of the file
%
%        Output:
%           -extOut: the extension of the file
%
%     Copyright (C) 2011-2016  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

k = strfind(filename, '.');

%there is no extension
if(isempty(k))
    extOut = [];
else
    %get the real extension
    k = k(end);
    extOut = filename((k + 1):max(size(filename)));
end

end

author: banterle


text:
%% ODF Estimation from EBSD data
%
%%
% In order to discuss ODF estimation from individual orientation data we
% start by loading an EBSD data set

mtexdata copper

plot(ebsd,ebsd.orientations)

%%
% of copper. The orientation distribution function can now be computed by

odf = calcDensity(ebsd('copper').orientations)

plotSection(odf,'contourf')
mtexColorMap LaboTeX
mtexColorbar

%%
% The function <orientation.calcODF.html calcODF> implements the ODF
% estimation from EBSD data in MTEX. The underlying statistical method is
% called kernel density estimation, which can be seen as a generalized
% histogram. To be more precise, let's $\psi : SO(3) \to R$ be a radially
% symmetric, unimodal model ODF. Then the kernel density estimator for the
% individual orientation data $o_1,o_2,\ldots,o_M$ is defined as
%
% $$f(o) = \frac{1}{M} \sum_{i=1}^{M} \psi(o o_i^{-1})$$
%
% The choice of the model ODF $\psi$ and in particular its halfwidth has a
% great impact in the resulting ODF. If no halfwidth is specified the
% default halfwidth of 10 degrees is selected.
%
%% Automatic halfwidth selection
%
% MTEX includes an automatic halfwidth selection algorithm which is called
% by the command <orientation.calcKernel.html calcKernel>. To work
% properly, this algorithm needs spatially independent EBSD data as in the
% case of this dataset of very rough EBSD measurements (only one
% measurement per grain).

% try to compute an optimal kernel
psi = calcKernel(ebsd.orientations)

%%
% In the above example, the EBSD measurements are spatial dependent and the
% resulting halfwidth is too small. To avoid this problem we have to perform
% grain reconstruction first and then estimate the halfwidth from the
% grains.

% grains reconstruction
grains = calcGrains(ebsd);

% correct for to small grains
grains = grains(grains.grainSize>5);

% compute optimal halfwidth from the meanorientations of grains
psi = calcKernel(grains('co').meanOrientation)

% compute the ODF with the kernel psi
odf = calcDensity(ebsd('co').orientations,'kernel',psi)

%%
% Once an ODF is estimated all the functionality MTEX offers for
% <ODFCharacteristics.html ODF analysis> and <ODFPlot.html ODF
% visualization> is available.

h = [Miller(1,0,0,odf.CS),Miller(1,1,0,odf.CS),Miller(1,1,1,odf.CS)];
plotPDF(odf,h,'antipodal','silent')

%% Effect of halfwidth selection
%
% As mentioned above a proper halfwidth selection is crucial for ODF
% estimation. The following simple numerical experiment illustrates the
% dependency between the kernel halfwidth and the estimated error.
%
% Let's start with a model ODF and simulate some individual orientation data.

modelODF = fibreODF(Miller(1,1,1,crystalSymmetry('cubic')),xvector);
ori = discreteSample(modelODF,10000)

%%
% Next we define a list of kernel halfwidth ,

hw = [1*degree, 2*degree, 4*degree, 8*degree, 16*degree, 32*degree];

%%
% estimate for each halfwidth an ODF and compare it to the original ODF.

e = zeros(size(hw));
for i = 1:length(hw)

  odf = calcDensity(ori,'halfwidth',hw(i),'silent');
  e(i) = calcError(modelODF, odf);

end

%%
% After visualizing the estimation error we observe that its value is large
% either if we choose a very small or a very large halfwidth.
% In this specific example, the optimal halfwidth seems to be about 4
% degrees.

close all
plot(hw/degree,e)
xlabel('halfwidth in degree')
ylabel('esimation error')

author: mtex-toolbox

repo: mtex

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


text:
function dist = lines_exp_dist_3d_2 ( p1, p2, q1, q2 )

%*****************************************************************************80
%
%% LINES_EXP_DIST_3D_2 computes the distance between two explicit lines in 3D.
%
%  Discussion:
%
%    The explicit form of a line in 3D is:
%
%      the line through the points P1 and P2.
%
%    This routine uses a method that is essentially independent of dimension.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 August 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(3,1), P2(3,1), two points on the first line.
%
%    Input, real Q1(3,1), Q2(3,1), two points on the second line.
%
%    Output, real DIST, the distance between the lines.
%
  dim_num = 3;
%
%  Let U = (P2-P1) and V = (Q2-Q1) be the direction vectors on
%  the two lines.
%
  u(1:dim_num,1) = p2(1:dim_num,1) - p1(1:dim_num,1);
  v(1:dim_num,1) = q2(1:dim_num,1) - q1(1:dim_num,1);
%
%  Let SN be the unknown coordinate of the nearest point PN on line 1,
%  so that PN = P(SN) = P1 + SN * (P2-P1).
%
%  Let TN be the unknown coordinate of the nearest point QN on line 2,
%  so that QN = Q(TN) = Q1 + TN * (Q2-Q1).
%
%  Let W0 = (P1-Q1).
%
  w0(1:dim_num,1) = p1(1:dim_num,1) - q1(1:dim_num,1);
%
%  The vector direction WC = P(SN) - Q(TC) is unique (among directions)
%  perpendicular to both U and V, so
%
%    U dot WC = 0
%    V dot WC = 0
%
%  or, equivalently:
%
%    U dot ( P1 + SN * (P2 - P1) - Q1 - TN * (Q2 - Q1) ) = 0
%    V dot ( P1 + SN * (P2 - P1) - Q1 - TN * (Q2 - Q1) ) = 0
%
%  or, equivalently:
%
%    (u dot u ) * sn - (u dot v ) tc = -u * w0
%    (v dot u ) * sn - (v dot v ) tc = -v * w0
%
%  or, equivalently:
%
%   ( a  -b ) * ( sn ) = ( -d )
%   ( b  -c )   ( tc )   ( -e )
%
  a = u' * u;
  b = u' * v;
  c = v' * v;
  d = u' * w0;
  e = v' * w0;
%
%  Check the determinant.
%
  det = - a * c + b * b;

  if ( det == 0.0 )
    sn = 0.0;
    if ( abs ( b ) < abs ( c ) )
      tn = e / c;
    else
      tn = d / b;
    end
  else
    sn = ( c * d - b * e ) / det;
    tn = ( b * d - a * e ) / det;
  end

  pn(1:dim_num,1) = p1(1:dim_num,1) + sn * ( p2(1:dim_num,1) - p1(1:dim_num,1) );
  qn(1:dim_num,1) = q1(1:dim_num,1) + tn * ( q2(1:dim_num,1) - q1(1:dim_num,1) );

  dist = sqrt ( sum ( ( pn(1:dim_num,1) - qn(1:dim_num,1) ).^2 ) );

  return
end


text:
function v = haar_1d_inverse ( n, u )

%*****************************************************************************80
%
%% HAAR_1D_INVERSE inverts the Haar transform of a vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 March 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the length of the vector.
%    N must be a power of 2.
%
%    Input, real U(N,1), the vector to be transformed.
%
%    Output, real V(N,1), the transformed vector.
%
  v = u(:);

  s = sqrt ( 2.0 );

  w = zeros ( n, 1 );

  m = 1;

  while ( m * 2 <= n )

    w(1:2:2*m-1) = ( v(1:m) + v(1+m:m+m) ) / s;
    w(2:2:2*m)   = ( v(1:m) - v(1+m:m+m) ) / s;

    v(1:2*m) = w(1:2*m);
    m = m * 2;

  end

  return
end