text:
classdef prtDataSetBigClass < prtDataSetBig & prtDataInterfaceCategoricalTargetsBig
    % prtDataSetBigClass is a class for prtDataSetBig that are for
    % classification. It is currently a placeholder for future
    % classification specific helper methods.

    methods (Hidden)

        function self = summaryClear(self)
            % self = summaryClear(self)
            self.summaryCache = [];
            self.targetCacheInitialized = false;
        end

        function self = summaryBuild(self,force)
            % self = summaryBuild(self,force)
            %
            if nargin < 2
                force = false;
            end
            if force || isempty(self.summaryCache)
                mrSummary = prtMapReduceSummarizeDataSetClass;
                self.summaryCache = mrSummary.run(self);
                self.targetCache = self.summaryCache.targetCache;
                self.classNames = self.targetCache.classNames;
                self.targetCacheInitialized = true;
            end
        end
    end

    methods

        function self = prtDataSetBigClass(varargin)
            self = prtUtilAssignStringValuePairs(self, varargin{:});
        end

        function n = getNumFeatures(self)
            summary = self.summarize;
            n = summary.nFeatures;
        end

        function n = getNumObservations(self)
            summary = self.summarize;
            n = summary.nObservations;
        end

        function n = getNumTargetDimensions(self)
            summary = self.summarize;
            n = summary.nTargetDimensions;
        end

        function summary = summarize(self)
            summary = self.summaryCache;
            if isempty(summary)
                summary = getDefaultSummaryStruct(self);
            end
        end

        function varargout = plot(self)
            plotHandles = plot(self.getRandomBlock());
            varargout = {};
            if nargout
                varargout = {plotHandles};
            end
        end
    end
    methods (Hidden)
        function summary = getDefaultSummaryStruct(self)

            warning('prtDataSetBigClass:defaultSummary','prtDataSetBig* classes provide default values (nan) for many fields, since these take time to calculate; use ds = ds.summaryBuild to build the summary');
            warning off prtDataSetBigClass:defaultSummary

            summary = struct('upperBounds',nan,'lowerBounds',nan,'nFeatures',nan,...
                'nTargetDimensions',nan,'uniqueClasses',nan,'nClasses',nan,...
                'isMary',nan,'nObservations',nan,'nBlocks',nan,'targetCache',nan,...
                'isSummaryValid',false);
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function nsm3d_test01 ( prefix )

%*****************************************************************************80
%
%% NSM3D_TEST01 examines a mesh from the NAVIER_STOKES_MESH3D data.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 May 2013
%
%  Author:
%
%    John Burkardt
%
  mesh_filename = sprintf ( '%s.mat', prefix );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'NSM3D_TEST01:\n' );
  fprintf ( 1, '  "%s"\n', mesh_filename );

  [ nodeco, elnode, bdynde ] = mesh3d_extract ( mesh_filename );
%
%  NODECO has an unused 4th dimension.  Suppress it.
%
  nodeco = nodeco(:,1:3);
  [ node_num, node_dim ] = size ( nodeco );
  [ element_num, element_order ] = size ( elnode );
  [ bdynde_num, ~ ] = size ( bdynde );

  fprintf ( 1, '  NODECO = %d nodes x %d spatial dimensions,\n', ...
    node_num, node_dim );
  fprintf ( 1, '  ELNODE = %d elements x %d nodes per element,\n', ...
    element_num, element_order );
  fprintf ( 1, '  BDYNDE = %d boundary nodes\n', bdynde_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Dim       Min          Max\n' );
  fprintf ( 1, '\n' );
  for j = 1 : node_dim
    xmin = min ( nodeco(:,j) );
    xmax = max ( nodeco(:,j) );
    fprintf ( 1, '  %d  %12f  %12f\n', j, xmin, xmax );
  end
%
%  Write the nodes to a file.
%
  node_filename = sprintf ( '%s_nodes.txt', prefix );
  nodeco = nodeco';
  r8mat_write ( node_filename, node_dim, node_num, nodeco );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Created the node file "%s".\n', node_filename );

  element_filename = sprintf ( '%s_elements.txt', prefix );
  elnode = elnode';
  i4mat_write ( element_filename, element_order, element_num, elnode );
  fprintf ( 1, '  Created the element file "%s".\n', element_filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function x = sq(x, varargin)

if numel(varargin)>0
    idims = varargin{1};
    xsiz = size(x);
    idims = idims(xsiz(idims)==1);
    xsiz(idims) = [];
    x = reshape(x, xsiz);
else
    x = squeeze(x);
end

author: cortex-lab


text:
function [sMap,D] = som_make(D, varargin)

%SOM_MAKE Create, initialize and train Self-Organizing Map.
%
% sMap = som_make(D, [[argID,] value, ...])
%
%  sMap = som_make(D);
%  sMap = som_make(D, 'munits', 20);
%  sMap = som_make(D, 'munits', 20, 'hexa', 'sheet');
%  sMap = som_make(D, 'msize', [4 6 7], 'lattice', 'rect');
%
%  Input and output arguments ([]'s are optional):
%   D        (matrix) training data, size dlen x dim
%            (struct) data struct
%   [argID,  (string) See below. The values which are unambiguous can
%    value]  (varies) be given without the preceeding argID.
%
%   sMap     (struct) map struct
%
% Here are the valid argument IDs and corresponding values. The values
% which are unambiguous (marked with '*') can be given without the
% preceeding argID.
%   'init'       *(string) initialization: 'randinit' or 'lininit' (default)
%   'algorithm'  *(string) training: 'seq' or 'batch' (default), 'imp' or 'sompak'
%   'munits'      (scalar) the preferred number of map units
%   'msize'       (vector) map grid size
%   'mapsize'    *(string) do you want a 'small', 'normal' or 'big' map
%                          Any explicit settings of munits or msize override this.
%   'lattice'    *(string) map lattice, 'hexa' or 'rect'
%   'shape'      *(string) map shape, 'sheet', 'cyl' or 'toroid'
%   'neigh'      *(string) neighborhood function, 'gaussian', 'cutgauss',
%                          'ep' or 'bubble'
%   'topol'      *(struct) topology struct
%   'som_topol','sTopol' = 'topol'
%   'mask'        (vector) BMU search mask, size dim x 1
%   'name'        (string) map name
%   'comp_names'  (string array / cellstr) component names, size dim x 1
%   'tracking'    (scalar) how much to report, default = 1
%   'training'    (string) 'short', 'default', 'long'
%                 (vector) size 1 x 2, first length of rough training in epochs,
%                          and then length of finetuning in epochs
%
%
% For more help, try 'type som_make' or check out online documentation.
% See also SOM_MAP_STRUCT, SOM_TOPOL_STRUCT, SOM_TRAIN_STRUCT,
%          SOM_RANDINIT, SOM_LININIT, SOM_SEQTRAIN, SOM_BATCHTRAIN.

%%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% som_make
%
% PURPOSE
%
% Creates, initializes and trains a SOM using default parameters.
%
% SYNTAX
%
%  sMap = som_make(D);
%  sMap = som_make(...,'argID',value,...);
%  sMap = som_make(...,value,...);
%
% DESCRIPTION
%
% Creates, initializes and trains a SOM with default parameters. Uses functions
% SOM_TOPOL_STRUCT, SOM_TRAIN_STRUCT, SOM_DATA_STRUCT and SOM_MAP_STRUCT to come
% up with the default values.
%
% First, the number of map units is determined. Unless they are
% explicitly defined, function SOM_TOPOL_STRUCT is used to determine this.
% It uses a heuristic formula of 'munits = 5*dlen^0.54321'. The 'mapsize'
% argument influences the final number of map units: a 'big' map has
% x4 the default number of map units and a 'small' map has x0.25 the
% default number of map units.
%
% After the number of map units has been determined, the map size is
% determined. Basically, the two biggest eigenvalues of the training
% data are calculated and the ratio between sidelengths of the map grid
% is set to this ratio. The actual sidelengths are then set so that
% their product is as close to the desired number of map units as
% possible.
%
% Then the SOM is initialized. First, linear initialization along two
% greatest eigenvectors is tried, but if this can't be done (the
% eigenvectors cannot be calculated), random initialization is used
% instead.  After initialization, the SOM is trained in two phases:
% first rough training and then fine-tuning. If the 'tracking'
% argument is greater than zero, the average quantization error and
% topographic error of the final map are calculated. If there are missing
% values (NaNs) in the data, it may be wise to use the Imputation SOM
% algorithm ('imp' agrument).
%
%
% REQUIRED INPUT ARGUMENTS
%
%  D           The data to use in the training.
%     (struct) A data struct. If a struct is given, '.comp_names' field as
%              well as '.comp_norm' field is copied to the map struct.
%     (matrix) A data matrix, size dlen x dim. The data matrix may
%              contain unknown values, indicated by NaNs.
%
% OPTIONAL INPUT ARGUMENTS
%
%  argID (string) Argument identifier string (see below).
%  value (varies) Value for the argument (see below).
%
% Here are the valid argument IDs and corresponding values. The values
% which are unambiguous (marked with '*') can be given without the
% preceeding argID.
%   'init'       *(string) initialization: 'randinit' or 'lininit' (default)
%   'munits'      (scalar) the preferred number of map units
%   'msize'       (vector) map grid size
%   'mapsize'    *(string) do you want a 'small', 'normal' or 'big' map
%                          Any explicit settings of munits or msize override this.
%   'lattice'    *(string) map lattice, 'hexa' or 'rect'
%   'shape'      *(string) map shape, 'sheet', 'cyl' or 'toroid'
%   'neigh'      *(string) neighborhood function, 'gaussian', 'cutgauss',
%                          'ep' or 'bubble'
%   'topol'      *(struct) topology struct
%   'som_topol','sTopol' = 'topol'
%   'mask'        (vector) BMU search mask, size dim x 1
%   'name'        (string) map name
%   'comp_names'  (string array / cellstr) component names, size dim x 1
%   'tracking'    (scalar) how much to report, default = 1
%   'training'    (string) 'short', 'default' or 'long'
%                 (vector) size 1 x 2, first length of rough training in epochs,
%                          and then length of finetuning in epochs
%
% OUTPUT ARGUMENTS
%
%  sMap (struct) the trained map struct
%
% EXAMPLES
%
%  To simply train a map with default parameters:
%
%   sMap = som_make(D);
%
%  With the optional arguments, the initialization and training can be
%  influenced. To change map size, use 'msize', 'munits' or 'mapsize'
%  arguments:
%
%   sMap = som_make(D,'mapsize','big'); or sMap=som_make(D,'big');
%   sMap = som_make(D,'munits', 100);
%   sMap = som_make(D,'msize', [20 10]);
%
%  Argument 'algorithm' can be used to switch between 'seq' and 'batch'
%  algorithms. 'batch' is the default, so to use 'seq' algorithm:
%
%   sMap = som_make(D,'algorithm','seq'); or sMap = som_make(D,'seq');
%
%  The 'tracking' argument can be used to control the amout of reporting
%  during training. The argument is used in this function, and it is
%  passed to the training functions. To make the function work silently
%  set it to 0.
%
%   sMap = som_make(D,'tracking',0);
%
% SEE ALSO
%
%  som_map_struct   Create a map struct.
%  som_topol_struct Default values for SOM topology.
%  som_train_struct Default values for SOM training parameters.
%  som_randinint    Random initialization algorithm.
%  som_lininit      Linear initialization algorithm.
%  som_seqtrain     Sequential training algorithm.
%  som_batchtrain   Batch training algorithm.
%  som_impbatch     The imputation SOM batch training algorithm.

% Copyright (c) 1999-2012 by the SOM toolbox programming team.
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 2.0beta juuso 111199
% Version 2.1beta tommi 121212

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% check arguments

% D
if isstruct(D)
  data_name = D.name;
  comp_names = D.comp_names;
  comp_norm = D.comp_norm;
  D = D.data;
else
  data_name = inputname(1);
  sDummy = som_data_struct(D(1,:));
  comp_names = sDummy.comp_names;
  comp_norm = sDummy.comp_norm;
end
[dlen dim] = size(D);

% defaults
mapsize = '';
sM = som_map_struct(dim);
sTopol = sM.topol;
munits = prod(sTopol.msize); % should be zero
mask = sM.mask;
name = sM.name;
neigh = sM.neigh;
tracking = 1;
algorithm = 'batch';
initalg = 'lininit';
training = 'default';

% varargin
i=1;
while i<=length(varargin),
  argok = 1;
  if ischar(varargin{i}),
    switch varargin{i},
      % argument IDs
     case 'mask',       i=i+1; mask = varargin{i};
     case 'munits',     i=i+1; munits = varargin{i};
     case 'msize',      i=i+1; sTopol.msize = varargin{i};
                        munits = prod(sTopol.msize);
     case 'mapsize',    i=i+1; mapsize = varargin{i};
     case 'name',       i=i+1; name = varargin{i};
     case 'comp_names', i=i+1; comp_names = varargin{i};
     case 'lattice',    i=i+1; sTopol.lattice = varargin{i};
     case 'shape',      i=i+1; sTopol.shape = varargin{i};
     case {'topol','som_topol','sTopol'},
                        i=i+1; sTopol = varargin{i}; munits = prod(sTopol.msize);
     case 'neigh',      i=i+1; neigh = varargin{i};
     case 'tracking',   i=i+1; tracking = varargin{i};
     case 'init',       i=i+1; initalg = varargin{i};
     case 'training',   i=i+1; training = varargin{i};
     case 'mweight',    i=i+1; sTopol.mweight = varargin{i};
      % unambiguous values
     case {'hexa','rect'}, sTopol.lattice = varargin{i};
     case {'sheet','cyl','toroid'}, sTopol.shape = varargin{i};
     case {'gaussian','cutgauss','ep','bubble'}, neigh = varargin{i};
     case {'seq','batch','sompak','imp'}, algorithm = varargin{i};
     case {'small','normal','big'}, mapsize = varargin{i};
     case {'randinit','lininit'}, initalg = varargin{i};
     case {'short','default','long'}, training = varargin{i};
      otherwise, argok=0;
    end
  elseif isstruct(varargin{i}) && isfield(varargin{i},'type'),
    switch varargin{i}(1).type,
     case 'som_topol', sTopol = varargin{i};
       otherwise, argok=0;
    end
  else
    argok = 0;
  end
  if ~argok,
    disp(['(som_make) Ignoring invalid argument #' num2str(i+1)]);
  end
  i = i+1;
end

if any(isnan(D(:))) && ~strcmp('imp',algorithm)
  disp('Missing values in data. Consider using the Imputation SOM (see help)')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% make the map struct

%% map size
if isempty(sTopol.msize) || ~prod(sTopol.msize),
  if tracking>0, fprintf(1,'Determining map size...\n'); end
  if ~munits,
    sTemp = som_topol_struct('dlen',dlen);
    munits = prod(sTemp.msize);
    switch mapsize,
     case 'small', munits = max(9,ceil(munits/4));
     case 'big',   munits = munits*4;
     otherwise % nil
    end
  end
  sTemp = som_topol_struct('data',D,'munits',munits);
  sTopol.msize = sTemp.msize;
  if tracking>0,
    fprintf(1,' map size [%d, %d]\n',sTopol.msize(1), sTopol.msize(2));
  end
end

% map struct
sMap = som_map_struct(dim,sTopol,neigh,'mask',mask,'name',name, ...
                      'comp_names', comp_names, 'comp_norm', comp_norm);

% function
if strcmp(algorithm,'sompak'),
  algorithm = 'seq';
  func = 'sompak';
else
  func = algorithm;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% initialization

if tracking>0, fprintf(1,'Initialization...\n'); end

switch initalg,
 case 'randinit', sMap = som_randinit(D, sMap);
 case 'lininit', sMap = som_lininit(D, sMap);
end
sMap.trainhist(1) = som_set(sMap.trainhist(1),'data_name',data_name);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% training
if tracking>0, fprintf(1,'Training using %s algorithm...\n',algorithm); end

% rough train
if tracking>0, fprintf(1,'Rough training phase...\n'); end
sTrain = som_train_struct(sMap,'dlen',dlen,'algorithm',algorithm,'phase','rough');
sTrain = som_set(sTrain,'data_name',data_name);
if isnumeric(training), sTrain.trainlen = training(1);
else
  switch training,
   case 'short', sTrain.trainlen = max(1,sTrain.trainlen/4);
   case 'long',  sTrain.trainlen = sTrain.trainlen*4;
  end
end
switch func,
 case 'seq',    sMap = som_seqtrain(sMap,D,sTrain,'tracking',tracking,'mask',mask);
 case 'sompak', sMap = som_sompaktrain(sMap,D,sTrain,'tracking',tracking,'mask',mask);
 case 'batch',  sMap = som_batchtrain(sMap,D,sTrain,'tracking',tracking,'mask',mask);
 case 'imp',    sMap = som_impbatch(sMap,D,sTrain,'tracking',tracking,'mask',mask);
end

% finetune
if tracking>0, fprintf(1,'Finetuning phase...\n'); end
sTrain = som_train_struct(sMap,'dlen',dlen,'phase','finetune');
sTrain = som_set(sTrain,'data_name',data_name,'algorithm',algorithm);
if isnumeric(training), sTrain.trainlen = training(2);
else
  switch training,
   case 'short', sTrain.trainlen = max(1,sTrain.trainlen/4);
   case 'long',  sTrain.trainlen = sTrain.trainlen*4;
  end
end
switch func,
 case 'seq',    sMap = som_seqtrain(sMap,D,sTrain,'tracking',tracking,'mask',mask);
 case 'sompak', sMap = som_sompaktrain(sMap,D,sTrain,'tracking',tracking,'mask',mask);
 case 'batch',  sMap = som_batchtrain(sMap,D,sTrain,'tracking',tracking,'mask',mask);
 case 'imp',    sMap = som_impbatch(sMap,D,sTrain,'tracking',tracking,'mask',mask);
end

% quality
if tracking>0,
  [mqe,tge,cbe] = som_quality(sMap,D);
  fprintf('Final quantization error: %5.3f\n',mqe)
  fprintf('Final topographic error:  %5.3f\n',tge)
  fprintf('Final combined error:  %5.3f\n',cbe)
end

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

author: ilarinieminen

repo: SOM-Toolbox


text:
function fs = get_seq ( d, norm )

%*****************************************************************************80
%
%% GET_SEQ generates all positive integer D-vectors that sum to NORM.
%
%  Discussion:
%
%    This function computes a list, in reverse dictionary order, of
%    all D-vectors of positive values that sum to NORM.
%
%    For example, fs = get_seq ( 3, 5 ) returns
%
%      3  1  1
%      2  2  1
%      2  1  2
%      1  3  1
%      1  2  2
%      1  1  3
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 May 2010
%
%  Author:
%
%    Original MATLAB version by Florian Heiss, Viktor Winschel.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Florian Heiss, Viktor Winschel,
%    Likelihood approximation by numerical integration on sparse grids,
%    Journal of Econometrics,
%    Volume 144, 2008, pages 62-80.
%
%  Parameters:
%
%    Input, integer D, the dimension.
%
%    Input, integer NORM, the value that each row must sum to.
%    NORM must be at least D.
%
%    Output, integer FS(*,D).  Each row of FS represents one vector
%    with all elements positive and summing to NORM.
%
  if ( norm < d )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GET_SEQ - Fatal error!\n' );
    fprintf ( 1, '  NORM = %d < D = %d.\n', norm, d );
    error ( 'GET_SEQ - Fatal error!' );
  end

  seq = zeros ( 1, d );
%
%  The algorithm is written to work with vectors whose minimum value is
%  allowed to be zero.  So we subtract D from NORM at the beginning and
%  then increment the result vectors by 1 at the end!
%
  a = norm - d;
  seq(1) = a;
  fs = seq;
  c = 1;

  while ( seq (d) < a )

    if ( c == d )
      for i = (c-1) : -1 : 1
        c = i;
        if ( seq(i) ~= 0 )
          break
        end
      end
    end

    seq(c) = seq(c) - 1;
    c = c + 1;
    seq(c) = a - sum ( seq(1:(c-1)) );

    if ( c < d )
      seq((c+1):d) = zeros ( 1, d - c );
    end

    fs = [ fs; seq ];

  end
%
%  Increment all entries by 1.
%
  fs = fs + 1;

  return
end


text:
function sphere_triangle_test01 ( )

%*****************************************************************************80
%
%% SPHERE_TRIANGLE_TEST01 tests SPHERE01_TRIANGLE_QUAD_01, 02, 03.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2014
%
%  Author:
%
%    John Burkardt
%
  e_test = [ ...
    0, 0, 0; ...
    1, 0, 0; ...
    0, 1, 0; ...
    0, 0, 1; ...
    2, 0, 0; ...
    0, 2, 2; ...
    2, 2, 2; ...
    0, 2, 4; ...
    0, 0, 6; ...
    1, 2, 4; ...
    2, 4, 2; ...
    6, 2, 0; ...
    0, 0, 8; ...
    6, 0, 4; ...
    4, 6, 2; ...
    2, 4, 8; ...
   16, 0, 0 ]';

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  Approximate an integral on a random spherical triangle.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  QUAD_01 uses centroids of spherical triangles.\n' );
  fprintf ( 1, '  QUAD_02 uses vertices of spherical triangles.\n' );
  fprintf ( 1, '  QUAD_03 uses midsides of spherical triangles.\n' );
%
%  Choose three points at random to define a spherical triangle.
%
  [ v1, seed ] = sphere01_sample ( 1, seed );
  [ v2, seed ] = sphere01_sample ( 1, seed );
  [ v3, seed ] = sphere01_sample ( 1, seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Vertices of random spherical triangle:\n' );
  fprintf ( 1, '\n' );
  r8vec_transpose_print ( 3, v1, '  V1:' );
  r8vec_transpose_print ( 3, v2, '  V2:' );
  r8vec_transpose_print ( 3, v3, '  V3:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  QUAD_01      QUAD_02      QUAD_03\n' );

  for j = 1 : 17

    e(1:3) = e_test(1:3,j);

    e = polyterm_exponent ( 'SET', e );

    polyterm_exponent ( 'PRINT', e );

    result_01 = sphere01_triangle_quad_01 ( v1, v2, v3, @polyterm_value_3d );

    result_02 = sphere01_triangle_quad_02 ( v1, v2, v3, @polyterm_value_3d );

    result_03 = sphere01_triangle_quad_03 ( v1, v2, v3, @polyterm_value_3d );

    fprintf ( 1, '  %14.6g  %14.6g  %14.6g\n', result_01, result_02, result_03 );

  end

  return
end


text:
%GM_PHD_Predict_Existing
%Last modified 27th August 2013
%Matlab code by Bryan Clarke b.clarke@acfr.usyd.edu.au

%This file performs prediction for existing targets
s = sprintf('Step 2: Prediction for existing targets.');
disp(s);

mk_k_minus_1_before_prediction = mk_minus_1;

for j = 1:size(mk_minus_1,2)
    wk_minus_1(j) = prob_survival * wk_minus_1(j);
    mk_minus_1(:,j) = F * mk_minus_1(:,j); %Assume constant velocity.

    P_range = calculateDataRange4(j);
    P_i = Q + F * Pk_minus_1(:,P_range) * F';

    prevState = mk_k_minus_1_before_prediction(:,j);
    newState = mk_minus_1(:,j);

    Pk_minus_1(:,P_range) = P_i;

    if(VERBOSE == 1)
        s = sprintf('\t\tExisting target %d. Previously at %3.4f %3.4f, now at %3.4f %3.4f.', j, prevState(1), prevState(2), newState(1), newState(2));
        disp(s);

        s = sprintf('\t\tP was %3.4f %3.4f, NOW %3.4f %3.4f', Pk_minus_1(1,P_range(1)), Pk_minus_1(2,P_range(2)), P_i(1,1), P_i(2,2));
        disp(s);
    end
end

%% Now we combine the birthed targets with the existing ones.
%Append newly birthed targets (in m_k_minus_1) to back of old ones
wk_k_minus_1 = [wk_minus_1, w_birth, w_spawn];
mk_k_minus_1 = [mk_minus_1, m_birth, m_spawn];
Pk_k_minus_1 = [Pk_minus_1, P_birth, P_spawn];
numTargets_Jk_k_minus_1 = numTargets_Jk_minus_1 + numBirthedTargets + numSpawnedTargets;
%Create a backup to allow for augmenting the measurement in the update
mk_k_minus_1_before_prediction = [mk_k_minus_1_before_prediction, m_birth_before_prediction];%m_birth_before_prediction also contains the spawned targets before prediction

if(VERBOSE == 1)
    s = sprintf('\tPerformed prediction for %d birthed and existing targets in total.', numTargets_Jk_k_minus_1);
    disp(s);
end

author: Sable

repo: mcbench-benchmarks


text:
function S=nowStr(f)

switch nargin
    case 0
        S = datestr(now,'yyyy_mm_dd');
    case 1
        S = datestr(now,f);
end

%%
% _*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 [mappedX, mapping] = landmark_isomap(X, no_dims, k, percentage)
%ISOMAP Runs the Isomap algorithm
%
%   [mappedX, mapping] = landmark_isomap(X, no_dims, k, percentage);
%
% The functions runs the Landmark Isomap algorithm on dataset X to reduce the
% dimensionality of the dataset to no_dims. The number of neighbors used in
% the compuations is set by k (default = 12). The variable percentage has to be
% between 0 and 1, and determines the number of landmarks that is used (default = 0.1).
%
% If the neighborhood graph that is constructed is not completely
% connected, only the largest connected component is embedded. The indices
% of this component are returned in mapping.conn_comp.
%
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

    if ~exist('no_dims', 'var')
        no_dims = 2;
    end
    if ~exist('k', 'var')
        k = 12;
    end
	if ~exist('percentage', 'var')
		percentage = 0.2;
    end

    % Construct neighborhood graph
    disp('Constructing neighborhood graph...');
    D = find_nn(X, k);

    % Select largest connected component
    blocks = components(D)';
    count = zeros(1, max(blocks));
    for i=1:max(blocks)
        count(i) = length(find(blocks == i));
    end
    [count, block_no] = max(count);
    conn_comp = find(blocks == block_no);
    D = D(conn_comp, conn_comp);
    mapping.D = D;
    n = size(D, 1);

    % Compute shortest paths
    disp('Computing shortest paths...');
    landmarks = randperm(n);
	landmarks = landmarks(1:round(percentage * n));
    nl = length(landmarks);
    D = dijkstra(D, landmarks);
	D = full(D)';
    mapping.DD = D;
    mapping.landmarks = landmarks;

	% Do not embed in more dimensions than (nl - 1)
	if no_dims > nl - 1
		no_dims = nl - 1;
		warning(['Target dimensionality reduced to ' num2str(no_dims) '...']);
	end

    % Performing MDS using eigenvector implementation
    disp('Constructing low-dimensional embedding...');
    D = D .^ 2;
    subB = -.5 .* (bsxfun(@minus, bsxfun(@minus, D, sum(D, 2) ./ nl), sum(D, 1) ./ n) + sum(D(:)) ./ (n * nl));
    subB2 = subB' * subB;
	subB2(isnan(subB2)) = 0;
	subB2(isinf(subB2)) = 0;
    [alpha, beta] = eig(subB2);
    val = beta .^ (1 / 2);
    mapping.invVal = inv(val);
    vec = subB * alpha * mapping.invVal;
	if size(vec, 2) < no_dims
		no_dims = size(vec, 2);
		warning(['Target dimensionality reduced to ' num2str(no_dims) '...']);
	end

    % Computing final embedding
    [val, ind] = sort(real(diag(val)), 'descend');
    vec = vec(:,ind(1:no_dims));
    val = val(1:no_dims);
    mappedX = real(bsxfun(@times, vec, sqrt(val)'));

    % Store data for out-of-sample extension
    mapping.conn_comp = conn_comp;
    mapping.k = k;
    mapping.X = X(conn_comp,:);
    mapping.alpha = alpha;
    mapping.beta = beta;
    mapping.no_dims = no_dims;

repo: cluster


text:
function [HDR] = save2bkr(arg1,arg2,arg3);
% SAVE2BKR loads EEG data and saves it in BKR format
% The following data formats are supported:
%	CNT, EDF, BKR, MAT, etc. format
%
%       HDR = save2bkr(sourcefile [, destfile [, option]]);
%
%       HDR = eegchkhdr();
%	HDR = save2bkr(HDR,data);
%
%   sourcefile	sourcefile wildcards are allowed
%   destfile	destination file in BKR format
%	if destfile is empty or a directory, sourcefile but with extension .bkr is used.
%   options
%       gain            Gain factor for unscaled EEG data (e.g. old Matlab files)
%       'removeDC'      removes mean
%       'regressEOG k:l,m:n'     removes EOG (channels m:n) from EEG (channels k:l)
%       'autoscale k:l'	uses only channels from k to l for scaling
%       'detrend k:l'	channels from k to l are detrended with an FIR-highpass filter.
%       'PhysMax=XXX'	uses a fixed scaling factor; might be important for concanating BKR files
%			+XXX and -XXX correspond to the maximum and minimum physical value, resp.
% 		You can concanate several options by separating with space, komma or semicolon
%
%   HDR		Header, HDR.FileName must contain target filename
%   data	data samples
%
% Examples:
%   save2bkr('/tmp/*.cnt',[],'autoscale 5:30');
%	converts all CNT-files from subdir /tmp/ into BKR files
%       and saves them in the current directory
%   save2bkr('/tmp/*.cnt','/tmp2/','autoscale 5:30, PhysMax=200');
%	converts all CNT-files from subdir /tmp/ into BKR files
%       and saves them in the directory /tmp2/
%
%
%
% see also: EEGCHKHDR, REGRESS_EOG, SLOAD

%	$Revision: 1.23 $
% 	$Id: save2bkr.m 2205 2009-10-27 12:18:15Z schloegl $
%	Copyright (C) 2002-2003 by Alois Schloegl <a.schloegl@ieee.org>
%    	This is part of the BIOSIG-toolbox http://biosig.sf.net/

% This library is free software; you can redistribute it and/or
% modify it under the terms of the GNU Library General Public
% License as published by the Free Software Foundation; either
% Version 2 of the License, or (at your option) any later version.
%
% This library 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
% Library General Public License for more details.
%
% You should have received a copy of the GNU Library General Public
% License along with this library; if not, write to the
% Free Software Foundation, Inc., 59 Temple Place - Suite 330,
% Boston, MA  02111-1307, USA.

FLAG_REGRESS_EOG = 0;
FLAG_REMOVE_DC = 0;
FLAG_AUTOSCALE = 0;
FLAG_DETREND = 0;
FLAG_PHYSMAX = 0;
FLAG_removeDrift = 0;

chansel = 0;

if nargin<2, arg2=[]; end;

if nargin<3,
        cali = NaN;
elseif isnumeric(arg3)
        cali = arg3;
else
        FLAG_REMOVE_DC = findstr(lower(arg3),'removedc');

        tmp = findstr(arg3,'autoscale');
        if ~isempty(tmp);
                [chansel,tmp] = strtok(arg3(tmp+9:length(arg3)),' ;+');
                tmp = str2num(chansel);
                if isempty(tmp),
                        fprintf(2,'invalid autoscale argument %s',chansel);
                        return;
                else
                        FLAG_AUTOSCALE = 1;
                        chansel = tmp;
                end;
        end;

        tmp = findstr(lower(arg3),'detrend');
        if ~isempty(tmp);
                [chansel_dt,tmp] = strtok(arg3(tmp+7:length(arg3)),' ;,+');
                tmp = str2num(chansel_dt);
                if isempty(tmp),
                        fprintf(2,'invalid detrend argument %s',chansel_dt);
                        return;
                else
                        FLAG_DETREND = 1;
                        chansel_dt = tmp;
                end;
        end;

        tmp = findstr(lower(arg3),'removedrift');
        if ~isempty(tmp);
                [chansel_dt2,tmp] = strtok(arg3(tmp+11:length(arg3)),' ;+');
                tmp = str2num(chansel_dt2);
                if isempty(tmp),
                        fprintf(2,'invalid RemoveDrift argument %s',chansel_dt2);
                        return;
                else
                        FLAG_removeDrift = 1;
                        chansel_dt2 = tmp;
                end;
        end;

        tmp = findstr(lower(arg3),'regresseog');
        if ~isempty(tmp);
                [chansel_dt3,tmp] = strtok(arg3(tmp+11:length(arg3)),' ;,+');
                [chansel_dt4,tmp] = strtok(tmp,' ;,+');
                tmp = str2num(chansel_dt3);
                FLAG_REGRESS_EOG = ~isempty(tmp);
                if isempty(tmp),
                        fprintf(2,'invalid REGRESSEOG argument %s',chansel_dt3);
                        return;
                else
                        FLAG_REGRESS_EOG = 1;
                        chansel_dt3 = tmp;
                end;
                tmp = str2num(chansel_dt4);
                FLAG_REGRESS_EOG = FLAG_REGRESS_EOG * ~isempty(tmp);
                if isempty(tmp),
                        fprintf(2,'invalid REGRESSEOG argument %s',chansel_dt4);
                        return;
                else
                        chansel_dt4 = tmp;
                end;
        end;
        tmp = findstr(lower(arg3),'physmax=');
        if ~isempty(tmp);
                [tmp,tmp1] = strtok(arg3(tmp+8:length(arg3)),' ;,');
                PHYSMAX = str2num(tmp);
                if isempty(PHYSMAX ),
                        fprintf(2,'invalid PhysMax argument %s',tmp);
                        return;
                else
                        FLAG_PHYSMAX = 1;
                end;
        end;
end;

if 0,
elseif exist(arg1,'file')
	inpath='';
	infile.name=arg1;
        outfile = arg2;
elseif exist(arg1,'dir')
        inpath = fileparts(arg1);
        infile = dir(arg1);	% input  file
        if isempty(infile)
                fprintf(2,'ERROR SAVE2BKR: file %s not found.\n',arg1);
                return;
        end;
        outfile = arg2;
elseif isstruct(arg1) & isnumeric(arg2),
        HDR  = arg1;
        data = arg2;
else  %if isstruct(arg1) & isnumeric(arg2),
        fprintf(2,'Error SAVE2BKR: invalid input arguments\n');
        return;
end;

if isstruct(arg1),
        %HDR.FileName 	= destfile;	% Assign Filename
        if isfield(HDR,'NS')
                if HDR.NS==size(data,2),
                        % It's ok.
                elseif HDR.NS==size(data,1),
                        warning('data is transposed\n');
                        data = data';
                else
                        fprintf(2,'HDR.NS=%i is not equal number of data columns %i\n',HDR.NS,size(data,2));
                        return;
                end;
        else
                HDR.NS = size(data,2);	% number of channels
        end;
        if ~isfield(HDR,'NRec'),
                HDR.NRec = 1;		% number of trials (1 for continous data)
        end;
        HDR.SPR 	= size(data,1)/HDR.NRec;	% number of samples per trial
        %HDR.SampleRate	= 100;		% Sampling rate
        %HDR.Label  	= hdr.Label; % Labels,

        %HDR.PhysMax 	= max(abs(data(:)));	% Physical maximum
        %HDR.DigMax 	= max(2^15-1);		% Digital maximum
        % --- !!! Previous scaling gave an error up to 6% and more !!! ---

        %HDR.Filter.LowPass  = 30;	% upper cutoff frequency
        %HDR.Filter.HighPass = .5;	% lower cutoff frequency
        HDR.FLAG.REFERENCE   = ' ';	% reference '', 'LOC', 'COM', 'LAP', 'WGT'
        %HDR.FLAG.REFERENCE  = HDR.Recording;
        HDR.FLAG.TRIGGERED   = HDR.NRec>1;	% Trigger Flag

        if FLAG_REMOVE_DC,
                %y = center(y,1);
                data = data - repmat(mean(data,1),size(data,1),1);
        end;
        if chansel == 0;
                chansel=1:HDR.NS;
        end;
        tmp = data(:,chansel);
        HDR.PhysMax = max(abs(tmp(:))); %gives max of the whole matrix
        HDR.DigMax = 2^15-1;            % maximum resulution
        for k = 1:HDR.NS,
                if any(k==chansel),
                        data(:,k) = data(:,k)*HDR.DigMax/HDR.PhysMax;
                else
                        mm = max(abs(data(:,k)));
                        data(:,k) = data(:,k)*HDR.DigMax/mm;
                end;
        end;
        HDR.FLAG.UCAL = 1;              % data is de-calibrated, no rescaling within SWRITE
        %HDR = eegchkhdr(HDR);
        HDR.TYPE = 'BKR';

        HDR = sopen (HDR,'w',0);     	% OPEN BKR FILE
        HDR = swrite(HDR,data);  	% WRITE BKR FILE
        %fwrite(HDR.FILE.FID,data','int16');  	% WRITE BKR FILE
        HDR = sclose(HDR);            % CLOSE BKR FILE

        % save Classlabels
        if isfield(HDR,'Classlabel'),
                fid = fopen([HDR.FileName(1:length(HDR.FileName)-4) '.par'],'wt');
                fprintf(fid, '%i\n', HDR.Classlabel);
                fclose(fid);
        end;
        if isfield(HDR,'ArtifactSelection'),
                fid = fopen([HDR.FileName(1:length(HDR.FileName)-4) '.sel'],'w');
                fprintf(fid, '%i\r\n', HDR.ArtifactSelection);
                fclose(fid);
        end;

        % final test
        try
                HDR = sopen(HDR.FileName,'r');
                HDR = sclose(HDR);
        catch
                fprintf(2,'Error SAVE2BKR: saving file %s failed\n',HDR.FileName);
        end;
        return;
end;

for k=1:length(infile);
        filename = fullfile(inpath,infile(k).name);
        [pf,fn,ext] = fileparts(filename);

        % load eeg data
        [y,HDR] = sload(filename);

        % load classlabels if the exist
        tmp = fullfile(HDR.FILE.Path,[HDR.FILE.Name,'.mat']);
        if exist(tmp),
                tmp=load(tmp);
                if isfield(tmp,'classlabel');
                        HDR.Classlabel = tmp.classlabel;
                end;
        end;

        if isempty(y),
                fprintf(2,'Error SAVE2BKR: file %s not found\n',filename);
                return;
        end;

        if ~isfield(HDR,'NS'),
                warning(['number of channels undefined in ',filename]);
                HDR.NS = size(y,2);
        end;
        if ~HDR.FLAG.TRIGGERED,
                HDR.NRec = 1;
                HDR.SPR = size(y,1);
        end;
        if ~isfield(HDR,'NRec'),
                HDR.NRec = 1;
        end;
        if ~isfield(HDR,'SPR'),
                HDR.SPR = size(y,1)/HDR.NRec;
        elseif length(HDR.SPR)>1,       % use just one sampling rate
                HDR.SPR = HDR.AS.MAXSPR;
                HDR.SampleRate = HDR.AS.MAXSPR/HDR.Dur;
                FLAG_PHYSMAX = 1;
                PHYSMAX = max(abs(y(:)));
                HDR.DigMax  = 2^15-1;
        end;

        if FLAG_REGRESS_EOG,
                fprintf(1,'\tREGRESS_EOG \n');
                [R,y] = regress_eog(y,chansel_dt3,chansel_dt4);
        end;

        if FLAG_REMOVE_DC,
                fprintf(1,'\tREMOVE_DC \n');
                y = y - repmat(mean(y,1),size(y,1),1);
        end;
        if FLAG_DETREND,
                B = -ones(1,HDR.SampleRate)/HDR.SampleRate;
                B(HDR.SampleRate/2) = B(HDR.SampleRate/2)+1;
                HDR.Filter.B = B;
                HDR.Filter.A = 1;
                %HDR.Filter.B=B;%conv(-B, HDR.Filter.B);
                Delay = length(B)/2;
                HDR.FLAG.FILT = 1;
                HDR.Filter.HighPass = .5;

                for k = chansel_dt,
                        tmp = filter(B,1,[y(:,k);zeros(length(B),1)]);
                        y(:,k) = tmp(Delay+1:size(y,1)+Delay);
                end;
        end;

        if FLAG_removeDrift,
                B = .5*(1 - cos(2*pi*(1:4*HDR.SampleRate+1)'/(4*HDR.SampleRate+2)));
                B = -B/sum(B);
                B(2*HDR.SampleRate) = B(HDR.SampleRate)+1;

                B = -ones(1,HDR.SampleRate)/HDR.SampleRate;
                B(HDR.SampleRate/2) = B(HDR.SampleRate/2)+1;
                %B(1) = B(1)+1;

                HDR.Filter.B = B;
                HDR.Filter.A = 1;
                %HDR.Filter.B=B;%conv(-B, HDR.Filter.B);
                Delay = (length(B)-1)/2;
                HDR.FLAG.FILT = 1;
                HDR.Filter.HighPass = .5;

                for k = chansel_dt2,
                        y(:,k) = filtfilt(B,1,y(:,k));
                        %y(:,k) = tmp(Delay+1:size(y,1)+Delay);
                end;
        end;

        if chansel == 0;
                chansel=1:HDR.NS;
        end;

        % add event channel
        if isfield(HDR,'EVENT')
                if ~length(HDR.EVENT.TYP),
                elseif 0,
                        % TypeList = unique(HDR.EVENT.TYP); but ignores NaN's
                        [sY ,idx] = sort(HDR.EVENT.TYP(:));
                        TypeList  = sY([1;find(diff(sY,1)>0)+1]);

                        event = zeros(size(y,1),length(TypeList));
                        for k2 = 1:length(TypeList),
                                tmp = (HDR.EVENT.TYP==TypeList(k2));
                                event(HDR.EVENT.POS(tmp),k2) = HDR.EVENT.TYP(tmp);
                        end;
                        HDR.NS = HDR.NS + size(event,2);
                        y = [y, event];
                elseif all(HDR.EVENT.TYP < 256),  % only NeuroScan Events are converted into separate channels
                        K = 0; event = [];
                        for k1 = 0:7,
                                tmp = bitand(HDR.EVENT.TYP,2^k1);
                                if any(tmp),
                                        K = K+1;
                                        event(size(y,1),K) = 0;
                                        event(HDR.EVENT.POS(tmp>0),K) = 1;
                                end;
                        end;
                        if any(sum(event,2)>1),
                                fprintf(2,'Warning SAVE2BKR: simulateneous events occur. \n');
                        end;
                        HDR.NS = HDR.NS + size(event,2);
                        y = [y, event];
                end;
        end;

        % re-scale data to account for the scaling factor in the header
        HDR.DigMax = 2^15-1;
        if FLAG_PHYSMAX,
                HDR.PhysMax = PHYSMAX;
        else
                tmp = y(:,chansel);
                HDR.PhysMax = max(abs(tmp(:))); %gives max of the whole matrix
        end;
        for k = 1:HDR.NS,
                if any(k==chansel),
                        y(:,k) = y(:,k)*HDR.DigMax/HDR.PhysMax; % keep correct scaling factor
                else
                        mm = max(abs(y(:,k)));
                        y(:,k) = y(:,k)*HDR.DigMax/mm;          % scale to maximum resolution
                end;
        end;
        HDR.FLAG.UCAL = 1;              % data is de-calibrated, no rescaling within SWRITE

        tmp = round(HDR.PhysMax);
        fprintf(1,'Rounding of PhysMax yields %f%% error.\n',abs((HDR.PhysMax-tmp)/tmp)*100);
        HDR.PhysMax = tmp;
        HDR.TYPE = 'BKR';
        HDR.FLAG.REFERENCE = ' ';
        HDR.FLAG.TRIGGERED = (HDR.NRec>1);

        if isempty(outfile), 	% default destination directory
                ix = max(find(filename=='.'));
                %HDR.FileName = [filename(1:ix-1),'.bkr'];  % destination directory is same as source directory
                HDR.FileName  = [HDR.FILE.Name,'.bkr'];     % destination directory is current working directory
        elseif isdir(outfile),	% output file
                HDR.FILE.Path = outfile;
                HDR.FileName  = fullfile(outfile,[HDR.FILE.Name,'.bkr']);
        else
                [HDR.FILE.Path,HDR.FILE.Name,Ext] = fileparts(outfile);
                HDR.FileName = fullfile(HDR.FILE.Path,[HDR.FILE.Name,Ext]);
        end;
        %HDR = eegchkhdr(HDR);

        HDR = sopen(HDR,'w');
        if HDR.FILE.FID < 0,
                fprintf(1,'Error SAVE2BKR: couldnot open file %s.\n',HDR.FileName);
                return;
        end;
        % writes data
        HDR = swrite(HDR,y(:,1:HDR.NS));  	% WRITE BKR FILE
        %count = fwrite(HDR.FILE.FID,y','short');
        HDR = sclose(HDR);

        % save classlabels
        if isfield(HDR,'Classlabel'),
                if ~isempty(HDR.Classlabel),
                        fid = fopen([HDR.FileName(1:length(HDR.FileName)-4) '.par'],'w');
                        fprintf(fid, '%i\r\n', HDR.Classlabel);
                        fclose(fid);
                end;
        end;
        if isfield(HDR,'ArtifactSelection'),
                fid = fopen([HDR.FileName(1:length(HDR.FileName)-4) '.sel'],'w');
                fprintf(fid, '%i\r\n', HDR.ArtifactSelection);
                fclose(fid);
        end;

        % final test
        try
                HDR = sopen(HDR.FileName,'r');
                HDR = sclose(HDR);
        catch
                fprintf(2,'Error SAVE2BKR: saving file %s failed\n',HDR.FileName);
        end;
end;

author: PatternRecognition

repo: OpenBMI


text:
clear all; close all;
I=imread('rice.png');
I=im2double(I);
M=2*size(I,1);
N=2*size(I,2);
u=-M/2:(M/2-1);
v=-N/2:(N/2-1);
[U,V]=meshgrid(u, v);
D=sqrt(U.^2+V.^2);
D0=30;
n=6;
H=1./(1+(D0./D).^(2*n));
J=fftshift(fft2(I, size(H, 1), size(H, 2)));
K=J.*H;
L=ifft2(ifftshift(K));
L=L(1:size(I,1), 1:size(I, 2));
figure;
subplot(121);
imshow(I);
subplot(122);
imshow(L);

author: HuangCongQing


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
% function [ph,th] = plotLM(LM,varargin);
%
% a convenient tool for plotting landmarks
%
% Input
%   LM            landmark positions
%   varargin    optional parameter, see below
%
% Output:
%  ph           handle to plot
%  th           handle to text
%
% see also E5_Hands_TPS for an example
%==============================================================================

function [ph,th] = plotLM(LM,varargin)

if nargin==0
    help(mfilename)
    runMinimalExample;
    return;
end

% setup default parameter
dim       = size(LM,2);
numbering = 'off';
dx        = 0.1;
fontsize  = 20;
color     = 'y';
for k=1:2:length(varargin), % overwrites default parameter
  eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

ph = []; th = []; % create empty handles

% handle options 'numbering', 'fontsize' and 'dx'
J = strcmp(varargin,'numbering') ...
  | strcmp(varargin,'dx') | strcmp(varargin,'fontsize');
if any(J),
  I = find(J);
  K = setdiff(1:length(varargin),[I,I+1]);
  varargin = {varargin{K}};
end;

% use plot in 2D or 3D
varargin = {varargin{:},'color',color};
switch dim,
  case 2, ph = plot(LM(:,1),LM(:,2),'x',varargin{:});
  case 3, ph = plot3(LM(:,1),LM(:,2),LM(:,3),'x',varargin{:});
end;

if strcmp(numbering,'off'), return;  end;

% add labels to the landmarks
dx = [dx,zeros(1,dim-1)];
pos = @(j) LM(j,:) + dx;
for j=1:size(LM,1),
  th(j) = text('position',pos(j),'string',sprintf('%d',j));
end;
set(th,'fontsize',fontsize,'color',color);
%------------------------------------------------------------------------------
function runMinimalExample
help(mfilename);
clear
close all
setup2DhandData
figure;
subplot(1,2,1);
viewImage2D(dataT,omega,m,'colormap','bone(265)');
hold on; ph =plotLM(LM(:,1:2),'color','g'); hold off;
set(ph,'marker','s')
subplot(1,2,2);
viewImage2D(dataR,omega,m,'colormap','bone(265)');
hold on; ph = plotLM(LM(:,3:4),'color','r'); hold off;
set(ph,'marker','o')
%------------------------------------------------------------------------------


text:
function err = hdf5write_safe(fname, name, val, options)

    err = -1;
    if isempty(val)
        return;
    end
    if ~exist('options','var')
        options = '';
    end

    force_scalar = false;
    force_array = false;
    if any(strcmp(options, 'array'))
       force_array = true;
    elseif any(strcmp(options, 'scalar'))
        force_scalar = true;
    end
    % Identify type of val and use SNIRF v1.1-compliant write function

    if exist(fname,'file')
        fid = H5F.open(fname, 'H5F_ACC_RDWR', 'H5P_DEFAULT');
    else
        fid = H5F.create(fname, 'H5F_ACC_TRUNC', 'H5P_DEFAULT', 'H5P_DEFAULT');
    end
    if fid < 0
        err = -1;
        return;
    end

    % Create group where dataset is located
    group = filesepStandard(fileparts(name), 'nameonly');
    gid = HDF5_CreateGroup(fid, group);
    if gid < 0
        err = -1;
        return;
    end

    if iscell(val) || isstring(val)
        if length(val) > 1 && ~force_scalar || force_array  % Returns true for single strings, believe it or not
            write_string_array(fid, fname, name, val);
        else
            write_string(fid, fname, name, val);
        end
        return
    elseif ischar(val)
        write_string(fid, fname, name, val);
        return
    elseif isfloat(val)
        if length(val) > 1 && ~force_scalar || force_array
            write_numeric_array(fname, name, val);
        else
            write_numeric(fid, fname, name, val);
        end
        return
    elseif isinteger(val)
        if length(val) > 1 && ~force_scalar || force_array
            write_numeric_array(fid, fname, name, val);  % As of now, no integer arrays exist
        else
            write_integer(fid, fname, name, val);
        end
        return
    else
        warning(['An unrecognized variable was saved to ', name, ' in ', fname])
    end

    H5F.close(fid);

end

function err = write_string(fid, fname, name, val)
    sid = H5S.create('H5S_SCALAR');
    tid = H5T.copy('H5T_C_S1');
    H5T.set_size(tid, 'H5T_VARIABLE');
    did = H5D.create(fid, name, tid, sid, 'H5P_DEFAULT');
    H5D.write(did, tid, 'H5S_ALL', 'H5S_ALL', 'H5P_DEFAULT', {val});
    err = 0;
end

function err = write_string_array(fid, fname, name, val)
    val = HDF5_Transpose(val);
    sid = H5S.create_simple(1, numel(val), H5ML.get_constant_value('H5S_UNLIMITED'));
    tid = H5T.copy('H5T_C_S1');
    H5T.set_size(tid, 'H5T_VARIABLE');
    pid = H5P.create('H5P_DATASET_CREATE');
    H5P.set_chunk(pid, 2);
    did = H5D.create(fid, name, tid, sid, pid);
    H5D.write(did, tid, 'H5S_ALL', 'H5S_ALL', 'H5P_DEFAULT', val);
    err = 0;
end

function err = write_numeric(fid, fname, name, val)
    fid = H5F.open(fname, 'H5F_ACC_RDWR', 'H5P_DEFAULT');
    tid = H5T.copy('H5T_NATIVE_DOUBLE');
    sid = H5S.create('H5S_SCALAR');
    H5D.create(fid, name, tid, sid, 'H5P_DEFAULT');
    h5write(fname, name, val);
    err = 0;
end

function err = write_numeric_array(fname, name, val)
    val = HDF5_Transpose(val);
    sizeval = size(val);
    if sizeval(1) == 1 || sizeval(2) == 1
        n = length(val);
    else
        n = sizeval;
    end
    h5create(fname, name, n, 'Datatype', 'double');
    h5write(fname, name, val);
    err = 0;
end

function err = write_integer(fid, fname, name, val)
    warning off;  % Suppress the int truncation warning
    tid = H5T.copy('H5T_NATIVE_INT');
    sid = H5S.create('H5S_SCALAR');
    H5D.create(fid, name, tid, sid, 'H5P_DEFAULT');
    h5write(fname, name, val);
    err = 0;
    warning on;
end

author: BUNPC


text:
%% Reads ECG recordings in Mortara format.
% Reads ECG recordings in Mortara format.
%
% Arguments:
%   + filename: recording to be read or folder containing HourXXRawData.bin
%   files.
%   + start_sample: (opt) start sample to read. Default 1.
%   + end_sample: (opt) end sample to read. Default min(All recording, ECG block of 200 Mbytes)
%
% Output:
%   + ECG: the ECG block
%   + heasig: header with the ECG properties.
%   + ann: annotations for the ECG recordings.
%
% Limits:
% This routine is limited to read blocks smaller than 200 Mbytes for
% performance reasons. You can disable this limit by doing:
% MaxIOread = Inf; %megabytes
%
% See also read_Mortara_header, read_ECG, ECGwrapper
%
% Author: Mariano Llamedo Soria
% <matlab:web('mailto:llamedom@electron.frba.utn.edu.ar','-browser') (email)>
% Version: 0.1 beta
% Birthdate: 29/05/2014
% Last update: 19/11/2014
% Copyright 2008-2015
%
function [ECG heasig end_sample] = read_Mortara_format( filename, start_sample, end_sample )

heasig.freq = 1000; %Hz
heasig.desc = char({'I','II','V1','V2','V3','V4','V5','V6'});
leads_in_mortara_files = 8;

ECG = [];
last_sample = [];

%No leer bloques mas grandes de 200 megabytes
MaxIOread = 200; %megabytes

if( nargin < 2 || isempty( start_sample ) )
    start_sample = 1;
else
    start_sample = max(1,start_sample);
end

if( isdir(filename) )
    rec_path = filename;
else
    [rec_path, rec_name] = fileparts(filename);
end

hours_files = dir([rec_path filesep 'Hour*.bin' ]);

hours_available = length(hours_files);

samples_per_hour = zeros(hours_available,1);

expected_samples = 60*60*heasig.freq;

for ii = 1:hours_available

    hour_filename = [rec_path filesep 'Hour' num2str(ii) 'RawData.bin' ];
    fid = fopen(hour_filename);

    if( fid > 0 )
        fseek(fid, 0, 'eof');
        bytes_totales = ftell(fid);
        samples_per_hour(ii) = bytes_totales / 2 / leads_in_mortara_files;
        fclose(fid);
        if( ii ~= expected_samples && samples_per_hour(ii) ~= expected_samples)
            error('read_Mortara_format:Impossible2read', 'Not enough samples in %s. Expected %d - found %d. Problem copying ?\n', hour_filename, expected_samples, samples_per_hour(ii));
        end
    else
        error('read_Mortara_format:Impossible2read', 'Problem opening %s.\n', hour_filename);
    end

end

total_samples = sum(samples_per_hour);

if( nargin < 3 || isempty( end_sample ) )
    %Intento la lectura total por defecto
    end_sample = total_samples;
end

samples2read = end_sample - start_sample + 1;

if( (samples2read*leads_in_mortara_files*2) > (MaxIOread * 1024^2) )
    samples2read = (MaxIOread * 1024^2) / leads_in_mortara_files / 2;
    end_sample = start_sample + samples2read - 1;
    warning(['No es recomendable leer mas de ' num2str(MaxIOread) ' Mb. Realice varias lecturas.'])
end

end_samples = cumsum(samples_per_hour);
start_samples = [1; end_samples(1:end-1)+1 ];
files2read = find(start_sample <= end_samples & end_sample >= start_samples );

for ii = rowvec(files2read)

    hour_filename = [rec_path filesep 'Hour' num2str(ii) 'RawData.bin' ];
    fidECG = fopen( hour_filename, 'r');

    if( fidECG > 0 )

        aux_offset = 0;

        if(ii == files2read(1))

            aux_offset = start_sample - start_samples(ii);
            fseek(fid, aux_offset * 2 * leads_in_mortara_files, 'cof');

        end

        if(ii == files2read(end))

            if( length(files2read) == 1 )
                % read samples from a single file
                this_samples2read = samples2read;
            else
                % read samples from several files
                this_samples2read = end_sample - start_samples(ii) + 1;
            end

        else

            if( aux_offset == 0)
                this_samples2read = samples_per_hour(ii);
            else
                this_samples2read = samples_per_hour(ii) - aux_offset;
            end

        end

        ECG = [ECG; fread(fid, [leads_in_mortara_files this_samples2read], '*int16')'];

        fclose(fid);
    else
        error('read_Mortara_format:Impossible2read', 'Problem opening %s.\n', hour_filename);
    end

end

[heasig.nsamp, heasig.nsig] = size(ECG);

if( abs(heasig.nsamp - samples2read) > 1 )
    warning('read_Mortara_format:ReadError', 'Problem reading %s. Required %d samples, and only read %d.\n', filename, samples2read, heasig.nsamp);
end

heasig.adczero = zeros(heasig.nsig,1);
heasig.gain = ones(heasig.nsig,1);
heasig.units = repmat('uV',heasig.nsig,1);

if( rec_path(end) == filesep )
    rec_path = rec_path(1:end-1);
end
[~, rec_name] = fileparts(rec_path);

heasig.recname = rec_name;
heasig.btime = calc_btime( start_sample, heasig.freq );
heasig.bdate = '01/01/2000';

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function test_ft_megrealign

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_megrealign ft_prepare_neighbours ft_transform_geometry

datainfo = ref_datasets;
sel      = match_str({datainfo.datatype},{'bti148' 'bti248' 'ctf151' 'ctf275' 'itab153' 'yokogawa160'}');
datainfo = datainfo(sel);

% as of yet, the code does not explicitly test anything, it just checks
% whether all runs through smoothly
for k = 1:numel(datainfo)
  fname = fullfile(datainfo(k).origdir,'latest/raw',datainfo(k).type,['preproc_',datainfo(k).datatype]);
  load(fname);

  cfg = [];
  cfg.channel = 'MEG';
  cfg.demean = 'yes';
  data = ft_preprocessing(cfg, data);

  vol      = [];
  vol.o    = [0 0 4];
  vol.r    = 10;
  vol.unit = 'cm';

  % ensure units in the gradiometer array and volume conductor to be equal
  data.grad = ft_convert_units(data.grad, 'cm');
  vol       = ft_convert_units(vol, data.grad.unit);

  % make 2 copies
  data2 = data;
  data3 = data;

  data2.grad = ft_transform_geometry([[eye(3) [0 0 -2]'];[0 0 0 1]], data.grad);
  data3.grad = ft_transform_geometry([1 0 0 0;0 cos(0.1) -sin(0.1) 0;0 sin(0.1) cos(0.1) 0;0 0 0 1], data.grad);

  cfg = [];
  cfg.template{1} = data.grad;
  cfg.template{2} = data2.grad;
  cfg.template{3} = data3.grad;
  cfg.inwardshift = 1;
  cfg.headmodel   = vol;

  interp = ft_megrealign(cfg, data);
  interp2 = ft_megrealign(cfg, data2);
  interp3 = ft_megrealign(cfg, data3);

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function q = times(q1,q2,takeRight)
% quaternion .* quaternion and quaternion .* vector3d
%
% Syntax
%   q = q1 .* q2
%
%   q = times(q1,q2)
%   q = times(q1,q2,takeRight)
%
% Input
%  q1 - @quaternion
%  q2 - @quaternion
%  takeRight - logical, use as output left or right input
%
% Output
%  q  - @quaternion

if isa(q1,'quaternion') && isa(q2,'quaternion')

  % which input will become the output?
  if nargin == 3
    if takeRight, q = q2; else, q = q1; end
  elseif isa(q1,'rotation')
    q = q1;
  else
    q = q2;
  end

  a1 = q1.a; b1 = q1.b; c1 = q1.c; d1 = q1.d;
  a2 = q2.a; b2 = q2.b; c2 = q2.c; d2 = q2.d;

  %standard algorithm
  q.a = a1 .* a2 - b1 .* b2 - c1 .* c2 - d1 .* d2;
  q.b = b1 .* a2 + a1 .* b2 - d1 .* c2 + c1 .* d2;
  q.c = c1 .* a2 + d1 .* b2 + a1 .* c2 - b1 .* d2;
  q.d = d1 .* a2 - c1 .* b2 + b1 .* c2 + a1 .* d2;

elseif isa(q1,'quaternion') && isa(q2,'double')

  q1.a = q1.a .* q2;
  q1.b = q1.b .* q2;
  q1.c = q1.c .* q2;
  q1.d = q1.d .* q2;
  q = q1;

elseif isa(q2,'quaternion') && isa(q1,'double')

  q2.a = q1 .* q2.a;
  q2.b = q1 .* q2.b;
  q2.c = q1 .* q2.c;
  q2.d = q1 .* q2.d;
  q = q2;

else

  q = rotate(q2,q1);

end

author: mtex-toolbox

repo: mtex

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


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

%%
clear; close all; clc;

%% Syntax
% |[D]=logic2levelset(varargin);|

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


text:
function montage_clusters_text(ovl,clusters,varargin)
% :Inputs:
%
%   **varargin:**
%        additional clusters structures
%        this function puts text cluster numbers on cluster centers
%
%   **color:**
%        cell array of text strings indicating colors {'r' 'g'} etc...
%
% ..
%    Tor Wager
% ..

if isempty(ovl), ovl = which('scalped_single_subj_T1.img');, end

myc = {'b' 'r' 'y' 'g' 'm' 'c'};

XYZmm = cat(2,clusters.XYZmm);

if length(varargin) > 0
    for i = 1:length(varargin)
        cl = varargin{i};
        clXYZmm = cat(2,cl.XYZmm);
        XYZmm = [XYZmm clXYZmm];
    end
end

% ----------------------------------------
% * overlay image
% ----------------------------------------
V = spm_vol(ovl);
oimg = spm_read_vols(V);
V.M = V.mat;

textx = size(oimg,2) - 50;
texty = size(oimg,1) - 6;

%[array,hdr,h,whichslices,rows,cols,figh] = readim2(ovl,'p');

% how many slices, which ones

XYZ = mm2voxel(XYZmm,V,2)';
whsl = unique(XYZ(3,:));
nsl = length(whsl) + 1;
rc = ceil(sqrt(nsl));
h = [];

% get centers for text plotting
cen = getcenters(clusters);

figure; colormap gray; set(gcf,'Color','w')

% plot first cluster structure
for i = 1:length(clusters)
    index = plot_cluster(clusters(i),oimg,rc,V,whsl,myc{1},textx,texty,1,size(oimg),cen(i,:),num2str(i));
    drawnow
end

% plot additional cluster structures

if length(varargin) > 0
    for i = 1:length(varargin)
        cl = varargin{i};
        cen = getcenters(cl);
        index = plot_cluster(cl,[],rc,V,whsl,myc{i+1},textx,texty,i+1,size(oimg),cen);

    end
end

return

% sub-functions
%

function index = plot_cluster(clusters,oimg,rc,V,whsl,myc,textx,texty,clind,odims,varargin)
% var argument is list of centers, x y z are rows, cols are coords

XYZmm = cat(2,clusters.XYZmm);
XYZ = mm2voxel(XYZmm,V)';

% text numbers, if plotting text at xyz coords of centers (varargin)
if length(varargin) > 0
    cencoo = varargin{1};
    cencoo = mm2voxel(cencoo,V)';
    centext = varargin{2};
    %cenind = 1;
end

% surface patch method
% ----------------------------------------------------------------------------------------
vol = voxel2mask(XYZ',odims);
vol = smooth3(vol);

index = 1;
for z = whsl

    subplot(rc,rc,index);

    if ~isempty(oimg)
        set(gca,'YDir','reverse');
        imagesc(oimg(:,:,z))
        hold on; axis image; axis off
        zmm = voxel2mm([1 1 z]',V.mat);
        text(textx,texty,['z = ' num2str(zmm(3))],'Color','w')
    else
        hold on
    end

    if z>1,FVC = isocaps(vol(:,:,z-1:z),0,'zmax');
    else FVC = isocaps(vol(:,:,z:z+1),0,'zmax');
    end
    try
	patch(FVC,'EdgeColor','none','FaceColor',myc,'FaceAlpha',.7)
    catch
	patch(FVC,'EdgeColor','none','FaceColor',myc)
    end

    % text numbers, if plotting text at xyz coords of centers (varargin)
    if length(varargin) > 0
        cencooz = cencoo(:,cencoo(3,:) == z);
        for i = 1:size(cencooz,2)
            text(cencooz(2,i),cencooz(1,i),centext,'Color','k'); %num2str(cenind),'Color','k')
            %cenind = cenind + 1;
        end
    end

    % plot method

    %myXYZ = XYZ(1:2,XYZ(3,:) == z);
    %plot(myXYZ(2,:),myXYZ(1,:),[myc 's'],'MarkerFaceColor',myc,'MarkerSize',3)

    index = index + 1;
    drawnow

end

subplot(rc,rc,index)
a = pwd; a = a(end-6:end);
b = num2str(clusters(1).threshold);

c = num2str(length(clusters));
text(0,clind-1,[myc ': ' clusters(1).title ' ' a ' u = ' b ', ' c ' clusters'])
axis off
axis([0 1 -1 clind])

return

function cen = getcenters(clusters)
    for i = 1:length(clusters)
        cen(i,:) = clusters(i).mm_center;
    end
    cen = cen';
return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function hmObj = hlp_validateHeadModelObject(hmObj)
% helper function to check the head model object, and
% load from base workspace if necessary
% Author: Tim Mullen, SCCN/INC/UCSD 2013

if isempty(hmObj)
    hmObj = 'hmObj'; % default name of variable in base wkspace
end
if ischar(hmObj)
    % load variable from base workspace
    try
       hmObj = evalin('base',hmObj);
    catch err1
        % try to open it as a file names
        try
            if ~exist(hmObj,'file'), hmObj = env_translatepath(hmObj); end
            hmObj = hlp_microcache('headmodels',@headModel.loadFromFile,hmObj);
        catch err2
            disp_once('BCILAB:HeadmodelMissing',['Note: The head model named "' hmObj '" was neither found in the base workspace nor could it be loaded as a file.']);
            hmObj = '';
            return;
        end
    end
end

% check that the object is a valid head model object
if ~isobject(hmObj)
    warning('You did not supply a valid head model object');
    hmObj = '';
    return;
end
if ~isprop(hmObj,'leadFieldFile') || isempty(hmObj.leadFieldFile)
    warning('Head model object does not appear to be valid. The lead field matrix does not exist.');
    hmObj = '';
    return;
end

repo: REST


text:
classdef BayesianfastF0NLS< handle

%     usuage

% obj= BayesianfastF0NLS(N, L, pitchBounds, A_var,voicingProb,varargin)
%  N is the frame length
%  L is the maximum allowed harmonic order
%  pitchBounds denotes the pitch range e.g.,[f_min/fs, f_max/fs], where fs is
%  the sampling frequency
%  A_var denotes the frequncy transition variance.
%  voicingProb denotes the initial voicing probability

    % Can only be set in the constructor
    properties (SetAccess=immutable)
        N
        L
        F
        pitchBoundsOuter = [0.0 0.5] % The outer bounds that are acceptable

        % default values
        epsilon = 0.0
        dcIsIncluded = false
        epsilon_ref = 0.0

        % Precomputed quantities
        crossCorrelationVectors
        Gamma1
        Gamma2
        fftShiftVector
    end

    properties (SetAccess=private)
        pitchBounds  % The current active bounds
        fullPitchGrid
        validFftIndices
        defaultRefinementTol
        refinementTol
        gPriorParam = 3
        logPitchPdfs % the pitch pdf (column) for every candidate model order
        logModelPmf % an (L+1)x1 vector
         logpi
        A
        scaled_alpha_buffer
        scaled_alpha_buffer2
        unvoicing_scaled_alpha_buffer
        logModelPrior;
        %% varSpeech is a regularization parameter and it should be chosen according to the input SNR
        % In high SNR, a larger value leads to better performance from our
        % experience.
        varSpeech=5e-6;
        %%
        B
        C
        norml_factor=0;
        cnt=1;
        prew_state
    end

    methods

        function obj = BayesianfastF0NLS(N, L, pitchBounds, A_var,voicingProb,varargin)
            switch nargin
                case 0
                     N=0.03*16000;
                     L=10;
                     pitchBounds=[50/16000,400/16000];
                     A_var=[2/16000];
                     voicingProb=0.7;
                case 1
                     L=10;
                     pitchBounds=[50/16000,400/16000];
                     A_var=[2/16000];
                     voicingProb=0.7;
                case 2
                    pitchBounds=[50/16000,400/16000];
                     A_var=[2/16000];
                     voicingProb=0.7;
                case 3
                     A_var=[2/16000];
                     voicingProb=0.7;
                case 4
                     voicingProb=0.7;
            end

            % validate input
            if length(varargin) >= 1
                if islogical(varargin{1})
                    obj.dcIsIncluded = varargin{1};
                else
                    error('Argument 4 is not of type logical (true/false)')
                end
            end
            if length(varargin) >= 2
                if isscalar(varargin{2})
                    obj.F = varargin{2};
                else
                    error('Argument 5 is not a scalar')
                end
            else
                obj.F = 2^14;
            end
            obj.defaultRefinementTol = 1/obj.F;

            if length(varargin) >= 3
                if isscalar(varargin{3})
                    obj.epsilon = varargin{3};
                else
                    error('Argument 6 is not a scalar')
                end
            end

            if length(varargin) > 4
                error('Too many input arguments')
            end

            if ~isscalar(N)
                error('Input argument N is not a scalar')
            else
                obj.N = N;
            end
            if ~isscalar(L)
                error('Input argument L is not a scalar')
            else
                obj.L = L;
            end

            if ~ (isvector(pitchBounds) && length(pitchBounds) == 2)
                error(['Input argument pitchBound must be a 2-' ...
                       'vector'])
            elseif pitchBounds(1) < obj.pitchBoundsOuter(1) || ...
                pitchBounds(2) > obj.pitchBoundsOuter(2)

                error(['Input argument pitchBounds must be within the ' ...
                       'bounds specified in the constructor (at ' ...
                       'least [0.0 0.5])'])
            else
                obj.pitchBounds = pitchBounds;
            end

            if pitchBounds(1) < 1/N
                warning(['The lower pitch bound is set lower than one '...
                    ' period/segment. Inaccurate results might be '...
                    'produced - especially if you do not set the'...
                    ' regularisation parameter.']);
            end

            % Init
            F = obj.F;
            minFftIndex = ceil(F*pitchBounds(1));
            maxFftIndex = floor(F*pitchBounds(2));
            obj.validFftIndices = (minFftIndex:maxFftIndex)';
            obj.fullPitchGrid = obj.validFftIndices/F;
            nPitches = length(obj.fullPitchGrid);
            obj.logModelPrior=[voicingProb,1-voicingProb];

            % cross-correlation vectors
            obj.crossCorrelationVectors = ...
                [N*ones(1, nPitches)/2 + N*obj.epsilon;...
                sin(pi*(1:2*L)'*obj.fullPitchGrid'*N)./...
                (2*sin(pi*(1:2*L)'*obj.fullPitchGrid'))];

            obj.fftShiftVector = ...
                exp(1i*2*pi*(0:ceil(F/2)-1)'*(N-1)/(2*F));

            % Compute Gamma for the T+H and T-H systems
            [obj.Gamma1, obj.Gamma2] = computeGamma(L, F, pitchBounds, ...
                obj.crossCorrelationVectors, nPitches, ...
                obj.validFftIndices, ...
                obj.dcIsIncluded);
            %%          added by Liming Shi
            obj.scaled_alpha_buffer=nan;
            grid_width=diff(obj.fullPitchGrid(1:2));
%             obj.logpi=(log(1/(nPitches-1))-log(grid_width))*ones(nPitches,1);
            obj.logpi=log((voicingProb)/(nPitches*L))*ones(nPitches,1);
            obj.A=zeros(nPitches,nPitches);

            for jj=1:nPitches
%                 f0=obj.fullPitchGrid(jj);
%                  obj.A(jj,:)= abs(sin(pi*(obj.fullPitchGrid-f0)/f0)./pi./(obj.fullPitchGrid-f0)*f0)+.1;
                 obj.A(jj,:)=(normpdf(obj.fullPitchGrid,obj.fullPitchGrid(jj),A_var));
%                  var_value=1e-1;
%                  obj.logA(jj,:)=log(lognpdf(obj.fullPitchGrid,log(obj.fullPitchGrid(jj)*exp(var_value)),var_value));
%                    obj.logA(jj,:)=log(normpdf(log(obj.fullPitchGrid),log(obj.fullPitchGrid(jj)),.1));
            end

            obj.A=obj.A./repmat(sum(obj.A,2),1,nPitches);

            obj.B=zeros(obj.L,obj.L);

            for jj=1:obj.L
                 obj.B(jj,:)=normpdf(jj,[1:obj.L],1);
            end

            obj.B=obj.B./repmat(sum(obj.B,2),1,obj.L);

            obj.C=[.7 .3;.4,.6];

            obj.scaled_alpha_buffer2=1/obj.L/nPitches*ones(nPitches,obj.L);

        end

        function varargout = computeCostFunctions(obj, x)
            % Compute and returns the cost functions for l=1,...,L

            % validate input
            if ~isvector(x) && length(x) == obj.N
                error(['First argument x must be vector of ' ...
                               'length N=', num2str(obj.N)]);
            end
            x = reshape(x, obj.N, 1); % make sure its a column
                                      % vector

            varargout{1} = computeAllCostFunctions(x, obj.L, ...
            	obj.fullPitchGrid, obj.fftShiftVector, ...
                obj.crossCorrelationVectors, obj.Gamma1, obj.Gamma2, ...
                obj.dcIsIncluded);

            if nargout == 2
                varargout{2} = obj.fullPitchGrid;
            end

        end

        function varargout = estimate(obj, x, flag,tinc,varargin)
            % Estimate fundamental frequency and order of the signal x

            % validate input
            if ~isvector(x) && length(x) == obj.N
                error(['First argument x must be vector of ' ...
                               'length N=', num2str(obj.N)]);
            end
            x = reshape(x, obj.N, 1); % make sure its a column
                                      % vector

            if length(varargin) == 1
                if isscalar(varargin{1})
                    obj.refinementTol = varargin{1};
                else
                    error('Argument 2 is not a scalar')
                end
            else
                obj.refinementTol = obj.defaultRefinementTol;
            end
            if flag==1
                spec_pow_all=(abs(fft(x,2048)).^2)';
                if isempty(obj.prew_state)
                    [psdouput,obj.prew_state]=estnoiseg_orig(spec_pow_all,tinc);
                else
                    [psdouput,obj.prew_state]=estnoiseg_orig(spec_pow_all,obj.prew_state);
                end
                mmse_covariance_all = (real(ifft(psdouput')));
                mmse_lpc_all= levinson(mmse_covariance_all,30);
                [xp] = filter(mmse_lpc_all,1,...
                    x);
                  x=xp/norm(xp)*norm(x);
            end

            % if the data consists of all zeros then return null model
            if  x'*x/obj.N<1e-10
                estimatedPitch = nan;
                estimatedOrder = 0;
                unvoicing_scaled_alpha=log(1);
            else
                % Step 2: compute the profile log-likelihood function
                % over only the fundamental frequency (the linear
                % parameters, noise variance and g have been integrated
                % out)
%                 x=x/norm(x);
                costs = obj.computeCostFunctions(x);

%                 for ii=1:size(costs,1)
%                     regvalue=1e-5;
%                     tempt=costs(ii,:)/(2*regvalue);
%                     pitchLogLikelihood(ii,:)=tempt+log(tempt)*(-3/2-ii+1)+...
%                         log(gammainc(costs(ii,:),3/2+ii-1));
%                 end

                cod = costs*(1/(x'*x+obj.varSpeech*obj.N));
% %                  cod = (costs)*(1/(x'*x));
% %                   cod = costs*(1/(x'*x));
                [~, pitchLogLikelihood] = ...
                    computePitchLogLikelihood(cod, obj.N, obj.gPriorParam);
% %                 pitchLogLikelihood=pitchLogLikelihood;
                null_modellike=1;

                    %% added/modified by Liming Shi (implements the hidden markov chain)
                if isnan(obj.scaled_alpha_buffer)
                    inx=isnan(pitchLogLikelihood);
                    pitchLogLikelihood(inx)=-inf;
                    bar_alpha=obj.logpi+pitchLogLikelihood';

                    unvoicing_bar_alpha=log(obj.logModelPrior(2)*null_modellike);
                    log_scale=log_sumsum_exp_ls([[unvoicing_bar_alpha,-inf*ones(1,obj.L-1)];bar_alpha]);
                    scaled_alpha=bar_alpha-log_scale;
                    unvoicing_scaled_alpha=unvoicing_bar_alpha-log_scale;

                else
                    inx=isnan(pitchLogLikelihood);
                    pitchLogLikelihood(inx)=-inf;
                    state_prior=obj.C(1,1)*obj.A'*obj.scaled_alpha_buffer*obj.B;

                    state_prior=state_prior+obj.scaled_alpha_buffer2*obj.C(2,1)*obj.unvoicing_scaled_alpha_buffer;

                    bar_alpha=log(state_prior)+pitchLogLikelihood';

                    temp=[1-obj.unvoicing_scaled_alpha_buffer;obj.unvoicing_scaled_alpha_buffer];

                    unvoicing_bar_alpha=log(obj.C(:,2)'*temp*null_modellike);
                    log_scale=log_sumsum_exp_ls([[unvoicing_bar_alpha,-inf*ones(1,obj.L-1)];bar_alpha]);
                    scaled_alpha=bar_alpha-log_scale;
                    unvoicing_scaled_alpha=unvoicing_bar_alpha-log_scale;

                end

                    [value, inx] = ...
                        max(scaled_alpha);
                    [~,inx_2]=max(value);
                    pitchIndex=inx(inx_2);
                    estimatedOrder=inx_2;
                    estimatedPitch = obj.fullPitchGrid(pitchIndex(1));
             %% added by Liming Shi (buffering the previous meaningful estimates of the posterior states)

                if unvoicing_scaled_alpha < log(0.5)
                    obj.scaled_alpha_buffer2=exp(scaled_alpha-log(1-exp(unvoicing_scaled_alpha)));
                else
%                     estimatedOrder=nan;
                end
                obj.scaled_alpha_buffer=exp(scaled_alpha);
                obj.unvoicing_scaled_alpha_buffer=exp(unvoicing_scaled_alpha);
            end

            %%

            if nargout >= 1
                varargout{1} = estimatedPitch;
            end
            if nargout >= 2
                varargout{2} = estimatedOrder;
            end
            if nargout >= 3
%                 [~, alpha] = objFunction(estimatedPitch, x, ...
%                     estimatedOrder, obj.dcIsIncluded);
                varargout{3} = 1-exp(unvoicing_scaled_alpha);
            end

        end
        function obj=destroy(obj)
        end
    end
end
function y=log_sumsum_exp_ls(x)
max_temp=max(max(x));
inf_ind=isinf(max_temp);
y=log(sum(sum(exp(x-max_temp))))+max_temp;
y(inf_ind)=-inf;
end

author: LimingShi


text:
function cascade_data(model, data_year, pca_dim)
% cascade_data(model, data_year, pca)
%
% Compute score statistics for filters and deformation models
% using the model and training data associated with the model
% and dataset year 'data_year'.  If PCA is given a value > 0, then
% the score statistics are computed using a PCA projection of the
% model and training data.
%
% The score statistics are written in to the file
% class_year_cascade_data_{no_pca,pca_<pca_dim>}_data_year.mat,
% which is used by cascade_model.m.
%
% model      Object model
% data_year  Dataset year as a string (e.g., '2007')
% pca_dim    Number of PCA components to use

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2009-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

conf = voc_config();
cscdir = conf.cascade.data_dir;

model.interval = conf.training.interval_fg;

% get training data
pos = pascal_data(model.class, model.year);

%[model, pca_model] = project_model_adapt(model, pca_dim);
% Loads <coeff>
load('pca.mat');
[model, pca_model] = project_model(model, coeff, pca_dim);

numpos  = length(pos);
pixels  = model.minsize * model.sbin / 2;
minsize = prod(pixels);
pars    = cell(1, numpos);

% compute latent filter locations and record target bounding boxes
parfor i = 1:numpos
  pars{i}.scores           = [];
  pars{i}.pca_scores       = [];
  pars{i}.total_scores     = [];
  pars{i}.total_pca_scores = [];
  fprintf('%s %s: cascade data: %d/%d\n', procid(), model.class, i, numpos);
  bbox = pos(i).boxes;
  % skip small examples
  if (bbox(3)-bbox(1)+1)*(bbox(4)-bbox(2)+1) < minsize
    continue;
  end
  % get example
  im = imreadx(pos(i));
  [im, bbox] = croppos(im, bbox);
  [pyra, model_dp] = gdetect_pos_prepare(im, model, bbox, 0.5);
  [ds, bs, trees] = gdetect_pos(pyra, model_dp, 1, 1, 0.5);
  if ~isempty(ds)
    % collect cascade score statistics
    pars{i}.scores = get_block_scores(pyra, model, trees);
    pars{i}.total_scores = ds(1,end);

    % Valid root location for the PCA model
    t = tree_mat_to_struct(trees{1});
    valid = [];
    valid.c = t(1).rule_index;
    valid.y = t(1).y;
    valid.x = t(1).x;
    valid.l = t(1).l;

    % Get detection with PCA model restricted to the root filter being placed
    % exactly where the original model's root filter was placed
    [pca_pyra, pca_model_dp] = gdetect_pos_prepare_c(im, pca_model, valid);
    [ds, bs, trees] = gdetect_pos_c(pca_pyra, pca_model_dp, valid);

    % Sanity check
    t = tree_mat_to_struct(trees{1});
    assert(valid.c == t(1).rule_index);
    assert(valid.y == t(1).y);
    assert(valid.x == t(1).x);
    assert(valid.l == t(1).l);

    pars{i}.pca_scores = get_block_scores(pca_pyra, pca_model, trees);
    pars{i}.total_pca_scores = ds(1,end);
  end
end
% Collate
pars = cell2mat(pars);
scores = cat(1, pars(:).scores);
pca_scores = cat(1, pars(:).pca_scores);

% Sanity check that detection scores match sum of block scores
% computed in get_block_scores
total_scores = cat(1, pars(:).total_scores);
total_pca_scores = cat(1, pars(:).total_pca_scores);
assert(max(abs(scores(:,1) - total_scores)) < 1e-5);
assert(max(abs(pca_scores(:,1) - total_pca_scores)) < 1e-5);

% Save cascade score statistics
class_year = [model.class '_' model.year];
file     = [cscdir class_year '_cascade_data_no_pca_' model.year];
pca_file = [cscdir class_year '_cascade_data_pca_' num2str(pca_dim) ...
            '_' model.year];
save(file, 'scores');
save(pca_file, 'pca_scores');

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function [Pj closest] = P_r2j(Pr, Pj, closest)
%P_R2J  Change reference frame of projection matrices
%
%   [Pk I] = P_r2j(Pr, Pj, closest)
%
% Tranforms the coordinate frame of a set of projection matrices such that
% a reference projection matrix becomes the identity matrix eye(3, 4). Can
% also select input matrices based on camera centre proximity to the
% reference frame.
%
%IN:
%   Pr - 3x4 projection matrix of the reference frame.
%   Pj - 3x4xN array of input projection matrices.
%   closest - 1xP indices of matrices of Pj to output, once arranged in
%             order of proximity of camera centre to Pr. Default: 1:N.
%
%OUT:
%   Pk - 3x4x(min(N,P)) array of transformed projection matrices.
%   I - 1x(min(N,P)) indices of matrices of Pj that make up Pk, such
%       that Pk = T(Pj(:,:,I)), T() being the coordinate transformation.

% $Id: P_r2j.m,v 1.2 2007/12/10 10:58:31 ojw Exp $

if nargin > 2
    % Determine which Pj matrices are required based on distance of optical
    % centres from that of Pr.
    T = zeros(3, size(Pj, 3));
    for a = 1:size(Pj, 3)
        T(:,a) = -Pj(:,1:3,a) \ Pj(:,4,a);
    end
    order = ojw_bsxfun(@minus, T, -Pr(:,1:3) \ Pr(:,4));
    [T order] = sort(sum(order .^ 2));
    closest = order(closest);
    Pj = Pj(:,:,closest);
end

% Calculate Pj relative to Pr
% Assume projection matrices make 0,0 the top left corner of the top left
% pixel. Convert to Matlab form, i.e. 0.5,0.5 is the top left corner of the
% top left pixel.
T = [1 0 0.5; 0 1 0.5; 0 0 1];
Pr = inv([T * Pr; 0 0 0 1]);
for a = 1:size(Pj, 3)
    Pj(:,:,a) = T * Pj(:,:,a) * Pr;
end
return

repo: SelfExSR


text:
function f = eval(SO3VF,rot,varargin)
%
% syntax
%   f = eval(SO3VF,rot)
%
% Input
%   rot - @rotation
%
% Output
%   f - @vector3d
%

% if isa(rot,'orientation')
%   ensureCompatibleSymmetries(SO3VF,rot)
% end

f = vector3d(SO3VF.SO3F.eval(rot)');
f = reshape(f',size(rot));

end


text:
function [Z_hat, E_hat] = exact_alm_lrr_l21v2(D, A, lambda, tol, maxIter,display)

% Aug 2013
% This matlab code implements the Exact ALM algorithm for
% min_{Z,E}  |Z|_* + lambda |E|_2,1  s.t.  D = AZ + E
%
% D - m x n matrix of observations/data (required input)
% A - m x k matrix of the dictionary (required input)

% lambda - weight on sparse error term in the cost function
%
% tol - tolerance for stopping criterion.
%     - DEFAULT 1e-7 if omitted or -1.
%
% maxIter - maximum number of iterations
%         - DEFAULT 1000, if omitted or -1.
%
[m n] = size(D);
k = size(A,2);

if nargin < 4 || isempty(tol)
    tol = 1e-7;
end

if nargin < 5 || isempty(maxIter)
    maxIter = 1000;
end

if nargin<6 || isempty(display)
    display = false;
end

maxIter_primal = 10000;
% initialize
Y = sign(D);
norm_two = norm(Y,2);
norm_inf = norm( Y(:), inf) / lambda;
dual_norm = max(norm_two, norm_inf);
Y = Y / dual_norm;

W = zeros(k,n);

Z_hat = zeros(k,n);
E_hat = zeros(m,n);
%parameters
dnorm = norm(D, 'fro');
tolProj1 = 1e-6 * dnorm;

anorm = norm(A,2);
tolProj2 = 1e-6 * dnorm/anorm;

mu = .5/norm_two; % this one can be tuned
rho = 6;          % this one can be tuned

%pre-computation
if m>=k
    inv_ata = inv(eye(k) + A'*A);
else
    inv_ata = eye(k) - A'/(eye(m)+A*A')*A;
end

iter = 0;
while iter < maxIter
    iter = iter + 1;

    % solve the primal problem by alternative projection
    primal_iter = 0;

    while primal_iter < maxIter_primal
        primal_iter = primal_iter + 1;
        temp_Z = Z_hat;
        temp_E = E_hat;

        %update J
        temp = temp_Z + W/mu;

        [U,S,V] = svd(temp, 'econ'); % stable
        %[U,S,V] = svdecon(temp); % fastest, but EIG must not contain NaN or Inf.

        diagS = diag(S);
        svp = length(find(diagS > 1/mu));
        diagS = max(0,diagS - 1/mu);

        if svp < 0.5 %svp = 0
            svp = 1;
        end
        J_hat = U(:,1:svp)*diag(diagS(1:svp))*V(:,1:svp)';

        % update Z
        temp = J_hat + A'*(D - temp_E) + (A'*Y-W)/mu;
        Z_hat = inv_ata*temp;

        %update E
        temp = D - A*Z_hat + Y/mu;
        E_hat =  solve_l1l2(temp, lambda/mu);

        if norm(E_hat - temp_E, 'fro') < tolProj1 && norm(Z_hat - temp_Z)<tolProj2
            break;
        end
    end

    H1 = D - A*Z_hat - E_hat;
    H2 = Z_hat - J_hat;
    Y = Y + mu*H1;
    W = W + mu*H2;
    mu = rho * mu;

    %% stop Criterion
    stopCriterion = max(norm(H1, 'fro')/dnorm, norm(H2,'fro')/dnorm*anorm);
    if display
        disp(['LRR: Iteration' num2str(iter) '(' num2str(primal_iter) '), mu ' num2str(mu) ', |E|_2,0 ' num2str(sum(sum(E_hat.^2,1)>0))...
        ', stopCriterion ' num2str(stopCriterion)]);
    end

    if stopCriterion < tol
        break;
    end

end

end

author: andrewssobral

repo: lrslibrary


text:
function ell = inertiaEllipsoid(points)
%INERTIAELLIPSOID Inertia ellipsoid of a set of 3D points.
%
%   Note: Deprecated! Use equivalentEllipsoid instead.
%
%
%   ELL = inertiaEllipsoid(PTS)
%   Compute the inertia ellipsoid of the set of points PTS. The result is
%   an ellipsoid defined by:
%   ELL = [XC YC ZC A B C PHI THETA PSI]
%   where [XC YC ZY] is the center, [A B C] are lengths of semi-axes (in
%   decreasing order), and [PHI THETA PSI] are euler angles representing
%   the ellipsoid orientation, in degrees.
%
%   Example
%     pts = randn(300, 3);
%     pts = transformPoint3d(pts, createScaling3d([6 4 2]));
%     pts = transformPoint3d(pts, createRotationOx(pi/6));
%     pts = transformPoint3d(pts, createRotationOy(pi/4));
%     pts = transformPoint3d(pts, createRotationOz(pi/3));
%     pts = transformPoint3d(pts, createTranslation3d([5 4 3]));
%     elli = inertiaEllipsoid(pts);
%     figure; drawPoint3d(pts); axis equal;
%     hold on; drawEllipsoid(elli, ...
%         'drawEllipses', true, 'EllipseColor', 'b', 'EllipseWidth', 3);
%
%   See also
%     equivalentEllipsoid
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2011-03-12, using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform

% deprecation warning
warning('geom3d:deprecated', ...
    [mfilename ' is deprecated, use ''equivalentEllipsoid'' instead']);

% number of points
n = size(points, 1);

% compute centroid
center = mean(points);

% compute the covariance matrix
covPts = cov(points)/n;

% perform a principal component analysis with 2 variables,
% to extract inertia axes
[U, S] = svd(covPts);

% extract length of each semi axis
radii = sqrt(5) * sqrt(diag(S)*n)';

% sort axes from greater to lower
[radii, ind] = sort(radii, 'descend');

% format U to ensure first axis points to positive x direction
U = U(ind, :);
if U(1,1) < 0
    U = -U;
    % keep matrix determinant positive
    U(:,3) = -U(:,3);
end

% convert axes rotation matrix to Euler angles
angles = rotation3dToEulerAngles(U);

% concatenate result to form an ellipsoid object
ell = [center, radii, angles];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
classdef CEC2010_F10 < PROBLEM
% <single> <real> <constrained>
% CEC'2010 constrained optimization benchmark problem

%------------------------------- Reference --------------------------------
% R. Mallipeddi and P. N. Suganthan, Problem definitions and evaluation
% criteria for the CEC 2010 competition on constrained real-parameter
% optimization, Nanyang Technological University, Singapore, 2010.
%------------------------------- 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".
%--------------------------------------------------------------------------

    properties
        O;      % Optimal decision vector
        Mat;	% Rotation matrix
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2010.mat'),'Data');
            obj.O = Data{10}.O;
            obj.M = 1;
            if isempty(obj.D) || obj.D < 30
                obj.D   = 10;
                obj.Mat = Data{10}.M_10;
            else
                obj.D   = 30;
                obj.Mat = Data{10}.M_30;
            end
            if isempty(obj.D); obj.D = 10; end
            obj.D = min(obj.D,length(obj.O));
            obj.lower    = zeros(1,obj.D) - 500;
            obj.upper    = zeros(1,obj.D) + 500;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            Z = 1 + PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopObj = sum(100*(Z(:,1:end-1).^2-Z(:,2:end)).^2+(Z(:,1:end-1)-1).^2,2);
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            Y = Z*obj.Mat;
            PopCon = abs(sum((Y.*sin(sqrt(abs(Y)))),2)) - 1e-4;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function sv=geph2pos(time,geph,sv)

TSTEP=60.0;
ERREPH_GLO=5.0;

t=timediff(time,geph.toe);
dts=-geph.taun+geph.gamn*t;

x(1:3)=geph.pos;
x(4:6)=geph.vel;

if t<0,tt=-TSTEP;else,tt=TSTEP;end

while abs(t)>1E-9
    if abs(t)<TSTEP,tt=t;end
    x=glorbit(tt,x,geph.acc);
    t=t-tt;
end

rs=x(1:3);
vars=ERREPH_GLO^2;

sv.pos=rs'; sv.dts=dts; sv.vars=vars;

return

repo: GINav


text:
function rezToPhy(rez, savePath, varargin)
% pull out results from kilosort's rez to either return to workspace or to
% save in the appropriate format for the phy GUI to run on. If you provide
% a savePath it should be a folder, and you will need to have npy-matlab
% available (https://github.com/kwikteam/npy-matlab)

[~, Nfilt, Nrank] = size(rez.W);
rez.Wphy = cat(1, zeros(1+rez.ops.nt0min, Nfilt, Nrank), rez.W); % for Phy, we need to pad the spikes with zeros so the spikes are aligned to the center of the window

% spikeTimes will be in samples, not seconds
rez.W = gather(single(rez.Wphy));
rez.U = gather(single(rez.U));
rez.mu = gather(single(rez.mu));

if size(rez.st3,2)>4
    rez.st3 = rez.st3(:,1:4);
end

[~, isort]   = sort(rez.st3(:,1), 'ascend');
rez.st3      = rez.st3(isort, :);
if ~isempty(rez.cProj)
    rez.cProj    = rez.cProj(isort, :);
    rez.cProjPC  = rez.cProjPC(isort, :, :);
end

% ix = rez.st3(:,4)>12;
% rez.st3 = rez.st3(ix, :);
% rez.cProj = rez.cProj(ix, :);
% rez.cProjPC = rez.cProjPC(ix, :,:);

fs = dir(fullfile(savePath, '*.npy'));
for i = 1:length(fs)
   delete(fullfile(savePath, fs(i).name));
end
if exist(fullfile(savePath, '.phy'), 'dir')
    rmdir(fullfile(savePath, '.phy'), 's');
end

spikeTimes = uint64(rez.st3(:,1));
% account for ops.trange(1) to accomodate real time
spikeTimes = spikeTimes - rez.ops.trange(1)*rez.ops.fs;
% [spikeTimes, ii] = sort(spikeTimes);
spikeTemplates = uint32(rez.st3(:,2));
if size(rez.st3,2)>4
    spikeClusters = uint32(1+rez.st3(:,5));
end
amplitudes = rez.st3(:,3);

Nchan = rez.ops.Nchan;

xcoords     = rez.xcoords(:);
ycoords     = rez.ycoords(:);
chanMap     = rez.ops.chanMap(:);
chanMap0ind = chanMap - 1;

nt0 = size(rez.W,1);
U = rez.U;
W = rez.W;

Nfilt = size(W,2);

templates = zeros(Nchan, nt0, Nfilt, 'single');
for iNN = 1:size(templates,3)
   templates(:,:,iNN) = rez.mu(iNN,1) * squeeze(U(:,iNN,:)) * squeeze(W(:,iNN,:))';
end
templates = permute(templates, [3 2 1]); % now it's nTemplates x nSamples x nChannels
templatesInds = repmat([0:size(templates,3)-1], size(templates,1), 1); % we include all channels so this is trivial

templateFeatures = rez.cProj;
templateFeatureInds = uint32(rez.iNeigh);
pcFeatures = rez.cProjPC;
pcFeatureInds = uint32(rez.iNeighPC);

% whiteningMatrix = rez.Wrot/rez.ops.scaleproc;
whiteningMatrix = eye(size(rez.Wrot)) / rez.ops.scaleproc;
whiteningMatrixInv = whiteningMatrix^-1;

% here we compute the amplitude of every template...

% unwhiten all the templates
tempsUnW = zeros(size(templates));
for t = 1:size(templates,1)
    tempsUnW(t,:,:) = squeeze(templates(t,:,:))*whiteningMatrixInv;
end

% The amplitude on each channel is the positive peak minus the negative
tempChanAmps = squeeze(max(tempsUnW,[],2))-squeeze(min(tempsUnW,[],2));

% The template amplitude is the amplitude of its largest channel
tempAmpsUnscaled = max(tempChanAmps,[],2);

% assign all spikes the amplitude of their template multiplied by their
% scaling amplitudes
spikeAmps = tempAmpsUnscaled(spikeTemplates).*amplitudes;

% take the average of all spike amps to get actual template amps (since
% tempScalingAmps are equal mean for all templates)
ta = clusterAverage(spikeTemplates, spikeAmps);
tids = unique(spikeTemplates);
tempAmps = zeros(numel(rez.mu),1);
tempAmps(tids) = ta; % because ta only has entries for templates that had at least one spike
gain = getOr(rez.ops, 'gain', 1);
tempAmps = gain*tempAmps'; % for consistency, make first dimension template number

templateFeatures = [];
if ~isempty(savePath)
    fileID = fopen(fullfile(savePath, 'cluster_KSLabel.tsv'),'w');
    fprintf(fileID, 'cluster_id%sKSLabel', char(9));
    fprintf(fileID, char([13 10]));

    fileIDCP = fopen(fullfile(savePath, 'cluster_ContamPct.tsv'),'w');
    fprintf(fileIDCP, 'cluster_id%sContamPct', char(9));
    fprintf(fileIDCP, char([13 10]));

    fileIDA = fopen(fullfile(savePath, 'cluster_Amplitude.tsv'),'w');
    fprintf(fileIDA, 'cluster_id%sAmplitude', char(9));
    fprintf(fileIDA, char([13 10]));

    rez.est_contam_rate(isnan(rez.est_contam_rate)) = 1;
    for j = 1:length(rez.good)
        if rez.good(j)
            fprintf(fileID, '%d%sgood', j-1, char(9));
        else
            fprintf(fileID, '%d%smua', j-1, char(9));
        end
        fprintf(fileID, char([13 10]));

        fprintf(fileIDCP, '%d%s%.1f', j-1, char(9), rez.est_contam_rate(j)*100);
        fprintf(fileIDCP, char([13 10]));

        fprintf(fileIDA, '%d%s%.1f', j-1, char(9), tempAmps(j));
        fprintf(fileIDA, char([13 10]));

    end
    fclose(fileID);
    fclose(fileIDCP);
    fclose(fileIDA);

    writeNPY(spikeTimes, fullfile(savePath, 'spike_times.npy'));
    writeNPY(uint32(spikeTemplates-1), fullfile(savePath, 'spike_templates.npy')); % -1 for zero indexing
    if size(rez.st3,2)>4
        writeNPY(uint32(spikeClusters-1), fullfile(savePath, 'spike_clusters.npy')); % -1 for zero indexing
    else
        writeNPY(uint32(spikeTemplates-1), fullfile(savePath, 'spike_clusters.npy')); % -1 for zero indexing
    end
    writeNPY(amplitudes, fullfile(savePath, 'amplitudes.npy'));
    writeNPY(templates, fullfile(savePath, 'templates.npy'));
    writeNPY(templatesInds, fullfile(savePath, 'templates_ind.npy'));

    %chanMap0ind = int32(chanMap0ind);
    chanMap0ind = int32([1:rez.ops.Nchan]-1);
    writeNPY(chanMap0ind, fullfile(savePath, 'channel_map.npy'));
    writeNPY([xcoords ycoords], fullfile(savePath, 'channel_positions.npy'));

    if ~isempty(templateFeatures)
        writeNPY(templateFeatures, fullfile(savePath, 'template_features.npy'));
        writeNPY(templateFeatureInds'-1, fullfile(savePath, 'template_feature_ind.npy'));% -1 for zero indexing
        writeNPY(pcFeatures, fullfile(savePath, 'pc_features.npy'));
        writeNPY(pcFeatureInds'-1, fullfile(savePath, 'pc_feature_ind.npy'));% -1 for zero indexing
    end

    writeNPY(whiteningMatrix, fullfile(savePath, 'whitening_mat.npy'));
    writeNPY(whiteningMatrixInv, fullfile(savePath, 'whitening_mat_inv.npy'));

    if isfield(rez, 'simScore')
        similarTemplates = rez.simScore;
        writeNPY(similarTemplates, fullfile(savePath, 'similar_templates.npy'));
    end

    % save a list of "good" clusters for Phy
%     fileID = fopen(fullfile(savePath, 'channel_names.tsv'), 'w');
%     fprintf(fileID, 'cluster_id%sKSLabel', char(9));
%     for j = 1:Nchan
%         fprintf(fileID, '%d%s%d', j-1,char(9),chanMap0ind(j));
%         fprintf(fileID, char([13 10]));
%     end
%     fclose(fileID);

    % Duplicate "KSLabel" as "group", a special metadata ID for Phy, so that
    % filtering works as expected in the cluster view
    KSLabelFilename = fullfile(savePath, 'cluster_KSLabel.tsv');
    copyfile(KSLabelFilename, fullfile(savePath, 'cluster_group.tsv'));

     %make params file
    if ~exist(fullfile(savePath,'params.py'),'file')
        fid = fopen(fullfile(savePath,'params.py'), 'w');

%        [~, fname, ext] = fileparts(rez.ops.fbinary);
%         fprintf(fid,['dat_path = ''',fname ext '''\n']);
%         fprintf(fid,'n_channels_dat = %i\n',rez.ops.NchanTOT);
        if ~isempty(varargin)
            [root, fname, ext] = fileparts(rez.ops.fbinary);
        else
            [root, fname, ext] = fileparts(rez.ops.fproc);
        end
%         fprintf(fid,['dat_path = ''',fname ext '''\n']);
        fprintf(fid,['dat_path = ''', strrep(rez.ops.fproc, '\', '/') '''\n']);

        fprintf(fid,'n_channels_dat = %i\n',rez.ops.Nchan);

        fprintf(fid,'dtype = ''int16''\n');
        fprintf(fid,'offset = 0\n');
        if mod(rez.ops.fs,1)
            fprintf(fid,'sample_rate = %i\n',rez.ops.fs);
        else
            fprintf(fid,'sample_rate = %i.\n',rez.ops.fs);
        end
%         fprintf(fid,'hp_filtered = False');
        fprintf(fid,'hp_filtered = True');

        fclose(fid);
    end
end

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort