text:
function polygon_triangulate_test ( )

%*****************************************************************************80
%
%% POLYGON_TRIANGULATE_TEST tests the POLYGON_TRIANGULATE library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 May 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYGON_TRIANGULATE_TEST\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the POLYGON_TRIANGULATE library.\n' );

  polygon_triangulate_test01 ( );

  polygon_triangulate_test02 ( 'comb' );
  polygon_triangulate_test02 ( 'hand' );
  polygon_triangulate_test02 ( 'i18' );

  polygon_triangulate_test03 ( 'comb' );
  polygon_triangulate_test03 ( 'hand' );
  polygon_triangulate_test03 ( 'i18' );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYGON_TRIANGULATE_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function polygon_triangulate_test01 ( )

%*****************************************************************************80
%
%% POLYGON_TRIANGULATE_TEST01 tests the "comb_10" polygon.
%
%  Discussion:
%
%    There are N-3 triangles in the triangulation.
%
%    For the first N-2 triangles, the first edge listed is always an
%    internal diagonal.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 May 2014
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYGON_TRIANGULATE_TEST01\n' );
  fprintf ( 1, '  Triangulate the comb_10 polygon.\n' );

  x = [ ...
    8; ...
    7; ...
    6; ...
    5; ...
    4; ...
    3; ...
    2; ...
    1; ...
    0; ...
    4];

  y = [ ...
    0; ...
   10; ...
    0; ...
   10; ...
    0; ...
   10; ...
    0; ...
   10; ...
    0; ...
   -2 ];

  n = 10;

  triangles = polygon_triangulate ( n, x, y );

  i4mat_print ( 3, n - 2, triangles, '  Triangles' );

  return
end
function polygon_triangulate_test02 ( prefix )

%*****************************************************************************80
%
%% POLYGON_TRIANGULATE_TEST02 triangulates a polygon described in a file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 May 2014
%
%  Author:
%
%    John Burkardt.
%

%
%  Create filenames.
%
  node_filename = strcat ( prefix, '_nodes.txt' );
  element_filename = strcat ( prefix, '_elements.txt' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYGON_TRIANGULATE_TEST02\n' );
  fprintf ( 1, '  Read polygon coordinates in "%s"\n', node_filename );
%
%  Read the node coordinates.
%
  [ ~, n ] = r8mat_header_read ( node_filename );

  xy = r8mat_data_read ( node_filename, 2, n );
%
%  Get the triangulation.
%
  triangle_num = n - 2;
  triangles = polygon_triangulate ( n, xy(1,:), xy(2,:) );
%
%  Write the triangulation to a file.
%
  i4mat_write ( element_filename, 3, triangle_num, triangles );

  fprintf ( 1, '  Write triangulation to "%s"\n', element_filename );

  return
end
function polygon_triangulate_test03 ( prefix )

%*****************************************************************************80
%
%% POLYGON_TRIANGULATE_TEST03 plots a triangulation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 May 2014
%
%  Author:
%
%    John Burkardt.
%

%
%  Create filenames.
%
  node_filename = strcat ( prefix, '_nodes.txt' );
  plot_filename = strcat ( prefix, '.png' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYGON_TRIANGULATE_TEST03\n' );
  fprintf ( 1, '  Read polygon coordinates in "%s"\n', node_filename );
%
%  Read the node coordinates.
%
  [ ~, n ] = r8mat_header_read ( node_filename );

  xy = r8mat_data_read ( node_filename, 2, n );
%
%  Get the triangulation.
%
  triangle_num = n - 2;
  triangles = polygon_triangulate ( n, xy(1,:), xy(2,:) );
%
%  Plot the triangulation.
%
  figure ( 1 )
  clf
  grid on
  hold on

  plot ( xy(1,1:n), xy(2,1:n), 'k.', 'MarkerSize', 25 )

  for i = 1 : n

    ip1 = mod ( i, n ) + 1;

    line ( [ xy(1,i), xy(1,ip1) ], ...
           [ xy(2,i), xy(2,ip1) ], ...
           'LineWidth', 3, 'Color', 'b' )

  end

  for i = 1 : n - 3
    i1 = triangles(1,i);
    i2 = triangles(2,i);
    line ( [ xy(1,i1), xy(1,i2) ], ...
           [ xy(2,i1), xy(2,i2) ], ...
           'LineWidth', 3, 'Color', 'r' )
  end

  xlabel ( '<-- X -->', 'Fontsize', 16 )
  ylabel ( '<-- Y -->', 'Fontsize', 16 )
  title ( prefix, 'Fontsize', 24 );
  axis equal

  hold off

  print ( '-dpng', plot_filename );
  fprintf ( 1, '  Plot triangulation, save to graphics file "%s"\n', plot_filename );

  return
end
function column_num = file_column_count ( input_file_name )

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

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

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( s_len_trim ( line ) == 0 )

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

    else
      got_one = TRUE;
      break;
    end

  end

  fclose ( input_unit );

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

  column_num = s_word_count ( line );

  return
end
function row_num = file_row_count ( input_file_name )

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

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

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

  record_num = 0;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

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

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

  end

  fclose ( input_unit );

  return
end
function i4mat_print ( m, n, a, title )

%*****************************************************************************80
%
%% I4MAT_PRINT prints an I4MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 September 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of the matrix.
%
%    Input, integer A(M,N), an M by N matrix to be printed.
%
%    Input, string TITLE, a title.
%
  i4mat_print_some ( m, n, a, 1, 1, m, n, title );

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

%*****************************************************************************80
%
%% I4MAT_PRINT_SOME prints out a portion of an I4MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of the matrix.
%
%    Input, integer A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, a title.
%
  incx = 10;

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

  if ( m <= 0 || n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (None)\n' );
    return
  end

  for j2lo = max ( jlo, 1 ): incx : min ( jhi, n )

    j2hi = j2lo + incx - 1;
    j2hi = min ( j2hi, n );
    j2hi = min ( j2hi, jhi );

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

    for j = j2lo : j2hi
      fprintf ( 1, '%7d  ', j );
    end

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

    i2lo = max ( ilo, 1 );
    i2hi = min ( ihi, m );

    for i = i2lo : i2hi

      fprintf ( 1, '%5d: ', i );
      for j = j2lo : j2hi
        fprintf ( 1, '%7d  ', a(i,j) );
      end

      fprintf ( 1, '\n' );

    end

  end

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

%*****************************************************************************80
%
%% I4MAT_WRITE writes an I4MAT file.
%
%  Discussion:
%
%    An I4MAT is an array of I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 June 2010
%
%  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.
%

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

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

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

%*****************************************************************************80
%
%% R8MAT_DATA_READ reads data from an R8MAT file.
%
%  Discussion:
%
%    An R8MAT is an array of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns of data.
%
%    Output, real TABLE(M,N), the point coordinates.
%

%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

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

  input_unit = fopen ( input_filename );

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

  table = zeros(m,n);

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      break;
    end

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

    elseif ( s_len_trim ( line ) == 0 )

    else

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

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

    end

  end

  fclose ( input_unit );

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

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

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

  n = file_row_count ( input_filename );

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

  return
end
function len = s_len_trim ( s )

%*****************************************************************************80
%
% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer 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 word_num = s_word_count ( s )

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

  word_num = 0;
  s_length = length ( s );

  if ( s_length <= 0 )
    return;
  end

  blank = TRUE;

  for i = 1 : s_length

    if ( s(i) == ' ' )
      blank = TRUE;
    elseif ( blank == TRUE )
      word_num = word_num + 1;
      blank = FALSE;
    end

  end

  return
end
function timestamp ( )

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

  return
end

author: johannesgerer

repo: jburkardt-m


text:
lambda = 5;
limits = [0 10; 0 10];

clutter_model = PoissonRateUniformPositionX('ClutterRate',lambda,'Limits',limits);

lambda_samples = clutter_model.random(1000,'cardinality');
spatial_samples = clutter_model.random();

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function varargout = transformPoint3d(varargin)
%TRANSFORMPOINT3D Transform a point with a 3D affine transform
%
%   PT2 = transformPoint3d(PT1, TRANS);
%   PT2 = transformPoint3d(X1, Y1, Z1, TRANS);
%   where PT1 has the form [xp yp zp], and TRANS is a [3x3], [3x4], [4x4]
%   matrix, return the point transformed according to the affine transform
%   specified by TRANS.
%
%   Format of TRANS is a 4-by-4 matrix.
%
%   The function accepts transforms given using the following formats:
%   [a b c]   ,   [a b c j] , or [a b c j]
%   [d e f]       [d e f k]      [d e f k]
%   [g h i]       [g h i l]      [g h i l]
%                                [0 0 0 1]
%
%   PT2 = transformPoint3d(PT1, TRANS)
%   also work when PT1 is a [Nx3xMxPxETC] array of double. In this case,
%   PT2 has the same size as PT1.
%
%   PT2 = transformPoint3d(X1, Y1, Z1, TRANS);
%   also work when X1, Y1 and Z1 are 3 arrays with the same size. In this
%   case, PT2 will be a 1-by-3 cell containing {X Y Z} outputs of size(X1).
%
%   [X2 Y2 Z2] = transformPoint3d(...);
%   returns the result in 3 different arrays the same size as the input.
%   This form can be useful when used with functions like meshgrid or warp.
%
%   MESH2 = transformPoint3d(MESH, TRANS) transforms the field 'vertices'
%   of the struct MESH and returns the same struct with the transformed
%   vertices
%
%   See also:
%   points3d, transforms3d, translation3d, meshgrid
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 10/02/2005.
%

%   23/03/2006 add support for non vector point data
%   26/10/2006 better support for large data handling: iterate on points
%       in the case of a memory lack.
%   20/04/2007 add link to rotationXX functions
%   29/09/2010 fix bug in catch case
%   12/03/2011 slightly reduce memory usage

% parse input arguments
meshStructSwitch = false;
if length(varargin) == 2
    if isstruct(varargin{1}) && isfield(varargin{1}, 'vertices')
        % If first argument is a struct with the field 'vertices' the
        % output will be the same struct, but with the transformed vertices
        meshStructSwitch = true;
        meshStruct = varargin{1};
        varargin{1}=varargin{1}.vertices;
    end
    % Point coordinates are given in a single N-by-3-by-M-by-etc argument.
    % Preallocate x, y, and z to size N-by-1-by-M-by-etc, then fill them in
    dim = size(varargin{1});
    dim(2) = 1;
    [x, y, z] = deal(zeros(dim, class(varargin{1})));
    x(:) = varargin{1}(:,1,:);
    y(:) = varargin{1}(:,2,:);
    z(:) = varargin{1}(:,3,:);
    trans  = varargin{2};
elseif length(varargin) == 4
    % Point coordinates are given in 3 different arrays
    x = varargin{1};
    y = varargin{2};
    z = varargin{3};
    dim = size(x);
    trans = varargin{4};
end

% eventually add null translation
if size(trans, 2) == 3
    trans = [trans zeros(size(trans, 1), 1)];
end

% eventually add normalization
if size(trans, 1) == 3
    trans = [trans;0 0 0 1];
end

% convert coordinates
NP  = numel(x);
try
    % vectorial processing, if there is enough memory
    %res = (trans*[x(:) y(:) z(:) ones(NP, 1)]')';
    %res = [x(:) y(:) z(:) ones(NP, 1)]*trans';
    res = [x(:) y(:) z(:) ones(NP,1,class(x))] * trans';

    % Back-fill x,y,z with new result (saves calling costly reshape())
    x(:) = res(:,1);
    y(:) = res(:,2);
    z(:) = res(:,3);
catch ME
    disp(ME.message)
    % process each point one by one, writing in existing array
    for i = 1:NP
        res = [x(i) y(i) z(i) 1] * trans';
        x(i) = res(1);
        y(i) = res(2);
        z(i) = res(3);
    end
end

% process output arguments
if nargout <= 1
    % results are stored in a unique array
    if length(dim) > 2 && dim(2) > 1
        warning('geom3d:shapeMismatch',...
            'Shape mismatch: Non-vector xyz input should have multiple x,y,z output arguments. Cell {x,y,z} returned instead.')
        varargout{1} = {x,y,z};
    else
        if meshStructSwitch
            meshStruct.vertices = [x y z];
            varargout{1}=meshStruct;
        else
            varargout{1} = [x y z];
        end
    end

elseif nargout == 3
    varargout{1} = x;
    varargout{2} = y;
    varargout{3} = z;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [neglogli, dL, H] = neglogli_poissGLM(prs,XX,YY,dtbin)
% [neglogli, dL, H] = Loss_GLM_logli_exp(prs,XX);
%
% Compute negative log-likelihood of data undr Poisson GLM model with
% exponential nonlinearity
%
% Inputs:
%   prs [d x 1] - parameter vector
%    XX [T x d] - design matrix
%    YY [T x 1] - response (spike count per time bin)
% dtbin [1 x 1] - time bin size used
%
% Outputs:
%   neglogli   = negative log likelihood of spike train
%   dL [d x 1] = gradient
%   H  [d x d] = Hessian (second deriv matrix)

% Compute GLM filter output and condititional intensity
vv = XX*prs; % filter output
rr = exp(vv)*dtbin; % conditional intensity (per bin)

% ---------  Compute log-likelihood -----------
Trm1 = -vv'*YY; % spike term from Poisson log-likelihood
Trm0 = sum(rr);  % non-spike term
neglogli = Trm1 + Trm0;

% ---------  Compute Gradient -----------------
if (nargout > 1)
    dL1 = -XX'*YY; % spiking term (the spike-triggered average)
    dL0 = XX'*rr; % non-spiking term
    dL = dL1+dL0;
end

% ---------  Compute Hessian -------------------
if nargout > 2
    H = XX'*bsxfun(@times,XX,rr); % non-spiking term
end

author: pillowlab

repo: GLMspiketraintutorial


text:
function [f_constr,label,auxdata] = getFrictionCone(obj, f, fric_coef)
    % returns the symbolic expression of the friction cone
    % constraints of the contact
    %
    % @note The friction cone constraints does not includes the
    % ZMP/CWS constraints.
    %
    % Parameters:
    %  f: the SymVariable of the constraint wrenches
    %  @type SymVariable
    %  fric_coef: the coefficient of friction. @type struct
    %
    % Return values:
    %  f_constr: symbolic function of the friction cone
    %  @type SymFunction
    %  label: the label of the constraints @type cellstr
    %  auxdata: constant data used in the constraint function
    %  @type cell
    %
    % Optional fields of fric_coef:
    %  mu: the (static) coefficient of friction. @type double
    %  gamma: the coefficient of torsional friction @type double

    assert(size(obj.WrenchBase,2) == length(f),...
        ['The dimension of the constraint wrenchs is incorrect.\n',...
        'Expected %d, instead %d'], size(obj.WrenchBase,2), length(f));

    mu = SymVariable('mu');
    gamma = SymVariable('gamma');

    fun_name = ['u_friction_cone_', obj.Name];
    switch obj.Type
        case 'PlanarLineContactWithFriction'
            % x, y, z
            constr = [f(2); % fz >= 0
                f(1) + (mu/sqrt(2))*f(2);
                -f(1) + (mu/sqrt(2))*f(2)];

            % create a symbolic function object
            f_constr = SymFunction(fun_name,...
                constr,{f},{mu});

            % create the label text
            label = {'normal_force';
                'friction_x_pos';
                'friction_x_neg';
                };

            % validate the provided static friction coefficient
            validateattributes(fric_coef.mu,{'double'},...
                {'scalar','real','>=',0},...
                'ContactFrame.getFrictionCone','mu');
            auxdata = {fric_coef.mu};
        case 'PlanarPointContactWithFriction'
            % x, y, z
            constr = [f(2); % fz >= 0
                f(1) + (mu/sqrt(2))*f(2);
                -f(1) + (mu/sqrt(2))*f(2)];

            % create a symbolic function object
            f_constr = SymFunction(fun_name,...
                constr,{f},{mu});

            % create the label text
            label = {'normal_force';
                'friction_x_pos';
                'friction_x_neg';
                };

            % validate the provided static friction coefficient
            validateattributes(fric_coef.mu,{'double'},...
                {'scalar','real','>=',0},...
                'ContactFrame.getFrictionCone','mu');
            auxdata = {fric_coef.mu};
        case 'PointContactWithFriction'
            % x, y, z
            constr = [f(3); % fz >= 0
                f(1) + (mu/sqrt(2))*f(3);
                -f(1) + (mu/sqrt(2))*f(3);
                f(2) + (mu/sqrt(2))*f(3);
                -f(2) + (mu/sqrt(2))*f(3)];
%                 f(4) + gamma * f(3);       % -gamma * fz < wz
%                 -f(4) + gamma * f(3)];     % wz < gamma * fz

            % create a symbolic function object
            f_constr = SymFunction(fun_name,...
                constr,{f},{[mu]});

            % create the label text
            label = {'normal_force';
                'friction_x_pos';
                'friction_x_neg';
                'friction_y_pos';
                'friction_y_neg'};
%                 'tor_firction_neg';
%                 'tor_firction_pos';
%                 };

            % validate the provided static friction coefficient
            validateattributes(fric_coef.mu,{'double'},...
                {'scalar','real','>=',0},...
                'ContactFrame.getFrictionCone','mu');
            % validate the provided torsional friction coefficient
            %             validateattributes(fric_coef.gamma,{'double'},...
            %                 {'scalar','real','>=',0},...
            %                 'ContactFrame.getFrictionCone','gamma');
            auxdata = [fric_coef.mu];
        case 'PointContactWithoutFriction'
            % z
            constr = f;

            % create a symbolic function object
            f_constr = SymFunction(fun_name,...
                constr,{f});

            % create the label text
            label = {'normal_force';
                };
            auxdata = [];
        case 'LineContactWithFriction'
            % x, y, z, roll, yaw
            constr = [f(3); % fz >= 0
                f(1) + (mu/sqrt(2))*f(3);  % -mu/sqrt(2) * fz < fx
                -f(1) + (mu/sqrt(2))*f(3); % fx < mu/sqrt(2) * fz
                f(2) + (mu/sqrt(2))*f(3);  % -mu/sqrt(2) * fz < fu
                -f(2) + (mu/sqrt(2))*f(3); % fy < mu/sqrt(2) * fz
                f(5) + gamma * f(3);       % -gamma * fz < wy
                -f(5) + gamma * f(3)];     % wy < gamma * fz

            % create a symbolic function object
            f_constr = SymFunction(fun_name,...
                constr,{f},{[mu;gamma]});

            % create the label text
            label = {'normal_force';
                'friction_x_pos';
                'friction_x_neg';
                'friction_y_pos';
                'friction_y_neg';
                'tor_firction_neg';
                'tor_firction_pos';
                };

            % validate the provided static friction coefficient
            validateattributes(fric_coef.mu,{'double'},...
                {'scalar','real','>=',0},...
                'ContactFrame.getFrictionCone','mu');

            % validate the provided torsional friction coefficient
            validateattributes(fric_coef.gamma,{'double'},...
                {'scalar','real','>=',0},...
                'ContactFrame.getFrictionCone','gamma');
            auxdata = [fric_coef.mu; fric_coef.gamma];
        case 'LineContactWithoutFriction'
            % z, roll
            constr = f;
            % create a symbolic function object
            f_constr = SymFunction(fun_name,...
                constr,{f});

            % create the label text
            label = {'normal_force';
                };
            auxdata = [];
        case 'PlanarContactWithFriction'
            % x, y, z, roll, pitch, yaw
            constr = [f(3); % fz >= 0
                f(1) + (mu/sqrt(2))*f(3);  % -mu/sqrt(2) * fz < fx
                -f(1) + (mu/sqrt(2))*f(3); % fx < mu/sqrt(2) * fz
                f(2) + (mu/sqrt(2))*f(3);  % -mu/sqrt(2) * fz < fu
                -f(2) + (mu/sqrt(2))*f(3); % fy < mu/sqrt(2) * fz
                f(6) + gamma * f(3);       % -gamma * fz < wz
                -f(6) + gamma * f(3)];     % wz < gamma * fz

            % create a symbolic function object
            f_constr = SymFunction(fun_name,...
                constr,{f},{[mu;gamma]});

            % create the label text
            label = {'normal_force';
                'friction_x_pos';
                'friction_x_neg';
                'friction_y_pos';
                'friction_y_neg';
                'tor_firction_neg';
                'tor_firction_pos';
                };

            % validate the provided static friction coefficient
            validateattributes(fric_coef.mu,{'double'},...
                {'scalar','real','>=',0},...
                'ContactFrame.getFrictionCone','mu');

            % validate the provided torsional friction coefficient
            validateattributes(fric_coef.gamma,{'double'},...
                {'scalar','real','>=',0},...
                'ContactFrame.getFrictionCone','gamma');
            auxdata = [fric_coef.mu; fric_coef.gamma];
        case 'PlanarContactWithoutFriction'
            % z, roll, pitch,
            constr = f;
            % create a symbolic function object
            f_constr = SymFunction(fun_name,...
                constr,{f});

            % create the label text
            label = {'normal_force';
                };
            auxdata = [];
    end

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function oev = getoe(ioev)

% interactive request of classical orbital elements

% NOTE: all angular elements are returned in radians

% input

%  ioev = request array (1 = yes, 0 = no)

% output

%  oev(1) = semimajor
%  oev(2) = orbital eccentricity
%  oev(3) = orbital inclination
%  oev(4) = argument of perigee
%  oev(5) = right ascension of the ascending node (ioev(5) = 1)
%  oev(5) = east longitude of the ascending node (ioev(5) = 2)
%  oev(6) = true anomaly

% Orbital Mechanics with MATLAB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dtr = pi / 180;

oev = zeros(6,1);

if (ioev(1) == 1)
    % semimajor axis (kilometers)
    % (semimajor axis > req)

    for itry = 1:1:5
        fprintf('\nplease input the semimajor axis (kilometers)');
        fprintf('\n(semimajor axis > 0)\n');

        sma = input('? ');

        if (sma > 0)
            break;
        end
    end
    oev(1) = sma;
end

if (ioev(2) == 1)
    % orbital eccentricity (non-dimensional)
    % (0 <= eccentricity < 1)

    for itry = 1:1:5
        fprintf('\nplease input the orbital eccentricity (non-dimensional)');
        fprintf('\n(0 <= eccentricity < 1)\n');

        ecc = input('? ');

        if (ecc >= 0 && ecc < 1)
            break;
        end
    end
    oev(2) = ecc;
end

if (ioev(3) == 1)
    % orbital inclination (degrees)
    % (0 <= inclination <= 180)

    for itry = 1:1:5
        fprintf('\nplease input the orbital inclination (degrees)');
        fprintf('\n(0 <= inclination <= 180)\n');

        orbinc = input('? ');

        if (orbinc >= 0 && orbinc <= 180)
            break;
        end
    end
    oev(3) = dtr * orbinc;
end

if (ioev(4) == 1)
    % argument of perigee (degrees)
    % (0 <= argument of perigee <= 360)

    if (ecc > 0)

        for itry = 1:1:5

            fprintf('\nplease input the argument of perigee (degrees)');
            fprintf('\n(0 <= argument of perigee <= 360)\n');

            argper = input('? ');

            if (argper >= 0 && argper <= 360)
                break;
            end
        end
    else
        argper = 0;
    end
    oev(4) = dtr * argper;
end

if (ioev(5) >= 1)

    if (ioev(5) == 1)
        % raan (degrees)
        % (0 <= raan <= 360)

        if (oev(3) > 0)
            for itry = 1:1:5
                fprintf('\nplease input the right ascension of the ascending node (degrees)');
                fprintf('\n(0 <= raan <= 360)\n');

                raan = input('? ');

                if (raan >= 0 && raan <= 360)
                    break;
                end
            end
            oev(5) = dtr * raan;
        else
            oev(5) = 0;
        end
    end

    if (ioev(5) == 2)
        % elan (degrees)
        % (0 <= elan <= 360)

        if (oev(3) > 0)
            for itry = 1:1:5
                fprintf('\nplease input the east longitude of the ascending node (degrees)');
                fprintf('\n(0 <= elan <= 360)\n');

                elan = input('? ');

                if (elan >= 0 && elan <= 360)
                    break;
                end
            end
            oev(5) = dtr * elan;
        else
            oev(5) = 0;
        end
    end
end

if (ioev(6) == 1)
    % true anomaly (degrees)
    % (0 <= true anomaly <= 360)

    for itry = 1:1:5
        fprintf('\nplease input the true anomaly (degrees)');
        fprintf('\n(0 <= true anomaly <= 360)\n');

        tanom = input('? ');

        if (tanom >= 0 && tanom <= 360)
            break;
        end
    end
    oev(6) = dtr * tanom;
end

author: Sable

repo: mcbench-benchmarks


text:
% simple demo to convert a polygon mesh into a voxel representation

% code is 99.99% based on
% [1] http://www.mathworks.com/matlabcentral/fileexchange/24086-polygon2voxel
% [2] http://www.mathworks.com/matlabcentral/fileexchange/21044-3d-voxelizer

%   % Compile the c-coded function
%   mex polygon2voxel_double.c -v

load model;

vertices = vertices - repmat(mean(vertices,1),size(vertices,1),1);

FV.faces = faces;
FV.vertices = vertices;

Volume=polygon2voxel(FV,[100 100 100],'auto');

%% visualization 1
figure
[X,Y,Z]=ind2sub(size(Volume),find(Volume(:)));
plot3(X,Y,Z,'.');
axis equal;
xlabel('x');
ylabel('y');
zlabel('z');

%% visualization 2
% 3d pirnter style visualization to add layer by layer
plot3D(Volume,'timed', 0.1)

figure
%% visualization 3
for i=1:size(Volume,1)
    imagesc(squeeze(Volume(i,:,:)));
    axis equal;
    axis tight;
    axis off
    title(i);
    pause(0.1);
end

for i=1:size(Volume,2)
    imagesc(squeeze(Volume(:,i,:)));
    axis equal;
    axis tight;
    axis off
    title(i);
    pause(0.1);
end

for i=1:size(Volume,3)
    imagesc(squeeze(Volume(:,:,i)));
    axis equal;
    axis tight;
    axis off
    title(i);
    pause(0.1);
end

author: jianxiongxiao

repo: ProfXkit


text:
function words = ExtractWordsFromString(inputString)

inWord = 0;
words = {};

for k=1:length(inputString)
    if inputString(k) == ' '
        isSpace = 1;
    elseif ~isempty( regexp(inputString(k), '\s', 'once') )
        isSpace = 1;
    else
        isSpace = 0;
    end

    if inWord == 1
        if isSpace == 1
            words{end+1} = inputString(word_start:k-1);
            inWord = 0;
        end

    else
        if isSpace == 0
            inWord = 1;
            word_start = k;
        end
    end
end

% if the last symbol is not space and we are still in a word.
if inWord==1
    words{end+1} = inputString(word_start:end);
end

author: singaxiong

repo: SignalGraph


text:
function [ s ] = kpath_length( k, A, mode )
%KPATH_LENGTH Find the distance along a path in k-space.
%   [ s ] = kpath_length( k, A, mode ) calculates the legnth along a series
%   of N points in k-space whose fractional coordinates are given by the
%   3xN array k. This is useful for creating band structure plots. The 3x3
%   matrix A gives the lattice vectors (mode = 'direct') or the reciprocal
%   lattice vectors (mode = 'reciprocal'). If mode is not specified, the
%   default is reciprocal.

    if nargin ==3
        switch mode
            case 'direct'
                B = reciprocal_lattice(A);
            case 'reciprocal'
                B = A;
        end
    end

    k = k*B; % convert to Cartesian coordinates

    s = zeros(1,size(k,1));

    for i = 1:size(k,1)-1
        d = sqrt(sum((k(i,:)-k(i+1,:)).^2));
        s(i+1) = s(i) + d;
    end

end


text:
function f = times(f, g, varargin)
%.*   TRIGTECH multiplication.
%   F.*G multiplies TRIGTECH objects F and G or a TRIGTECH by a scalar if either
%   F or G is a scalar.
%
%   If F is an array-valued TRIGTECH, then F.*C is supported if C is a row
%   vector of doubles with the same number of columns as F.
%
% See also MTIMES, RDIVIDE.

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

% TRIGTECH * [] = []:
if ( isempty(f) || isempty(g) )
    f = [];
    return
end

if ( ~isa(f, 'trigtech') )      % Ensure F is a TRIGTECH.

    f = times(g, f, varargin{:});
    return

elseif ( isa(g, 'double') )     % TRIGTECH .* double.

    % Do the multiplication:
    if ( size(g, 2) > 1 )
        f.values = bsxfun(@times, f.values, g);
        f.coeffs = f.vals2coeffs(f.values);
    else
        f.values = f.values*g;
        f.coeffs = f.coeffs*g;
    end
    f.isReal = f.isReal & isreal(g);
    return

elseif ( ~isa(f, 'trigtech') || ~isa(g, 'trigtech') )
    % Don't know how to do the operation.

    error('CHEBFUN:TRIGTECH:times:typeMismatch', ...
        ['Incompatible operation between objects.\n', ...
         'Make sure functions are of the same type.']);

elseif ( size(f.values, 1) == 1 )
    % If we have (constant TRIGTECH).*TRIGTECH, reverse the order and call TIMES
    % again:
    f = times(g, f.values);
    return

elseif ( size(g.values, 1) == 1)
    % If we have TRIGTECH.*(constant TRIGTECH), convert the (constant TRIGTECH)
    % to a scalar and call TIMES again:
    f = times(f, g.values);
    return
end

% Get the size of each TRIGTECH:
[fn, fm] = size(f.values);
[gn, gm] = size(g.values);

% The length of the product is known:
fNew = prolong(f, fn + gn - 1);

% Check dimensions:
if ( fm ~= gm )
    if ( fm == 1 )
        % Allow [Inf x 1] .* [Inf x m].
        fNew.values = repmat(fNew.values, 1, gm);
        fNew.coeffs = repmat(fNew.coeffs, 1, gm);
    elseif ( gm == 1 )
        % Allow [Inf x m] .* [Inf x 1].
        g.values = repmat(g.values, 1, fm);
        g.coeffs = repmat(g.coeffs, 1, fm);
    else
        error('CHEBFUN:TRIGTECH:times:dim2', ...
            'Inner matrix dimensions must agree.');
    end
end

% Check for two cases where the output is known in advance to be positive,
% namely F == conj(G) or F == G and isreal(F).
pos = false;

% Multiply values:
if ( isequal(f, g) )
   values = fNew.values.^2;
   if ( isreal(f) )
       pos = true;
   end
elseif ( isequal(conj(f), g) )
   values = conj(fNew.values).*fNew.values;
   pos = true;
else
   gNew = prolong(g, fn + gn - 1);
   values = fNew.values.*gNew.values;
end

% Assign values and coefficients back to f:
f.values = values;
f.coeffs = f.vals2coeffs(values);

% Update ishappy:
f.ishappy = f.ishappy && g.ishappy;

% Simplify.
f = simplify(f);

if ( pos )
    % Here we know that the product of F and G should be positive. However,
    % SIMPLIFY may have destroyed this property, so we enforce it.
    f.values = abs(f.values);
    f.coeffs = f.vals2coeffs(f.values);
    f.isReal = true(1, size(f.coeffs, 2));
else
    f.isReal = f.isReal & g.isReal;
end

% If f and g are real then make the result real.
f.values(:,f.isReal) = real(f.values(:,f.isReal));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [vhat,iteration]=decode_ldpc_matlab(rx_waveform,No,h,rows,cols,ind,r,c,max_iter)

% bit to check
% check to bit
% update APP LLR
% hard decision
% if valid codework, then exit

vhat(1,1:cols)=0;   %length of decoded word = number of colums of H

s=struct('a',sparse(rows,cols,0),'b',sparse(rows,cols,0),'g',sparse(rows,cols,0));
% a:check to bit messages of previous iteration, b: current check to bit messages,
% g:updated codebit APP LLR

gamma_n=(4/No)*rx_waveform;

s.g(ind)=gamma_n(c);

bitmessage=sparse(rows,cols,0);
bitmessage_2=sparse(rows,cols,0);

for iteration=1:max_iter       %iteration loop

    bitmessage(ind)=tanh( ( -(s.g(ind))+ s.a(ind))  /2);    %bit to check messages

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %loop for computing check-to-bit messages
    for i=1:rows
        cc=c(find(r==i));   %
        rowtemp=bitmessage(i,:);    %ith row of bitmessage

        rowtemp_int=zeros(1,length(cc));
        rowtemp_int=rowtemp(cc);    %interesting elements corresponding to 1's in H
        ze=find(rowtemp_int==0);    %number of zeros in interesting elements

        if length(ze)>=2    %if #of zeros >=2, all interesting elements are zero
            bitmessage_2(i,cc)=0;

        elseif isempty(ze)==1
            rowtemp_int1=zeros(1,length(cc));
            rowtemp_int1(1:length(cc))=prod(rowtemp_int);
            bitmessage_2(i,cc)=rowtemp_int1./rowtemp_int;

        else  %one zero in interesting elements
            zuss=find(rowtemp_int==0);
            rowtemp_int2=rowtemp_int;
            rowtemp_int2(zuss)=1;
            rowtemp_int3=zeros(1,length(cc));
            rowtemp_int3(zuss)=prod(rowtemp_int2);
            bitmessage_2(i,cc)=rowtemp_int3;
        end

    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    s.b(ind)=-2*atanh(bitmessage_2(ind));   %check to bit messages

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    s.b(find(s.b==inf))=1000;    % limiting the checknode to bitnode LLR values to +-50
    s.b(find(s.b==-inf))=-1000;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    s.a(ind)=s.b(ind);                      %current check to bit messages for use in next iteration

    sum_of_b=sum(s.b,1);

    tempsum=sum_of_b+gamma_n;               % updating APP LLR
    s.g(ind)=tempsum(c);

    vhat(find(tempsum>=0))=1;
    vhat(find(tempsum<0))=0;

    if mod(vhat*h',2)==0
        break
    end

end      %iteration loop -end

%   Copyright (c) 2005 by Shaikh Faisal Zaheer, faisal.zaheer@gmail.com
%   $Revision: 1.0 $  $Date: 10/11/2005 $


text:
function trans = createScaling(varargin)
%CREATESCALING Create the 3*3 matrix of a scaling in 2 dimensions.
%
%   TRANS = createScaling(SX, SY);
%   return the matrix corresponding to scaling by SX and SY in the 2
%   main directions.
%   The returned matrix has the form:
%   [SX  0  0]
%   [0  SY  0]
%   [0   0  1]
%
%   TRANS = createScaling(SX);
%   Assume SX and SY are equals.
%
%   TRANS = createScaling(CENTER, ...);
%   Specifies the center of the scaling transform. The argument CENTER
%   should be a 1-by-2 array representing coordinates of center.
%
%   See also
%   transforms2d, transformPoint, createTranslation, createRotation
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2004-04-07
% Copyright 2004-2022 INRA - TPV URPOI - BIA IMASTE

% defined default arguments
sx = 1;
sy = 1;
cx = 0;
cy = 0;

% process input arguments
if nargin == 1
    % the argument is either the scaling factor in both direction,
    % or a 1x2 array containing scaling factor in each direction
    var = varargin{1};
    sx = var(1);
    sy = var(1);
    if length(var)>1
        sy = var(2);
    end
elseif nargin == 2
    % the 2 arguments are the scaling factors in each dimension
    sx = varargin{1};
    sy = varargin{2};
elseif nargin == 3
    % first argument is center, 2nd and 3d are scaling factors
    center = varargin{1};
    cx = center(1);
    cy = center(2);
    sx = varargin{2};
    sy = varargin{3};
end

% concatenate results in a 3-by-3 matrix
trans = [sx 0 cx*(1-sx); 0 sy cy*(1-sy); 0 0 1];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function diagnostic = solvesdp_multiple(varargin);

yalmiptime = clock;
h = varargin{2};
varargin{2} = sum(recover(depends(h)));

if ~is(h,'linear')
    error('Parts of your matrix objective is not linear (multiple solutions can currently only be obtained for linear objectives)');
end

if is(h,'complex')
    error('Parts of your matrix objective is complex-valued (which makes no sense since complex numbers have no natural ordering');
end

if nargin<3
    ops = sdpsettings('saveyalmipmodel',1);
    varargin{3} = ops;
else
    ops = varargin{3};
    varargin{3}.saveyalmipmodel = 1;
end
varargin{3}.pureexport = 1;
yalmiptime = clock;
model = solvesdp(varargin{:});
diagnostic.yalmiptime = etime(clock,yalmiptime);
yalmiptime = clock;

h_variables = getvariables(h);
h_base = getbase(h);
for i = 1:length(h)
    model.f = h_base(i,1);
    model.c = mapObjective(h_variables,model.used_variables,h_base(i,2:end));

    % *************************************************************************
    % TRY TO SOLVE PROBLEM
    % *************************************************************************
    if ops.debug
        eval(['output = ' model.solver.call '(model);']);
    else
        try
            eval(['output = ' model.solver.call '(model);']);
        catch
            output.Primal = zeros(length(model.c),1)+NaN;
            output.Dual  = [];
            output.Slack = [];
            output.solvertime   = nan;
            output.solverinput  = [];
            output.solveroutput = [];
            output.problem = 9;
            output.infostr = yalmiperror(output.problem,lasterr);
        end
    end

    if ops.dimacs
        try
            b = -model.c;
            c = model.F_struc(:,1);
            A = -model.F_struc(:,2:end)';
            x = output.Dual;
            y = output.Primal;
            % FIX this nonlinear crap (return variable type in
            % compilemodel)
            if options.relax == 0 & any(full(sum(model.monomtable,2)~=0))
                if ~isempty(find(sum(model.monomtable | model.monomtable,2)>1))
                    z=real(exp(model.monomtable*log(y+eps)));
                    y = z;
                end
            end

            if isfield(output,'Slack')
                s = output.Slack;
            else
                s = [];
            end

            dimacs = computedimacs(b,c,A,x,y,s,model.K);
        catch
            dimacs = [nan nan nan nan nan nan];
        end
    else
        dimacs = [nan nan nan nan nan nan];
    end

    % ********************************
    % ORIGINAL COORDINATES
    % ********************************
    if isempty(output.Primal)
       output.Primal = zeros(size(model.recoverdata.H,2),1);
    end
    output.Primal(:,i) = model.recoverdata.x_equ+model.recoverdata.H*output.Primal;

    % ********************************
    % OUTPUT
    % ********************************
    diagnostic.yalmiptime = diagnostic.yalmiptime + etime(clock,yalmiptime)-output.solvertime;
    diagnostic.solvertime(:,i) = output.solvertime;
    try
        diagnostic.info = output.infostr;
    catch
        diagnostic.info = yalmiperror(output.problem,model.solver.tag);
    end
    diagnostic.problem(:,i) = output.problem;
    diagnostic.dimacs(:,i) = dimacs;

    % Some more info is saved internally
    solution_internal = diagnostic;
    solution_internal.variables = model.recoverdata.used_variables(:);
    solution_internal.optvar = output.Primal(:,i);

    if ~isempty(model.parametric_variables)
        diagnostic.mpsol = output.solveroutput;
        options.savesolveroutput = actually_save_output;
    end;

    if model.options.savesolveroutput
        diagnostic.solveroutput = output.solveroutput;
    end
    if model.options.savesolverinput
        diagnostic.solverinput = output.solverinput;
    end
    if model.options.saveyalmipmodel
        diagnostic.yalmipmodel = model;
    end

    if ops.warning && warningon && isempty(strfind(output.infostr,'No problems detected'))
        disp(['Warning: ' output.infostr]);
    end

    if ismember(output.problem,ops.beeponproblem)
        try
            beep; % does not exist on all ML versions
        catch
        end
    end

    % And we are done! Save the result
    if ~isempty(output.Primal)
        yalmip('setsolution',solution_internal,i);
    end
    if model.options.saveduals & model.solver.dual
        if isempty(model.Fremoved) | (nnz(model.Q)>0)
            try
                setduals(F,output.Dual,model.K);
            catch
            end
        else
            try
                % Duals related to equality constraints/free variables
                % have to be recovered b-A*x-Ht == 0
                b = -model.oldc;
                A = -model.oldF_struc(1+model.oldK.f:end,2:end)';
                H = -model.oldF_struc(1:model.oldK.f,2:end)';
                x = output.Dual;
                b_equ = b-A*x;
                newdual =  H\b_equ;
                setduals(model.Fremoved + F,[newdual;output.Dual],model.oldK);
            catch
                % this is a new feature...
                disp('Dual recovery failed. Please report this issue.');
            end
        end
    end
end

function newbase = mapObjective(local_vars,global_vars,base)
newbase = spalloc(length(global_vars),1,nnz(base));
for i = 1:length(local_vars)
    j = find(local_vars(i)==global_vars);
    newbase(j) = base(i);
end

function yesno = warningon

s = warning;
yesno = isequal(s,'on');

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%% STUDY_01_wrist_scaffold_single_fiber_iFEA_elastic
% Below is a demonstration for:
%
% * Inverse FEA based material parameter optimisation
% * febio_spec version 3.0
% * febio, FEBio
% * hexahedral elements, hex8
% * static, solid
% * isotropic elastic, Ogden
% * displacement logfile
% * stress logfile

clear; close all; clc;

%%
% Plot settings
fontSize=40;
markerSize=25;
lineWidth=3;

%% Control parameters
% Path names
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'2021_Nataliya_Perevoshchikova_wrist_project','data','temp');
loadPathMatProp=fullfile(defaultFolder,'2021_Nataliya_Perevoshchikova_wrist_project','data','mat_prop','PCL');

%%
% Building a quadrilateral circular mesh
pointSpacingHeight=0.1;
d=0.35;
h=10;
xSpacing=0.5;
ySpacing=0.4;
nCopies_x=8;
nCopies_y=8;

r=d/2;
ne=2; %Elements in radius
f=0.6; %Fraction (with respect to outer radius) where central square appears
NumStrands=4*7;
%Create the mesh
[Fq,Vq]=discQuadMesh(ne,r,f);

%%
% Lofting to a fiber
Vq(:,3)=0;
Vq1=Vq;
Vq1(:,3)=Vq1(:,3)-h/2;
Vq2=Vq;
Vq2(:,3)=Vq2(:,3)+h/2;

numStepsSweep=ceil(h/pointSpacingHeight);
xc=zeros(numStepsSweep,1);
yc=zeros(numStepsSweep,1);
zc=linspace(-h/2,h/2,numStepsSweep)';
Vc=[xc yc zc];

[~,~,~,S]=sweepLoft(Vq1,Vq2,[0 0 1],[0 0 1],Vc,numStepsSweep,0,0);

%Compose hexahedral elements
X=S.X'; Y=S.Y'; Z=S.Z'; %Coordinate matrices
V=[X(:) Y(:) Z(:)]; %Create node list

I=size(Vq,1)*((1:1:numStepsSweep)-1);
I=I(ones(size(Fq,1),1),:);
I=I(:);

FQ=repmat(Fq,numStepsSweep,1)+I(:,ones(size(Fq,2),1));
E=[FQ(1:end-size(Fq,1),:) FQ(size(Fq,1)+1:end,:)]; %The hexahedral elements
F_start=FQ(1:size(Fq,1),:);
F_end=FQ(end-size(Fq,1)+1:end,:);

[E,V,ind1,ind2]=mergeVertices(E,V); %Merge nodes (start and end are not shared yet)

[F]=element2patch(E); %Element faces
F_start=ind2(F_start); %Start faces for BC's
F_end=ind2(F_end); %End faces for BC's
ind=tesBoundary(F,V); %Indices of boundary faces
Fb=F(ind,:); %Boundary faces (for plotting)

%%
% Visualizing mesh

cFigure; hold on;
title('Hexahedral mesh one fiber','fontSize',fontSize);
gpatch(Fb,V,'gw','k',0.5);
gpatch(F_start,V,'rw','k',1);
gpatch(F_end,V,'bw','k',1);
axisGeom;
camlight headlight;
set(gca,'FontSize',fontSize);
drawnow;

% Defining file names
febioFebFileNamePart='tempModel';
febioFebFileName=fullfile(savePath,[febioFebFileNamePart,'.feb']); %FEB file name
febioLogFileName=[febioFebFileNamePart,'.txt']; %FEBio log file name
febioLogFileName_disp=[febioFebFileNamePart,'_disp_out.txt']; %Log file name for exporting displacement
febioLogFileName_force=[febioFebFileNamePart,'_force_out.txt']; %Log file name for exporting force
febioLogFileName_stress=[febioFebFileNamePart,'_stress_out.txt']; %Log file name for exporting stress

% FEA control settings
timeEnd=10;
numTimeSteps=10; %Number of time steps desired to reach yield strength
max_refs=50; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=6; %Optimum number of iterations
max_retries=60; %Maximum number of retires
dtmin=(timeEnd/numTimeSteps)/100; %Minimum time step size
dtmax=timeEnd/numTimeSteps; %Maximum time step size
tStep=timeEnd/numTimeSteps;%s
min_residual=1e-20;
symmetric_stiffness=1;
lstol=0.9;

%Initial material parameter set
k_factor=100;

%Material model is Ogden 2nd order with initial estimates
c1_ini=57;
m1_ini=2;
c2_ini=57;
m2_ini=-2;
k_ini=(c1_ini+c2_ini)*k_factor;
P=[c1_ini m1_ini c2_ini m2_ini];

% Experimental elastic properties
matname = fullfile(loadPathMatProp,'Elastic_properties.mat');
load(matname);

cFigure;
title('Experimental elastic region','fontSize',fontSize);
hold on;
xlabel('Average displacement, mm','FontSize',fontSize);
ylabel('Average load, N','FontSize',fontSize);
plot(disp_exp,load_exp,'k.-','LineWidth',3);
set(gca,'FontSize',fontSize);
camlight headlight;
displacementMagnitude=max(disp_exp);

%% DEFINE BC's

%Supported nodes
bcRigidList=unique(F_start(:));
%Prescribed force nodes
bcPrescribeList=unique(F_end(:));
bcPrescribeMagnitudes=displacementMagnitude(ones(1,numel(bcPrescribeList)),:);

%% Find side face sets
[Nb]=patchNormal(F,V);
zVec=[0 0 1];
d=dot(Nb,zVec(ones(size(Nb,1),1),:),2);
Z=V(:,3);
ZF=mean(Z(F),2);
pointSpacing=mean(patchEdgeLengths(F,V));

logicSides_F1=(d>-1) & ZF>=(max(V(:,3))-2*pointSpacing);
logicSides_F2=(d<1) & ZF<=(min(V(:,3))+2*pointSpacing);

logicSides_FS=logicSides_F1;
logicSides_FL=logicSides_F2;

faceBoundaryMarker(logicSides_FS)=1;
faceBoundaryMarker(logicSides_FL)=2;

%Visualize BC's
cFigure; hold on;
patch('faces',F,'vertices',V,'FaceColor',[0, 0.5, 0],'FaceAlpha',0.1,'EdgeColor','None');

hl(1)=plotV(V(bcRigidList,:),'k.','MarkerSize',markerSize);
hl(2)=plotV(V(bcPrescribeList,:),'.','Color',[0.6350, 0.0780, 0.1840],'MarkerSize',markerSize);

legend(hl,{'BC full support','BC prescribed Z displacement'})
axisGeom;
camlight headlight;
set(gca,'FontSize',fontSize);

axes('position',[.2 .6 .2 .2], 'NextPlot', 'add');
plotV(V(bcPrescribeList,:),'.','Color',[0.6350, 0.0780, 0.1840],'MarkerSize',markerSize*2);
patch('faces',F(faceBoundaryMarker==1,:),'vertices',V,'FaceColor',[0, 0.5, 0],'FaceAlpha',0.5,'EdgeColor','k');
view(3);
axis equal; axis vis3d; axis tight;
set(gca,'XColor','none','YColor','none','ZColor','none');

axes('position',[.2 .3 .2 .2], 'NextPlot', 'add');
plotV(V(bcRigidList,:),'k.','MarkerSize',markerSize*2);
patch('faces',F(faceBoundaryMarker==2,:),'vertices',V,'FaceColor',[0, 0.5, 0],'FaceAlpha',0.5,'EdgeColor','k');
view(3);
axis equal; axis vis3d; axis tight;
set(gca,'XColor','none','YColor','none','ZColor','none');

drawnow;

%% Defining the FEBio input structure
% See also |febioStructTemplate| and |febioStruct2xml| and the FEBio user
% manual.

%Get a template with default settings
[febio_spec]=febioStructTemplate;

%febio_spec version
febio_spec.ATTR.version='3.0';

%Module section
febio_spec.Module.ATTR.type='solid';

%Control section
febio_spec.Control.analysis='STATIC';
febio_spec.Control.time_steps=numTimeSteps;
febio_spec.Control.step_size=tStep;
febio_spec.Control.time_stepper.dtmin=dtmin;
febio_spec.Control.time_stepper.dtmax=dtmax;
febio_spec.Control.time_stepper.max_retries=max_retries;
febio_spec.Control.time_stepper.opt_iter=opt_iter;
febio_spec.Control.solver.max_refs=max_refs;
febio_spec.Control.solver.max_ups=max_ups;

%Material section
materialName1='Material1';
febio_spec.Material.material{1}.ATTR.name=materialName1;
febio_spec.Material.material{1}.ATTR.type='Ogden';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.c1=c1_ini;
febio_spec.Material.material{1}.m1=m1_ini;
febio_spec.Material.material{1}.c2=c2_ini;
febio_spec.Material.material{1}.m2=m2_ini;
febio_spec.Material.material{1}.k=k_ini;

%Geometry section
% -> Nodes
febio_spec.Mesh.Nodes{1}.ATTR.name='Object1'; %The node set name
febio_spec.Mesh.Nodes{1}.node.ATTR.id=(1:size(V,1))'; %The node id's
febio_spec.Mesh.Nodes{1}.node.VAL=V; %The nodel coordinates

% -> Elements
partName1='Part1';
febio_spec.Mesh.Elements{1}.ATTR.name=partName1; %Name of this part
febio_spec.Mesh.Elements{1}.ATTR.type='hex8'; %Element type of this set
febio_spec.Mesh.Elements{1}.elem.ATTR.id=(1:1:size(E,1))'; %Element id's
febio_spec.Mesh.Elements{1}.elem.VAL=E;

% -> NodeSets
nodeSetName1='bcSupportList';
nodeSetName2='bcPrescribeList_Z';

febio_spec.Mesh.NodeSet{1}.ATTR.name=nodeSetName1;
febio_spec.Mesh.NodeSet{1}.node.ATTR.id=bcRigidList(:);

febio_spec.Mesh.NodeSet{2}.ATTR.name=nodeSetName2;
febio_spec.Mesh.NodeSet{2}.node.ATTR.id=bcPrescribeList(:);

%MeshDomains section
febio_spec.MeshDomains.SolidDomain.ATTR.name=partName1;
febio_spec.MeshDomains.SolidDomain.ATTR.mat=materialName1;

%Boundary condition section
% -> Fix boundary conditions
febio_spec.Boundary.bc{1}.ATTR.type='fix';
febio_spec.Boundary.bc{1}.ATTR.node_set=nodeSetName1;
febio_spec.Boundary.bc{1}.dofs='x';

febio_spec.Boundary.bc{2}.ATTR.type='fix';
febio_spec.Boundary.bc{2}.ATTR.node_set=nodeSetName1;
febio_spec.Boundary.bc{2}.dofs='y';

febio_spec.Boundary.bc{3}.ATTR.type='fix';
febio_spec.Boundary.bc{3}.ATTR.node_set=nodeSetName1;
febio_spec.Boundary.bc{3}.dofs='z';

febio_spec.Boundary.bc{4}.ATTR.type='fix';
febio_spec.Boundary.bc{4}.ATTR.node_set=nodeSetName2;
febio_spec.Boundary.bc{4}.dofs='x';

febio_spec.Boundary.bc{5}.ATTR.type='fix';
febio_spec.Boundary.bc{5}.ATTR.node_set=nodeSetName2;
febio_spec.Boundary.bc{5}.dofs='y';

% -> Prescribe boundary conditions
febio_spec.Boundary.bc{6}.ATTR.type='prescribe';
febio_spec.Boundary.bc{6}.ATTR.node_set=nodeSetName2;
febio_spec.Boundary.bc{6}.dof='z';
febio_spec.Boundary.bc{6}.scale.ATTR.lc=1;
febio_spec.Boundary.bc{6}.scale.VAL=1;
febio_spec.Boundary.bc{6}.relative=0;

%LoadData section
% -> load_controller
febio_spec.LoadData.load_controller{1}.ATTR.id=1;
febio_spec.LoadData.load_controller{1}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{1}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{1}.points.point.VAL=[0 0; timeEnd displacementMagnitude];

%Output section
% -> log file
febio_spec.Output.logfile.ATTR.file=febioLogFileName;
febio_spec.Output.logfile.node_data{1}.ATTR.file=febioLogFileName_disp;
febio_spec.Output.logfile.node_data{1}.ATTR.data='ux;uy;uz';
febio_spec.Output.logfile.node_data{1}.ATTR.delim=',';
febio_spec.Output.logfile.node_data{1}.VAL=1:size(V,1);

febio_spec.Output.logfile.node_data{2}.ATTR.file=febioLogFileName_force;
febio_spec.Output.logfile.node_data{2}.ATTR.data='Rx;Ry;Rz';
febio_spec.Output.logfile.node_data{2}.ATTR.delim=',';
febio_spec.Output.logfile.node_data{2}.VAL=1:size(V,1);

febio_spec.Output.logfile.element_data{1}.ATTR.file=febioLogFileName_stress;
febio_spec.Output.logfile.element_data{1}.ATTR.data='sz';
febio_spec.Output.logfile.element_data{1}.ATTR.delim=',';
febio_spec.Output.logfile.element_data{1}.VAL=1:size(E,1);

%% Quick viewing of the FEBio input file structure
% The |febView| function can be used to view the xml structure in a MATLAB
% figure window.

%%
% febView(febio_spec); %Viewing the febio file|

%% Exporting the FEBio input file
% Exporting the febio_spec structure to an FEBio input file is done using
% the |febioStruct2xml| function.

febioStruct2xml(febio_spec,febioFebFileName); %Exporting to file and domNode

%% Running the FEBio analysis
% To run the analysis defined by the created FEBio input file the
% |runMonitorFEBio| function is used. The input for this function is a
% structure defining job settings e.g. the FEBio input file name. The
% optional output runFlag informs the user if the analysis was run
% succesfully.

febioAnalysis.run_filename=febioFebFileName; %The input file name
febioAnalysis.run_logname=febioLogFileName; %The name for the log file
febioAnalysis.disp_on=1; %Display information on the command window
febioAnalysis.runMode='internal';%external';

[runFlag]=runMonitorFEBio(febioAnalysis);%START FEBio NOW!!!!!!!!

%% Import FEBio results

if runFlag==1 %i.e. a succesful run

    % Importing nodal displacements from a log file
    [~, N_disp_mat,~]=importFEBio_logfile(fullfile(savePath,febioLogFileName_disp));

    N_disp_mat=N_disp_mat(:,2:end,:);
    sizImport=size(N_disp_mat);
    sizImport(3)=sizImport(3)+1;
    N_disp_mat_n=zeros(sizImport);
    N_disp_mat_n(:,:,2:end)=N_disp_mat;
    N_disp_mat=N_disp_mat_n;
    DN=N_disp_mat(:,:,end);
    DN_magnitude=sqrt(sum(DN(:,1).^2,2)); %X
    V_def=V+DN;
    [CF]=vertexToFaceMeasure(F,DN_magnitude);

    %Importaing force from log file
    [time_mat, Force_mat,~]=importFEBio_logfile(fullfile(savePath,febioLogFileName_force));
    time_mat=[0; time_mat(:)]; %Time
    Force_mat=Force_mat(:,2:end,:);
    sizImport=size(Force_mat);
    sizImport(3)=sizImport(3)+1;
    Force_mat_n=zeros(sizImport);
    Force_mat_n(:,:,2:end)=Force_mat;
    Force_mat=Force_mat_n;

    %%
    % Plotting the simulated results using |anim8| to visualize and animate
    % deformations

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    hp=gpatch(F,V_def,CF,'k',1); %Add graphics object to animate
    gpatch(F,V,0.5*ones(1,3),'none',0.25); %A static graphics object

    axisGeom(gca,fontSize);
    colormap(gjet(250)); colorbar;
    caxis([0 max(DN_magnitude)]);
    axis([min(V_def(:,1)) max(V_def(:,1)) min(V_def(:,2)) max(V_def(:,2)) min(V_def(:,3)) max(V_def(:,3))]); %Set axis limits statically
    view(130,25); %Set view direction
    camlight headlight;

    % Set up animation features
    animStruct.Time=time_mat; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments
        DN=N_disp_mat(:,:,qt); %Current displacement
        DN_magnitude=sqrt(sum(DN.^2,2)); %Current displacement magnitude
        V_def=V+DN; %Current nodal coordinates
        [CF]=vertexToFaceMeasure(F,DN_magnitude); %Current color data to use

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp hp]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData'}; %Properties of objects to animate
        animStruct.Set{qt}={V_def,CF}; %Property values for to set in order to animate
    end
    anim8(hf,animStruct); %Initiate animation feature
    drawnow;

    %%
    % Calculate the simulated applied uniaxial stretch
    DZ_set=N_disp_mat(bcPrescribeList,3,:); %Z displacements of the prescribed set
    DZ_set=mean(DZ_set,1); %Calculate mean Z displacements across nodes
    DZ_set=squeeze(DZ_set);

    %Calculated the simulated force
    FZ_set=Force_mat(bcPrescribeList,3,:);
    FZ_set=sum(FZ_set,1); %Calculate mean Z force across nodes
    FZ_set=squeeze(FZ_set);

    disp_sim=DZ_set;
    load_sim=FZ_set*NumStrands;

    [~, indUnique] = unique(disp_sim);
    disp_sim=disp_sim(indUnique);
    load_sim=load_sim(indUnique);
    load_sim_exp = interp1(disp_sim,load_sim,disp_exp,'pchip');

    % Visualize curves
    hf=cFigure;

    hold on;
    title('Uniaxial load-displacement curve','FontSize',fontSize);
    xlabel('Displacement','FontSize',fontSize); ylabel('Load','FontSize',fontSize);

    hl(1)=plot(disp_sim,load_sim,'r.-','lineWidth',lineWidth*1.1,'markerSize',markerSize*2);
    hl(2)=plot(disp_exp,load_exp,'k-','lineWidth',lineWidth);

    legend(hl,{'Simulation','Experiment'},'Location','northwest');

    view(2); axis tight;  grid on; axis square; box on;
    set(gca,'FontSize',fontSize);
    drawnow;

end

%% Create structures for optimization

% Material structure
mat_struct.id=1; %Material id
mat_struct.par_names={'c1','m1','c2','m2','k'}; %Parameter names
mat_struct.par_values={c1_ini m1_ini c2_ini m2_ini k_ini}; %Parameter values

objectiveStruct.parNormFactors=abs(P); %This will normalize the parameters to ones(size(P))
objectiveStruct.Pb_struct.xx_c=P; %Parameter constraining centre
objectiveStruct.Pb_struct.xxlim=[[P(1)/10000 -10000 P(3)/10000 -10000]' [P(1)*10000 10000 P(3)*10000 10000]']; %Parameter bounds

febioAnalysis.disp_on=0;
febioAnalysis.disp_log_on=0;

%What should be known to the objective function:
objectiveStruct.h=hl(1);
objectiveStruct.bcPrescribeList=bcPrescribeList;
objectiveStruct.disp_exp=disp_exp;
objectiveStruct.load_exp=load_exp;
objectiveStruct.febioAnalysis=febioAnalysis;
objectiveStruct.febio_spec=febio_spec;
objectiveStruct.febioFebFileName=febioFebFileName;
objectiveStruct.mat_struct=mat_struct;
objectiveStruct.k_factor=k_factor;
objectiveStruct.NumStrands=NumStrands;

%Optimisation settings
maxNumberIterations=100; %Maximum number of optimization iterations
maxNumberFunctionEvaluations=maxNumberIterations*10; %Maximum number of function evaluations, N.B. multiple evaluations are used per iteration
functionTolerance=1e-25; %Tolerance on objective function value
parameterTolerance=0.001;%Tolerance on parameter variation
displayTypeIterations='iter';

objectiveStruct.method=2;

%File names of output files
output_names.displacement=fullfile(savePath,febioLogFileName_disp);
output_names.load=fullfile(savePath,febioLogFileName_force);

objectiveStruct.run_output_names=output_names;

%% start optimization

Pn=P./objectiveStruct.parNormFactors;

switch objectiveStruct.method
    case 1 %fminsearch and Nelder-Mead
        OPT_options=optimset('fminsearch'); % 'Nelder-Mead simplex direct search'
        OPT_options = optimset(OPT_options,'MaxFunEvals',maxNumberFunctionEvaluations,...
            'MaxIter',maxNumberIterations,...
            'TolFun',functionTolerance,...
            'TolX',parameterTolerance,...
            'Display',displayTypeIterations,...
            'FinDiffRelStep',1e-2,...
            'DiffMaxChange',0.5);
        [Pn_opt,OPT_out.fval,OPT_out.exitflag,OPT_out.output]= fminsearch(@(Pn) objectiveFunctionIFEA(Pn,objectiveStruct),Pn,OPT_options);
    case 2 %lsqnonlin and Levenberg-Marquardt
        OPT_options = optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt');
        OPT_options = optimoptions(OPT_options,'MaxFunEvals',maxNumberFunctionEvaluations,...
            'MaxIter',maxNumberIterations,...
            'TolFun',functionTolerance,...
            'TolX',parameterTolerance,...
            'Display',displayTypeIterations,...
            'FinDiffRelStep',1e-2,...
            'DiffMaxChange',0.5);
        [Pn_opt,OPT_out.resnorm,OPT_out.residual]= lsqnonlin(@(Pn) objectiveFunctionIFEA(Pn,objectiveStruct),Pn,[],[],OPT_options);
end

%%
[Fopt,OPT_stats_out]=objectiveFunctionIFEA(Pn_opt,objectiveStruct);

%% Unnormalize and constrain parameters

P_opt=Pn_opt.*objectiveStruct.parNormFactors; %Scale back, undo normalization

%Constraining parameters
for q=1:1:numel(P)
    [P(q)]=boxconstrain(P(q),objectiveStruct.Pb_struct.xxlim(q,1),objectiveStruct.Pb_struct.xxlim(q,2),objectiveStruct.Pb_struct.xx_c(q));
end

disp_text=sprintf('%6.16e,',P_opt); disp_text=disp_text(1:end-1);
disp(['P_opt=',disp_text]);

%%
hf1=cFigure;
%title('displacement load curves optimised','FontSize',fontSize);
xlabel('Displacement (mm)','FontSize',fontSize); ylabel('Force (N)','FontSize',fontSize); zlabel('Z','FontSize',fontSize); hold on;

Hn(1)=plot(OPT_stats_out.disp_sim,OPT_stats_out.load_sim,'r.-','lineWidth',lineWidth*1.1,'markerSize',markerSize*2);
Hn(2)=plot(disp_exp,load_exp,'k-','lineWidth',lineWidth);

legend(Hn,{'Simulation','Experiment'},'Location','northwest');
view(2); axis tight;  grid on;
set(gca,'FontSize',fontSize);
drawnow;

function [Fopt,OPT_stats_out]=objectiveFunctionIFEA(Pn,objectiveStruct)

%%

febioFebFileName=objectiveStruct.febioFebFileName;
febio_spec=objectiveStruct.febio_spec;

%% Unnormalize and constrain parameters

P=Pn.*objectiveStruct.parNormFactors; %Scale back, undo normalization
P_in=P; %Proposed P

%Constraining parameters
for q=1:1:numel(P)
    [P(q)]=boxconstrain(P(q),objectiveStruct.Pb_struct.xxlim(q,1),objectiveStruct.Pb_struct.xxlim(q,2),objectiveStruct.Pb_struct.xx_c(q));
end

%% Setting material parameters

%Acces material parameters
mat_struct=objectiveStruct.mat_struct;

mat_struct.par_values={P(1) P(2) P(3) P(4) (P(1)+P(3))*objectiveStruct.k_factor};

disp('SETTING MATERIAL PARAMETERS...');
disp(['Proposed (norm.): ',sprintf(repmat('%6.6e ',[1,numel(Pn)]),Pn)]);
disp(['Proposed        : ',sprintf(repmat('%6.6e ',[1,numel(P_in)]),P_in)]);
disp(['Set (constr.)   : ',sprintf(repmat('%6.6e ',[1,numel(P)]),P)]);

%Assign material parameters
matId=mat_struct.id;
for q=1:1:numel(mat_struct.par_names)
    parNameNow=mat_struct.par_names{q};
    parValuesNow=mat_struct.par_values{q};
    febio_spec.Material.material{matId}.(parNameNow)=mat2strIntDouble(parValuesNow);
end
febioStruct2xml(febio_spec,febioFebFileName); %Exporting to file and domNode

disp('Done')

%% START FEBio

[runFlag]=runMonitorFEBio(objectiveStruct.febioAnalysis);

pause(0.1);

bcPrescribeList=objectiveStruct.bcPrescribeList;
disp_exp=objectiveStruct.disp_exp;
load_exp=objectiveStruct.load_exp;

if runFlag==1

    %Importing displacement
    [~,N_disp_mat,~]=importFEBio_logfile(objectiveStruct.run_output_names.displacement);
    N_disp_mat=N_disp_mat(:,2:end,:); %Remove index
    sizImport=size(N_disp_mat);
    sizImport(3)=sizImport(3)+1;
    N_disp_mat_n=zeros(sizImport);
    N_disp_mat_n(:,:,2:end)=N_disp_mat;
    N_disp_mat=N_disp_mat_n; %With added zeros for t=0

    %Derive applied stretch
    DZ_set=N_disp_mat(bcPrescribeList,3,:); %Final nodal displacements
    DZ_set=mean(DZ_set,1);
    DZ_set=squeeze(DZ_set);

    %Importing force from log file
    [time_mat, Force_mat,~]=importFEBio_logfile(objectiveStruct.run_output_names.load);
    time_mat=[0; time_mat(:)]; %Time
    Force_mat=Force_mat(:,2:end,:); %Remove index
    sizImport=size(Force_mat);
    sizImport(3)=sizImport(3)+1;
    Force_mat_n=zeros(sizImport);
    Force_mat_n(:,:,2:end)=Force_mat;
    Force_mat=Force_mat_n; %With added zeros for t=0

    %Calculated the simulated force
    FZ_set=Force_mat(bcPrescribeList,3,:);
    FZ_set=sum(FZ_set,1); %Calculate mean Z force across nodes
    FZ_set=squeeze(FZ_set);

    disp_sim=DZ_set;
    load_sim=FZ_set*objectiveStruct.NumStrands; %Multiply force by number of strands

    %Remove double entries
    [~, indUnique] = unique(disp_sim);
    disp_sim=disp_sim(indUnique);
    load_sim=load_sim(indUnique);

    if ~isempty(objectiveStruct.h)
        objectiveStruct.h.XData=disp_sim;
        objectiveStruct.h.YData=load_sim;
        drawnow;
    end

    %Interpolate experiment onto simulated points
    load_sim_exp = interp1(disp_sim,load_sim,disp_exp,'pchip');

    %Derive Fopt
    loadDev=load_exp-load_sim_exp;

    switch objectiveStruct.method
        case 1
            Fopt=sum((loadDev).^2); %Sum of squared differences
        case 2
            Fopt=sum((loadDev).^2)/length(loadDev);%loadDev(:); %Squared differences
    end

    OPT_stats_out.load_sim=load_sim;
    OPT_stats_out.disp_sim=disp_sim;
    OPT_stats_out.loadDev=loadDev;
    OPT_stats_out.Fopt=Fopt;

else %Output NaN
    switch objectiveStruct.method
        case 1
            Fopt=NaN;
        case 2
            Fopt=NaN(size(load_exp)); %Squared differences
    end
    OPT_stats_out=[];
end

end

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function run_writesnap(varargin)
% Send real-time BCI outputs to a SNAP-based stimulus presentation system.
% run_writetcp(Model,SourceStream,OutputHost,OutputPort,OutputForm,MessageFormat,UpdateFrequency,PredictorName,ConnTimeout)
%
% This function runs in the background and processes data from some MATLAB stream (created with some
% other background input plugin, e.g., the BioSemi reader). The processed data is periodically
% forwarded to a program running the SNAP framework.
%
% In:
%   Model : A model data structure (as obtained from bci_train) based on which the predictor shall be
%           created; typically this is a model struct, but for convenience it can be a file name,
%           variable name in the base workspace, or a cell array of {file name, variable name} to
%           refer to a variable inside a .mat file. The model is not modified by this function.
%           (default: 'lastmodel')
%
%   SourceStreamNames : Optional names of stream data structures in the MATLAB base workspace to
%                       consider as possible data sources (previously created with onl_newstream);
%                       if a stream contains all channels that are needed by the predictor, or
%                       alternatively has the right number and type of channels it will be considered
%                       as a potential source stream unless ambiguous. (default: 'laststream')
%
%   OutputHost : destination host name to send results to (computer name, URL or IP address)
%                (default: 'localhost')
%
%   OutputPort : destination port to which the results are sent (default: 7897)
%
%   OutputForm : output data form, see onl_predict (default: 'distribution')
%
%   TargetVariable : name of the variable in the SNAP module which should carry the BCI predictions
%                    (default: 'bci')
%
%   UpdateFrequency : The rate at which new outputs will be computed. (default: 10)
%
%   PredictAt : Predict at markers. If nonempty, this is a cell array of online target markers relative
%               to which predictions shall be made. If empty, predictions are always made on the most recently
%               added sample. (default: {})
%
%   PredictorName : Name of the predictor to be created; a variable of this name will be created in
%                   the MATLAB base workspace to hold the predictor's state. If a variable with this
%                   name already exists it will be overridden. (default: 'lastpredictor')
%
%   Verbose : whether to display verbose outputs (e.g. connection failure) (default: false)
%
% Examples:
%   % open a new BCILAB processing stream, using the previously learned predictive model 'mymodel',
%   % and reading from a previously opened input stream named 'mystream'; send outputs to some IP:port
%   run_writesnap('mymodel','mystream','192.168.1.5')
%
%   % as before, but pass arguments by name
%   run_writetcp('Model','mymodel','SourceStream','mystream','OutputHost','192.168.1.5')
%
%   % as before, but transmit the data in a raw binary encoding
%   run_writetcp('Model','mymodel','SourceStream','mystream','OutputHost','192.168.1.5','TargetVariable','alpha')
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-11-19

% declare the name of this component (shown in the menu)
declare_properties('name','SNAP');

% define arguments
opts = arg_define(varargin, ...
    arg({'pred_model','Model'}, 'lastmodel', [], 'Predictive model. As obtained via bci_train or the Model Calibration dialog.','type','expression'), ...
    arg({'in_stream','SourceStreamNames','SourceStream'}, 'laststream',[],'Input Matlab stream name(s). Optional names of stream data structures in the MATLAB base workspace to consider as possible data sources (previously created with onl_newstream); if a stream contains all channels that are needed by the predictor, or alternatively has the right number and type of channels it will be considered as a potential source stream unless ambiguous.','typecheck',false,'shapecheck',false), ...
    arg({'out_hostname','OutputHost','Host'}, 'localhost',[],'Destination TCP hostname. Can be a computer name, URL, or IP address.'), ...
    arg({'out_port','OutputPort','Port'}, 7897, [],'Destination TCP port. Matches the SNAP default.'), ...
    arg({'out_form','OutputForm','Form'},'expectation',{'expectation','distribution','mode','raw'},'Output form. Can be the expected value (posterior mean) of the target variable, or the distribution over possible target values (probabilities for each outcome, or parametric distribution), or the mode (most likely value) of the target variable.'), ...
    arg({'target_variable','TargetVariable','Target'},'bci',[],'Target Variable. This is the variable name in the SNAP module that should receive the BCI value'), ...
    arg({'update_freq','UpdateFrequency'},10,[0 Inf],'Update frequency. The rate at which new outputs will be computed.'), ...
    arg({'predict_at','PredictAt'}, {},[],'Predict at markers. If nonempty, this is a cell array of online target markers relative to which predictions shall be made. If empty, predictions are always made on the most recently added sample.','type','expression'), ...
    arg({'pred_name','PredictorName'}, 'lastpredictor',[],'Name of new predictor. A variable of this name will be created in the MATLAB base workspace to hold the predictor''s state. If a variable with this name already exists it will be overridden.'), ...
    arg({'verbose_output','Verbose'}, false,[],'Verbose output. Whether to display verbose outputs (e.g., connection failure).'));

% convert format strings to formatting functions
if ~isvarname(opts.target_variable)
    disp('Note: The given target variable is likely not a valid variable name.'); end

% get a fresh connection id
id = sprintf('c%.0f',fresh_id('snapsocket'));

% start background writer job
onl_write_background( ...
    'ResultWriter',@(y)send_message(y,opts.verbose_output,id,opts.out_hostname,opts.out_port,opts.target_variable),...
    'MatlabStream',opts.in_stream, ...
    'Model',opts.pred_model, ...
    'OutputFormat',opts.out_form, ...
    'UpdateFrequency',opts.update_freq, ...
    'PredictorName',opts.pred_name, ...
    'PredictAt',opts.predict_at, ...
    'Verbose',opts.verbose_output, ...
    'StartDelay',0,...
    'EmptyResultValue',[]);

disp('Now writing...');

% background message sending function
function send_message(yy,verbose_output,id,host,port,varname)
persistent conns;
try
    % try to connect if necessary
    if ~isfield(conns,id)
        conns.(id) = connect(host,port); end
    % for each prediction...
    for k=1:size(yy,1)
        y = yy(k,:);
        try
            % send it off
            strm = conns.(id).strm;
            if isscalar(y)
                strm.writeBytes(char([sprintf('setup %s=%.5f',varname,y) 10]));
            elseif ~isempty(y)
                strm.writeBytes(char(['setup ' varname '=(' sprintf('%.5f,',y) ')' 10]));
            end
            strm.flush();
        catch e1
            if strcmp(e1.identifier, 'MATLAB:Java:GenericException')
                % failed to send: try to re-connect...
                conns.(id) = connect(host,port);
            else
                rethrow(e1);
            end
        end
    end
catch e2
    if strcmp(e2.identifier, 'MATLAB:Java:GenericException')
        if verbose_output
            fprintf('Could not connect to %s:%.0f\n',host,port); end
    else
        hlp_handleerror(e2);
    end
end

function newconn = connect(host,port)
import java.io.*
import java.net.*
import java.lang.*
newconn.sock = Socket();
newconn.sock.setTcpNoDelay(1);
newconn.sock.setSoTimeout(5);
newconn.sock.connect(InetSocketAddress(host,port),1);
newconn.strm = DataOutputStream(newconn.sock.getOutputStream());
disp(['Connection to ' host ' established.']);

repo: REST


text:
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
clc; zoom on; clf; echo on;

%TFDEMO1 Introduction to the Time-Frequency Toolbox.
%	O. Lemoine - May 1996.
%	Copyright (c) CNRS.

% Welcome to the Time-Frequency Toolbox for MATLAB. This demonstration
% follows the plan of the tutorial, and consider exactly the same examples.
% Therefore, for further information about these illustrations, we advise
% you to refer to the corresponding chapter of the tutorial.
%
% Press any key to continue...

pause; clc;

% First let's create an analytic linear frequency modulated signal, whose
% normalized frequency is changing from 0 to 0.5 :

sig1=fmlin(128,0,0.5);
plot(real(sig1)); axis([1 128 -1 1]);
xlabel('Time'); ylabel('Real part');
title('Linear frequency modulation'); grid

% From this time-domain representation, it is difficult to say what kind
% of modulation is contained in this signal.
%
% Press any key to continue...

pause; clc;

% Now let's consider its energy spectrum :

dsp1=fftshift(abs(fft(sig1)).^2);
plot((-64:63)/128,dsp1);
xlabel('Normalized frequency'); ylabel('Squared modulus');
title('Spectrum'); grid

% We still can not say, from this plot, anything about the evolution in
% time of the frequency content.
%
% In order to have a more informative description of such a signal, it
% would be better to directly represent their frequency content while
% still keeping the time description parameter : this is precisely the
% aim of time-frequency analysis. To illustrate this, let's try the
% Wigner-Ville distribution on this signal :
%
% press any key to continue...

pause;

tfrwv(sig1);

% We can see on this representation that the linear progression of the
% frequency with time, from 0 to 0.5, is clearly shown.
%
% Press any key to continue...

pause; clc;

% If we now add some complex white Gaussian noise on this signal, with
% a 0 dB signal to noise ratio,

sig2=sigmerge(sig1,noisecg(128),0);
Min=min(real(sig2)); Max=max(real(sig2));
clf; plot(real(sig2)); axis([1 128 Min Max]);
xlabel('Time'); ylabel('Real part');
title('Linear frequency modulation plus noise'); grid

% press any key to continue...

pause;

% and consider the spectrum of it :

dsp2=fftshift(abs(fft(sig2)).^2);
plot((-64:63)/128,dsp2);
xlabel('Normalized frequency'); ylabel('Squared modulus');
title('Spectrum'); grid

% It is worse than before to interpret these plots. On the other hand, the
% Wigner-Ville distribution still show quite clearly the linear progression
% of the frequency with time :
%
% press any key to continue...

pause;

tfrwv(sig2);

% press any key to continue...

pause; clc;

% The second example we consider is a bat sonar signal, recorded with a
% sampling frequency of 230.4 kHz and an effective bandwidth equal to
% [8 kHz, 80 kHz].
%  First, load the signal from the MAT-file bat.mat :

echo off
DirectoryStr='../data/';
while (exist([DirectoryStr 'bat.mat'])==0),
 fprintf('I can''t find %s\n', [DirectoryStr 'bat.mat']);
 DirectoryStr=input('name of the directory where bat.mat is : ','s');
end;
eval(['load ' DirectoryStr 'bat.mat']);
echo on

Fe=230.4; t0=(0:2047).'/Fe;
clf; plot(t0,bat); xlabel('Time [ms]');
axis([t0(1) t0(2048) -900 800]); grid;

% From this plot, we can not say precisely what is the frequency content
% at each time instant t ; similarly, if we look at its spectrum,
%
% press any key to continue...

pause;
dsp=fftshift(abs(fft(bat)).^2); dsp=dsp/max(dsp);
f0=(-1024:1023)*Fe/2048;
plot(f0,10*log10(dsp)); xlabel('Frequency [kHz]'); ylabel('Squared modulus  (dB)');
axis([0 Fe/2,-30 0])
title('Spectrum  (dB)'); grid; pause

% we can not say at what time the signal is located around 38 kHz, and at
% what time around 40 kHz. Let us now consider a representation called
% the pseudo Wigner-Ville distribution, applied on the most interesting
% part of this signal :
%
% press any key to continue...

tfrpwv(hilbert(bat),8:8:2048,256,tftb_window(127,'hanning'),1);

% We then have a nice description of its spectral content varying with
% time : it is narrow-band signal, whose frequency content is decreasing
% from around 55 kHz to 38kHz, with a non-linear frequency modulation
% (approximately of hyperbolic shape).
%
% press any key to continue...

pause; clc;

% The last introductory example presented here is a transient signal
% embedded in a -5 dB white Gaussian noise. This transient signal is a
% constant frequency modulated by a one-sided exponential amplitude :

trans=amexpo1s(92,1,64).*fmconst(92);
sig=[zeros(100,1) ; trans ; zeros(64,1)];
sign=sigmerge(sig,noisecg(256),-5);
Min=min(real(sign)); Max=max(real(sign));
subplot(211); plot(real(sign)); axis([1 256 Min Max]);
xlabel('Time'); title('Noisy transient signal'); grid
dsp=fftshift(abs(fft(sign)).^2);
subplot(212); plot((-128:127)/256,dsp); grid
xlabel('Normalized frequency'); title('Energy spectrum');

% From these representations, it is difficult to localize precisely the
% signal in the time-domain as well as in the frequency domain.  Now let us
% have a look at the spectrogram of this signal :
%
% press any key to continue...

pause;
tfrsp(sign);

% the transient signal appears distinctly around the normalized frequency
% 0.25, and between time points 125 and 160.
%
% Press any key to return to the main menu.

echo off

author: HeLiangHIT


text:
function indx_max = r8vec_indexed_heap_d_max ( n, a, indx )

%*****************************************************************************80
%
%% R8VEC_INDEXED_HEAP_D_MAX: maximum value in heap descending indexed R8VEC.
%
%  Discussion:
%
%    An R8VEC is a vector of R8's.
%
%    An indexed R8VEC is an R8VEC of data values, and an R8VEC of N indices,
%    each referencing an entry of the data vector.
%
%    This is one of three functions needed to model a priority queue.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 August 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Thomas Cormen, Charles Leiserson, Ronald Rivest,
%    Introduction to Algorithms,
%    MIT Press, 2001,
%    ISBN: 0262032937,
%    LC: QA76.C662.
%
%  Parameters:
%
%    Input, integer N, the number of items in the index vector.
%
%    Input, real A(*), the data vector.
%
%    Input, integer INDX(N), the index vector.
%
%    Output, integer INDX_MAX, the index in A of the maximum value
%    in the heap.
%
  indx_max = indx(1);

  return
end


text:
function dx = HIVsys_ZURAKOWSKI(t,x,u,p)

try
    run_HIV_params

    dx = zeros(5,1);
    dx(1) = lambda1 - d*x(1) - alpha1*(1-eta*u(1))*x(1)*x(2);
    dx(2) = alpha1*(1-eta*u(1))*x(2)*x(1) - a*x(2) - p1*x(4)*x(2) - p2*x(5)*x(2);
    dx(3) = c2*x(1)*x(2)*x(3) - c2*q*x(2)*x(3) - b2*x(3);
    dx(4) = c1*x(2)*x(4) - b1*x(4);
    dx(5) = c2*q*x(2)*x(3) - h*x(5);
catch
    keyboard
end

author: eurika-kaiser

repo: SINDY-MPC


text:
function obj = misc_history(obj)
%MISC_HISTORY - Records the history of BTB toolbox function calls that
%               affected the object.
%
%Synopsis:
%  MISC_HISTORY(OBJ)
%
%Arguments:
%  OBJ:     struct with neuro data, or montage, or marker
%
%Returns:
%  OBJ:     struct with extra history entry in the substruct .history with
%           the following fields:
%
%  fcn:           function handle
%  fcn_params:    names of named function parameters.
%  varargin:      optional arguments
%  date:          date the function was executed
%
%  All named function arguments also appear as fields except for the object
%  that is modified (e.g. the EEG struct). The latter appears in the
%  parameter list as 'obj'. (note: to avoid confusion proc_ functions
%  should not have another parameter called 'obj')
%
% Note: misc_history saves all arguments passed to functions unless they
% exceed a certain byte size.
%
%See also: misc_applyHistory

% 08-2012 Matthias Treder
global BTB

if ~isempty(BTB.History) && BTB.History==0 ||isempty(obj)
  return
end

misc_checkType(obj,'STRUCT');

if isfield(obj,'history')
  ht = obj.history;
  N = numel(ht)+1;
else
  ht = cell(1,1);
  N = 1;
end
ht{N} = struct();

objname = inputname(1);  % Name of object in caller's workspace

maxSize = 10*1000*1000;   % max byte size of arguments saved in history
% Get function name (of caller)
ST = dbstack('-completenames',1);
ht{N}.fcn = eval(['@' ST(1).name]);

% Find labels of function parameters (arguments)
code = (evalc(['type ' ST(1).name])); % get function code
expr = ['^\s*function s*\[*\s*(\w+|\s*|,)*\]*\s*=\s*' ST(1).name '\s*\((\w+|\s+|,)*\)'];
token = regexpi(code,expr,'tokens');
token = strrep(token{1}{end},' ','');  % remove whitespace
token = regexp(token,',','split');   % split arguments by commas

% Remove params that have not been used (nargin < number of params)
nActualArguments= evalin('caller', 'nargin');
if nActualArguments < numel(token)
  token(nActualArguments+1:end) = [];
end

ht{N}.fcn_params = token(~ismember(token,'varargin','legacy'));

% Get argument values for named arguments
for ii=1:numel(ht{N}.fcn_params)
  if strcmp(ht{N}.fcn_params{ii},objname)
    ht{N}.fcn_params{ii} = 'obj';  % always call the dat/epo/cnt struct obj
  else
    s =evalin('caller',['whos(''' token{ii} ''')']);  % Get size of variable
    if s.bytes > maxSize
       ht{N}.(ht{N}.fcn_params{ii}) = sprintf('variable exceeds maximum size %d kb',maxSize/1000);
    else
       ht{N}.(ht{N}.fcn_params{ii}) = evalin('caller',ht{N}.fcn_params{ii});
    end
  end
end

% Get optional arguments (varargin)
if any(ismember(token,'varargin','legacy'))
  va = evalin('caller','varargin');
  % Check size of variable
  s=evalin('caller','whos(''varargin'')');
  if s.bytes > maxSize
     ht{N}.varargin_note = sprintf('varargin exceeds maximum size %d kb',maxSize/1000);
  else
    for ii=1:numel(va)
        ht{N}.(sprintf('varargin%d',ii)) = va{ii};
    end
  end
end

ht{N}.date= datestr(now,0);
obj.history = ht;

author: bbci


text:
function [x,fval,exitflag,info] = opti_hybrj(fun,grad,x0,opts)
%OPTI_HYBRJ Solve a SNLE using HYBRJ (Powell-Hybrid MINPACK Routine)
%
%   F(x) = 0
%
%   x = opti_hybrj(fun,grad,x0) solves a SNLE where fun is a vector of
%   functions to solve. grad is an optional gradient of the nonlinear
%   equations and x0 is a starting guess.
%
%   x = opti_hybrj(fun,grad,x0,opts) uses opts to pass optiset options to
%   the solver.
%
%   [x,fval,exitflag,info] = opti_hybrj(...) returns the objective value at
%   the solution, together with the solver exitflag, and an information
%   structure.
%
%   THIS IS A WRAPPER FOR HYBRJ + HYBRD
%   See supplied License

%   Copyright (C) 2012 Jonathan Currie (I2C2)

if(nargin < 4), opts = optiset; end
if(nargin < 3), error('HYBRJ requires at least 3 arguments'); end

%Setup display level
opts.display = dispLevel(opts.display);

%Check we have a valid x0
if(isempty(x0) || any(isnan(x0)))
    error('HYBRJ requires an initial guess, x0!');
end

t = tic;
% Run HYBRJ
[x, fval, exitflag, iter] = hybrj(fun,grad,x0,opts);

%Collect Results
info.FuncEvals = iter;
info.Time = toc(t);
info.Algorithm = 'HYBRJ: MINPACK Powell Hybrid';

switch(exitflag)
    case 1
        info.Status = 'Optimal';
    case 0
        info.Status = 'Exceeded Function Evaluations';
    case -1
        info.Status = 'Infeasible / Could not Converge';
    case -2
        info.Status = 'HYBRJ Error';
    case -5
        info.Status = 'User Exited';
    otherwise
        info.Status = 'HYBRJ Error';
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function F = calc_expdecay2(vX, vXdata)
%
% Compute exponential decay combined with x^2 function
%F = vX(1)*exp(-vX(2)*vXdata)+vX(3)*exp(vX(4)*vXdata)
% F = vX(1)+vX(2).*vXdata+vX(3).*(vXdata).^2;
F = vX(1).*exp(vX(2).*vXdata);

author: CelsoReyes


text:
function [out] = coeff2func(alpha,rho,x)
% out = coeff2func(alpha,rho,x)
% Returns the values at x of the function
% \sum_{p} alpha_p e_p(r)
% Where e_p(r) = C_p J_0(rho_p r) (see function Cp.m for a definition of
% C_p).

y = x(:)';
C = Cp(rho);

J0vals = besselj(0,rho(:)*y);
out = J0vals'*(C(:).*alpha(:));

reshape(out,size(x,1),size(x,2));

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [J, grad] = lrCostFunction(theta, X, y, lambda)
%LRCOSTFUNCTION Compute cost and gradient for logistic regression with
%regularization
%   J = LRCOSTFUNCTION(theta, X, y, lambda) computes the cost of using
%   theta as the parameter for regularized logistic regression and the
%   gradient of the cost w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Hint: The computation of the cost function and gradients can be
%       efficiently vectorized. For example, consider the computation
%
%           sigmoid(X * theta)
%
%       Each row of the resulting matrix will contain the value of the
%       prediction for that example. You can make use of this to vectorize
%       the cost function and gradient computations.
%
% Hint: When computing the gradient of the regularized cost function,
%       there're many possible vectorized solutions, but one solution
%       looks like:
%           grad = (unregularized gradient for logistic regression)
%           temp = theta;
%           temp(1) = 0;   % because we don't add anything for j = 0
%           grad = grad + YOUR_CODE_HERE (using the temp variable)
%

addpath('../ex2');
[J, grad] = costFunctionReg(theta, X, y, lambda);

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

end

author: SaveTheRbtz

repo: ml-class


text:
function I = localSum( I, dims, shape, op )
% Fast routine for box, max and min filtering.
%
% Same effect as calling 'C=convn( I, ones(dims), shape)', except more
% efficient. Computes local sums by using running sums. To get sum in
% non-overlapping windows, use shape='block'. Equivalent to doing localSum,
% and then subsampling (except more efficient). If operation op is set to
% 'max' or 'min', computes local maxes or mins instead of sums. Note, that
% when applicable convBox and convMax are significantly faster.
%
% USAGE
%  I = localSum( I, dims, [shape], [op] )
%
% INPUTS
%  I       - matrix to compute sum over
%  dims    - size of volume to compute sum over, can be scalar
%  shape   - ['full'] 'valid', 'full', 'same', or 'block'
%  op      - ['sum'] 'max', or 'min'
%
% OUTPUTS
%  C       - matrix of sums
%
% EXAMPLE - 1
%  A=rand(500,500,1); dim=25; f=ones(dim,1); shape='same'; r=20;
%  tic, for i=1:r, B = localSum(A,dim,shape); end; toc
%  tic, for i=1:r, C = conv2(conv2(A,f,shape),f',shape); end; toc
%  diff=B-C; im(diff), sum(abs(diff(:)))
%
% EXAMPLE - 2
%  load trees; I=ind2gray(X,map); figure(1); im(I);
%  I1=localSum(I,3,'block','sum'); figure(2); im(I1); title('sum')
%  I2=localSum(I,3,'block','max'); figure(3); im(I2); title('max')
%  I3=localSum(I,3,'block','min'); figure(4); im(I3); title('min')
%
% See also convBox, convMax, imShrink
%
% Piotr's Computer Vision Matlab Toolbox      Version 3.22
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<3 || isempty(shape)), shape='full'; end
if( nargin<4 || isempty(op)), op='sum'; end
assert(any(strcmp(shape,{'same','valid','full','block'})));
assert(any(strcmp(op,{'sum','max','min'})));

% standard case (use nlfiltersep)
if( nargin>=3 && strcmp(shape,'block') )
  switch op
    case 'sum', fun = @sum;
    case 'max', fun = @(x,d) max(x,[],d);
    case 'min', fun = @(x,d) min(x,[],d);
  end
else
  switch op
    case 'sum', fun = @nlfiltersep_sum;
    case 'max', fun = @nlfiltersep_max;
    case 'min', error('min not implemented');
  end
end
I = nlfiltersep( I, dims, shape, fun );

end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function update_editfield_value(src,~)
    %update_editfield_value general function to update x.Value when x.String is updated, where x is an editfield
    % use as the callback
    %
    %  uicontrol('Style','edit','string','123','Value',num2str(123),'Callback',@update_editfield_value);
    %
    %  whenever the string (text) is changed, then the value updates, too.
    %
    %  only works for single values, not arrays, or non-numeric text.
    try
        src.Value = str2double(src.String);
    catch
        warning('unable to convert "%s" to a number',src.String);
        src.Value = nan;
    end
end


text:
function obj2 = unitaryopp(f, obj1, varargin);

    % make new memory mapped data file (blank)
    % ----------------------------------------
    newFileName = mmo.getnewfilename;
    data = memmapfile(obj1.dataFile, 'writable', obj1.writable, 'format', { 'single' obj1.dimensions 'x' });
    fid = fopen(newFileName, 'w');
    if obj1.dimensions(end) == 1, obj1.dimensions(end) = []; end;
    ss.type = '()';
    ss.subs(1:length(obj1.dimensions)-1) = { ':' };
    for index = 1:obj1.dimensions(end)
        ss.subs{length(obj1.dimensions)} = index;
        tmpdata = f(subsref(data.Data.x, ss), varargin{:});
        fwrite(fid, tmpdata, 'float');
    end;
    fclose(fid);

    % copy the data
    % -------------
    obj2 = mmo(newFileName, obj1.dimensions, true, obj1.transposed);
    obj2 = updateWorkspace(obj2);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
classdef MimCopyUID < MimGuiPlugin
    % MimCopyUID. Copy the UID of the current dataset to the clipboard
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Copy UID'
        SelectedText = 'Copy UID'
        ToolTip = 'Copy the UID of the current dataset to the clipboard'
        Category = 'Developer tools'
        Visibility = 'Dataset'
        Mode = 'Segment'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 1

        Location = 103
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            image_info = gui_app.GetImageInfo;
            uid = image_info.ImageUid;
            disp(['Current dataset UID is: ' uid]);
            clipboard('copy', uid);
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.DeveloperMode && gui_app.IsDatasetLoaded;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


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

clear all;
%close all;
clf

num=[1 1];
den=conv([1 0],[0.1 -1]);
w=logspace(-1,2,100);
[m,p]=bode(num,den,w);
subplot(2,1,1),loglog(w,m,'-',w,ones(size(w)),'-');
ylabel('Magnitude');
title('Fig. 6.30 Bode Plot for G=(s+1)/[s(s/10 -1)] (a) magnitude');
bodegrid;
subplot(2,1,2)
if p>0, p=p-360; end      % Matlab quadrant control varies in different versions
semilogx(w,p,'-',w,-180*ones(size(w)),'-');
axis([.1 100 -270 -90])

xlabel('\omega (rad/sec)');
ylabel('Phase (deg)');
title('(b) phase');
bodegrid;

author: Sable