text:
function box_plot ( point_filename, header )

%*****************************************************************************80
%
%% BOX_PLOT plots integer pairs as filled boxes.
%
%  Discussion:
%
%    The input file contains pairs of integer coordinates.
%
%  Usage:
%
%    box_plot ( 'point_file_name' )
%
%    reads 5 values per line (x,y,r,g,b) from the file,
%    and plots a box of color rgb at position xy.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string POINT_FILENAME, the name of the file containing the
%    coordinates of the points.
%
%    Input, string HEADER, an optional title for the plot.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BOX_PLOT\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Read a file of POINT_NUM points in 2 dimensions;\n' );
  fprintf ( 1, '  Display the points in a MATLAB graphics window.\n' );
%
%  First argument is the point file.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BOX_PLOT:\n' );
    point_filename = input ( '  Enter the name of the point file: ' );
  end
%
%  Second item is the optional header.
%
  if ( nargin < 2 )
    header = s_escape_tex ( point_filename );
  end
%
%  Read the data.
%
  [ dim_num, point_num ] = r8mat_header_read ( point_filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the header of "%s".\n', point_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Data dimension DIM_NUM   = %d\n', dim_num );
  fprintf ( 1, '  Number of points POINT_NUM  = %d\n', point_num );

  xyrgb = r8mat_data_read ( point_filename, dim_num, point_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the data in "%s".\n', point_filename );
%
%  Extract the XY data and guarantee that the data is integral.
%
  xy(1:2,1:point_num) = round ( xyrgb(1:2,1:point_num) );

  i4mat_transpose_print_some ( 2, point_num, xy, 1, 1, ...
    2, 5, '  A portion of the XY data:' );
%
%  Extract the RGB data.
%
  rgb(1:3,1:point_num) = xyrgb(3:5,1:point_num);

  r8mat_transpose_print_some ( 3, point_num, rgb, 1, 1, ...
    3, 5, '  A portion of the RGB data:' );
%
%  Make the plot.
%
  x_data_min = min ( xy(1,:) );
  x_data_max = max ( xy(1,:) );
  y_data_min = min ( xy(2,:) );
  y_data_max = max ( xy(2,:) );

  x_plot_min = x_data_min - 0.5;
  x_plot_max = x_data_max + 0.5;

  y_plot_min = y_data_min - 0.5;
  y_plot_max = y_data_max + 0.5;

  x_plot_range = x_plot_max - x_plot_min;
  y_plot_range = y_plot_max - y_plot_min;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Data Min:  %12f  %12f\n', x_data_min, y_data_min );
  fprintf ( 1, '  Data Max:  %12f  %12f\n', x_data_max, y_data_max );

  margin = 0.025 * max ( x_plot_range, y_plot_range );

  x_axes_min = x_plot_min - margin;
  x_axes_max = x_plot_max + margin;
  y_axes_min = y_plot_min - margin;
  y_axes_max = y_plot_max + margin;

  clf
%
%  Every box starts out WHITE.
%
  for x = x_data_min : x_data_max

    for y = y_data_min : y_data_max

      a = x - 0.44;
      b = x + 0.44;
      c = y - 0.44;
      d = y + 0.44;

      patch ( [ a, b, b, a ], [ c, c, d, d ], 'w', 'EdgeColor', 'none' );

    end

  end
%
%  Now draw boxes of the user's selected color.
%
  for point = 1 : point_num

    a = xy(1,point) - 0.44;
    b = xy(1,point) + 0.44;
    c = xy(2,point) - 0.44;
    d = xy(2,point) + 0.44;
    patch ( [ a, b, b, a ], [ c, c, d, d ], rgb(1:3,point)', 'EdgeColor', 'none');

  end
%
%  TEMPORARY!
%  Draw a line to indicate desired accuracy level for sparse grids.
%
  if ( 0 )
    line ( [-0.5, 7.5 ], [ 7.5, -0.50 ], 'LineWidth', 3, 'Color', 'w' )
  end
  if ( 0 )
    line ( [-0.5, 9.5 ], [ 9.5, -0.50 ], 'LineWidth', 3, 'Color', 'w' )
  end
  if ( 0 )
    line ( [-0.5, 11.5 ], [ 11.5, -0.50 ], 'LineWidth', 3, 'Color', 'w' )
  end
  if ( 1 )
    line ( [-0.5, 13.5 ], [ 13.5, -0.50 ], 'LineWidth', 3, 'Color', 'w' )
  end
%
%  The TITLE function will interpret underscores in the title.
%  We need to unescape such escape sequences!
%
  if ( 0 )
    xlabel ( '<--- I --->', 'Fontsize', 16 )
    ylabel ( '<--- J --->', 'Fontsize', 16 )
    title ( header, 'Fontsize', 24 )
  end

  axis ( [ x_axes_min, x_axes_max, y_axes_min, y_axes_max ] );
  axis equal
  axis square
  axis tight

  plot_filename = filename_ext_swap ( point_filename, 'png' );
  print ( '-dpng', plot_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A copy of this image was saved as "%s"\n', plot_filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BOX_PLOT\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function column_num = file_column_count ( input_file_name )

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

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

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( s_len_trim ( line ) == 0 )

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

    else
      got_one = TRUE;
      break;
    end

  end

  fclose ( input_unit );

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

  column_num = s_word_count ( line );

  return
end
function filename_new = filename_ext_swap ( filename, ext )

%*****************************************************************************80
%
%% FILENAME_EXT_SWAP replaces the current "extension" of a file name.
%
%  Discussion:
%
%    The "extension" of a filename is the string of characters
%    that appears after the LAST period in the name.  A file
%    with no period, or with a period as the last character
%    in the name, has a "null" extension.
%
%  Example:
%
%          Input           Output
%    ================     =============
%    FILENAME     EXT     FILENAME_NEW
%
%    bob.for      obj     bob.obj
%    bob.bob.bob  txt     bob.bob.txt
%    bob          yak     bob.yak
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 August 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILENAME, a file name.
%    On output, the extension of the file has been changed.
%
%    Input, string EXT, the extension to be used on the output
%    copy of FILE_NAME, replacing the current extension if any.
%
%    Output, string FILENAME_NEW, a copy of the input file name,
%    with the new extension.
%
  filename_len = length ( filename );

  ext_len = length ( ext );

  period = filename_len + 1;

  for i = filename_len : -1 : 1
    if ( filename(i:i) == '.' )
      period = i;
      break
    end
  end

  filename_new(1:period-1) = filename(1:period-1);
  filename_new(period) = '.';
  filename_new(period+1:period+ext_len) = ext(1:ext_len);

  return
end
function row_num = file_row_count ( input_file_name )

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

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

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

  record_num = 0;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

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

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

  end

  fclose ( input_unit );

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

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

  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d  ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );
    fprintf ( 1, '\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d  ', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%7d  ', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

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

%*****************************************************************************80
%
%% R8MAT_DATA_READ reads data from an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns of data.
%
%    Output, real TABLE(M,N), the point coordinates.
%
  table = zeros ( m, n );
%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

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

  input_unit = fopen ( input_filename );

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

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      break;
    end

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

    elseif ( s_len_trim ( line ) == 0 )

    else

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

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

    end

  end

  fclose ( input_unit );

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

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

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

  n = file_row_count ( input_filename );

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

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

%*****************************************************************************80
%
%% R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, an optional title.
%
  incx = 5;

  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d       ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d ', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%12f', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

  return
end
function s2 = s_escape_tex ( s1 )

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

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

  while ( s1_pos < s1_length )

    s1_pos = s1_pos + 1;

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

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

  end

  return
end
function len = s_len_trim ( s )

%*****************************************************************************80
%
%% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LEN, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end

  return
end
function 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


text:
function M = spm_mesh_isosurface(V, t, s)
% Compute isosurface geometry from volume data
% FORMAT M = spm_mesh_isosurface(V, t, s)
% V        - volume data
%            spm_vol struct, nifti object or 3D array
% t        - isosurface value
% s        - Gaussian filter width (FWHM) in {edges} [Default: 0]
%
% M        - patch structure
%
% This is merely a wrapper around isosurface.
%__________________________________________________________________________
% Copyright (C) 2019 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_mesh_isosurface.m 7677 2019-10-24 09:55:24Z guillaume $

if ischar(V) || isstruct(V)
    V = spm_vol(V);
    V = struct('dat',spm_read_vols(V),'mat',V.mat);
elseif isa(V,'nifti')
    V = struct('dat',full(V.dat),'mat',V.mat);
elseif isnumeric(V) || islogical(V)
    V = struct('dat',V,'mat',eye(4));
else
    error('Invalid volume data type.');
end

if nargin < 3, s = 0; end

if any(s)
    spm_smooth(V.dat, V.dat, s);
end

for i=1:numel(t)

    [faces,vertices] = isosurface(V.dat, t(i));

    if isempty(vertices)
        faces    = zeros(0,3);
        vertices = zeros(0,3);
    end

    % Swap around x and y because isosurface uses meshgrid and not ndgrid
    mat      = V(1).mat(1:3,:) * [0 1 0 0;1 0 0 0;0 0 1 0; 0 0 0 1];
    vertices = (mat * [vertices'; ones(1,size(vertices,1))])';

    M(i) = struct('faces',faces, 'vertices',vertices);

end

author: spm


text:
classdef CNSDEDVC < ALGORITHM
% <multi> <real/integer> <robust>
% Constrained nondominated sorting differential evolution based on decision variable classification
% SN     ---     4 --- Number of perturbed solutions
% PN     ---     6 --- Number of perturbations
% TN     ---    15 --- Number of repeated times of perturbation
% theta  --- 0.001 --- Threshold for DVC operation
% eta    --- 0.001 --- Desired level of robustness

%------------------------------- Reference --------------------------------
% W. Du, W. Zhong, Y. Tang, W. Du, and Y. Jin, High-dimensional robust
% multi-objective optimization for order scheduling: A decision variable
% classification approach, IEEE Transactions on Industrial Informatics,
% 15(1): 293-304.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            [SN,PN,TN,theta,eta] = Algorithm.ParameterSet(4,6,15,0.001,0.001);

            %% Generate random population
            Population = Problem.Initialization();
            [HR,LR]    = DVC(Problem,Population,SN,PN,TN,theta);
            [~,FrontNo,CrowdDis] = EnvironmentalSelection(Problem,Population,Problem.N,false,eta);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                if ~isempty(HR)
                    for subgen = 1 : 10
                        OffDec       = Population(TournamentSelection(2,end,FrontNo,-CrowdDis)).decs;
                        NewDec       = OperatorDE(Problem,Population.decs,Population(randi(end,1,end)).decs,Population(randi(end,1,end)).decs,{0.9,0.5,1,20});
                        OffDec(:,HR) = NewDec(:,HR);
                        Offspring    = Problem.Evaluation(OffDec);
                        [Population,FrontNo,CrowdDis] = EnvironmentalSelection(Problem,[Population,Offspring],Problem.N,true,eta);
                    end
                end
                if ~isempty(LR)
                    for subgen = 1 : 2
                        OffDec       = Population(TournamentSelection(2,end,FrontNo,-CrowdDis)).decs;
                        NewDec       = OperatorDE(Problem,Population.decs,Population(randi(end,1,end)).decs,Population(randi(end,1,end)).decs,{0.9,0.5,1,20});
                        OffDec(:,LR) = NewDec(:,LR);
                        Offspring    = Problem.Evaluation(OffDec);
                        [Population,FrontNo,CrowdDis] = EnvironmentalSelection(Problem,[Population,Offspring],Problem.N,false,eta);
                    end
                end
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%*****************************************************************************************************
%NAME: eplot.m
%AUTHOR: Andri M. Gretarsson
%DATE: 01/21/98
%
%SYNTAX: eplot(X,Y, 'colour')
%
%Note that 'colour' is NOT an optional argument.
%
%This function acts just like the built-in function 'plot' but plots error-bars. The error-bars
%are plotted in the colour given by 'colour'.  'X' and 'Y' are Nx2 matrixes, the first column
%representing the values of the coordinate (x or y), the second column representing the uncertainty
%('error') of those values.  'colour' is a one-letter string which must be one of the letters allowed
%in the built-in matlab function 'plot', to specify the plot colour.  Note that the function exits
%with "hold" set to "off".
%
%This function does not print points in addition to the error bars.  Where the error bars cross, is
%the coordinate point.  This means that if both error bars are exceedingly small complared to the
%coordinate values, the mark will be correspondingly small.  In such situations, it may be better to
%use 'plot' directly and specify that the error is smaller than the size of the mark.
%
%EXAMPLE:
%
%X=[1.0	0.2
%   2.0	0.2]
%Y=[1.0 	0.25
%   2.0	0.25]
%eplot(X,Y,'g')
%
%plots a green cross of width 0.2 and height 0.25 at coordinate (1.0,1.0), and a cross of width 0.2
%and height 0.25 at coordinate( (2.0,2.0).
%
%LAST MODIFIED:  01/21/98
%*****************************************************************************************************

function eplot=linearplot(x,y,colourstring)

    xvalue=x(:,1);												%For clarity
    xerror=x(:,2);
    yvalue=y(:,1);
    yerror=y(:,2);

    plot(xvalue-xerror,yvalue-yerror,'w-',xvalue+xerror,yvalue+yerror,'w-'); hold on;
    %Sets appropriate axes but otherwise invisible on a white background

    for i=1:length(xvalue)
        plot([xvalue(i)-xerror(i) xvalue(i)+xerror(i)], [yvalue(i) yvalue(i)], colourstring);
        plot([xvalue(i) xvalue(i)], [yvalue(i)-yerror(i) yvalue(i)+yerror(i)], colourstring);
    end
    hold off;

author: CelsoReyes


text:
function [ w, x ] = lyness_rule ( rule, order )

%*****************************************************************************80
%
%% LYNESS_RULE returns the points and weights of a Lyness quadrature rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 October 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    James Lyness, Dennis Jespersen,
%    Moderate Degree Symmetric Quadrature Rules for the Triangle,
%    Journal of the Institute of Mathematics and its Applications,
%    Volume 15, Number 1, February 1975, pages 19-32.
%
%  Parameters:
%
%    Input, integer RULE, the index of the rule.
%
%    Input, integer ORDER, the order of the rule.
%
%    Output, real W(ORDER), the weights.
%
%    Output, real X(2,ORDER), the points.
%

%
%  Get the suborder information.
%
  suborder_num = lyness_suborder_num ( rule );

  suborder = lyness_suborder ( rule, suborder_num );

  [ sub_xyz, sub_w ] = lyness_subrule ( rule, suborder_num );
%
%  Expand the suborder information to a full order rule.
%
  x = zeros ( 2, order );
  w = zeros ( 1, order );

  o = 0;

  for s = 1 : suborder_num

    if ( suborder(s) == 1 )

      o = o + 1;
      x(1:2,o) = sub_xyz(1:2,s);
      w(o) = sub_w(s);

    elseif ( suborder(s) == 3 )

      for k = 1 : 3
        o = o + 1;
        x(1,o) = sub_xyz ( i4_wrap(k,  1,3), s );
        x(2,o) = sub_xyz ( i4_wrap(k+1,1,3), s );
        w(o) = sub_w(s) / 3.0;
      end

    elseif ( suborder(s) == 6 )

      for k = 1 : 3
        o = o + 1;
        x(1,o) = sub_xyz ( i4_wrap(k,  1,3), s );
        x(2,o) = sub_xyz ( i4_wrap(k+1,1,3), s );
        w(o) = sub_w(s) / 6.0;
      end

      for k = 1 : 3
        o = o + 1;
        x(1,o) = sub_xyz ( i4_wrap(k+1,1,3), s );
        x(2,o) = sub_xyz ( i4_wrap(k,  1,3), s );
        w(o) = sub_w(s) / 6.0;
      end

    else

      fprintf ( 1, '\n' );
      fprintf ( 1, 'LYNESS_RULE - Fatal error!\n' );
      fprintf ( 1, '  Illegal SUBORDER(%d) = %d\n', s, suborder(s) );
      error ( 'LYNESS_RULE - Fatal error!' );

    end

  end

  return
end

repo: jburkardt-m


text:
function [imgs, ds, pix_resolutions] = ptv_create_pyramid(img, nlevels, pix_resolution, varargin)
%    create_iso_pyramid(img, nlevels, pix_resolution, varargin, {force_isotropic, k_down, interp_type, fine_rescale})
%
    if isempty(img)
        imgs = [];
        ds = [];
        pix_resolutions = [];
        return;
    end
    interp_type = 0;
    force_isotropic = false;
    if nargin >= 4
        force_isotropic = varargin{1};
    end
    if force_isotropic
        error('not implemented');
    end
    k_down = 0.5;
    if nargin >= 5
        k_down = varargin{2};
    end
    if nargin >= 6
        interp_type = varargin{3};
    end
    fine_rescale = false;
    if nargin >= 7
        fine_rescale = varargin{4};
    end

    ndims = 2;
    if size(img, 3) > 1
        ndims = 3;
    end

    Nch = size(img, 4);
    Nimgs = size(img, 5);
    pix_resolutions = zeros(nlevels, 3);

    pix_resolutions(1, :) = pix_resolution;
    imgs = cell(nlevels, 1);
    imgs{1} = img;
    ds = ones(nlevels, ndims);
    for i = 2 : nlevels
        if ~force_isotropic
            d = ones(1, ndims);
        end

        idxs = cell(ndims, 1);
        newsz = zeros(ndims, 1);
        if fine_rescale
            sigmas = 0.5 * ones(ndims, 1) * (0.5/k_down * (i-1));
        else
            sigmas = 0.4 * ones(ndims, 1) * (0.5/k_down);
        end

        newsz = ceil(k_down * size(imgs{i-1}(:,:,:, 1)));
        if ndims == 3
            imgs{i} = zeros([newsz, Nch, Nimgs]);
        elseif ndims == 2
            imgs{i} = zeros([newsz, 1, Nch, Nimgs]);
        end
        pix_resolutions(i, :) = pix_resolutions(i - 1, :) / k_down;

        for jj = 1 : Nch
            for ii = 1 : Nimgs
                if ndims == 2
                    if fine_rescale
                        t = imgaussfilt(imgs{1}(:,:, 1, jj, ii), sigmas);
                    else
                        t = imgaussfilt(imgs{i-1}(:,:, 1, jj, ii), sigmas);
                    end
                    imgs{i}(:,:, :, jj, ii) = imresize(t, size(imgs{i}(:,:, 1, jj, ii)), 'bilinear');
                elseif ndims == 3
                    if fine_rescale
                        t = imgaussfilt3(imgs{1}(:,:,:, jj, ii), sigmas);
                    else
                        t = imgaussfilt3(imgs{i-1}(:,:,:, jj, ii), sigmas);
                    end
                    imgs{i}(:,:,:, jj, ii) = volresize(t, size(imgs{i}(:,:,:, jj, ii)), interp_type);
                end
            end
        end
    end
end

repo: pTVreg


text:
function [EI,Pth,Qthth] = lme_mass_RgEI1(X,Zcols,W,CBhat,L,phi,ni,G,GDa,GDb)
% [EI,Pth,Qthth] = lme_mass_RgEI1(X,Zcols,W,CBhat,L,phi,ni,G,GDa,GDb)
%
% Expected information matrix of the restricted log-likelihood of a whole
% region. This is less computationally efficient than lme_mass_RgEI
% but more numerically stable.
%
% Input
% X: Ordered design Matrix (according to time for each subject).
% Zcols: Vector with the indices of the colums of X that will be considered
% as random effects.
% W: Inverses of the estimated temporal covariance matrices for each
% subject stacked in W.
% CBhat: Asymptotic covariance matrix of the fixed effects.
% L: Cholesky factor of the covariance matrix of the random effects (D).
% phi: Within-subject standard deviation of the errors.
% ni: Vector whose entries are the number of repeated measures for each
% subject in the study (ordered according to X).
% G: Spatial covariance matrix.
% GDa: Derivative of the spatial covariance matrix for the first spatial
% parameter.
% GDb: Derivative of the spatial covariance matrix for the second spatial
% parameter (empty for spatial models with a single parameter).
%
% Output
% EI: Expected information matrix.
% Pth,Qthth: Matrices that are useful for inferences on the fixed effects.
%
% $Revision: 1.2 $  $Date: 2015/01/06 17:14:55 $
% Original Author: Jorge Luis Bernal Rusiel
% CVS Revision Info:
%    $Author: mreuter $
%    $Date: 2015/01/06 17:14:55 $
%    $Revision: 1.2 $
% References: Bernal-Rusiel J.L., Greve D.N., Reuter M., Fischl B., Sabuncu
% M.R., 2012. Statistical Analysis of Longitudinal Neuroimage Data with Linear
% Mixed Effects Models, NeuroImage, doi:10.1016/j.neuroimage.2012.10.065.
%
m = length(ni);
n = sum(ni);
nv = size(G,1);
q = length(Zcols);
nth = q*(q+1)/2+1;
p = size(X,2);
Pth = zeros(nth,p,p);
Qthth = zeros(nth,nth,p,p);
Der = zeros(nth,n,max(ni));
%Computation of the first order derivatives of the temporal covariance matrix
jk = 0;
for k=1:q
    for j=1:k
        jk = jk + 1;
        posi = 1;
        for i=1:m
            posf = posi+ni(i)-1;
            Zi = X(posi:posf, Zcols);
            Zki = Zi(:,k);
            Mjki = Zki*L(j,:)*Zi';
            Mjki = Mjki + Mjki';
            Der(jk,posi:posf,1:ni(i)) = Mjki;
            posi = posf+1;
        end;
    end;
end;
posi = 1;
for i=1:m
    posf = posi+ni(i)-1;
    Der(nth,posi:posf,1:ni(i)) = 2*phi*eye(ni(i));
    posi = posf+1;
end;
%Computation of Pis,Qijs and the expected information matrix EI.
for j=1:nth
    posi = 1; Pj = 0;
    Bj = squeeze(Der(j,:,:));
    for i=1:m
        posf = posi+ni(i)-1;
        Wi = W(posi:posf,1:ni(i));
        Pj = Pj - X(posi:posf,:)'*Wi*Bj(posi:posf,1:ni(i))*Wi*X(posi:posf,:);
        posi = posf+1;
    end;
    Pth(j,:,:) = Pj;
end;
EI = zeros(nth+2,nth+2);
%Expected information among Lijs (including phi)
for k=1:nth
    Bk = squeeze(Der(k,:,:));
    Pk = squeeze(Pth(k,:,:));
    for j=1:k
        Bj = squeeze(Der(j,:,:));
        Pj = squeeze(Pth(j,:,:));
        posi = 1; Qkj = 0;
        traceBkj = 0;
        for i=1:m
            posf = posi+ni(i)-1;
            Wi = W(posi:posf,1:ni(i));
            Bkji = Wi*Bk(posi:posf,1:ni(i))*Wi*Bj(posi:posf,1:ni(i));
            traceBkj = traceBkj + trace(Bkji);
            Bkji = Bkji*Wi;
            Qkji = X(posi:posf,:)'*Bkji*X(posi:posf,:);
            Qkj = Qkj + Qkji;
            posi = posf+1;
        end;
        Qthth(k,j,:,:) = Qkj;
        Qthth(j,k,:,:) = Qkj;
        EI(k,j) = nv*(traceBkj - trace(CBhat*(2*Qkj-Pk*CBhat*Pj)));
        EI(j,k) = EI(k,j);
    end;
end;
%Expected information between Lijs (including phi) and a (first spatial parameter)
Mauxa = G\GDa;
trMauxa = trace(Mauxa);
for k=1:nth
    Bk = squeeze(Der(k,:,:));
    Pk = squeeze(Pth(k,:,:));
    posi = 1;
    traceBk = 0;
    for i=1:m
        posf = posi+ni(i)-1;
        Wi = W(posi:posf,1:ni(i));
        Bki = Wi*Bk(posi:posf,1:ni(i));
        traceBk = traceBk + trace(Bki);
        posi = posf+1;
    end;
    EI(k,nth+1) = trMauxa*(traceBk + trace(CBhat*Pk));
    EI(nth+1,k) = EI(k,nth+1);
end
%Expected information between a and a
EI(nth+1,nth+1) = (n-p)*trace(Mauxa*Mauxa);
if ~isempty(GDb)
    %Expected information between Lijs (including phi) and b (second spatial parameter)
    Mauxb = G\GDb;
    trMauxb = trace(Mauxb);
    for k=1:nth
        Bk = squeeze(Der(k,:,:));
        Pk = squeeze(Pth(k,:,:));
        posi = 1;
        traceBk = 0;
        for i=1:m
            posf = posi+ni(i)-1;
            Wi = W(posi:posf,1:ni(i));
            Bki = Wi*Bk(posi:posf,1:ni(i));
            traceBk = traceBk + trace(Bki);
            posi = posf+1;
        end;
        EI(k,nth+2) = trMauxb*(traceBk + trace(CBhat*Pk));
        EI(nth+2,k) = EI(k,nth+2);
    end
    %Expected information between b and b
    EI(nth+2,nth+2) = (n-p)*trace(Mauxb*Mauxb);
    %Expected information between a and b
    EI(nth+1,nth+2) = (n-p)*trace(Mauxa*Mauxb);
    EI(nth+2,nth+1) = EI(nth+1,nth+2);
else
    EI = EI(1:nth+1,1:nth+1);
end;
EI = 0.5*EI;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [X_poly] = polyFeatures(X, p)
%POLYFEATURES Maps X (1D vector) into the p-th power
%   [X_poly] = POLYFEATURES(X, p) takes a data matrix X (size m x 1) and
%   maps each example into its polynomial features where
%   X_poly(i, :) = [X(i) X(i).^2 X(i).^3 ...  X(i).^p];
%

% You need to return the following variables correctly.
X_poly = zeros(numel(X), p);

% ====================== YOUR CODE HERE ======================
% Instructions: Given a vector X, return a matrix X_poly where the p-th
%               column of X contains the values of X to the p-th power.
%
%

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

end

author: worldveil

repo: coursera-ml


text:
function [ t, rank ] = subset_next ( n, t, rank )

%*****************************************************************************80
%
%% SUBSET_NEXT computes the subset lexicographic successor.
%
%  Discussion:
%
%    This is a lightly modified version of "subset_lex_successor()" from COMBO.
%
%  Example:
%
%    On initial call, N is 5 and the input value of RANK is -1.
%    Then here are the successive outputs from the program:
%
%   Rank   T1   T2   T3   T4   T5
%   ----   --   --   --   --   --
%      0    0    0    0    0    0
%      1    0    0    0    0    1
%      2    0    0    0    1    0
%      3    0    0    0    1    1
%     ..   ..   ..   ..   ..   ..
%     30    1    1    1    1    0
%     31    1    1    1    1    1
%     -1    0    0    0    0    0  <-- Reached end of cycle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 May 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer N, the number of elements in the master set.
%    N must be positive.
%
%    Input/output, integer T(N), describes a subset.  T(I) is 0 if
%    the I-th element of the master set is not in the subset, and is
%    1 if the I-th element is part of the subset.
%    On input, T describes a subset.
%    On output, T describes the next subset in the ordering.
%
%    Input/output, integer RANK, the rank.
%    If RANK = -1 on input, then the routine understands that this is
%    the first call, and that the user wishes the routine to supply
%    the first element in the ordering, which has RANK = 0.
%    In general, the input value of RANK is increased by 1 for output,
%    unless the very last element of the ordering was input, in which
%    case the output value of RANK is -1.
%

%
%  Return the first element.
%
  if ( rank == -1 )
    t = zeros ( n, 1 );
    rank = 0;
    return
  end

  for i = n : -1 : 1

    if ( t(i) == 0 )
      t(i) = 1;
      rank = rank + 1;
      return
    else
      t(i) = 0;
    end

  end

  rank = -1;

  return
end


text:
classdef ListSelector < wt.abstract.BaseWidget & wt.mixin.Enableable &...
        wt.mixin.FontStyled & wt.mixin.ButtonColorable &...
        wt.mixin.FieldColorable
    % Select from an array of items and add them to a list

    % Copyright 2020-2021 The MathWorks Inc.

    %% Events
    events (HasCallbackProperty, NotifyAccess = protected)

        % Triggered when a button is pushed
        ButtonPushed

        % Triggered when the value of the list selection changes
        ValueChanged

        % Triggered when the highlighted value changes
        HighlightedValueChanged

    end %events

    %% Public properties
    properties (AbortSet)

        % List of items to add to the list
        Items (1,:) string = ["Item 1" "Item 2" "Item 3" "Item 4"]

        % Data associated with  items (optional)
        ItemsData (1,:)

        % Indicates whether to allow duplicate entries in the list
        AllowDuplicates  (1,1) matlab.lang.OnOffSwitchState = false

        % Indicates whether to allow sort controls %RAJ - Future feature
        %Sortable  (1,1) matlab.lang.OnOffSwitchState = true

        % Inidicates what to do when add button is pressed (select from
        % Items or custom using ButtonPushed event or ButtonPushedFcn)
        AddSource (1,1) wt.enum.ListAddSource = wt.enum.ListAddSource.Items

    end %properties

    properties (AbortSet, Dependent)

        % Indices of displayed items that are currently added to the list
        SelectedIndex (1,:)

        % The current selection
        Value (1,:)

        % The current highlighted selection
        HighlightedValue (1,:)

    end %properties

    properties (AbortSet, Dependent, UsedInUpdate = false)

        % Width of the buttons
        ButtonWidth

    end %properties

    %% Read-Only properties
    properties (SetAccess = private)

        % Additional user buttons
        UserButtons wt.ButtonGrid

    end %properties

    %% Internal Properties
    properties ( Transient, NonCopyable, ...
            Access = {?wt.abstract.BaseWidget, ?wt.test.BaseWidgetTest} )

        % The ListBox control
        ListBox (1,1) matlab.ui.control.ListBox

        % The list sorting buttons
        ListButtons wt.ButtonGrid

        % Listen to button pushes in sections
        ButtonPushedListener event.listener

    end %properties

    %% Protected methods
    methods (Access = protected)

        function setup(obj)

            % Call superclass setup to establish the main grid
            obj.setup@wt.abstract.BaseWidget();

            % Set default size
            obj.Position(3:4) = [120 130];

            % Configure grid
            obj.Grid.Padding = 3;
            obj.Grid.ColumnWidth = {'1x',25};
            obj.Grid.RowHeight = {106,'1x'};

            % Create the list buttons
            obj.ListButtons = wt.ButtonGrid(obj.Grid);
            obj.ListButtons.Icon = ["add_24.png", "delete_24.png", "up_24.png", "down_24.png"];
            obj.ListButtons.ButtonTag = ["Add", "Remove", "Up", "Down"];
            obj.ListButtons.Layout.Column = 2;
            obj.ListButtons.Layout.Row = 1;
            obj.ListButtons.Orientation = "vertical";
            obj.ListButtons.ButtonHeight(:) = {25};

            % Create an additional button grid for custom buttons
            obj.UserButtons = wt.ButtonGrid(obj.Grid,"Icon",[]);
            obj.UserButtons.Layout.Column = 2;
            obj.UserButtons.Layout.Row = 2;
            obj.UserButtons.Orientation = "vertical";

            % Create the ListBox
            obj.ListBox = uilistbox(obj.Grid);
            obj.ListBox.Multiselect = true;
            obj.ListBox.ValueChangedFcn = @(h,e)obj.onSelectionChanged(e);
            obj.ListBox.Layout.Column = 1;
            obj.ListBox.Layout.Row = [1 2];

            % Update listeners
            obj.ButtonPushedListener = event.listener(...
                [obj.ListButtons obj.UserButtons],...
                'ButtonPushed',@(h,e)obj.onButtonPushed(e) );

            % Update the internal component lists
            obj.BackgroundColorableComponents = [obj.ListButtons, obj.UserButtons];
            obj.FontStyledComponents = [obj.ListBox, obj.UserButtons, obj.ListButtons];
            obj.EnableableComponents = [obj.ListBox, obj.UserButtons, obj.ListButtons];
            obj.ButtonColorableComponents = [obj.UserButtons obj.ListButtons];
            obj.FieldColorableComponents = [obj.ListBox];

        end %function

        function update(obj)

            % What is selected?
            selIdx = obj.SelectedIndex;

            % Update the list
            obj.ListBox.Items = obj.Items(selIdx);
            obj.ListBox.ItemsData = selIdx;

            % Update button enable states
            obj.updateEnables();

        end %function

        function updateEnables(obj)

            % Button enables
            if obj.Enable

                % What is selected?
                selIdx = obj.SelectedIndex;

                % Highlighted selection in list?
                hiliteIdx = obj.getListBoxSelectedIndex();

                % How many items selected into list
                numRows = numel(selIdx);
                numHilite = numel(hiliteIdx);

                % Toggle button enables
                obj.ListButtons.ButtonEnable = [
                    obj.AllowDuplicates || ( numel(selIdx) < numel(obj.Items) ) %Add Button
                    ~isempty(hiliteIdx) % Delete Button
                    numHilite && ( hiliteIdx(end) > numHilite ) %Up Button
                    numHilite && ( hiliteIdx(1) <= (numRows - numHilite) ) %Down Button
                    ];

            end %if obj.Enable

        end %function

        function onSelectionChanged(obj,evt)

            % Get the new and old values
            if isempty(obj.ItemsData)
                itemsData = obj.Items;
            else
                itemsData = obj.ItemsData;
            end

            if isempty(evt.PreviousValue)
                oldValue = itemsData([]);
            else
                oldValue = itemsData(evt.PreviousValue);
            end

            if isempty(evt.Value)
                newValue = itemsData([]);
            else
                newValue = itemsData(evt.Value);
            end

            % Update button enable states
            obj.updateEnables();

            % Trigger event
            evtOut = wt.eventdata.ValueChangedData(newValue, oldValue);
            notify(obj,"HighlightedValueChanged",evtOut);

        end %function

        function onButtonPushed(obj,evt)

            % Which button?
            switch evt.Tag

                case 'Add'

                    switch obj.AddSource

                        case wt.enum.ListAddSource.Items
                            obj.promptToAddListItems()

                        case wt.enum.ListAddSource.ButtonPushedFcn
                            notify(obj,"ButtonPushed",evt);

                    end %switch obj.AddSource

                case 'Remove'
                    obj.removeListBoxSelection();

                case 'Up'
                    obj.shiftListBoxIndex(-1);

                case 'Down'
                    obj.shiftListBoxIndex(1);

                otherwise
                    % Trigger event for user buttons
                    notify(obj,"ButtonPushed",evt);

            end %switch

            % Request update
            obj.requestUpdate();

        end %function

        function promptToAddListItems(obj)
            % Prompt a dialog to add items to the listbox

            % Prompt for stuff to add
            if obj.AllowDuplicates
                newSelIdx = listdlg("ListString",obj.Items);
            else
                newSelIdx = listdlg(...
                    "ListString",obj.Items,...
                    "InitialValue",obj.ListBox.ItemsData);
            end

            if isempty(newSelIdx)
                % User cancelled
                return
            elseif obj.AllowDuplicates
                newSelIdx = [obj.SelectedIndex newSelIdx];
            end

            % Was a change made?
            if ~isequal(obj.SelectedIndex, newSelIdx)

                % Get the original value
                oldValue = obj.Value;

                % Make the update
                obj.SelectedIndex = newSelIdx;

                % Trigger event
                evtOut = wt.eventdata.ValueChangedData(obj.Value, oldValue);
                notify(obj,"ValueChanged",evtOut);

            end %if

        end %function

        function selIdx = getListBoxSelectedIndex(obj)
            % Get the current selected row indices in the listbox

            warnState = warning('off','MATLAB:structOnObject');
            s = struct(obj.ListBox);
            warning(warnState);
            selIdx = s.SelectedIndex;
            if isequal(selIdx, -1)
                selIdx = [];
            end

        end %function

        function removeListBoxSelection(obj)
            % Removes the currently selected items from the listbox

            % What's currently selected?
            idxSel = obj.getListBoxSelectedIndex();

            % Is there something to remove?
            if ~isempty(idxSel)

                % Get the original value
                oldValue = obj.Value;

                % Remove it
                obj.ListBox.Items(idxSel) = [];
                obj.ListBox.ItemsData(idxSel) = [];

                % Trigger event
                evtOut = wt.eventdata.ValueChangedData(obj.Value, oldValue);
                notify(obj,"ValueChanged",evtOut);

            end %if

        end %function

        function shiftListBoxIndex(obj, shift)
            % Shift selected items up/down within a listbox
            % This assumes ItemsData contains unique values

            % What is the current order?
            selIdx = obj.getListBoxSelectedIndex();

            % Make indices to all items as they are now
            idxNew = 1:numel(obj.ListBox.Items);
            idxOld = idxNew;

            % Find the last stable item that doesn't move
            [~,idxStable] = setdiff(idxNew, selIdx, 'stable');
            if ~isempty(idxStable)
                idxFirstStable = idxStable(1);
                idxLastStable = idxStable(end);
            else
                idxFirstStable = inf;
                idxLastStable = 0;
            end

            % Which way do we loop?
            if shift > 0 %Shift to end

                for idxToMove = numel(selIdx):-1:1

                    % Calculate if there's room to move this item
                    idxThisBefore = selIdx(idxToMove);
                    thisShift = max( min(idxLastStable-idxThisBefore, shift), 0 );

                    % Where does this item move from/to
                    idxThisAfter = idxThisBefore + thisShift;

                    % Where do other items move from/to
                    idxOthersBefore = selIdx(idxToMove)+1:1:idxThisAfter;
                    idxOthersAfter = idxOthersBefore - thisShift;

                    % Move the items
                    idxNew([idxThisAfter idxOthersAfter]) = idxNew([idxThisBefore idxOthersBefore]);

                end

            elseif shift < 0 %Shift to start

                for idxToMove = 1:numel(selIdx)

                    % Calculate if there's room to move this item
                    idxThisBefore = selIdx(idxToMove);
                    thisShift = min( max(idxFirstStable-idxThisBefore, shift), 0 );

                    % Where does this item move from/to
                    idxThisAfter = idxThisBefore + thisShift;

                    % Where do other items move from/to
                    idxOthersBefore = idxThisAfter:1:selIdx(idxToMove)-1;
                    idxOthersAfter = idxOthersBefore - thisShift;

                    % Move the items
                    idxNew([idxThisAfter idxOthersAfter]) = idxNew([idxThisBefore idxOthersBefore]);

                end

            end %if shift > 0

            % Was a change made?
            if ~isequal(idxOld, idxNew)

                % Get the original value
                oldValue = obj.Value;

                % Make the shift
                obj.ListBox.Items = obj.ListBox.Items(idxNew);
                obj.ListBox.ItemsData = obj.ListBox.ItemsData(idxNew);

                % Trigger event
                evtOut = wt.eventdata.ValueChangedData(obj.Value, oldValue);
                notify(obj,"ValueChanged",evtOut);

            end %if

        end %function

    end %methods

    %% Accessors
    methods

        function value = get.SelectedIndex(obj)
            value = obj.ListBox.ItemsData;
        end
        function set.SelectedIndex(obj,value)
            obj.ListBox.Items = obj.Items(value);
            obj.ListBox.ItemsData = value;
        end

        function value = get.Value(obj)
            if isempty(obj.ItemsData)
                value = obj.Items(:,obj.ListBox.ItemsData);
            else
                value = obj.ItemsData(:,obj.ListBox.ItemsData);
            end
        end
        function set.Value(obj,value)
            if isempty(value)
                obj.SelectedIndex = [];
            else
                if isempty(obj.ItemsData)
                    [tf, selIdx] = ismember(value, obj.Items);
                else
                    [tf, selIdx] = ismember(value, obj.ItemsData);
                end
                if ~all(tf)
                    warning("widgets:ListSelector:InvalidValue",...
                        "Attempt to set an invalid Value to the list.")
                    selIdx(~tf) = [];
                end
                obj.SelectedIndex = selIdx;
            end
        end

        function value = get.HighlightedValue(obj)
            selIdx = obj.ListBox.Value;
            if isempty(selIdx) || ~isnumeric(selIdx)
                selIdx = [];
            end
            if isempty(obj.ItemsData)
                value = obj.Items(:,selIdx);
            else
                value = obj.ItemsData(:,selIdx);
            end
        end
        function set.HighlightedValue(obj,value)
            if isempty(obj.ItemsData)
                [~, obj.ListBox.Value] = ismember(value, obj.Items);
            else
                [~, obj.ListBox.Value] = ismember(value, obj.ItemsData);
            end
        end

        function value = get.ButtonWidth(obj)
            value = obj.Grid.ColumnWidth{2};
        end
        function set.ButtonWidth(obj,value)
            obj.Grid.ColumnWidth{2} = value;
        end

    end %methods

end % classdef

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [xi,w]=fourthOrderPyramidCubPoints()
%%FOURTHORDERPYRAMIDCUBPOINTS Generate fourth-order cubature points for
%  integration over a 3-dimensional pyramid with a square base with the
%  peak vertex at (0,0,1) and the base vertices at (1,-1,-1), (-1,-1,-1),
%  (-1,1,-1), and (1,1,-1).
%
%INPUTS: None
%
%OUTPUTS: xi This is a 3XnumCubPoints set of points for the standard
%            square pyramid.
%          w A 1XnumCubPoints set of cubature weights. This sums to the
%            volume of the standard square pyramid (8/3).
%
%This function implements the points given in [1] (10 points).
%
%EXAMPLE:
%We compare a 4th-order moment computed using these cubature points
%to one computed using monomialIntPyramid. The results are the same within
%typical finite precision limits.
% [xi,w]=fourthOrderPyramidCubPoints();
% alpha=[2;0;2];
% theMoment=findMomentFromSamp(alpha,xi,w);
% intVal=monomialIntPyramid(alpha);
% RelErr=(theMoment-intVal)/intVal
%
%REFERENCES:
%[1] F. D. Witherden and P. E. Vincent, "On the identification of symmetric
%    quadrature rules for finite element methods," Computer and Mathematics
%    with Applications, vol. 69, no. 10, pp. 1232-1241, May 2015.
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

M=[                                     0,                                          0,   0.35446557777227471722730849524904581806,   0.30331168845504517111391728481208001144;
                                        0,                                          0,  -0.74972609378250711033655604388057044149,   0.55168907357213937275730716433358729608;
  0.6505815563982325146829577797417295398,                                          0,  -0.35523170084357268589593075201816127231,   0.28353223437153468006819777082540613962;
                                        0,    0.6505815563982325146829577797417295398,  -0.35523170084357268589593075201816127231,   0.28353223437153468006819777082540613962;
 -0.6505815563982325146829577797417295398,                                          0,  -0.35523170084357268589593075201816127231,   0.28353223437153468006819777082540613962;
                                        0,   -0.6505815563982325146829577797417295398,  -0.35523170084357268589593075201816127231,   0.28353223437153468006819777082540613962;
 0.65796699712169008954533549931479427127,   0.65796699712169008954533549931479427127,  -0.92150343220236930457646242598412224897,   0.16938424178833585063066278355484370017;
 0.65796699712169008954533549931479427127,  -0.65796699712169008954533549931479427127,  -0.92150343220236930457646242598412224897,   0.16938424178833585063066278355484370017;
-0.65796699712169008954533549931479427127,   0.65796699712169008954533549931479427127,  -0.92150343220236930457646242598412224897,   0.16938424178833585063066278355484370017;
-0.65796699712169008954533549931479427127,  -0.65796699712169008954533549931479427127,  -0.92150343220236930457646242598412224897,   0.16938424178833585063066278355484370017];

w=M(:,4);
xi=M(:,1:3)';

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function test_approx_test06 ( )

%*****************************************************************************80
%
%% TEST_APPROX_TEST06 uses cubic spline interpolation on all problems.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 August 2011
%
%  Author:
%
%    John Burkardt
%
  ibcbeg = 0;
  ibcend = 0;
  ybcbeg = 0.0;
  ybcend = 0.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_APPROX_TEST06\n' );
  fprintf ( 1, '  Cubic spline interpolation.\n' );

  prob_num = p00_prob_num ( );

  for prob = 1 : prob_num

    title = p00_title ( prob );

    data_num = p00_data_num ( prob );

    [ xdata, ydata ] = p00_dat ( prob, data_num );

    a = xdata(1);
    b = xdata(data_num);
%
%  Set up the interpolation function.
%
    ypp = spline_cubic_set ( data_num, xdata, ydata, ibcbeg, ybcbeg, ...
      ibcend, ybcend );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Problem %d\n', prob );
    fprintf ( 1, '  %s\n', title );
    fprintf ( 1, '\n' );
    fprintf ( 1, '    X   Y\n' );
    fprintf ( 1, '\n' );
%
%  Evaluate the interpolation function.
%
    for i = 1 : data_num - 1

      jmax = 3;

      if ( i == data_num - 1 )
        jhi = jmax;
      else
        jhi = jmax - 1;
      end

      for j = 1 : jhi

        xval = ( ( jmax - j     ) * xdata(i)     ...
               + (        j - 1 ) * xdata(i+1) ) ...
               / ( jmax     - 1 );

        [ yval, ypval, yppval ] = spline_cubic_val ( data_num, xdata, ydata, ...
          ypp, xval );

        if ( j == 1 || j == 3 )
          mark = '*';
        else
          mark = ' ';
        end

        fprintf ( 1, '  %c  %14g  %14g\n', mark, xval, yval );

      end

    end

  end

  return
end


text:
function A = mesh3d2 (n)
% create an n-by-n-by-n 3D mesh for the 2nd difference operator
% Example:
%   A = mesh3d2 (10) ;  % a 10-by-10-by-10 mesh
% See also: cs_demo

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

nn = 1:n^3 ;
ii = [nn-n^2 ; nn-n ; nn-1 ; nn ; nn+1 ; nn+n ; nn+n^2] ;
jj = repmat (nn, 7, 1) ;
xx = repmat ([-1 -1 -1 6 -1 -1 -1]', 1, n^3) ;
keep = find (ii >= 1 & ii <= n^3 & jj >= 1 & jj <= n^3) ;
A = sparse (ii (keep), jj (keep), xx (keep)) ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% inverse dynamics (recursive Newton-Euler) using spatial vector notation
function  tau = ID( model, q, qd, qdd )

    a_grav = SpatialAcceleration([0;0;-9.81;0;0;0]);

    for i = 1:model.NB
        [ XJ, S(:,i) ] = jcalc( model.pitch(i), q(i) );
        XJ
        S(:,i)
        vJ = SpatialVelocity(S(:,i)*qd(i));
        Xup(i) = XJ*model.Xtree(i);
        if model.parent(i) == 0
            v(i) = vJ;
            a(i) = Xup(i)*(-a_grav) + SpatialAcceleration(S(:,i)*qdd(i));
            a
        else
            v(i) = Xup(i)*v(model.parent(i)) + vJ;
            a(i) = Xup(i)*a(model.parent(i)) ...
                + SpatialAcceleration(S(:,i)*qdd(i)) ...
                + cross(v(i),vJ);
        end
        f(i) = model.I(i)*a(i) + cross( v(i), model.I(i)*v(i) );
    end

    v
    a
    f

    for i = model.NB:-1:1
        tau(i,1) = S(:,i)' * double(f(i));
        if model.parent(i) ~= 0
            f(model.parent(i)) = f(model.parent(i)) + Xup(i)*f(i);
        end
    end
end

function  [Xj,S] = jcalc( pitch, q )  %FIXED VW ORDER

    % jcalc  Calculate joint transform and motion subspace.
    % [Xj,S]=jcalc(pitch,q) calculates the joint transform and motion subspace
    % matrices for a revolute (pitch==0), prismatic (pitch==inf) or helical
    % (pitch==any other value) joint.  For revolute and helical joints, q is
    % the joint angle.  For prismatic joints, q is the linear displacement.

    if pitch == 0				% revolute joint
        Xj = Twist(SE3.Rz(q));
        S = [0;0;0;0;0;1];
    elseif pitch == inf			% prismatic joint
        Xj = Twist(SE3([0 0 q]));
        S = [0;0;1;0;0;0];
    else					% helical joint
        Xj = Twist(SE3.Rz(q) * SE3([0 0 q*pitch]));
        S = [0;0;pitch0;0;1;];
    end
end

author: petercorke

repo: spatialmath-matlab


text:
function f = isint(m)

% function f = isint(m)
%
% <m> is a matrix
%
% return a logical matrix the same size as <m>.
% an element is 1 iff it is a float and finite and exactly equal to an integer.
% specifically:
%   f = isfloat(m) & isfinite(m) & m==round(m);
%
% example:
% isequal(isint([1 1.5 NaN Inf]),[1 0 0 0])

% do it
f = isfloat(m) & isfinite(m) & m==round(m);

author: cvnlab

repo: GLMsingle

save_path: github-repos/MATLAB/cvnlab-GLMsingle


text:
function rv = radvl (pos, vel, velobs, star, dist)

% this function predicts the radial velocity of the observed
% object as it would be measured by spectroscopic means.  radial
% velocity is here defined as the radial velocity measure (z)
% times the speed of light.  for a solar system body, it applies
% to a fictitious emitter at the center of the observed object,
% assumed massless (no gravitational red shift), and does not
% in general apply to reflected light.  for stars, it includes
% all effects, such as gravitational red shift, contained
% in the catalog barycentric radial velocity measure, a scalar
% derived from spectroscopy.  nearby stars with a known kinematic
% velocity vector (obtained independently of spectroscopy) can be
% treated like solar system objects.  see lindegren & dravins
% (2003), astronomy & astrophysics 401, 1185-1201.
%
%      pos    = geometric position vector of object with respect to
%               observer, corrected for light-time, in au (in)
%      vel    = velocity vector of object with respect to solar
%               system barycenter, components in au/day (in)
%      velobs = velocity vector of observer with respect to solar
%               system barycenter, components in au/day (in)
%      star   = 3-element array of catalog data for a star, to be
%               non-zero if observed object is a star for which the
%               catalog radial velocity is consistent with
%               the iau definition of barycentric radial velocity
%               measure (otherwise all elements should be set to
%               0.d0 exactly) (in)
%               star(1) = catalog ra in hours
%               star(2) = catalog dec in degrees
%               star(3) = z*c, the catalog barycentric radial
%                         velocity measure times the speed of light,
%                         in kilometers/second
%               all three data elements must apply to the same
%               epoch (usually j2000.0 = jd 2451545.0 tt)
%      dist   = 3-element array of distances in au (in)
%               dist(1) = distance of observer from the geocenter
%               dist(2) = distance of observer from the sun
%               dist(3) = distance of object from the sun
%      rv     = the observed radial velocity measure times
%               the speed of light, in kilometers/second (out)

% note 1:  all the input arguments are bcrs quantities, expressed
% with respect to the icrs axes.  vel and velobs are kinematic
% velocities -- derived from geometry or dynamics, not spectroscopy.

% note 2:  if any element of array star is non-zero, the algorithm
% used will be consistent with the iau definition of stellar
% radial velocity, specifically, the barycentric radial velocity
% measure, which is derived from spectroscopy.  in that case,
% the vector vel can be very approximate -- or, for distant stars
% or galaxies, zero -- as it will be used only for a small geometric
% correction that is proportional to proper motion.

% note 3:  any of the distances in array dist can be set to zero
% (0.d0) if the corresponding general relativistic gravitational
% potential term is not to be evaluated.  these terms
% generally are important only at the meter/second level.  if
% the first two distances are both zero, an average value
% will be used for the relativistic term for the observer,
% appropriate for an observer on the surface of the earth.  the
% third distance, if given, is used only for solar system objects.

% ported from NOVAS 3.1

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

radcon = pi / 180.0d0;

% get au, length of astronomical unit in meters

au = astcon ('au', 1.0d0);

% get c, the speed of light in meters/second

c = astcon ('c', 1.0d0);

% get gs, heliocentric gravitational constant

gs = astcon ('gs', 1.0d0);

% get ge, geocentric gravitational constant

ge = astcon ('ge', 1.0d0);

% (gs and ge are in meters^3/second^2)

c2 = c^2;

toms = au / 86400.0d0;

toms2 = toms^2;

rv = 0.0d0;

% compute length of position vector = distance to object in au

posmag = sqrt(pos(1)^2 + pos(2)^2 + pos(3)^2);

if (posmag < 1.0d-8)

    return

end

% determine how object is to be processed

dostar = star(1) ~= 0.0d0 || star(2) ~= 0.0d0 || star(3) ~= 0.0d0;

% compute unit vector toward object

for j = 1:3

    uk(j) = pos(j) / posmag;

end

% compute velocity-squared factors

v2 = (vel(1)^2 + vel(2)^2 + vel(3)^2) * toms2;

vo2 = (velobs(1)^2 + velobs(2)^2 + velobs(3)^2) * toms2;

% compute geopotential at observer, unless observer is geocentric

r = dist(1) * au;

phigeo = 0.0d0;

if (r > 1.0d6)

    phigeo = ge / r;

end

% compute solar potential at observer

r = dist(2) * au;

phisun = 0.0d0;

if (r > 1.0d8)

    phisun = gs / r;

end

% compute relativistic potential and velocity factor for observer

if (dist(1) ~= 0.0d0 || dist(2) ~= 0.0d0)

    % lindegren & dravins eq. (41), second factor in parentheses

    rel = 1.0d0 - (phigeo + phisun) / c2 - 0.5d0 * vo2 / c2;

else

    % lindegren & dravins eq. (42), inverse

    rel = 1.0d0 - 1.550d-8;

end

if (dostar == 1)

    % for stars, update barycentric radial velocity measure for change in view angle

    ra = star(1) * 15.0d0 * radcon;

    dc = star(2) * radcon;

    du(1) = uk(1) - (cos (dc) * cos (ra));

    du(2) = uk(2) - (cos (dc) * sin (ra));

    du(3) = uk(3) - (sin (dc));

    zc = star(3) * 1.0d3 + (vel(1) * du(1) + vel(2) * du(2) + vel(3) * du(3)) * toms;

    % compute observed radial velocity measure of a star (inverse of
    % lindegren & dravins eq. (41))

    zb1 = 1.d0 + zc / c;

    kvobs = (uk(1) * velobs(1) + uk(2) * velobs(2) + uk(3) * velobs(3)) * toms;

    zobs1 = zb1 * rel / (1.0d0 + kvobs / c);

else

    % compute solar potential at object, if within solar system

    r = dist(3) * au;

    phisun = 0.0d0;

    if (r > 1.0d8 && r < 1.0d16)

        phisun = gs / r;

    end

    % compute observed radial velocity measure of a planet or other
    % object -- including a nearby star -- where kinematic
    % barycentric velocity vector is known and gravitational
    % red shift is negligible (lindegren & dravins eq. (40),
    % applied as per s. klioner private communication (2006))

    kv = (uk(1) * vel(1) + uk(2) * vel(2) + uk(3) * vel(3)) * toms;

    zb1 = (1.0d0 + kv / c) / (10.0d0 - phisun / c2 - 0.5d0 * v2 / c2);

    kvobs = (uk(1) * velobs(1) + uk(2) * velobs(2) + uk(3) * velobs(3)) * toms;

    zobs1 = zb1 * rel / (1.0d0 + kvobs / c);

end

% convert observed radial velocity measure to kilometers/second

rv = (zobs1 - 1.0d0) * c / 1000.0d0;

author: goGPS-Project


text:
function Dict = edl_buildiidxdict(S, idxfn)
%EDL_BUILDIIDXDICT Builds a dictionary using internal index
%
% $ Syntax $
%   - Dict = edl_buildiidxdict(S, idxfn)
%
% $ Arguments $
%   - S:        a struct array on which the dictionary is built
%   - idxfn:    the field name of the internal index
%               (default = 'internal_index')
%   - Dict:     the built dictionary
%               a cell array, indexed by internal_index
%
% $ Remarks $
%   - The length of the cell array (dict) will be equal to the maximum
%     internal index in S.
%   - In S, different elements should have different internal indices.
%   - For the index corresponding to no element, the corresponding cells
%     are empty.
%   - The dictionary will not automatically keep track of the change of
%     S. If S is changed, the dictonary should be rebuilt in order to
%     reflect the updates.
%
% $ History $
%   - Created by Dahua Lin, on Aug 15, 2006
%

%% parse and verify input arguments

if ~isstruct(S)
    error('sltoolbox:invalidarg', ...
        'The S should be a struct');
end

if nargin < 2 || isempty(idxfn)
    idxfn = 'internal_index';
end

if ~isfield(S, idxfn)
    error('sltoolbox:invalidarg', ...
        'The S have no specified internal index field');
end

%% Build dictionary

% gather the internal indices
n = numel(S);
inds = zeros(n, 1);
for i = 1 : n
    idx = S(i).(idxfn);
    if ischar(idx)
        idx = str2double(idx);
    end
    inds(i) = idx;
end

% construct dictionary
m = max(inds);

Dict = cell(m, 1);

for i = 1 : n
    idx = inds(i);
    if ~isempty(Dict{idx})
        error('sltoolbox:rterror', ...
            'Repeated index %d is encountered', idx);
    end
    Dict{idx} = S(i);
end

author: lmthang


text:
    function xsol = fzero_data(x,y,y0)
        N=length(x);
        vectsign=sign(y-y0);
        pos=zeros(1,N);
        for i=1:(N-1)
            if vectsign(i)~=vectsign(i+1)
                pos(i)=1;
                pos(i+1)=1;
            end
        end

        indices=find(pos);
        Nmax=2*(length(indices)-1);
        vectint=zeros(1,Nmax);
        vectint(1)=indices(1);
        vectint(Nmax)=indices(length(indices));
        for j=2:2:(Nmax-2)
            vectint(j)=indices(j/2+1);
            vectint(j+1)=indices(j/2+1);
        end
        Nmaxsol=Nmax/2;
        indsol=zeros(1,2*Nmaxsol);
        for k=1:Nmaxsol
            if (y(vectint(2*k-1))-y0>=0 && y(vectint(2*k))-y0<=0)||(y(vectint(2*k-1))-y0<=0 && y(vectint(2*k))-y0>=0)
                indsol(2*k-1)=1;
                indsol(2*k)=1;
            else
                indsol(2*k-1)=0;
                indsol(2*k)=0;
            end
        end
        vectpos=vectint(find(indsol));
        Nsol=length(vectpos)/2;
        xsol=zeros(1,Nsol);
        for n=1:Nsol
            xsol(n)=interp1(y(vectpos(2*n-1:2*n))-y0,x(vectpos(2*n-1:2*n)),0);
        end

        clear vectsign pos indices Nmax vectint Nmaxsol indsol vectpos Nsol

    end

author: Sable

repo: mcbench-benchmarks


text:
function [r, c] = kurtosis_ijkl_to_rc(i, j, k, l)
%KURTOSIS_IJKL_TO_RC converts indices in a kurtosis tensor from 3x3x3x3
%format to 6x6 format

% Author: Darryl McClymont <darryl.mcclymont@gmail.com>
% Copyright © 2014-2016 University of Oxford
% Version: 0.1.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

switch i
    case 1
        r1 = [1 2 3];
    case 2
        r1 = [2 4 5];
    case 3
        r1 = [3 5 6];
end
switch j
    case 1
        r2 = [1 2 3];
    case 2
        r2 = [2 4 5];
    case 3
        r2 = [3 5 6];
end

if i == j
    switch i
        case 1
            r =1;
        case 2
            r = 4;
        case 3
            r = 6;
    end
else
    r = intersect(r1, r2);
end

switch k
    case 1
        c1 = [1 2 3];
    case 2
        c1 = [2 4 5];
    case 3
        c1 = [3 5 6];
end
switch l
    case 1
        c2 = [1 2 3];
    case 2
        c2 = [2 4 5];
    case 3
        c2 = [3 5 6];
end

if k == l
    switch k
        case 1
            c =1;
        case 2
            c = 4;
        case 3
            c = 6;
    end
else
    c = intersect(c1, c2);
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% MSH_EVALUATE_ELEMENT_LIST: evaluate the parameterization in a given list of elements.
%
%     msh_elems = msh_evaluate_element_list (msh, elements)
%
% INPUTS:
%
%    msh:          mesh object (see msh_cartesian)
%    element_list: numbering of the elements where the evaluations are performed.
%
% OUTPUT:
%
%     msh_elems: structure containing the quadrature rule in the given elements of the physical domain, which contains the following fields
%
%     FIELD_NAME         (SIZE)                  DESCRIPTION
%     npatch             (scalar)                number of patches
%     ndim               (scalar)                dimension of the parametric space
%     rdim               (scalar)                dimension of the physical space
%     nel                (scalar)                number of elements in the list
%     elem_list          (1 x nel)               numbering of the elements in the list
%     nqn                (scalar)                number of quadrature points per element (must be the same for every patch)
%     nqn_dir            (1 x ndim)              number of quadrature points in each direction (must be the same for every patch)
%     nel_per_patch      (1 x npatch)            number of selected elements on each patch
%     elem_list_of_patch (1 x npatch cell-array) selected elements on the patch, with local numbering
%     nel_dir_of_patch   (1 x npatch cell-array) the total number of elements in each direction, for each patch
%     quad_weights, geo_map, geo_map_jac, deo_map_der2, jacdet, element_size (see msh_evaluate_col for details)
%
% The function only works if the number of quadrature points is the same for all the patches and all directions.
%
% Copyright (C) 2015 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function msh_col = msh_evaluate_element_list (msh, elem_list, varargin)

  elem_list = elem_list(:)';

  msh_col.npatch = msh.npatch;
  msh_col.ndim = msh.ndim;
  msh_col.rdim = msh.rdim;

  msh_col.nel  = numel (elem_list);
  msh_col.elem_list = elem_list;

  if (isempty (elem_list)), return, end

  fields = {'quad_weights', 'geo_map', 'geo_map_jac', 'geo_map_der2', 'jacdet', 'element_size'};
  cat_position = [2 3 4 5 2 2];
  for ii = 1:numel (fields)
    msh_col.(fields{ii}) = [];
  end

  Nelem = cumsum ([0, msh.nel_per_patch]);

  indices = cell (1, msh.npatch);
  for iptc = 1:msh.npatch
    [~,indices{iptc},~] = intersect ((Nelem(iptc)+1):Nelem(iptc+1), elem_list);
  end
  active_patches = find (~cellfun (@isempty, indices));

  for iptc = active_patches
    msh_patch = msh_evaluate_element_list (msh.msh_patch{iptc}, indices{iptc});

    if (msh_patch.nqn_dir ~= msh.msh_patch{active_patches(1)}.nqn_dir)
      error ('msh_evaluate_element_list: for multipatch geometries, the number of quadrature points should be the same on each patch')
    end

    for ii = 1:numel (fields)
      if (isfield (msh_patch, fields{ii}))
        msh_col.(fields{ii}) = cat (cat_position(ii), msh_col.(fields{ii}), msh_patch.(fields{ii}));
      end
    end
  end

  msh_col.nel_per_patch = cellfun (@numel, indices);
  msh_col.elem_list_of_patch = indices;
  msh_col.nel_dir_of_patch = cell (1, msh.npatch);
  for iptc = 1:msh.npatch
    msh_col.nel_dir_of_patch{iptc} = msh.msh_patch{iptc}.nel_dir;
  end
  msh_col.nqn = msh.msh_patch{active_patches(1)}.nqn;
  msh_col.nqn_dir = msh.msh_patch{active_patches(1)}.nqn_dir;

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
% rf_scaleGrad.m
% Jamie Near, McGill University 2020.
%
% USAGE:
% RF_out=rf_scaleGrad(RF_in,scale);
%
% DESCRIPTION:
% Scale the gradient amplitude of a gradient modulated pulse by a fixed factor.
%
% INPUTS:
% RF_in     = Input RF pulse definition structure.  Must be a gradient
%             modulated pulse.
% scale     = Scaling factor by which you would like to multiply the
%             existing gradient waveform.
%
% OUTPUTS:
% RF_out    = Output rf waveform following the scaling of gradient
%             waveform.

function RF_out=rf_scaleGrad(RF_in,scale);

if ~isstruct(RF_in)
    error('ERROR:  the input RF pulse must be in structure format.  Try using rf_readwaveform to convert it!  Aborting.  ');
end

if ~RF_in.isGM
    error('ERROR:  the input RF pulse must be a gradient modulated pulse.  ABORTING!');
end

newWaveform=RF_in.waveform;
newWaveform(:,4)=newWaveform(:,4)*scale;

RF_out=RF_in;
RF_out.waveform=newWaveform;
RF_out.tthk=RF_in.tthk/scale;

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function [BC,side,r] = voxel_grid(V,side,varargin)
  % VOXEL_GRID Prepare a voxel grid around a set of points V
  %
  % [BC,side,r] = voxel_grid(V,side,varargin)
  %
  % Inputs:
  %   V  #V by 3 list of input point positions
  %   side  either:
  %     scalar specifying how many steps along the x-coordinate of V's bounding
  %     box
  %          or
  %     3 list specifying how many steps in x, y, and z coordinates of V's
  %       bounding box
  %   Optional:
  %     'Pad' followed by the number of "extra cells" to add on all six sides
  %       of the grid (pad_count).
  % Outputs:
  %   BC  prod(side+2*pad_count) by 3 list of cell centers
  %   side  number of cells on each side: side+2*pad_count
  %   r  size of step in each direciton
  %
  % See also: voxelize

  pad_count = 0;
  % default values
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Pad'}, ...
    {'pad_count'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  dim = size(V,2);

  assert(all(side>(pad_count*2+1)),'side should be > 2*pad_count+1');
  side = side-pad_count*2;
  switch numel(side)
  case 3
    side(1) = side(1);
    side(2) = side(2);
    side(3) = side(3);
    NV = min(V);
    XV = max(V);
    r = [XV-NV]./([side(1) side(2) side(3)] - 1);
  case 1
    % Enclose bounding box in regular mesh
    side(1) = side(1);
    NV = min(V);
    XV = max(V);
    for d = 2:dim
      side(d) = ceil(side(1) * (XV(d)-NV(d))/(XV(1)-NV(1)));
    end
    r = max((XV-NV)./(side-1));
    % recenter
    old_cen = 0.5*(XV + NV);
    XV = NV + r*(side-1);
    cen = 0.5*(XV + NV);
    XV = XV + old_cen - cen;
    NV = NV + old_cen - cen;
  otherwise
    error('side should be scalar or triplet');
  end
  assert(all(side==ceil(side)),'All side values should be integer');

  side = side+pad_count*2;
  old_cen = 0.5*(XV + NV);
  XV = NV+r.*(side-1);
  cen = 0.5*(XV + NV);
  XV = XV + old_cen - cen;
  NV = NV + old_cen - cen;

  r = (XV-NV)./(side-1);

  switch dim
  case 3
    [X,Y,Z] = meshgrid( ...
      NV(1)+linspace(0,1,side(1))*(XV(1)-NV(1)), ...
      NV(2)+linspace(0,1,side(2))*(XV(2)-NV(2)), ...
      NV(3)+linspace(0,1,side(3))*(XV(3)-NV(3)));
    % barycenters of cells
    BC = [X(:) Y(:) Z(:)];
  case 2
    [X,Y] = meshgrid( ...
      NV(1)+linspace(0,1,side(1))*(XV(1)-NV(1)), ...
      NV(2)+linspace(0,1,side(2))*(XV(2)-NV(2)));
    % barycenters of cells
    BC = [X(:) Y(:)];
  end

end

author: alecjacobson

repo: gptoolbox


text:
function h=plotgauss2d(mu, Sigma, plot_cross)
% PLOTGAUSS2D Plot a 2D Gaussian as an ellipse with optional cross hairs
% h=plotgauss2(mu, Sigma)
%
% h=plotgauss2(mu, Sigma, 1) also plots the major and minor axes
%
% Example
% clf; S=[2 1; 1 2]; plotgauss2d([0;0], S, 1); axis equal

h = plotcov2(mu, Sigma);
return;

%%%%%%%%%%%%%%%%%%%%%%%%
function old

if nargin < 3, plot_cross = 0; end
[V,D]=eig(Sigma);
lam1 = D(1,1);
lam2 = D(2,2);
v1 = V(:,1);
v2 = V(:,2);
%assert(approxeq(v1' * v2, 0))
if v1(1)==0
  theta = 0; % horizontal
else
  theta = atan(v1(2)/v1(1));
end
a = sqrt(lam1);
b = sqrt(lam2);
h=plot_ellipse(mu(1), mu(2), theta, a,b);

if plot_cross
  mu = mu(:);
  held = ishold;
  hold on
  minor1 = mu-a*v1; minor2 = mu+a*v1;
  hminor = line([minor1(1) minor2(1)], [minor1(2) minor2(2)]);

  major1 = mu-b*v2; major2 = mu+b*v2;
  hmajor = line([major1(1) major2(1)], [major1(2) major2(2)]);
  %set(hmajor,'color','r')
  if ~held
    hold off
  end
end


% path planner
%
% Modified:
%   - 4/06/2010 - RWB
%
% output is a vector containing P.num_waypoints waypoints
%
% input is the map of the environment and desired goal location
function out = path_planner1(in,P,map)
load  'uavW.mat'
i1 = uavW;
clear uavW;

 NN = 0;
  pn        = in(1+NN);
  pe        = in(2+NN);
  h         = in(3+NN);
  % Va      = in(4+NN);
  % alpha   = in(5+NN);
  % beta    = in(6+NN);
  % phi     = in(7+NN);
  % theta   = in(8+NN);
  chi     = in(9+NN);
  % p       = in(10+NN);
  % q       = in(11+NN);
  % r       = in(12+NN);
  % Vg      = in(13+NN);
  % wn      = in(14+NN);
  % we      = in(15+NN);
  % psi     = in(16+NN);
  flag_new_waypoints =  in(17+NN);
  NN = NN + 17;
  t         =  in(1+NN);

  persistent num_waypoints
  persistent wpp

  if (t==0) || flag_new_waypoints
    % format for each point is [pn, pe, pd, chi, Va^d] where the position
    % of the waypoint is (pn, pe, pd), the desired course at the waypoint
    % is chi, and the desired airspeed between waypoints is Va
    % if chi!=-9999, then Dubins paths will be used between waypoints.
    switch 1
          case 1

        [num_waypoints , wpp] = getWpp(P,i1);

          case 2  % Dubins
            num_waypoints = 4;
            wpp = [...
                    0, 0, -100, 0, P.Va0;...
                    300, 0, -100, 45*pi/180, P.Va0;...
                    0, 300, -100, 45*pi/180, P.Va0;...
                    300, 300, -100, -135*pi/180, P.Va0;...
                  ];
          case 3 % path through city using Dubin's paths
               % current configuration
              wpp_start = [pn, pe, -h, chi, P.Va0];
              % desired end waypoint
              if norm([pn; pe; -h]-[map.width; map.width; -h])<map.width/2,
                  wpp_end = [0, 0, -h, 0, P.Va0];
              else
                  wpp_end = [map.width, map.width, P.pd0, 0, P.Va0];
              end
              waypoints = planRRTDubins(wpp_start, wpp_end, P.R_min, map);
              num_waypoints = size(waypoints,1);
              wpp = [];
              for i=1:num_waypoints
                  wpp = [...
                            wpp;...
                            waypoints(i,1), waypoints(i,2), waypoints(i,3), waypoints(i,4), P.Va0;...
                        ];
              end
          case 4  % path through city using straight-line RRT
               % current configuration
              wpp_start = [pn, pe, -h, chi, P.Va0];
              % desired end waypoint
              if norm([pn; pe; -h]-[map.width; map.width; -h])<map.width/2,
                  wpp_end = [0, 0, -h, 0, P.Va0];
              else
                  wpp_end = [map.width, map.width, P.pd0, 0, P.Va0];
              end
              waypoints = planRRT(wpp_start, wpp_end, map);
              num_waypoints = size(waypoints,1);
              wpp = [];
              for i=1:num_waypoints
                  wpp = [...
                            wpp;...
                            waypoints(i,1), waypoints(i,2), waypoints(i,3), waypoints(i,4), P.Va0;...
                        ];
              end
          case 5  % cover path through city using Dubin's paths
               % current configuration
              wpp_start = [pn, pe, -h, chi, P.Va0];
              waypoints = planCoverRRTDubins(wpp_start, P.R_min, map);
              num_waypoints = size(waypoints,1);
              wpp = [];
              for i=1:num_waypoints
                  wpp = [...
                            wpp;...
                            waypoints(i,1), waypoints(i,2), waypoints(i,3), waypoints(i,4), P.Va0;...
                        ];
              end
    end
    for i=num_waypoints+1:P.size_waypoint_array
          wpp = [...
                    wpp;...
                    -9999, -9999, -9999, -9999, -9999;...
                ];
    end

  end

  out = [num_waypoints; reshape(wpp', 5*P.size_waypoint_array, 1)];

end


text:
function [ a, q, lambda, seed ] = r8symm_test ( n, lambda_mean, lambda_dev, ...
  seed )

%*****************************************************************************80
%
%% R8SYMM_TEST determines a symmetric matrix with a certain eigenstructure.
%
%  Discussion:
%
%    An R8SYMM is a real symmetric matrix stored using full storage, and
%    using R8 arithmetic.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real LAMBDA_MEAN, the mean value of the normal
%    distribution from which the eigenvalues will be chosen.
%
%    Input, real LAMBDA_DEV, the standard deviation of the normal
%    distribution from which the eigenvalues will be chosen.
%
%    Input/output, integer SEED, a seed for the random
%    number generator.
%
%    Output, real A(N,N), the test matrix.
%
%    Output, real Q(N,N), the eigenvector matrix.
%
%    Output, real LAMBDA(N), the eigenvalue vector.
%

%
%  Choose the eigenvalues LAMBDA.
%
  [ lambda, seed ] = r8vec_normal ( n, lambda_mean, lambda_dev, seed );
%
%  Get a random orthogonal matrix Q.
%
  [ q, seed ] = r8mat_orth_uniform ( n, seed );
%
%  Set A = Q * Lambda*I * Q'.
%
  a = q * diag ( lambda ) * q';

  return
end

author: johannesgerer


text:
function OUT = cell2num(datacell)
% =======================================================================
% Convert a cell array (of numbers) into a double array. If the cell arrqay
% does not include a number returns NaN
% =======================================================================
% OUT = cell2num(datacell)
% -----------------------------------------------------------------------
% INPUTS
%	- datacell: a (m x n) cell array with numbers
% -----------------------------------------------------------------------
% OUTPUT
%	- OUT: a (m x n) matrix with numbers
% =========================================================================
% Ambrogio Cesa Bianchi, March 2015
% ambrogiocesabianchi@gmail.com
%-------------------------------------------------------------------------

if ~iscell(datacell)
    error('Input has to be a cell array')
end

% Initialize
OUT = zeros(size(datacell));

% Convert
for c=1:size(datacell,2)
    for r=1:size(datacell,1)
        if isnumeric(datacell{r,c})
            OUT(r,c)=datacell{r,c};
        else
            OUT(r,c)=NaN;
        end
    end
end

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
%% runMonitorFEBio
% Below is a demonstration of the features of the |runMonitorFEBio| function

%%
clear; close all; clc;

%% Syntax
% |[runFlag]=runMonitorFEBio(FEBioRunStruct);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function sR=icassoCluster(sR,varargin)
%function sR=icassoCluster(sR,['identifier1',value1,'identifier2',value2,...])
%
%PURPOSE
%
%To cluster the ICA estimates and to compute a relative clustering validity
%index (R-index)
%
%EXAMPLE OF BASIC USAGE
%
% sR=icassoCluster(sR);
%
%where sR is an Icasso result structure. This applies hierarchical
%clustering using group-average linkage agglomeration strategy and
%stores the results back into workspace variable sR.
%
%INPUT
%
% sR (struct) Icasso result data structure
%
%Optional input arguments are given as argument identifier - value
%pairs: 'identifier1', value1, 'identifier2', value2,...
%(case insensitive)
%
% 'simfcn' (string) 'abscorr' (default) | (matrix)
%   Indicates how to compute similarities S between estimates i,j
%   'abscorr' S(i,j) is the absolute value of the linear
%     correlation coefficient
%   (matrix) explicitly given similarity matrix S
%     (elements should be 0...1)
% 'strategy' (string) 'AL' (default) | 'SL' | 'CL'
%  Sets the clustering strategy:
%  'AL' hierarchical group average linkage
%  'SL' hierarchical single linkage (nearest neighbor)
%  'CL' hierarchical complete linkage (furthest neighbor)
% 'L' (integer) | (string) 'rdim' (default)
%   computes a relative clustering validity index for 2...L
%   clusters. Default string 'rdim' sets L same as the (reduced)
%   data dimension
%'s2d' (string) (default is 'sim2dis')
%  the name of function that is used to make the transformation
%  from similarities between IC components to dissimilarities
%  D. Default 'sim2dis' makes simply D=1-S;
%
%OUTPUT
%
% sR (struct) updated Icasso result data structure
%
%The function updates the fields sR.cluster.* only.
%
%DETAILS
%1. The function computes similarities between the
%estimates. See Note 1
%2. stores the similarity matrix S into field sR.cluster.similarity
%and the method into field sR.cluster.simfcn
%2. transforms the similarities into dissimilarities (distances)
%D. See Note 2.
%3. applies the selected clustering strategy on dissimilarities. D
%The results is a partition matrix P (of size MxM). See explanation
%in function hcluster The function stores P in field
%sR.cluster.partition (Stores also other outputs of function
%hcluster into sR.cluster.dendrogram).
%4. computes a relative clustering validity index (see function
%rindex) for dissimilarities D and partitions P from 2 to L
%clusters. If not explicitly given (or string 'rdim' is given, L
%will be the (reduced data) dimension. Stores the result in field
%sR.cluster.index.R.
%
%NOTE 1
%
%By default, the similarities between estimates are computed as
%linear correlation coefficients (using the demixing matrix rows and
%dewhitening matrix; see function corrw). You can also explicitly
%give any MxM similarity matrix as input. The similarities
%should be between 0...1. (M is the total number of estimates.)
%
%NOTE 2
%The clustering and clustering validity indices are computed for
%dissimilarities. Therefore, the similarities must be transformed
%into dissimilarities. By default, icassoCluster uses subfunction
%sim2dis to do this; sim2dis performs simply D=1-S. If you wish to
%make the transformation otherwise, you can set a different
%function by using input identifier - valuepair
%'s2d','myfunctionname'. icassoCluster then computes
%dissimilarities by calling
%   D=feval('myfunctionname',sR.cluster.similarity);
%A possible function is sqrtsim2dis
%
%SEE ALSO
% hcluster
% som_linkage
% sqrtsim2dis
% rindex
% icassoDendrogram

%COPYRIGHT NOTICE
%This function is a part of Icasso software library
%Copyright (C) 2003-2005 Johan Himberg
%
%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 any later version.
%
%This program is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.
%
%You should have received a copy of the GNU General Public License
%along with this program; if not, write to the Free Software
%Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

% ver 1.21 070305 johan

% Init some variables

% total number of estimates
M=icassoGet(sR,'M');

%reduced data dimension
rdim=icassoGet(sR,'rdim');

% Set default parameters
default={'simfcn','abscorr','s2d','sim2dis','strategy','AL','L','rdim'};

%% Check optional arguments and add defaults
clusterparameters=processvarargin(varargin,default);
num_of_args=length(clusterparameters);

%% check arguments
for i=1:2:num_of_args;
  switch lower(clusterparameters{i})

   case 'simfcn'
    simfcn=clusterparameters{i+1};

    % Explicit similarity matrix?
    if isnumeric(simfcn),
      if size(simfcn,1)==M & size(simfcn,2)==M,
        sR.cluster.similarity=simfcn;
        sR.cluster.simfcn='<similarities given explicitly>';
      else
        error('Explicitly given similarity matrix has wrong size!');
      end
    else
      % should be a string
      switch lower(simfcn)
       case 'abscorr'
        % ok
        sR.cluster.simfcn=lower(simfcn);
       otherwise
        error('''simfcn'' must be string ''abscorr'' or an MxM similarity matrix');
      end
    end
   case 's2d'
    s2dfcn=lower(clusterparameters{i+1});
    if ~ischar(s2dfcn),
      error('''s2d'' must be a string (name of a function)');
    end
    sR.cluster.s2d=s2dfcn;
   case 'l'
    L=clusterparameters{i+1};
    if isnumeric(L),
      % The user has specified max number for clusters

      % Check L
      if fix(L)~=L,
        error('''L'' must be an integer.');
      elseif L<2,
        error('''L'' must be at least 2.');
      elseif L>M,
        error('''L'' cannot be more than the number of estimates.');
      end
    else
      if ~strcmp(lower(L),'rdim'),
        error('''L'' expects an integer value or ''rdim''.');
      end
      % set (reduced) data dimension
      L=icassoGet(sR,'rdim');
    end

    if L>100,
      warning(['R-index requested for more that 100 clusters: this can' ...
               ' be heavy...']);
    end

   case 'strategy'
    strategy=clusterparameters{i+1};
    if ~ischar(strategy),
      error('''strategy'' must be a string');
    end

    % we are case insensitive
    strategy=upper(strategy);
    sR.cluster.strategy=strategy;

    switch sR.cluster.strategy
     case {'AL','CL','SL'}
      ; % hierarchical clustering
     otherwise
      error(['Strategy ' strategy ' not implemented.']);
    end
   otherwise
    error(['Indentifier ' clusterparameters{i} ' not recognized.']);
  end
end

%%%% Compute similarities %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch lower(sR.cluster.simfcn)
 case '<similarities given explicitly>'
  % already handled
 case 'abscorr'
  sR.cluster.similarity=abs(corrw(icassoGet(sR,'W'),icassoGet(sR,'dewhitemat')));
  %just to make sure
  sR.cluster.similarity(sR.cluster.similarity>1)=1;
  sR.cluster.similarity(sR.cluster.similarity<0)=0;
end

%%%%% Convert to dissimilarities using .s2d

D=feval(sR.cluster.s2d, sR.cluster.similarity);

%%%% Make partition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[sR.cluster.partition,sR.cluster.dendrogram.Z,sR.cluster.dendrogram.order]=...
    hcluster(D,sR.cluster.strategy);

%%%%% Compute cluster validity %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% init R
sR.cluster.index.R=ones(M,1)*NaN;
% compute

sR.cluster.index.R(1:L,1)=rindex(D,sR.cluster.partition(1:L,:));