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

%SOM_SEQTRAIN  Use sequential algorithm to train the Self-Organizing Map.
%
% [sM,sT] = som_seqtrain(sM, D, [[argID,] value, ...])
%
%  sM     = som_seqtrain(sM,D);
%  sM     = som_seqtrain(sM,sD,'alpha_type','power','tracking',3);
%  [M,sT] = som_seqtrain(M,D,'ep','trainlen',10,'inv','hexa');
%
%  Input and output arguments ([]'s are optional):
%   sM      (struct) map struct, the trained and updated map is returned
%           (matrix) codebook matrix of a self-organizing map
%                    size munits x dim or  msize(1) x ... x msize(k) x dim
%                    The trained map codebook is returned.
%   D       (struct) training data; data struct
%           (matrix) training data, size dlen x dim
%   [argID, (string) See below. The values which are unambiguous can
%    value] (varies) be given without the preceeding argID.
%
%   sT      (struct) learning parameters used during the training
%
% Here are the valid argument IDs and corresponding values. The values which
% are unambiguous (marked with '*') can be given without the preceeding argID.
%   'mask'        (vector) BMU search mask, size dim x 1
%   'msize'       (vector) map size
%   'radius'      (vector) neighborhood radiuses, length 1, 2 or trainlen
%   'radius_ini'  (scalar) initial training radius
%   'radius_fin'  (scalar) final training radius
%   'alpha'       (vector) learning rates, length trainlen
%   'alpha_ini'   (scalar) initial learning rate
%   'tracking'    (scalar) tracking level, 0-3
%   'trainlen'    (scalar) training length
%   'trainlen_type' *(string) is the given trainlen 'samples' or 'epochs'
%   'train'      *(struct) train struct, parameters for training
%   'sTrain','som_train '  = 'train'
%   'alpha_type' *(string) learning rate function, 'inv', 'linear' or 'power'
%   'sample_order'*(string) order of samples: 'random' or 'ordered'
%   'neigh'      *(string) neighborhood function, 'gaussian', 'cutgauss',
%                          'ep' or 'bubble'
%   'topol'      *(struct) topology struct
%   'som_topol','sTopo l'  = 'topol'
%   'lattice'    *(string) map lattice, 'hexa' or 'rect'
%   'shape'      *(string) map shape, 'sheet', 'cyl' or 'toroid'
%
% For more help, try 'type som_seqtrain' or check out online documentation.
% See also  SOM_MAKE, SOM_BATCHTRAIN, SOM_TRAIN_STRUCT.

%%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% som_seqtrain
%
% PURPOSE
%
% Trains a Self-Organizing Map using the sequential algorithm.
%
% SYNTAX
%
%  sM = som_seqtrain(sM,D);
%  sM = som_seqtrain(sM,sD);
%  sM = som_seqtrain(...,'argID',value,...);
%  sM = som_seqtrain(...,value,...);
%  [sM,sT] = som_seqtrain(M,D,...);
%
% DESCRIPTION
%
% Trains the given SOM (sM or M above) with the given training data
% (sD or D) using sequential SOM training algorithm. If no optional
% arguments (argID, value) are given, a default training is done, the
% parameters are obtained from SOM_TRAIN_STRUCT function. Using
% optional arguments the training parameters can be specified. Returns
% the trained and updated SOM and a train struct which contains
% information on the training.
%
% REFERENCES
%
% Kohonen, T., "Self-Organizing Map", 2nd ed., Springer-Verlag,
%    Berlin, 1995, pp. 78-82.
% Kohonen, T., "Clustering, Taxonomy, and Topological Maps of
%    Patterns", International Conference on Pattern Recognition
%    (ICPR), 1982, pp. 114-128.
% Kohonen, T., "Self-Organized formation of topologically correct
%    feature maps", Biological Cybernetics 43, 1982, pp. 59-69.
%
% REQUIRED INPUT ARGUMENTS
%
%  sM          The map to be trained.
%     (struct) map struct
%     (matrix) codebook matrix (field .data of map struct)
%              Size is either [munits dim], in which case the map grid
%              dimensions (msize) should be specified with optional arguments,
%              or [msize(1) ... msize(k) dim] in which case the map
%              grid dimensions are taken from the size of the matrix.
%              Lattice, by default, is 'rect' and shape 'sheet'.
%  D           Training data.
%     (struct) data struct
%     (matrix) data matrix, size [dlen dim]
%
% OPTIONAL INPUT ARGUMENTS
%
%  argID (string) Argument identifier string (see below).
%  value (varies) Value for the argument (see below).
%
%  The optional arguments can be given as 'argID',value -pairs. If an
%  argument is given value multiple times, the last one is
%  used. The valid IDs and corresponding values are listed below. The values
%  which are unambiguous (marked with '*') can be given without the
%  preceeding argID.
%
%   'mask'       (vector) BMU search mask, size dim x 1. Default is
%                         the one in sM (field '.mask') or a vector of
%                         ones if only a codebook matrix was given.
%   'msize'      (vector) map grid dimensions. Default is the one
%                         in sM (field sM.topol.msize) or
%                         'si = size(sM); msize = si(1:end-1);'
%                         if only a codebook matrix was given.
%   'radius'     (vector) neighborhood radius
%                         length = 1: radius_ini = radius
%                         length = 2: [radius_ini radius_fin] = radius
%                         length > 2: the vector given neighborhood
%                                     radius for each step separately
%                                     trainlen = length(radius)
%   'radius_ini' (scalar) initial training radius
%   'radius_fin' (scalar) final training radius
%   'alpha'      (vector) learning rate
%                         length = 1: alpha_ini = alpha
%                         length > 1: the vector gives learning rate
%                                     for each step separately
%                                     trainlen is set to length(alpha)
%                                     alpha_type is set to 'user defined'
%   'alpha_ini'  (scalar) initial learning rate
%   'tracking'   (scalar) tracking level: 0, 1 (default), 2 or 3
%                         0 - estimate time
%                         1 - track time and quantization error
%                         2 - plot quantization error
%                         3 - plot quantization error and two first
%                             components
%   'trainlen'   (scalar) training length (see also 'tlen_type')
%   'trainlen_type' *(string) is the trainlen argument given in 'epochs'
%                         or in 'samples'. Default is 'epochs'.
%   'sample_order'*(string) is the sample order 'random' (which is the
%                         the default) or 'ordered' in which case
%                         samples are taken in the order in which they
%                         appear in the data set
%   'train'     *(struct) train struct, parameters for training.
%                         Default parameters, unless specified,
%                         are acquired using SOM_TRAIN_STRUCT (this
%                         also applies for 'trainlen', 'alpha_type',
%                         'alpha_ini', 'radius_ini' and 'radius_fin').
%   'sTrain', 'som_train' (struct) = 'train'
%   'neigh'     *(string) The used neighborhood function. Default is
%                         the one in sM (field '.neigh') or 'gaussian'
%                         if only a codebook matrix was given. Other
%                         possible values is 'cutgauss', 'ep' and 'bubble'.
%   'topol'     *(struct) topology of the map. Default is the one
%                         in sM (field '.topol').
%   'sTopol', 'som_topol' (struct) = 'topol'
%   'alpha_type'*(string) learning rate function, 'inv', 'linear' or 'power'
%   'lattice'   *(string) map lattice. Default is the one in sM
%                         (field sM.topol.lattice) or 'rect'
%                         if only a codebook matrix was given.
%   'shape'     *(string) map shape. Default is the one in sM
%                         (field sM.topol.shape) or 'sheet'
%                         if only a codebook matrix was given.
%
% OUTPUT ARGUMENTS
%
%  sM          the trained map
%     (struct) if a map struct was given as input argument, a
%              map struct is also returned. The current training
%              is added to the training history (sM.trainhist).
%              The 'neigh' and 'mask' fields of the map struct
%              are updated to match those of the training.
%     (matrix) if a matrix was given as input argument, a matrix
%              is also returned with the same size as the input
%              argument.
%  sT (struct) train struct; information of the accomplished training
%
% EXAMPLES
%
% Simplest case:
%  sM = som_seqtrain(sM,D);
%  sM = som_seqtrain(sM,sD);
%
% To change the tracking level, 'tracking' argument is specified:
%  sM = som_seqtrain(sM,D,'tracking',3);
%
% The change training parameters, the optional arguments 'train',
% 'neigh','mask','trainlen','radius','radius_ini', 'radius_fin',
% 'alpha', 'alpha_type' and 'alpha_ini' are used.
%  sM = som_seqtrain(sM,D,'neigh','cutgauss','trainlen',10,'radius_fin',0);
%
% Another way to specify training parameters is to create a train struct:
%  sTrain = som_train_struct(sM,'dlen',size(D,1),'algorithm','seq');
%  sTrain = som_set(sTrain,'neigh','cutgauss');
%  sM = som_seqtrain(sM,D,sTrain);
%
% By default the neighborhood radius goes linearly from radius_ini to
% radius_fin. If you want to change this, you can use the 'radius' argument
% to specify the neighborhood radius for each step separately:
%  sM = som_seqtrain(sM,D,'radius',[5 3 1 1 1 1 0.5 0.5 0.5]);
%
% By default the learning rate (alpha) goes from the alpha_ini to 0
% along the function defined by alpha_type. If you want to change this,
% you can use the 'alpha' argument to specify the learning rate
% for each step separately:
%  alpha = 0.2*(1 - log([1:100]));
%  sM = som_seqtrain(sM,D,'alpha',alpha);
%
% You don't necessarily have to use the map struct, but you can operate
% directly with codebook matrices. However, in this case you have to
% specify the topology of the map in the optional arguments. The
% following commads are identical (M is originally a 200 x dim sized matrix):
%  M = som_seqtrain(M,D,'msize',[20 10],'lattice','hexa','shape','cyl');
%
%  M = som_seqtrain(M,D,'msize',[20 10],'hexa','cyl');
%
%  sT= som_set('som_topol','msize',[20 10],'lattice','hexa','shape','cyl');
%  M = som_seqtrain(M,D,sT);
%
%  M = reshape(M,[20 10 dim]);
%  M = som_seqtrain(M,D,'hexa','cyl');
%
% The som_seqtrain also returns a train struct with information on the
% accomplished training. This is the same one as is added to the end of the
% trainhist field of map struct, in case a map struct is given.
%  [M,sTrain] = som_seqtrain(M,D,'msize',[20 10]);
%
%  [sM,sTrain] = som_seqtrain(sM,D); % sM.trainhist{end}==sTrain
%
% SEE ALSO
%
%  som_make         Initialize and train a SOM using default parameters.
%  som_batchtrain   Train SOM with batch algorithm.
%  som_train_struct Determine default training parameters.

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

% Version 1.0beta juuso 220997
% Version 2.0beta juuso 101199

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Check arguments

error(nargchk(2, Inf, nargin));  % check the number of input arguments

% map
struct_mode = isstruct(sMap);
if struct_mode,
  sTopol = sMap.topol;
else
  orig_size = size(sMap);
  if ndims(sMap) > 2,
    si = size(sMap); dim = si(end); msize = si(1:end-1);
    M = reshape(sMap,[prod(msize) dim]);
  else
    msize = [orig_size(1) 1];
    dim = orig_size(2);
  end
  sMap   = som_map_struct(dim,'msize',msize);
  sTopol = sMap.topol;
end
[munits dim] = size(sMap.codebook);

% data
if isstruct(D),
  data_name = D.name;
  D = D.data;
else
  data_name = inputname(2);
end
D = D(find(sum(isnan(D),2) < dim),:); % remove empty vectors from the data
[dlen ddim] = size(D);                % check input dimension
if dim ~= ddim, error('Map and data input space dimensions disagree.'); end

% varargin
sTrain = som_set('som_train','algorithm','seq','neigh', ...
		 sMap.neigh,'mask',sMap.mask,'data_name',data_name);
radius     = [];
alpha      = [];
tracking   = 1;
sample_order_type = 'random';
tlen_type  = 'epochs';

i=1;
while i<=length(varargin),
  argok = 1;
  if ischar(varargin{i}),
    switch varargin{i},
     % argument IDs
     case 'msize', i=i+1; sTopol.msize = varargin{i};
     case 'lattice', i=i+1; sTopol.lattice = varargin{i};
     case 'shape', i=i+1; sTopol.shape = varargin{i};
     case 'mask', i=i+1; sTrain.mask = varargin{i};
     case 'neigh', i=i+1; sTrain.neigh = varargin{i};
     case 'trainlen', i=i+1; sTrain.trainlen = varargin{i};
     case 'trainlen_type', i=i+1; tlen_type = varargin{i};
     case 'tracking', i=i+1; tracking = varargin{i};
     case 'sample_order', i=i+1; sample_order_type = varargin{i};
     case 'radius_ini', i=i+1; sTrain.radius_ini = varargin{i};
     case 'radius_fin', i=i+1; sTrain.radius_fin = varargin{i};
     case 'radius',
      i=i+1;
      l = length(varargin{i});
      if l==1,
        sTrain.radius_ini = varargin{i};
      else
        sTrain.radius_ini = varargin{i}(1);
        sTrain.radius_fin = varargin{i}(end);
        if l>2, radius = varargin{i}; tlen_type = 'samples'; end
      end
     case 'alpha_type', i=i+1; sTrain.alpha_type = varargin{i};
     case 'alpha_ini', i=i+1; sTrain.alpha_ini = varargin{i};
     case 'alpha',
      i=i+1;
      sTrain.alpha_ini = varargin{i}(1);
      if length(varargin{i})>1,
        alpha = varargin{i}; tlen_type = 'samples';
        sTrain.alpha_type = 'user defined';
      end
     case {'sTrain','train','som_train'}, i=i+1; sTrain = varargin{i};
     case {'topol','sTopol','som_topol'},
      i=i+1;
      sTopol = varargin{i};
      if prod(sTopol.msize) ~= munits,
        error('Given map grid size does not match the codebook size.');
      end
      % unambiguous values
     case {'inv','linear','power'}, sTrain.alpha_type = varargin{i};
     case {'hexa','rect'}, sTopol.lattice = varargin{i};
     case {'sheet','cyl','toroid'}, sTopol.shape = varargin{i};
     case {'gaussian','cutgauss','ep','bubble'}, sTrain.neigh = varargin{i};
     case {'epochs','samples'}, tlen_type = varargin{i};
     case {'random', 'ordered'}, sample_order_type = 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};
      if prod(sTopol.msize) ~= munits,
        error('Given map grid size does not match the codebook size.');
      end
     case 'som_train', sTrain = varargin{i};
     otherwise argok=0;
    end
  else
    argok = 0;
  end
  if ~argok,
    disp(['(som_seqtrain) Ignoring invalid argument #' num2str(i+2)]);
  end
  i = i+1;
end

% training length
if ~isempty(radius) || ~isempty(alpha),
  lr = length(radius);
  la = length(alpha);
  if lr>2 || la>1,
    tlen_type = 'samples';
    if     lr> 2 && la<=1, sTrain.trainlen = lr;
    elseif lr<=2 && la> 1, sTrain.trainlen = la;
    elseif lr==la,        sTrain.trainlen = la;
    else
      error('Mismatch between radius and learning rate vector lengths.')
    end
  end
end
if strcmp(tlen_type,'samples'), sTrain.trainlen = sTrain.trainlen/dlen; end

% check topology
if struct_mode,
  if ~strcmp(sTopol.lattice,sMap.topol.lattice) || ...
	~strcmp(sTopol.shape,sMap.topol.shape) || ...
	any(sTopol.msize ~= sMap.topol.msize),
    warning('Changing the original map topology.');
  end
end
sMap.topol = sTopol;

% complement the training struct
sTrain = som_train_struct(sTrain,sMap,'dlen',dlen);
if isempty(sTrain.mask), sTrain.mask = ones(dim,1); end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% initialize

M        = sMap.codebook;
mask     = sTrain.mask;
trainlen = sTrain.trainlen*dlen;

% neighborhood radius
if length(radius)>2,
  radius_type = 'user defined';
else
  radius = [sTrain.radius_ini sTrain.radius_fin];
  rini = radius(1);
  rstep = (radius(end)-radius(1))/(trainlen-1);
  radius_type = 'linear';
end

% learning rate
if length(alpha)>1,
  sTrain.alpha_type ='user defined';
  if ~(abs(length(alpha)-trainlen) < ...
          1e4*eps(min(abs(length(alpha)),abs(trainlen)))),
    error('Trainlen and length of neighborhood radius vector do not match.')
  end
  if any(isnan(alpha)),
    error('NaN is an illegal learning rate.')
  end
else
  if isempty(alpha), alpha = sTrain.alpha_ini; end
  if strcmp(sTrain.alpha_type,'inv'),
    % alpha(t) = a / (t+b), where a and b are chosen suitably
    % below, they are chosen so that alpha_fin = alpha_ini/100
    b = (trainlen - 1) / (100 - 1);
    a = b * alpha;
  end
end

% initialize random number generator
rand('state',sum(100*clock));

% distance between map units in the output space
%  Since in the case of gaussian and ep neighborhood functions, the
%  equations utilize squares of the unit distances and in bubble case
%  it doesn't matter which is used, the unitdistances and neighborhood
%  radiuses are squared.
Ud = som_unit_dists(sTopol).^2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Action

update_step = 100;
mu_x_1 = ones(munits,1);
samples = ones(update_step,1);
r = samples;
alfa = samples;

qe = 0;
start = clock;
if tracking >  0, % initialize tracking
  track_table = zeros(update_step,1);
  qe = zeros(floor(trainlen/update_step),1);
end

printedbytes = 0;

for t = 1:trainlen,

  % Every update_step, new values for sample indeces, neighborhood
  % radius and learning rate are calculated. This could be done
  % every step, but this way it is more efficient. Or this could
  % be done all at once outside the loop, but it would require much
  % more memory.
  ind = rem(t,update_step); if ind==0, ind = update_step; end
  if ind==1,
    steps = [t:min(trainlen,t+update_step-1)];
    % sample order
    switch sample_order_type,
     case 'ordered', samples = rem(steps,dlen)+1;
     case 'random',  samples = ceil(dlen*rand(update_step,1)+eps);
    end

    % neighborhood radius
    switch radius_type,
     case 'linear',       r = rini+(steps-1)*rstep;
     case 'user defined', r = radius(steps);
    end
    r=r.^2;        % squared radius (see notes about Ud above)
    r(r==0) = eps; % zero radius might cause div-by-zero error

    % learning rate
    switch sTrain.alpha_type,
     case 'linear',       alfa = (1-steps/trainlen)*alpha;
     case 'inv',          alfa = a ./ (b + steps-1);
     case 'power',        alfa = alpha * (0.005/alpha).^((steps-1)/trainlen);
     case 'user defined', alfa = alpha(steps);
    end
  end

  % find BMU
  x = D(samples(ind),:);                 % pick one sample vector
  known = ~isnan(x);                     % its known components
  Dx = M(:,known) - x(mu_x_1,known);     % each map unit minus the vector
  [qerr bmu] = min((Dx.^2)*mask(known)); % minimum distance(^2) and the BMU

  % tracking
  if tracking>0,
    track_table(ind) = sqrt(qerr);
    if ind==update_step,
      n = ceil(t/update_step);
      qe(n) = mean(track_table);
      printedbytes = trackplot(M,D,tracking,start,n,qe,printedbytes);
    end
  end

  % neighborhood & learning rate
  % notice that the elements Ud and radius have been squared!
  % (see notes about Ud above)
  switch sTrain.neigh,
  case 'bubble',   h = (Ud(:,bmu)<=r(ind));
  case 'gaussian', h = exp(-Ud(:,bmu)/(2*r(ind)));
  case 'cutgauss', h = exp(-Ud(:,bmu)/(2*r(ind))) .* (Ud(:,bmu)<=r(ind));
  case 'ep',       h = (1-Ud(:,bmu)/r(ind)) .* (Ud(:,bmu)<=r(ind));
  end
  h = h*alfa(ind);

  % update M
  M(:,known) = M(:,known) - h(:,ones(sum(known),1)).*Dx;

end; % for t = 1:trainlen

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Build / clean up the return arguments

if tracking, fprintf(1,'\n'); end

% update structures
sTrain = som_set(sTrain,'time',datestr(now,0));
if struct_mode,
  sMap = som_set(sMap,'codebook',M,'mask',sTrain.mask,'neigh',sTrain.neigh);
  tl = length(sMap.trainhist);
  sMap.trainhist(tl+1) = sTrain;
else
  sMap = reshape(M,orig_size);
end

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% subfunctions

%%%%%%%%
function [count] = trackplot(M,D,tracking,start,n,qe, printedbytes)

  l = length(qe);
  elap_t = etime(clock,start);
  tot_t = elap_t*l/n;
  % Carriage return does not work as it should (even on UNIX) when printing
  % to screen, so let's do this instead
  fprintf(1, repmat('\b', 1, printedbytes));
  count = fprintf(1,'Training: %3.0f/ %3.0f s',elap_t,tot_t);
  switch tracking
   case 1,
   case 2,
    plot(1:n,qe(1:n),(n+1):l,qe((n+1):l))
    title('Quantization errors for latest samples')
    drawnow
   otherwise,
    subplot(2,1,1), plot(1:n,qe(1:n),(n+1):l,qe((n+1):l))
    title('Quantization error for latest samples');
    subplot(2,1,2), plot(M(:,1),M(:,2),'ro',D(:,1),D(:,2),'b.');
    title('First two components of map units (o) and data vectors (+)');
    drawnow
  end
  % end of trackplot

author: ilarinieminen

repo: SOM-Toolbox


text:
function LU = getbounds(F)

K.f = 0;
K.l = 0;
L = [];
U = [];
LU =  yalmip('getbounds',1:yalmip('nvars'));
F = flatten(F);
for i = 1:length(F.clauses)
    if F.clauses{i}.type == 2
        X = F.clauses{i}.data;
        AB = getbase(X);
        K.l = prod(size(X));
        variables = getvariables(X);
        [lb,ub,cand_rows] = findulb(AB,K);
        LU(variables,1) = max([lb LU(variables,1)]')';
        LU(variables,2) = min([ub LU(variables,2)]')';
    elseif F.clauses{i}.type == 3
        % FIX : Extract from equalities and binary constraints
    end
end

binary = yalmip('binvariables');
LU(binary,1) = 0;
LU(binary,2) = 1;

% Try to bound some nonlinear terms
% FIX: complete code
[mt,variable_type] = yalmip('monomtable');
quadratic = find(variable_type == 2);
if ~isempty(quadratic)
    M = mt(quadratic,:);
    for i = 1:size(M,1)
        [ii,jj] = find(M(i,:));
        if length(ii) == 1
            LU(quadratic(i),1) = min([0 LU(jj,1)^2]);
            LU(quadratic(i),2) = max([LU(jj,1)^2 LU(jj,2)^2]);
        end
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [y]=axpx(a,x)
%[y]=AXPX(A,X)
%This is a technical function that multiplies
%Matrix A by a vector and projects on X itself;
%This is for time-stepping

d=x.d;
n=a.n;
m=a.m;
rx=x.r;
psx=x.ps;
att=a.tt;
corea=att.core;
psa=att.ps;
ra=att.r;
%Do right-to-left orthogonalization and computation of psi
%matrices
corex=x.core;
pos1=psx(d+1);
psi=cell(d+1,1); %Psi-matrices
psi{d+1}=1; psi{1}=1;
cr1=corex(psx(d):psx(d+1)-1);
y=x;
corey=y.core;
psy=y.ps;
ry=y.r;
for i=d:-1:2
   %  fprintf('i=%d \n');
   cr2=corey(psy(i-1):psy(i)-1);
   cr1=reshape(cr1,[ry(i),n(i)*ry(i+1)]);
   cr2=reshape(cr2,[ry(i-1)*n(i-1),ry(i)]);
   cr1=cr1.';
   [q,rm]=qr(cr1,0); rn=size(q,2); rm=rm.';
   q=q.';
   ry(i)=rn;
   corey(pos1-ry(i+1)*n(i)*ry(i):pos1-1)=q(:);
      %Convolution is now performed for psi(i) using psi(i+1) and corea, and
   %(new) core q
   cra=corea(psa(i):psa(i+1)-1); cra=reshape(cra,[ra(i),n(i),m(i),ra(i+1)]);
   cry=reshape(conj(q),[ry(i),n(i),ry(i+1)]);
   crx=corex(psx(i):psx(i+1)-1);
   crx=reshape(crx,[rx(i),m(i),rx(i+1)]);
   pscur=psi{i+1}; pscur=reshape(pscur,[ra(i+1),rx(i+1),ry(i+1)]); %ra,rx,ry
   %First, convolve over rx(i+1)
   crx=reshape(crx,[rx(i)*m(i),rx(i+1)]);
    pscur=permute(pscur,[2,1,3]); pscur=reshape(pscur,[rx(i+1),ra(i+1)*ry(i+1)]);
   pscur=crx*pscur; %pscur is now rx(i)*m(i)*ra(i+1)*ry(i+1)
   %Convolve over m(i),ra(i+1),n(i),ry(i+1)
    pscur=reshape(pscur,[rx(i),m(i)*ra(i+1),ry(i+1)]);
    pscur=permute(pscur,[1,3,2]);
    pscur=reshape(pscur,[rx(i)*ry(i+1),m(i)*ra(i+1)]);
    cra=reshape(cra,[ra(i)*n(i),m(i)*ra(i+1)]); cra=cra.';
    pscur=pscur*cra;
    %pscur is now rx(i)*ry(i+1)*ra(i)*n(i), it is left to convolve over
    %n(i)*ry(i+1)
    pscur=reshape(pscur,[rx(i),ry(i+1),ra(i),n(i)]);
    pscur=permute(pscur,[3,1,4,2]);
    pscur=reshape(pscur,[ra(i)*rx(i),n(i)*ry(i+1)]);
    cry=reshape(cry,[ry(i),n(i)*ry(i+1)]); cry=cry.';
    pscur=pscur*cry; %pscur is ra*rx*ry
    pscur=reshape(pscur,[rx(i),ra(i),ry(i)]);
    psi{i}=pscur;
   %End of psi-block
   pos1=pos1-ry(i+1)*n(i)*ry(i);
   cr1=cr2*rm;
end
corey(pos1-ry(2)*n(1)*ry(1):pos1-1)=cr1(:);
pos1=pos1-ry(2)*n(1)*ry(1);
corey=corey(pos1:end); %Truncate unused elements
%Now compute the projection itself while recomputing psi-matrices
 pos1=1;
 psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
 rm=1;
 %keyboard;
  for i=1:d
          %Our convolution is
     %ps1(ra(i),rx(i),ry(i))*cra(ra(i),n(i),m(i),ra(i+1))
     %*ps2(ra(i+1),rx(i+1),ry(i+1))*crx(rx(i),m(i)*rx(i+1))->
     %cry(ry(i),n(i),ry(i+1)

     ps1=psi{i}; ps2=psi{i+1};
     cra=corea(psa(i):psa(i+1)-1);
     crx=corex(psx(i):psx(i+1)-1);
%     corey(psy(i):psy(i+1)-1)=crx(:);
%     y.core=corey;
%          psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
%     y.ps=psy;
     ps1=reshape(ps1,[ra(i),rx(i),ry(i)]);
     ps1=permute(ps1,[2,3,1]);
     ps1=reshape(ps1,[rx(i)*ry(i),ra(i)]);
     cra=reshape(cra,[ra(i),n(i)*m(i)*ra(i+1)]);
     cry=ps1*cra;      %cry is rx(i)*ry(i)*n(i)*m(i)*ra(i+1)
     %convolve over m(i),rx(i) with crx
     cry=reshape(cry,[rx(i),ry(i),n(i),m(i),ra(i+1)]);
     cry=permute(cry,[2,3,5,1,4]);
     cry=reshape(cry,[ry(i)*n(i)*ra(i+1),rx(i)*m(i)]);
     crx=reshape(crx,[rx(i)*m(i),rx(i+1)]);
     cry=cry*crx; cry0=cry;
     %cry is ry(i)*n(i)*ra(i+1)*rx(i+1)
     cry=reshape(cry,[ry(i)*n(i),ra(i+1)*rx(i+1)]);
     ps2=reshape(ps2,[ra(i+1)*rx(i+1),ry(i+1)]);
     %cry0=cry; %cry0 is ry(i)*n(i)*ra(i+1)*rx(i+1), convolution over ry(i)*n(i) will bring new psi --- not true; new psi will be brought with and
     %orthogonalization
     cry=cry*ps2;
      corey(pos1:pos1+ry(i)*n(i)*ry(i+1)-1)=cry(:);
      y.core=corey;
      psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
      y.ps=psy;
       %keyboard;

     cry=reshape(cry,[ry(i)*n(i),ry(i+1)]);
     [q,rm]=qr(cry,0);
     rn=size(q,2);
     ry(i+1)=rn;
     %q=q*rm;
     corey(pos1:pos1+ry(i)*n(i)*ry(i+1)-1)=q(:);
     pos1=pos1+ry(i)*n(i)*ry(i+1);
     %And now --- orthogonalization! (of cry) (he-he)
     %Now compute "new" psi
     %We need to compute
          %cry is rx(i)*ry(i)*n(i)*m(i)*ra(i+1)
     %convolve over m(i),rx(i) with crx
     %cry0 is ry(i)*n(i)*ra(i+1)*rx(i+1);
     cry0=reshape(cry0,[ry(i),n(i),ra(i+1),rx(i+1)]);
     cry0=permute(cry0,[3,4,1,2]);
     cry0=reshape(cry0,[ra(i+1)*rx(i+1),ry(i)*n(i)]);
     q=reshape(conj(q),[ry(i)*n(i),ry(i+1)]);
     psi{i+1}=cry0*q;

  end
  psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
  corey(psy(i):psy(i+1)-1)=corey(psy(i):psy(i+1)-1)*rm;
  y.core=corey;
  y.r=ry;
  y.ps=psy;
  y.n=n;
  y.d=d;
  return
end

author: oseledets

repo: TT-Toolbox


text:
function misclassified = ml_get_misclassified(X,labels, f, options)
%ML_GET_MISCLASSIFIED
%
%   input -----------------------------------------------------------------
%
%       o  X        : (N x D), dataset of N datapoints of dimension D.
%
%       o labels    : (N x 1), ground truth class labels.
%
%       o f         : function handle, classifier f.
%
%               - y = f(X); y are predicted class labels.
%
%   output ----------------------------------------------------------------
%
%       o misclassified : struct,
%
%           - misclassified.index :     indicies of data points X which have
%                                       been misclassified.
%
%           - misclassified.predicted: the predicted class label (which was wrong).
%
%           - misclassified.truth:     class label which should have been
%                                      predicted.

class_id        = unique(labels);
num_classes     = length(class_id);
M               = zeros(num_classes,num_classes);
dim_swaped      = false; %false: (MXD) , true: (DXM)
if isfield(options,'dim_swaped'),       dim_swaped   = options.dim_swaped;   end

% predicted class label
if dim_swaped
    hc          = f(X');
else
    hc          = f(X);
end

% entries which are 0, mean they where correctly classified, otherwise not.
tmp         = hc(:) - labels(:);
idx         = find(tmp ~= 0);

misclassified.idx           = idx;
misclassified.predicted     = hc(idx);
misclassified.truth         = labels(idx);

end

author: epfl-lasa


text:
function res = smoothPolygon(poly, M)
%SMOOTHPOLYGON Smooth a polygon using local averaging.
%
%   RES = smoothPolygon(POLY, M)
%   POLY contains the polygon vertices, and M is the size of smoothing
%   (given as the length of the convolution window).
%
%
%   Example
%     img = imread('circles.png');
%     img = imfill(img, 'holes');
%     contours = bwboundaries(img');
%     contour = contours{1};
%     imshow(img); hold on; drawPolygon(contour, 'b');
%     contourf = smoothPolygon(contour, 11);
%     drawPolygon(contourf, 'm');
%
%   See also
%     polygons2d, smoothPolyline, simplifyPolygon, resamplePolygon

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2015-02-17, using Matlab 8.4.0.150421 (R2014b)
% Copyright 2015-2022 INRA - Cepia Software Platform

% compute the number of elements before and after
M1 = floor((M - 1) / 2);
M2 = ceil((M - 1) / 2);

% repeat beginning and end of contour
poly2 = [poly(end-M1+1:end, :) ; poly ; poly(1:M2,:)];

% create convolution vector
v2 = ones(M, 1) / M;

% apply contour filtering
res(:,1) = conv(poly2(:,1), v2, 'same');
res(:,2) = conv(poly2(:,2), v2, 'same');

% keep the interesting part
res = res(M1+1:end-M2, :);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [m,v,w,g,f,pp,gg]=v_gaussmix(x,c,l,m0,v0,w0,wx)
%V_GAUSSMIX fits a gaussian mixture pdf to a set of data observations [m,v,w,g,f]=(x,c,l,m0,v0,w0,wx)
%
% Usage:
%    (1) [m,v,w]=v_gaussmix(x,[],[],k);    % create GMM with k mixtures and diagonal covariances
%    (2) [m,v,w]=gaussmix(x,[],[],k,'v');    % create GMM with k mixtures and full covariances
%
% Inputs: n data values, k mixtures, p parameters, l loops
%
%     X(n,p)   Input data vectors, one per row.
%     C(1)     Minimum variance of normalized data (Use [] to take default value of 1/n^2)
%     L        The integer portion of l gives a maximum loop count. The fractional portion gives
%              an optional stopping threshold. Iteration will cease if the increase in
%              log likelihood density per data point is less than this value. Thus l=10.001 will
%              stop after 10 iterations or when the increase in log likelihood falls below
%              0.001.
%              As a special case, if L=0, then the first three outputs are omitted.
%              Use [] to take default value of 100.0001
%     M0       Number of mixtures required (or initial mixture means - see below)
%     V0       Initialization mode:
%                'f'    Initialize with K randomly selected data points [default]
%                'p'    Initialize with centroids and variances of random partitions
%                'k'    k-means algorithm ('kf' and 'kp' determine initialization)
%                'h'    k-harmonic means algorithm ('hf' and 'hp' determine initialization) [default]
%                's'    use unscaled data during initialization phase instead of scaling it first
%                'm'    M0 contains the initial centres
%                'v'    full covariance matrices
%              Mode 'hf' [the default] generally gives the best results but 'f' is faster and often OK
%     W0(n,1)  Data point weights
%
%     Alternatively, initial values for M0, V0 and W0 can be given  explicitly:
%
%     M0(k,p)  Initial mixture means, one row per mixture.
%     V0(k,p)  Initial mixture variances, one row per mixture.
%      or V0(p,p,k)  one full-covariance matrix per mixture
%     W0(k,1)  Initial mixture weights, one per mixture. The weights should sum to unity.
%     WX(n,1)  Data point weights
%
% Outputs: (Note that M, V and W are omitted if L==0)
%
%     M(k,p)   Mixture means, one row per mixture. (omitted if L==0)
%     V(k,p)   Mixture variances, one row per mixture. (omitted if L==0)
%       or V(p,p,k) if full covariance matrices in use (i.e. either 'v' option or V0(p,p,k) specified)
%     W(k,1)   Mixture weights, one per mixture. The weights will sum to unity. (omitted if L==0)
%     G       Average log probability of the input data points.
%     F        Fisher's Discriminant measures how well the data divides into classes.
%              It is the ratio of the between-mixture variance to the average mixture variance: a
%              high value means the classes (mixtures) are well separated.
%     PP(n,1)  Log probability of each data point
%     GG(l+1,1) Average log probabilities at the beginning of each iteration and at the end
%
% The fitting procedure uses one of several initialization methods to create an initial guess
% for the mixture centres and then uses the EM (expectation-maximization) algorithm to refine
% the guess. Although the EM algorithm is deterministic, the initialization procedures use
% random numbers and so the routine will not give identical answers if you call it multiple
% times with the same input data. See v_randvec() for generating GMM data vectors.

%  Bugs/Suggestions
%     (1) Allow processing in chunks by outputting/reinputting an array of sufficient statistics
%     (2) Other initialization options:
%              'l'    LBG algorithm
%              'm'    Move-means (dog-rabbit) algorithm
%     (3) Allow updating of weights-only, not means/variances
%     (4) Allow freezing of means and/or variances

%      Copyright (C) Mike Brookes 2000-2009
%      Version: $Id: v_gaussmix.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[n,p]=size(x); % n = number of training values, p = dimension of data vector
wn=ones(n,1);
memsize=v_voicebox('memsize');    % set memory size to use
if isempty(c)
    c=1/n^2;
else
    c=c(1);         % just to prevent legacy code failing
end
fulliv=0;           % initial variance is not full
if isempty(l)
    l=100+1e-4;         % max loop count + stopping threshold
end
if nargin<5 || isempty(v0) || ischar(v0)             % no initial values specified for m0, v0, w0
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %  No initialvalues given, so we must use k-means or equivalent
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if nargin<6
        if nargin<5 || isempty(v0)
            v0='hf';                 % default initialization mode: hf
        end
        wx=wn;                      % no data point weights
    else
        wx=w0(:);                   % data point weights
    end
    wx=wx/sum(wx);
    if any(v0=='m')
        k=size(m0,1);
    else
        k=m0;
    end
    fv=any(v0=='v');                % full covariance matrices requested
    mx0=wx'*x;         % calculate mean of input data in each dimension
    vx0=wx'*x.^2-mx0.^2; % calculate variance of input data in each dimension
    sx0=sqrt(vx0);
    sx0(sx0==0)=1;      % do not divide by zero when scaling
    if n<=k                         % each data point can have its own mixture
        xs=(x-mx0(wn,:))./sx0(wn,:);          % scale the data
        m=xs(mod((1:k)-1,n)+1,:);   % just include all points several times
        v=zeros(k,p);               % will be set to floor later
        w=zeros(k,1);
        w(1:n)=1/n;
        if l>0
            l=0.1;                  % no point in iterating
        end
    else                            % more points than mixtures
        if any(v0=='s')
            xs=x;                   % do not scale data during initialization
        else
            xs=(x-mx0(wn,:))./sx0(wn,:);  % else scale now
            if any(v0=='m')
                m=(m0-mx0(ones(k,1),:))./sx0(ones(k,1),:);  % scale specified means as well
            end
        end
        w=repmat(1/k,k,1);                      % all mixtures equally likely
        if any(v0=='k')                         % k-means initialization
            if any(v0=='m')
                [m,e,j]=v_kmeans(xs,k,m);
            elseif any(v0=='p')
                [m,e,j]=v_kmeans(xs,k,'p');
            else
                [m,e,j]=v_kmeans(xs,k,'f');
            end
        elseif any(v0=='h')                     % k-harmonic means initialization
            if any(v0=='m')
                [m,e,j]=v_kmeanhar(xs,k,[],4,m);
            else
                if any(v0=='p')
                    [m,e,j]=v_kmeanhar(xs,k,[],4,'p');
                else
                    [m,e,j]=v_kmeanhar(xs,k,[],4,'f');
                end
            end
        elseif any(v0=='p')                     % Initialize using a random partition
            j=ceil(rand(n,1)*k);                % allocate to random clusters
            j(v_rnsubset(k,n))=1:k;               % but force at least one point per cluster
            for i=1:k
                m(i,:)=mean(xs(j==i,:),1);
            end
        else
            if any(v0=='m')
                m=m0;                           % use specified centres
            else
                m=xs(v_rnsubset(k,n),:);          % Forgy initialization: sample k centres without replacement [default]
            end
            [e,j]=v_kmeans(xs,k,m,0);             % find out the cluster allocation
        end
        if any(v0=='s')
            xs=(x-mx0(wn,:))./sx0(wn,:);      % scale data now if not done previously
        end
        v=zeros(k,p);                   % diagonal covariances
        w=zeros(k,1);
        for i=1:k
            ni=sum(j==i);               % number assigned to this centre
            w(i)=(ni+1)/(n+k);          % weight of this mixture
            if ni
                v(i,:)=sum((xs(j==i,:)-repmat(m(i,:),ni,1)).^2,1)/ni;
            else
                v(i,:)=zeros(1,p);
            end
        end
    end
else
    %%%%%%%%%%%%%%%%%%%%%%%%
    % use initial values given as input parameters
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if nargin<7
        wx=wn;              % no data point weights
    end
    wx=wx(:)/sum(wx); % normalize weights and force a column vector
    mx0=wx'*x;         % calculate mean of input data in each dimension
    vx0=wx'*x.^2-mx0.^2; % calculate variance of input data in each dimension
    sx0=sqrt(vx0);
    sx0(sx0==0)=1;      % do not divide by zero when scaling
    [k,p]=size(m0);
    xs=(x-mx0(wn,:))./sx0(wn,:);          % scale the data
    m=(m0-mx0(ones(k,1),:))./sx0(ones(k,1),:);          % and the means
    v=v0;
    w=w0;
    fv=ndims(v)>2 || size(v,1)>k;                       % full covariance matrix is supplied
    if fv
        mk=eye(p)==0;                                    % off-diagonal elements
        fulliv=any(v(repmat(mk,[1 1 k]))~=0);            % check if any are non-zero
        if ~fulliv
            v=reshape(v(repmat(~mk,[1 1 k])),p,k)'./repmat(sx0.^2,k,1);   % just pick out and scale the diagonal elements for now
        else
            v=v./repmat(sx0'*sx0,[1 1 k]);              % scale the full covariance matrix
        end
    end
end
if length(wx)~=n
    error('%d datapoints but %d weights',n,length(wx));
end
lsx=sum(log(sx0));
xsw=xs.*repmat(wx,1,p); % weighted data points
if ~fulliv          % initializing with diagonal covariance
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Diagonal Covariance matrices  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    v=max(v,c);         % apply the lower bound
    xs2=xs.^2.*repmat(wx,1,p);          % square and weight the data for variance calculations

    % If data size is large then do calculations in chunks

    nb=min(n,max(1,floor(memsize/(8*p*k))));    % chunk size for testing data points
    nl=ceil(n/nb);                  % number of chunks
    jx0=n-(nl-1)*nb;                % size of first chunk

    im=repmat(1:k,1,nb); im=im(:);
    th=(l-floor(l))*n;
    sd=(nargout > 3*(l~=0)); % = 1 if we are outputting log likelihood values
    lp=floor(l)+sd;   % extra loop needed to calculate final G value

    lpx=zeros(1,n);             % log probability of each data point
    wk=ones(k,1);
    wp=ones(1,p);
    wnb=ones(1,nb);
    wnj=ones(1,jx0);

    % EM loop

    g=0;                            % dummy initial value for comparison
    gg=zeros(lp+1,1);
    ss=sd;                          % initialize stopping count (0 or 1)
    for j=1:lp
        g1=g;                       % save previous log likelihood (2*pi factor omitted)
        m1=m;                       % save previous means, variances and weights
        v1=v;
        w1=w;
        vi=-0.5*v.^(-1);                % data-independent scale factor in exponent
        lvm=log(w)-0.5*sum(log(v),2);   % log of external scale factor (excluding -0.5*p*log(2pi) term)

        % first do partial chunk (of length jx0)

        jx=jx0;
        ii=1:jx;                        % indices of data points in this chunk
        kk=repmat(ii,k,1);              % kk(jx,k): one row per data point, one column per mixture
        km=repmat(1:k,1,jx);            % km(jx,k): one row per data point, one column per mixture
        py=reshape(sum((xs(kk(:),:)-m(km(:),:)).^2.*vi(km(:),:),2),k,jx)+lvm(:,wnj); % py(k,jx) pdf of each point with each mixture
        mx=max(py,[],1);                % mx(1,jx) find normalizing factor for each data point to prevent underflow when using exp()
        px=exp(py-mx(wk,:));            % find normalized probability of each mixture for each datapoint
        ps=sum(px,1);                   % total normalized likelihood of each data point
        px=px./ps(wk,:);                % relative mixture probabilities for each data point (columns sum to 1)
        lpx(ii)=log(ps)+mx;
        pk=px*wx(ii);                   % pk(k,1) effective fraction of data points for each mixture (could be zero due to underflow)
        sx=px*xsw(ii,:);
        sx2=px*xs2(ii,:);
        for il=2:nl                     % process the data points in chunks
            ix=jx+1;
            jx=jx+nb;                   % increment upper limit
            ii=ix:jx;                   % indices of data points in this chunk
            kk=repmat(ii,k,1);
            py=reshape(sum((xs(kk(:),:)-m(im,:)).^2.*vi(im,:),2),k,nb)+lvm(:,wnb);
            mx=max(py,[],1);            % find normalizing factor for each data point to prevent underflow when using exp()
            px=exp(py-mx(wk,:));        % find normalized probability of each mixture for each datapoint
            ps=sum(px,1);               % total normalized likelihood of each data point
            px=px./ps(wk,:);            % relative mixture probabilities for each data point (columns sum to 1)
            lpx(ii)=log(ps)+mx;
            pk=pk+px*wx(ii);            % pk(k,1) effective fraction of data points for each mixture (could be zero due to underflow)
            sx=sx+px*xsw(ii,:);
            sx2=sx2+px*xs2(ii,:);
        end
        g=lpx*wx;                       % total log probability summed over all data points
        gg(j)=g;                        % save log prob at each iteration
        w=pk;                           % normalize to get the weights
        if pk                           % if all elements of pk are non-zero
            m=sx./pk(:,wp);             % calculate mixture means
            v=sx2./pk(:,wp);            % and variances
        else
            wm=pk==0;                   % mask indicating mixtures with zero weights
            nz=sum(wm);              	% number of zero-weight mixtures
            [vv,mk]=sort(lpx);          % find the lowest probability data points
            m=zeros(k,p);               % initialize means and variances to zero (variances are floored later)
            v=m;
            m(wm,:)=xs(mk(1:nz),:); 	% set zero-weight mixture means to worst-fitted data points
            w(wm)=1/n;               	% set these weights non-zero
            w=w*n/(n+nz);            	% normalize so the weights sum to unity
            wm=~wm;                 	% mask for non-zero weights
            m(wm,:)=sx(wm,:)./pk(wm,wp);  % recalculate means and variances for mixtures with a non-zero weight
            v(wm,:)=sx2(wm,:)./pk(wm,wp);
        end
        v=max(v-m.^2,c);                % apply floor to variances
        if g-g1<=th && j>1
            if ~ss, break; end  %  stop
            ss=ss-1;       % stop next time
        end
    end
    if sd && ~fv  % we need to calculate the final probabilities
        pp=lpx'-0.5*p*log(2*pi)-lsx;   % log of total probability of each data point
        gg=gg(1:j)-0.5*p*log(2*pi)-lsx;    % average log prob at each iteration
        g=gg(end);
        m=m1;       % back up to previous iteration
        v=v1;
        w=w1;
        mm=sum(m,1)/k;
        f=(m(:)'*m(:)-k*mm(:)'*mm(:))/sum(v(:));
    end
    if ~fv
        m=m.*sx0(ones(k,1),:)+mx0(ones(k,1),:);	% unscale means
        v=v.*repmat(sx0.^2,k,1);                % and variances
    else
        v1=v;
        v=zeros(p,p,k);
        mk=eye(p)==1;                           % mask for diagonal elements
        v(repmat(mk,[1 1 k]))=v1';              % set from v1
    end
end
if fv              % check if full covariance matrices were requested
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Full Covariance matrices  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    pl=p*(p+1)/2;
    lix=1:p^2;
    cix=repmat(1:p,p,1);
    rix=cix';
    lix(cix>rix)=[];                                        % index of lower triangular elements
    cix=cix(lix);                                           % index of lower triangular columns
    rix=rix(lix);                                           % index of lower triangular rows
    dix=find(rix==cix);
    lixi=zeros(p,p);
    lixi(lix)=1:pl;
    lixi=lixi';
    lixi(lix)=1:pl;                                        % reverse index to build full matrices
    v=reshape(v,p^2,k);
    v=v(lix,:)';                                            % lower triangular in rows

    % If data size is large then do calculations in chunks

    nb=min(n,max(1,floor(memsize/(24*p*k))));    % chunk size for testing data points
    nl=ceil(n/nb);                  % number of chunks
    jx0=n-(nl-1)*nb;                % size of first chunk
    %
    th=(l-floor(l))*n;
    sd=(nargout > 3*(l~=0)); % = 1 if we are outputting log likelihood values
    lp=floor(l)+sd;   % extra loop needed to calculate final G value
    %
    lpx=zeros(1,n);             % log probability of each data point
    wk=ones(k,1);
    wp=ones(1,p);
    wpl=ones(1,pl);             % 1 index for lower triangular matrix
    wnb=ones(1,nb);
    wnj=ones(1,jx0);

    % EM loop

    g=0;                        % dummy initial value for comparison
    gg=zeros(lp+1,1);
    ss=sd;                      % initialize stopping count (0 or 1)
    vi=zeros(p*k,p);            % stack of k inverse cov matrices each size p*p
    vim=zeros(p*k,1);       	% stack of k vectors of the form inv(v)*m
    mtk=vim;                  	% stack of k vectors of the form m
    lvm=zeros(k,1);
    wpk=repmat((1:p)',k,1);
    for j=1:lp
        g1=g;               	% save previous log likelihood (2*pi factor omitted)
        m1=m;                	% save previous means, variances and weights
        v1=v;
        w1=w;
        for ik=1:k

            % these lines added for debugging only
            %             vk=reshape(v(k,lixi),p,p);
            %             condk(ik)=cond(vk);
            %%%%%%%%%%%%%%%%%%%%
            [uvk,dvk]=eig(reshape(v(ik,lixi),p,p));	% convert lower triangular to full and find eigenvalues
            dvk=max(diag(dvk),c);                	% apply variance floor to eigenvalues
            vik=-0.5*uvk*diag(dvk.^(-1))*uvk';      % calculate inverse
            vi((ik-1)*p+(1:p),:)=vik;               % vi contains all mixture inverses stacked on top of each other
            vim((ik-1)*p+(1:p))=vik*m(ik,:)';       % vim contains vi*m for all mixtures stacked on top of each other
            mtk((ik-1)*p+(1:p))=m(ik,:)';           % mtk contains all mixture means stacked on top of each other
            lvm(ik)=log(w(ik))-0.5*sum(log(dvk));       % vm contains the weighted sqrt of det(vi) for each mixture
        end
        %
        %         % first do partial chunk
        %
        jx=jx0;
        ii=1:jx;
        xii=xs(ii,:).';
        py=reshape(sum(reshape((vi*xii-vim(:,wnj)).*(xii(wpk,:)-mtk(:,wnj)),p,jx*k),1),k,jx)+lvm(:,wnj);
        mx=max(py,[],1);                % find normalizing factor for each data point to prevent underflow when using exp()
        px=exp(py-mx(wk,:));            % find normalized probability of each mixture for each datapoint
        ps=sum(px,1);                   % total normalized likelihood of each data point
        px=px./ps(wk,:);                % relative mixture probabilities for each data point (columns sum to 1)
        lpx(ii)=log(ps)+mx;
        pk=px*wx(ii);                       % effective fraction of data points for each mixture (could be zero due to underflow)
        sx=px*xsw(ii,:);
        sx2=px*(xsw(ii,rix).*xs(ii,cix));	% accumulator for variance calculation (lower tri cov matrix as a row)
        for il=2:nl
            ix=jx+1;
            jx=jx+nb;        % increment upper limit
            ii=ix:jx;
            xii=xs(ii,:).';
            py=reshape(sum(reshape((vi*xii-vim(:,wnb)).*(xii(wpk,:)-mtk(:,wnb)),p,nb*k),1),k,nb)+lvm(:,wnb);
            mx=max(py,[],1);                % find normalizing factor for each data point to prevent underflow when using exp()
            px=exp(py-mx(wk,:));            % find normalized probability of each mixture for each datapoint
            ps=sum(px,1);                   % total normalized likelihood of each data point
            px=px./ps(wk,:);                % relative mixture probabilities for each data point (columns sum to 1)
            lpx(ii)=log(ps)+mx;
            pk=pk+px*wx(ii);                    % effective fraction of data points for each mixture (could be zero due to underflow)
            sx=sx+px*xsw(ii,:);             % accumulator for mean calculation
            sx2=sx2+px*(xsw(ii,rix).*xs(ii,cix));	% accumulator for variance calculation
        end
        g=lpx*wx;                   % total log probability summed over all data points
        gg(j)=g;                    % save convergence history
        w=pk;               	% w(k,1) normalize to get the column of weights
        if pk                       % if all elements of pk are non-zero
            m=sx./pk(:,wp);         % find mean and mean square
            v=sx2./pk(:,wpl);
        else
            wm=pk==0;                       % mask indicating mixtures with zero weights
            nz=sum(wm);                  % number of zero-weight mixtures
            [vv,mk]=sort(lpx);             % find the lowest probability data points
            m=zeros(k,p);                   % initialize means and variances to zero (variances are floored later)
            v=zeros(k,pl);
            m(wm,:)=xs(mk(1:nz),:);                % set zero-weight mixture means to worst-fitted data points
            w(wm)=1/n;                      % set these weights non-zero
            w=w*n/(n+nz);                   % normalize so the weights sum to unity
            wm=~wm;                         % mask for non-zero weights
            m(wm,:)=sx(wm,:)./pk(wm,wp);  % recalculate means and variances for mixtures with a non-zero weight
            v(wm,:)=sx2(wm,:)./pk(wm,wpl);
        end
        v=v-m(:,cix).*m(:,rix);                 % subtract off mean squared
        if g-g1<=th && j>1
            if ~ss, break; end  %  stop
            ss=ss-1;       % stop next time
        end
    end
    if sd  % we need to calculate the final probabilities
        pp=lpx'-0.5*p*log(2*pi)-lsx;   % log of total probability of each data point
        gg=gg(1:j)-0.5*p*log(2*pi)-lsx;    % average log prob at each iteration
        g=gg(end);
        %             gg' % *** DEBUG ONLY ***
        m=m1;                                           % back up to previous iteration
        v=zeros(p,p,k);                                 % reserve spave for k full covariance matrices
        trv=0;                                          % sum of variance matrix traces
        for ik=1:k                                      % loop for each mixture to apply variance floor
            [uvk,dvk]=eig(reshape(v1(ik,lixi),p,p));	% convert lower triangular to full and find eigenvectors
            dvk=max(diag(dvk),c);                       % apply variance floor to eigenvalues
            v(:,:,ik)=uvk*diag(dvk)*uvk';               % reconstitute full matrix
            trv=trv+sum(dvk);                           % add trace to the sum
        end
        w=w1;
        mm=sum(m,1)/k;
        f=(m(:)'*m(:)-k*mm(:)'*mm(:))/trv;
    else
        v1=v;                                           % lower triangular form
        v=zeros(p,p,k);                                 % reserve spave for k full covariance matrices
        for ik=1:k                                      % loop for each mixture to apply variance floor
            [uvk,dvk,]=eig(reshape(v1(ik,lixi),p,p));	% convert lower triangular to full and find eigenvectors
            dvk=max(diag(dvk),c);                       % apply variance floor
            v(:,:,ik)=uvk*diag(dvk)*uvk';               % reconstitute full matrix
        end
    end
    m=m.*sx0(ones(k,1),:)+mx0(ones(k,1),:);  % unscale means
    v=v.*repmat(sx0'*sx0,[1 1 k]);
end
if l==0         % suppress the first three output arguments if l==0
    m=g;
    v=f;
    w=pp;
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [doseHandle,cMap,window] = matRad_plotDoseSlice(axesHandle,doseCube,plane,slice,threshold,alpha,cMap,window)
% matRad function that generates a dose plot of a selected slice
% The function can also be used in personal matlab figures by passing the
% corresponding axes handle.
%
% call
%   [doseHandle,cMap,window] = matRad_plotDoseSlice(axesHandle, doseCube,plane,slice,threshold)
%   [doseHandle,cMap,window] = matRad_plotDoseSlice(axesHandle, doseCube,plane,slice,threshold,alpha)
%   [doseHandle,cMap,window] = matRad_plotDoseSlice(axesHandle, doseCube,plane,slice,threshold,cMap)
%   [doseHandle,cMap,window] = matRad_plotDoseSlice(axesHandle, doseCube,plane,slice,threshold,window)
%   [doseHandle,cMap,window] = matRad_plotDoseSlice(axesHandle, doseCube,plane,slice,threshold,alpha,cMap)
%   [doseHandle,cMap,window] = matRad_plotDoseSlice(axesHandle, doseCube,plane,slice,threshold,alpha,window)
%   [doseHandle,cMap,window] = matRad_plotDoseSlice(axesHandle, doseCube,plane,slice,threshold,cMap,window)
%   [doseHandle,cMap,window] = matRad_plotDoseSlice(axesHandle, doseCube,plane,slice,threshold,alpha,cMap,window)
%
% input
%   axesHandle  handle to axes the slice should be displayed in
%   doseCube    3D array of the dose to select the slice from
%   plane       plane view (coronal=1,sagittal=2,axial=3)
%   slice       slice in the selected plane of the 3D cube
%   threshold   threshold above which the dose shall be displayed
%               for negative values (i.e. difference maps), also the values
%               smaller than the negative threshold will be displayed
%               if empty, no threshold will be applied
%   alpha       optional argument defining the alpha value, default is 0.6.
%               To use the default when providing a custom culormap, put in
%               an empty array by [].
%   cMap        optional argument defining the colormap, default is jet
%               if you want to use the default map with the window argument
%               you can use an empty array []
%   window      optional argument defining the displayed range. default is
%               [min(doseCube(:)) max(doseCube(:))]
%
% output
%   doseHandle: handle of the plotted dose axes
%   cMap        used colormap (same as input if set)
%   window      used window (same as input if set)
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

matRad_cfg = MatRad_Config.instance();

%Use default colormap?
if nargin < 7 || isempty(cMap)
    cMap = jet(64);
end
if nargin < 6 || isempty(alpha)
    alpha = 0.6;
end
if nargin < 8 || isempty(window)
    window = [min(doseCube(:)) max(doseCube(:))];
end

maxDose = max(doseCube(:));

cMapScale = size(cMap,1) - 1;

if plane == 1  % Coronal plane
    dose_slice = squeeze(doseCube(slice,:,:));
elseif plane == 2 % sagittal plane
    dose_slice = squeeze(doseCube(:,slice,:));
elseif plane == 3 % Axial plane
    dose_slice = squeeze(doseCube(:,:,slice));
else
	matRad_cfg.dispError('Invalid plane ''%d'' selected for visualization!',plane);
end

if ~isempty(threshold)
    mask = alpha * (dose_slice < window(2) & dose_slice > window(1) & dose_slice > threshold*maxDose);
else
    mask = alpha * (dose_slice < window(2) & dose_slice > window(1));
end

dose_slice = uint8(cMapScale*(dose_slice - window(1))/(window(2)-window(1)));

%This circumenvents a bug in Octave when the index in the image hase the maximum value of uint8
if matRad_cfg.isOctave
	dose_slice(dose_slice == 255) = 254;
end

dose_rgb = ind2rgb(dose_slice,cMap);

% plot dose distribution
doseHandle = image('CData',dose_rgb,'Parent',axesHandle);

% alphadata for image objects is not yet implemented in Octave
set(doseHandle,'AlphaData',mask);

end

repo: matRad


text:
function [transformationMatrix] = ...
    optical_flow_linear_registration2d(moving, fixed, varargin)
% OPTICAL_FLOW_LINEAR_REGISTRATION2D Estimates a transformation matrix using optical flow
%
% INPUT ARGUMENTS
% moving                - Moving iamge
% fixed                 - Fixed image
%
% OPTIONAL INPUT ARGUMENTS
% 'transformationModel'     - Transformation model for estimating the
%                             displacement field
%                             'translation', 'affine' (default)
%
% 'multiModal'              - Set whether to perform multi-modal or
%                             uni-modal image registration
%                             false (default), true
%
% 'numberOfChannels'        - Number of channels to use in when computing
%                             the entropy (based on channel coding). This
%                             is only relevant if multiModal is set to
%                             true.
%                             Default value is 8
%
% OUTPUT ARGUMENTS
% transformationMatrix  - Estimated transformation matrix

% Copyright (c) 2012
% danne.forsberg@outlook.com
%
% 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/>.

%% Setup default parameters
% translation, affine
transformationModel = 'affine';

% multi-modal
multiModal = false;

% number of channels, only valid for multi-modal registration
numberOfChannels = 8;

% Overwrites default parameter
for k=1:2:length(varargin)
    eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

% Initialize transformation matrix
transformationMatrix = eye(3);

[b_moving(:,:,1) b_moving(:,:,2)] = gradient(moving);
[b_fixed(:,:,1) b_fixed(:,:,2)] = gradient(fixed);

if multiModal
    b = b_fixed(:,:,2);
    b(:,:,2) = b_fixed(:,:,1);

    [delta_c mask] = estimate_delta_c(fixed,moving,numberOfChannels);
    delta_c(mask ~= 1) = 0;
    mask = repmat(mask,[1 1 2]);
    b(mask ~= 1) = 0;
else
    b = (b_moving(:,:,2) + b_fixed(:,:,2))/2;
    b(:,:,2) = (b_moving(:,:,1) + b_fixed(:,:,1))/2;

    delta_c = fixed - moving;
end

[G, h] = build_G_h_linear2d(b, delta_c, transformationModel);

% Solve the equation system
switch transformationModel
    case 'translation'
        d = G \ h;
        transformationMatrix(1:2,3) = d;
    case {'rigid','affine'}
        p = G \ h;
        transformationMatrix(1:2,1:3) = [1+p(1) p(2) p(5);...
                                         p(3) 1+p(4) p(6)];
end

author: fordanic

repo: image-registration


text:
function r8mat_l_inverse_test ( )

%*****************************************************************************80
%
%% R8MAT_L_INVERSE_TEST tests R8MAT_L_INVERSE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%

  n = 4;
%
%  Each row of this definition is a COLUMN of the matrix.
%
  a = [ ...
    1.0, 2.0, 4.0,  7.0; ...
    0.0, 3.0, 5.0,  8.0; ...
    0.0, 0.0, 6.0,  9.0; ...
    0.0, 0.0, 0.0, 10.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8MAT_L_INVERSE_TEST\n' );
  fprintf ( 1, '  R8MAT_L_INVERSE inverts a lower triangular matrix.\n' );

  r8mat_print ( n, n, a, '  Matrix A to be inverted:' );

  b = r8mat_l_inverse ( n, a );

  r8mat_print ( n, n, b, '  Inverse matrix B:' );

  c = a * b;

  r8mat_print ( n, n, c, '  Product C = A * B:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function b = ispathvalid_startup(p, options)
% ispathvalid can replace matlab's o.isdir, o.isfile and exist functions
% all of which have flaws (e.g., isdir not fully backwards compatible,
% exist has bugs where it'll confuse files and folders)
%
b = false;
if ~exist('options','var')
    options = '';
end
if ~ischar(options)
    return;
end
if ~isempty(options)
    if optionExists_startup(options,'file')

        % Err check 1.
        if isdir_private_startup(p)
            return;
        end

        % Err check 2.
        if optionExists_startup(options, 'checkextension')
            [~, ~, ext] = fileparts(p);
            if isempty(ext)
                return;
            end
        end

    elseif optionExists_startup(options,'dir')
        if isfile_private_startup(p)
            return;
        end
    else
        return;
    end
end
b = ~isempty(dir(p));

author: BUNPC


text:
function [varargout] = fIdentity(varargin)
% FIDENTITY  The global identity function
varargout = varargin;

author: Sable

repo: mcbench-benchmarks


text:
% Partially reproduces figure 1 from "Online Prediction of Time Series Data
% With Kernels". (3 algorithms, 25 MC simulations.)
%
% Learning curves for KNLMS, NORMA, and KRLS on a nonlinear system.
% Execution time: 2.5 minutes (Intel Pentium Core2 Duo).
%
% C. Richard, J.C.M. Bermudez, and P. Honeine, "Online Prediction of Time
% Series Data With Kernels," IEEE Transactions on Signal Processing,
% vol. 57, no. 3, pp. 1058-1067, March 2009,
% http://dx.doi.org/10.1109/TSP.2008.2009895
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab
% https://github.com/steven2358/kafbox/

clear
close all

%% PARAMETERS

N = 3000; % number of training data points
ktype = 'gauss';
kpar = 1/sqrt(2*3.73);

setups{1} = norma(struct('lambda',0.98,'tau',38,'eta',1,'tcoeff',-1/2,'kerneltype',ktype,'kernelpar',kpar));
setups{2} = knlms(struct('eta',.09,'eps',0.03,'mu0',0.5,'kerneltype',ktype,'kernelpar',kpar));
setups{3} = krls(struct('nu',.6,'kerneltype',ktype,'kernelpar',kpar));

numsim = 25;

%% RUN ALGORITHMS
t1 = tic;
fprintf('Fig. 1 from "Online Prediction of Time Series Data With Kernels".\n');

num_setup = length(setups);
MSE = zeros(N,num_setup);
titles = cell(num_setup,1);

for sim_ind = 1:numsim
    fprintf('SIM %d:\n',sim_ind)

    % Generate the data
    [X,y,yref] = generate_doddbench(N);

    for setup_ind = 1:num_setup
        kaf = setups{setup_ind};
        titles{setup_ind} = upper(class(kaf));

        for n=1:N
            if ~mod(n,floor(N/10)), fprintf('.'); end % progress indicator, 10 dots

            y_est = kaf.evaluate(X(n,:)); % test on test set
            err = yref(n) - y_est;
            MSE(n,setup_ind) = MSE(n,setup_ind) + err.^2/numsim;

            kaf.train(X(n,:),y(n)); % train with one input-output pair
        end
        fprintf('\n');
    end
end

toc(t1)
%% OUTPUT

% MSE smoothing by moving average for visualization
MSE_smooth = filter(1/20*ones(20,1),1,MSE);

figure
plot(10*log10(MSE_smooth));
title('Learning curves')
grid on

xlabel('iteration')
ylabel('MSE (dB)')
legend(titles)

repo: kafbox


text:
% This function is borrowed from rastamat

function d = deltas(x, w)
% D = deltas(X,W)  Calculate the deltas (derivatives) of a sequence
%    Use a W-point window (W odd, default 9) to calculate deltas using a
%    simple linear slope.  This mirrors the delta calculation performed
%    in feacalc etc.  Each row of X is filtered separately.
% 2003-06-30 dpwe@ee.columbia.edu

if nargin < 2
  w = 9;
end

[nr,nc] = size(x);

if nc == 0
  % empty vector passed in; return empty vector
  d = x;

else
  % actually calculate deltas

  % Define window shape
  hlen = floor(w/2);
  w = 2*hlen + 1;
  win = hlen:-1:-hlen;

  % pad data by repeating first and last columns
  xx = [repmat(x(:,1),1,hlen),x,repmat(x(:,end),1,hlen)];

  % Apply the delta filter
  d = filter(win, 1, xx, [], 2);  % filter along dim 2 (rows)

  % Trim edges
  d = d(:,2*hlen + [1:nc]);

end

author: jtkim-kaist

repo: VAD

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


text:
%% This file will be used to compare identified PDE of the SINDy-PI and ground truth.
% Code By: K
% Last Updated: 2019/06/24
%% Clear all
close all;clear all;clc;
%% Add path
[fld_status, fld_msg, fld_msgID]=mkdir('Datas');
[fld_status, fld_msg, fld_msgID]=mkdir('Figures');
addpath('Datas')
addpath('Figures')
addpath('Functions')
%% Define parameters for the simulation
% Define diffusion rate
Dx=0.01;Dz=0.01;Ds=1;Du=1;

% Others
q=0.1;f=1.5;ksi=0.3;alpha=0.3;beta=0.26;gama=0.4;ksi2=1.5;ksi3=0.003;phi=0;

% Define the time horizon
dt=0.001;T=1;
tspan=0:dt:T;

% Define the spatial domain
L=20; % Total length of each dimension
n=128; % Discretization point of each dimension
N=n*n; % Total points used, n^2

x2=linspace(-L/2,L/2,n+1); x=x2(1:n);
y=x;dx=x(2)-x(1);

kx=(2*pi/L)*[0:(n/2-1) -n/2:-1];
ky=kx;

% Get n-dimensional grid
[X,Y]=meshgrid(x,y);
[KX,KY]=meshgrid(kx,ky);
K2=KX.^2+KY.^2; K22=reshape(K2,N,1);
Kx=reshape(KX,N,1);Ky=reshape(KY,N,1);Kxx=reshape(KX,N,1);Kyy=reshape(KY,N,1);

% Create a time matrix
r_t=zeros(n,n,length(tspan));r=zeros(n,n,length(tspan));r_x=zeros(n,n,length(tspan));r_xx=zeros(n,n,length(tspan));r_y=zeros(n,n,length(tspan));r_yy=zeros(n,n,length(tspan));
z_t=zeros(n,n,length(tspan));z=zeros(n,n,length(tspan));z_x=zeros(n,n,length(tspan));z_xx=zeros(n,n,length(tspan));z_y=zeros(n,n,length(tspan));z_yy=zeros(n,n,length(tspan));
s_t=zeros(n,n,length(tspan));s=zeros(n,n,length(tspan));s_x=zeros(n,n,length(tspan));s_xx=zeros(n,n,length(tspan));s_y=zeros(n,n,length(tspan));s_yy=zeros(n,n,length(tspan));
u_t=zeros(n,n,length(tspan));u=zeros(n,n,length(tspan));u_x=zeros(n,n,length(tspan));u_xx=zeros(n,n,length(tspan));u_y=zeros(n,n,length(tspan));u_yy=zeros(n,n,length(tspan));
r_lap=zeros(n,n,length(tspan));z_lap=zeros(n,n,length(tspan));s_lap=zeros(n,n,length(tspan));u_lap=zeros(n,n,length(tspan));
%% Set a Guassian for the initial condition
r0=cos(sqrt(X.^2+Y.^2)).^2+GaussianFilter(X,Y,0.05,0.05,0,0)+GaussianFilter(X,Y,1,0.01,0,0)+GaussianFilter(X,Y,0.01,1,0,0)+GaussianFilter(X,Y,0.25,0.25,5,5)+GaussianFilter(X,Y,0.25,0.25,5,-5)+GaussianFilter(X,Y,0.25,0.25,-5,5)+GaussianFilter(X,Y,0.25,0.25,-5,-5)+0.1;
z0=sin(0.1*sqrt(X.^2+Y.^2)).^2+GaussianFilter(X,Y,0.05,0.01,0,0)+GaussianFilter(X,Y,0.25,0.25,5,5)+GaussianFilter(X,Y,0.25,0.25,5,-5)+GaussianFilter(X,Y,0.25,0.25,-5,5)+GaussianFilter(X,Y,0.25,0.25,-5,-5)+0.1;
s0=sin(sqrt(X.^2+Y.^2)).^2+GaussianFilter(X,Y,0.5,0.5,0,5)+GaussianFilter(X,Y,0.5,0.5,5,0)+GaussianFilter(X,Y,0.5,0.5,0,-5)+GaussianFilter(X,Y,0.5,0.5,-5,0)+0.1;
u0=cos(sqrt(X.^2+Y.^2)).^2+GaussianFilter(X,Y,0.5,0.01,0,0)+GaussianFilter(X,Y,0.01,0.5,0,0)+GaussianFilter(X,Y,0.25,0.25,5,5)+GaussianFilter(X,Y,0.25,0.25,5,-5)+GaussianFilter(X,Y,0.25,0.25,-5,5)+GaussianFilter(X,Y,0.25,0.25,-5,-5)+0.01;

%% Plot the initial conditions
figure(1)
surf(x,y,r0)
set(gca,'FontSize',18);
set(gcf,'Position',[100 100 600 600]);
set(gcf,'PaperPositionMode','auto');
view(-18,35)
box('off')
axis('off')

figure(2)
surf(x,y,z0)
set(gca,'FontSize',18);
set(gcf,'Position',[100 100 600 600]);
set(gcf,'PaperPositionMode','auto');
view(-18,35)
box('off')
axis('off')

figure(3)
surf(x,y,s0)
set(gca,'FontSize',18);
set(gcf,'Position',[100 100 600 600]);
set(gcf,'PaperPositionMode','auto');
view(-18,35)
box('off')
axis('off')

figure(4)
surf(x,y,u0)
set(gca,'FontSize',18);
set(gcf,'Position',[100 100 600 600]);
set(gcf,'PaperPositionMode','auto');
view(-18,35)
box('off')
axis('off')
%%
close all
% Transfer the initial guess to fourier domain
xzsu0t=[reshape(fft2(r0),1,N) reshape(fft2(z0),1,N) reshape(fft2(s0),1,N) reshape(fft2(u0),1,N)].';

%% Simulate the original equation
tic
opts = odeset('RelTol',1e-12,'AbsTol',1e-13);
NeedDev=0;
[time,xzsusol]=ode45(@(time,xzsut)BZ_Reaction_PDE(time,xzsut,Kx,Kxx,Ky,Kyy,K22,n,N,Dx,Dz,Ds,Du,q,f,ksi,alpha,beta,gama,ksi2,ksi3,phi,NeedDev),tspan,xzsu0t,opts);
toc

% Simulate the DL-SINDy discovered equation
tic
opts = odeset('RelTol',1e-12,'AbsTol',1e-13);
NeedDev=0;
[time,xzsusol_DL]=ode45(@(time,xzsut)BZ_Reaction_DL_SINDy_PDE(time,xzsut,Kx,Kxx,Ky,Kyy,K22,n,N,NeedDev),tspan,xzsu0t,opts);
toc

%% Extract values
% After simulation, we extract the derivative
NeedDev=1;

for pin=1:length(tspan)
    % Get the derivative data
    [~,r_t(:,:,pin),z_t(:,:,pin),s_t(:,:,pin),u_t(:,:,pin),r(:,:,pin),z(:,:,pin),s(:,:,pin),u(:,:,pin),r_x(:,:,pin),z_x(:,:,pin),s_x(:,:,pin),u_x(:,:,pin),r_y(:,:,pin),z_y(:,:,pin),s_y(:,:,pin),u_y(:,:,pin),r_xx(:,:,pin),z_xx(:,:,pin),s_xx(:,:,pin),u_xx(:,:,pin),r_yy(:,:,pin),z_yy(:,:,pin),s_yy(:,:,pin),u_yy(:,:,pin),r_lap(:,:,pin),z_lap(:,:,pin),s_lap(:,:,pin),u_lap(:,:,pin)]=...
        BZ_Reaction_PDE(0,xzsusol(pin,:).',Kx,Kxx,Ky,Kyy,K22,n,N,Dx,Dz,Ds,Du,q,f,ksi,alpha,beta,gama,ksi2,ksi3,phi,NeedDev);

    % Do the same for the DL-SINDy
    [~,~,~,~,~,r_DL(:,:,pin),z_DL(:,:,pin),s_DL(:,:,pin),u_DL(:,:,pin),~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~]=...
        BZ_Reaction_SINDy_PI_PDE(0,xzsusol_DL(pin,:).',Kx,Kxx,Ky,Kyy,K22,n,N,NeedDev);
end

%% Plot the result for original data
for pin=2:4
    if pin==1
        Animate=r;
        limits=[0 1];
    elseif pin==2
        Animate=z;
        limits=[0 1];
    elseif pin==3
        Animate=s;
        limits=[0 1];
    elseif pin==4
        Animate=u;
        limits=[0 1];
    end
    % Normalize it
    Dummy=Animate(:,:,1);
    Max=max(max(Dummy))
    Animate=Animate/Max;
    %
    for j=1:200:length(tspan)
        clf
        surf(x,y,Animate(:,:,j))
        colormap(cool)
        %shading interp
        %caxis([0,1])
        %Make plot propotional
        %zlim([0,1.1])
        %set(gca,'FontSize',18);
        set(gcf,'Position',[100 100 600 600]);
        set(gcf,'PaperPositionMode','auto');
        view(-18,35)
        box('off')
        axis('off')
        drawnow limitrate
    end
end

%% Plot the result for DL-SINDy simulation
for pin=2:4
    if pin==1
        Animate=r_DL;
        limits=[0 1];
    elseif pin==2
        Animate=z_DL;
        limits=[0 1];
    elseif pin==3
        Animate=s_DL;
        limits=[0 1];
    elseif pin==4
        Animate=u_DL;
        limits=[0 1];
    end
    % Normalize it
    Dummy=Animate(:,:,1);
    Max=max(max(Dummy))
    Animate=Animate/Max;
    %
    for j=1:200:length(tspan)
        clf
        surf(x,y,Animate(:,:,j))
        colormap(cool)
        %caxis([0,1])
        %Make plot propotional
        %zlim([0,1.1])
        %set(gca,'FontSize',18);
        set(gcf,'Position',[100 100 600 600]);
        set(gcf,'PaperPositionMode','auto');
        view(-18,35)
        box('off')
        axis('off')
        drawnow limitrate
    end
end

author: dynamicslab

repo: SINDy-PI


text:
function [c, r, t, cpulse, verbose] = ...
    tapas_physio_read_physlogfiles_siemens(log_files, cardiac_modality, verbose, varargin)
% reads out physiological time series and timing vector for Siemens
% logfiles of peripheral cardiac monitoring (ECG/Breathing Belt or
% pulse oximetry)
%
%   [cpulse, rpulse, t, c] = tapas_physio_read_physlogfiles_siemens(...
%           logfile, vendor, cardiac_modality)
%
% IN    log_files
%       .log_cardiac        contains ECG or pulse oximeter time course
%                           for GE: ECGData...
%       .log_respiration    contains breathing belt amplitude time course
%                           for GE: RespData...
%       cardiac_modality    'ECG' or 'PULS'/'PPU'/'OXY' to determine
%                           which channel data to be returned
%                           if not given, will be read out from file name
%                           suffix
%       verbose
%       .level              debugging plots are created if level >=3
%       .fig_handles        appended by handle to output figure
%
%       varargin            propertyName/value pairs, as folloes
%           'ecgChannel'    'v1', 'v2', 'mean' (default)
%                           determines which ECG channel to use as
%                           output cardiac curve
%           'sqpar'         sqpar, needed for scan_align to 'last' volume
%                           of a run
%
% OUT
%   cpulse              time events of R-wave peak in cardiac time series (seconds)
%                       for GE: usually empty
%   r                   respiratory time series
%   t                   vector of time points (in seconds)
%                       NOTE: This assumes the default sampling rate of 40
%                       Hz
%   c                   cardiac time series (ECG or pulse oximetry)
%
% EXAMPLE
%   [ons_secs.cpulse, ons_secs.rpulse, ons_secs.t, ons_secs.c] =
%       tapas_physio_read_physlogfiles_siemens(logfile, vendor, cardiac_modality);
%
%   See also tapas_physio_main_create_regressors

% Author: Lars Kasper
%         file structure information from PhLeM Toolbox, T. Verstynen (November 2007);
%                and Deshpande and J. Grinstead, Siemens Medical Solutions (March 2009)
%         additional log information Miriam Sebold, Charite Berlin (2014)
%
% Created: 2014-07-08
% Copyright (C) 2014 Institute for Biomedical Engineering, ETH/Uni Zurich.
%
% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

%% read out values

if nargin < 3
    verbose.level = 0;
end
DEBUG = verbose.level >=2;

% process optional input parameters and overwrite defaults
defaults.sqpar = [];
defaults.endCropSeconds     = 1;
% used channel depends on cardiac modality
switch cardiac_modality
    case 'ECG'
        defaults.ecgChannel = 'mean'; %'mean'; 'v1'; 'v2'
    otherwise
        defaults.ecgChannel = 'v1';
end

args = tapas_physio_propval(varargin, defaults);
tapas_physio_strip_fields(args);

cpulse              = [];

dt                  = log_files.sampling_interval;

explicit_relative_start_acquisition = log_files.relative_start_acquisition;

if isempty(explicit_relative_start_acquisition)
    explicit_relative_start_acquisition = 0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Determine relative start of acquisition from dicom headers and
% logfile footers
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

hasScanTimingDicomImage = ~isempty(log_files.scan_timing);
hasCardiacData = ~isempty(log_files.cardiac);
hasRespData = ~isempty(log_files.respiration);

if hasScanTimingDicomImage
    % Times are in seconds since start of the day

    dicomHeader             = spm_dicom_headers(...
        fullfile(log_files.scan_timing));

    try % old example DICOM format
        tStartScanDicom    = dicomHeader{1}.AcquisitionTime;

        TR = dicomHeader{1}.RepetitionTime/1000;

        % TODO: Include AcquisitionNumber? InstanceNumber?

    catch % new XA30 DICOM export
        dc = dicomHeader{1};

        % tried different fields, Study/SeriesTime refer to 1st volume
        % ContentTime > InstanceCreationTime > dc.AcquisitionDateTime
        % difference is about one second each
        % (1.14 and 0.89s for vol 250, 1.35 and 1.97s for vol 001, TR was 1.2s)
        % Since naming was most similar to AcquisitionDateTime, we chose
        % that one...
        % parse Date/Time format from Dicom field, and reformat it as
        % seconds since start of day
        dateVector = datevec(dc.AcquisitionDateTime, 'yyyymmddHHMMSS.FFF');
        tStartScanDicom = dateVector(end-2:end)*[3600 60 1]';
        TR = dc.SharedFunctionalGroupsSequence{1}.MRTimingAndRelatedParametersSequence{1}.RepetitionTime/1000;

    end

    tStopScanDicom     = tStartScanDicom + TR;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Read in cardiac data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if hasCardiacData

    [lineData, logFooter] = tapas_physio_read_physlogfiles_siemens_raw(...
        log_files.cardiac);
    tLogTotal = logFooter.LogStopTimeSeconds - logFooter.LogStartTimeSeconds;

    if hasScanTimingDicomImage
        tStartScan = tStartScanDicom; % this is just the start of the selected DICOM volume
        tStopScan = tStopScanDicom;
    else
        % Just different time scale, gives bad scaling in plots, and not
        % needed...
        %     tStartScan = logFooter.ScanStartTimeSeconds;
        %     tStopScan = logFooter.ScanStopTimeSeconds;
        tStartScan = logFooter.LogStartTimeSeconds;
        tStopScan = logFooter.LogStopTimeSeconds;
    end

    switch log_files.align_scan
        case 'first'
            relative_start_acquisition = tStartScan ...
                - logFooter.LogStartTimeSeconds;
        case 'last'
            % shift onset of first scan by knowledge of run duration and
            % onset of last scan in run
            relative_start_acquisition = ...
                (tStopScan - sqpar.Nscans*sqpar.TR) ...
                - logFooter.LogStartTimeSeconds;
    end

    % add arbitrary offset specified by user
    relative_start_acquisition = relative_start_acquisition + ...
        explicit_relative_start_acquisition;

    data_table = tapas_physio_siemens_line2table(lineData, cardiac_modality);

    if isempty(dt)
        nSamplesC = size(data_table,1);
        dt_c = tLogTotal/(nSamplesC-1);
    else
        dt_c = dt(1);
    end

    dataCardiac = tapas_physio_siemens_table2cardiac(data_table, ecgChannel, dt_c, ...
        relative_start_acquisition, endCropSeconds);

    if DEBUG
        verbose.fig_handles(end+1) = ...
            tapas_physio_plot_raw_physdata_siemens(dataCardiac);
    end

    %% crop end of log file
    cpulse = dataCardiac.cpulse_on;
    c = dataCardiac.c;
    t_c = dataCardiac.t;
    stopSample = dataCardiac.stopSample;

    cpulse(cpulse > t_c(dataCardiac.stopSample)) = [];
    t_c(stopSample+1:end) = [];
    c(stopSample+1:end) = [];

else
    c = [];
    t_c = [];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Read in respiratory data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if hasRespData
    [lineData, logFooter] = tapas_physio_read_physlogfiles_siemens_raw(...
        log_files.respiration);
    tLogTotal = logFooter.LogStopTimeSeconds - logFooter.LogStartTimeSeconds;

    if hasScanTimingDicomImage
        tStartScan = tStartScanDicom;
        tStopScan = tStopScanDicom; % is incorrect, use tStartScan + TR!
    else
        % Just different time scale, gives bad scaling in plots, and not
        % needed...
        %     tStartScan = logFooter.ScanStartTimeSeconds;
        %     tStopScan = logFooter.ScanStopTimeSeconds;
        tStartScan = logFooter.LogStartTimeSeconds;
        tStopScan = logFooter.LogStopTimeSeconds;
    end

    switch log_files.align_scan
        case 'first'
            relative_start_acquisition = tStartScan - ...
                logFooter.LogStartTimeSeconds;
        case 'last'
            % shift onset of first scan by knowledge of run duration and
            % onset of last scan in run
            relative_start_acquisition = ...
                (tStartScan - (sqpar.Nscans-1)*sqpar.TR) ...
                - logFooter.LogStartTimeSeconds;
    end

    % add arbitrary offset specified by user
    relative_start_acquisition = relative_start_acquisition + ...
        explicit_relative_start_acquisition;

    data_table = tapas_physio_siemens_line2table(lineData, 'RESP');

    if isempty(dt)
        nSamplesR = size(data_table,1);
        dt_r = tLogTotal/(nSamplesR-1);
    else
        dt_r = dt(end);
    end

    dataResp = tapas_physio_siemens_table2cardiac(data_table, ecgChannel, ...
        dt_r, relative_start_acquisition, endCropSeconds);

    if DEBUG
        verbose.fig_handles(end+1) = ...
            tapas_physio_plot_raw_physdata_siemens(dataResp);
    end

    r = dataResp.c;
    t_r = dataResp.t;

    %
    %% crop end of log file???
    % stopSample = dataResp.stopSample;
    % t(stopSample+1:end) = [];
    % c(stopSample+1:end) = [];

else
    r = [];
    t_r = [];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Adapt time scales resp/cardiac, i.e. zero fill c or r
% to get equal length with max(c_t, r_t)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nSamplesR = numel(r);
nSamplesC = numel(c);

if  nSamplesR > nSamplesC
    t = t_r; % time scale defined by respiration now, since longer
    if nSamplesC > 0 % zero-fill cardiac data
        c(nSamplesC+1:nSamplesR) = 0;
    end
else
    t = t_c;
    if nSamplesR > 0
        r(nSamplesR+1:nSamplesC) = 0;
    end
end

end

author: translationalneuromodeling

repo: tapas


text:
function program_07 ( )

%*****************************************************************************80
%
%% PROGRAM_07: Monte Carlo integral estimate for arbitrary triangle.
%
%  Discussion:
%
%    The program
%    * reads a triangle T (defined by three points),
%    * reads a random number seed;
%    * reads integer exponents P and Q;
%    * reads N, the number of random values to generate;
%    * it then computes N random points in the triangle;
%    * evaluates X^P * Y^Q at each point, and averages
%      to estimate the integral over the triangle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 February 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PROGRAM_07 - Monte Carlo estimate of\n' );
  fprintf ( 1, '  Integral x^p y^q\n' );
  fprintf ( 1, '  over an arbitrary triangle.\n' );
%
%  Get the triangle vertices.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Define a triangle T:\n' );

  t_v1 = input ( '  Enter [ T.v1.x, T.v1.y]: ' );
  t_v2 = input ( '  Enter [ T.v2.x, T.v2.y]: ' );
  t_v3 = input ( '  Enter [ T.v3.x, T.v3.y]: ' );
%
%  Get the random number seed.
%
  seed = input ( 'Enter a random number seed:  ' );
  rand ( 'state', seed );
%
%  Get the powers P and Q.
%
  p = input ( 'Enter the power P for X^P:  ' );
  q = input ( 'Enter the power Q for Y^Q:  ' );
%
%  Get the number of values to generate.
%
  n = input ( 'Enter the number of samples to generate:  ' );
%
%  Compute the area.
%
  t_area = 0.5 * ( t_v1(1) * ( t_v2(2) - t_v3(2) ) ...
                 + t_v2(1) * ( t_v3(2) - t_v1(2) ) ...
                 + t_v3(1) * ( t_v1(2) - t_v2(2) ) );
%
%  Generate sample points in the unit triangle,
%  map them to points in the given triangle,
%  evaluate the integrand, add to QUAD.
%
  quad = 0.0;

  for i = 1 : n

    r = rand;
    s = rand;

    xi1 =   1.0       - sqrt ( s );
    xi2 = ( 1.0 - r ) * sqrt ( s );
    xi3 =         r   * sqrt ( s );

    x = xi1 * t_v1(1) + xi2 * t_v2(1) + xi3 * t_v3(1);
    y = xi1 * t_v1(2) + xi2 * t_v2(2) + xi3 * t_v3(2);

    quad = quad + x^p * y^q;

  end
%
%  Normalize QUAD by the area and the number of points.
%
  quad = quad * t_area / n;
%
%  We didn't work out the exact integral, so all we have is
%  this estimate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Estimated integral = %12.6f\n', quad );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PROGRAM_07\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end


text:
function results = PTKComputeVentilatedVolume(ventilation_mask, region_mask)
    % PTKComputeVentilatedVolume. Computes the percentage of a region showing
    %     ventilation in gas MRI images
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2014.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.

    ventilation_mask_raw = ventilation_mask.RawImage > 0;
    ventilation_mask_raw = ventilation_mask_raw & region_mask.RawImage > 0;

    sum_ventilated_voxels = sum(ventilation_mask_raw(:));
    sum_region_voxels = sum(region_mask.RawImage(:));

    percentage_ventilated = 100*sum_ventilated_voxels/sum_region_voxels;

    ventilated_volume_mm3 = sum_ventilated_voxels*prod(ventilation_mask.VoxelSize);

    results = PTKMetrics;
    results.AddMetric('VentilatedVolumePercent', percentage_ventilated, 'Ventilated volume (%)');
    results.AddMetric('VentilatedVolume', ventilated_volume_mm3/1000, 'ventilated Volume (cm^3)');
end

author: tomdoel

repo: pulmonarytoolkit


text:
function randi_test ( )

%*****************************************************************************80
%
%% RANDI_TEST shows how random integers are generated in MATLAB.
%
%  Discussion:
%
%    The RANDI function is meant to replace MATLAB's previous function
%    for computing random integers, called RANDINT().
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2013
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'RANDI_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the MATLAB RANDI function.\n' );

  randi_test01 ( );

  randi_test02 ( );

  seed = 123456789;
  randi_test03 ( seed );

  seed = 987654321;
  randi_test03 ( seed );

  seed = 123456789;
  randi_test03 ( seed );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'RANDI_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

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

  return
end
function randi_test01 ( )

%*****************************************************************************80
%
%% RANDI_TEST01 simply calls the random integer generator a few times.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'RANDI_TEST01:\n' );
  fprintf ( 1, '  In MATLAB, random integers are generated by calling RANDI:\n' );
  fprintf ( 1, '  The maximum value is specified as IMAX.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A = randi ( imax, 1 ) a random scalar value between 1 and IMAX.\n' );
  fprintf ( 1, '  B = randi ( imax, 5, 1 ) a random column vector of 5 entries.\n' );
  fprintf ( 1, '  C = randi ( imax, 1, 5 ) a random row vector of 5 entries.\n' );
  fprintf ( 1, '  D = randi ( imax, 3, 4 ) a 3 by 4 random matrix.\n' );
  fprintf ( 1, '  E = randi ( imax, 5 )    a 5 by 5 random matrix.\n' );

  a = randi ( 9, 1 )
  b = randi ( 9, 5, 1 )
  c = randi ( 9, 1, 5 )
  d = randi ( 4, 3, 4 )
  e = randi ( 4, 5 )

  return
end
function randi_test02 ( )

%*****************************************************************************80
%
%% RANDI_TEST02 specifies the lower and upper limits.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'RANDI_TEST02:\n' );
  fprintf ( 1, '  RANDI allows the user to specify the numeric range.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  A = randi ( [ 5,    10 ],  1, 1 ) a random scalar value.\n' );
  fprintf ( 1, '  B = randi ( [ 7,     8 ], 10, 1 ) a random column vector of 5 entries.\n' );
  fprintf ( 1, '  C = randi ( [ -1,   +1 ],  1, 5 ) a random row vector of 5 entries.\n' );
  fprintf ( 1, '  D = randi ( [ -5,   +5 ],  3, 4 ) a 3 by 4 random matrix.\n' );
  fprintf ( 1, '  E = randi ( [ 100, 200 ],  5, 5 ) a 5 by 5 random matrix.\n' );

  a = randi ( [   5,  10 ], 1, 1 )
  b = randi ( [   7,   8 ], 10, 1 )
  c = randi ( [  -1,  +1 ], 1, 5 )
  d = randi ( [  -5,  +5 ], 3, 4 )
  e = randi ( [ 100, 200 ], 5, 5 )

  return
end
function randi_test03 ( seed )

%*****************************************************************************80
%
%% RANDI_TEST03 sets the seed before calling RANDI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'RANDI_TEST03:\n' );
  fprintf ( 1, '  By setting the random number seed, you can control\n' );
  fprintf ( 1, '  how the random number sequence begins.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The command "rng ( 123456789 )" sets the seed to 123456789.\n' );

  rng ( seed );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Seed has been set to %d\n', seed );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now generate 5 random values.\n' );

  for i = 1 : 5
    a = randi (  [ 1, 100 ], 1, 1 );
    fprintf ( 1, '  RANDI([1,100],1,1) = %g\n', a );
  end

  rng ( seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Seed has been reset to %d\n', seed );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now generate 5 more random values.\n' );

  for i = 1 : 5
    a = randi ( [ 1, 100 ], 1, 1 );
    fprintf ( 1, '  RANDI([1,100],1,1) = %g\n', a );
  end

  return
end
function timestamp ( )

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

  return
end


text:
function [datfield, dimord] = getdatfield(data)

% GETDATFIELD
%
% Use as
%   [datfield, dimord] = getdatfield(data)
% where the output arguments are cell-arrays.
%
% See also GETDIMORD, GETDIMSIZ

datfield  = fieldnames(data);
datfield  = datfield(:)'; % it should be a row-vector throughout the subsequent code

% these descriptive fields are cell-arrays and not treated as data
% the descriptive fields such as time and freq can be treated as data
datfield  = setdiff(datfield, {'label' 'labelcmb'});
datfield  = setdiff(datfield, ignorefields('selectdata'));

% find substructure fields
structfield = {};
for i=1:numel(datfield)
  if isstruct(data.(datfield{i}))
    structfield(end+1) = datfield(i);
  end
end
datfield  = setdiff(datfield, structfield);

% find homogenous transformation matrices, i.e. a 4x4 matrix that is named xxx2yyy
transformfield = {};
for i=1:numel(datfield)
  if ~isempty(regexp(datfield, '^[a-zA-Z]*2[a-zA-Z]*$', 'once')) && isequal(size(data.(datfield{i})), [4 4])
    transformfield(end+1) = datfield(i);
  end
end
datfield  = setdiff(datfield, transformfield);

xxxlabel  = datfield(~cellfun(@isempty, regexp(datfield, 'label$'))); % xxxlabel
datfield  = setdiff(datfield, xxxlabel);

xxxdimord = datfield(~cellfun(@isempty, regexp(datfield, 'dimord$'))); % xxxdimord
datfield  = setdiff(datfield, xxxdimord);

dimord = cell(size(datfield));
for i=1:length(datfield)
  dimord{i} = getdimord(data, datfield{i});
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% ----------------------------------------------------------------------------
%
%         Convert Allen CCF indices to Franklin-Paxinos labels
%
%   Based on data from Chon et al. Enhanced and unified anatomical labeling
%   for a common mouse brain atlas (2020).
%
% ----------------------------------------------------------------------------

%% Start with coordinates within the Allen CCF mouse brain atlas
% in the form [AP1, DV1, ML1
%              AP2, DV2, ML2]

% This example is of points along a neuropixels probe track through cortex, SC, PAG
brain_points =    [893   189   475
                   890    57   454
                   891   114   468
                   922   156   482
                   923   114   472
                   920   215   492
                   920   268   505
                   954   199   479
                   948   256   491
                   943   312   508
                   938   369   525
                   932   421   537
                   928   471   551
                   951   281   507
                   944   342   522
                   937   400   538
                   931   450   551];

% directory of reference files
annotation_volume_location = 'C:\Drive\Histology\annotation_volume_10um_by_index.npy'; % from the allen inst (see readme)
structure_tree_location = 'C:\Drive\Histology\structure_tree_safe_2017.csv'; % located in github repo
CCF_to_FP_location =  'C:\Drive\Histology\CCF_to_FP.csv'; % located in github repo
FP_table_location = 'C:\Drive\Histology\FP_table_Chon_2020.csv'; % located in github repo
chon_images_loc = 'C:\Drive\Histology\Suppl_File1_Labels'; % from chon et al (supplementary data 4, https://www.nature.com/articles/s41467-019-13057-w)

% generate values for pixel-to-coordinate transformation
bregma = allenCCFbregma(); % estimated bregma position in reference data space
atlas_resolution = 0.010; % pixels to mm

% should the brain image be dark or light
black_brain = true;
brain_points_color = [.5 .5 1];

%% load the reference brain annotations
if ~exist('av','var') || ~exist('st','var')
    disp('loading reference atlas...')
    av = readNPY(annotation_volume_location);
    st = loadStructureTree(structure_tree_location);
end
if ~exist('CCFtoFPtable','var') || ~exist('FPtable','var')
    disp('loading CCF-FP lookup tables...')
    CCFtoFPtable = loadCCFtoFP(CCF_to_FP_location);
    FPtable = loadFPtable(FP_table_location);
end

% initialize array of region annotations
annotation_CCF = cell(size(brain_points,1),3);
annotation_FP = cell(size(brain_points,1),3);

%% process data

% loop through every point to get ROI locations and region annotations
for point = 1:size(brain_points,1)

    % find the annotation, name, and acronym of the current point from
    % Allen CCF data
    ann = av(brain_points(point,1),brain_points(point,2),brain_points(point,3));
    name = st.safe_name{ann};
    acr = st.acronym{ann};

    annotation_CCF{point,1} = ann;
    annotation_CCF{point,2} = name;
    annotation_CCF{point,3} = acr;

    % find the annotation, name, and acronym of the current ROI pixel
    % using Chon et al data synthesizing CCF and Franklin-Paxinos
    [ann_FP, name_FP, acr_FP] = CCF_to_FP(brain_points(point,1), brain_points(point,2), brain_points(point,3), ...
                                          CCFtoFPtable, FPtable, chon_images_loc);

    annotation_FP{point,1} = ann_FP;
    annotation_FP{point,2} = name_FP;
    annotation_FP{point,3} = acr_FP;

end

% get coordinates relative to bregm\
ap = -(brain_points(:,1)-bregma(1))*atlas_resolution;
dv = (brain_points(:,2)-bregma(2))*atlas_resolution;
ml = (brain_points(:,3)-bregma(3))*atlas_resolution;

% generate table
data_table = table(annotation_CCF(:,2),annotation_CCF(:,3), annotation_FP(:,2),annotation_FP(:,3),...
                        ap,dv,ml, annotation_CCF(:,1), annotation_FP(:,1),...
         'VariableNames', {'CCF_name', 'CCF_abbrv', 'FP_name', 'FP_abbrv', 'AP_location', 'DV_location', 'ML_location', 'CCF_index', 'FP_index'});

%% display results

% plot points on the wire frame brain
fwireframe = plotBrainGrid([], [], [], black_brain); hold on;
fwireframe.InvertHardcopy = 'off';
figure(fwireframe); hold on
hp = plot3(brain_points(:,1), brain_points(:,3), brain_points(:,2), '.','linewidth',2, 'color',brain_points_color,'markers',10);

% display table
disp(data_table)

author: cortex-lab

repo: allenCCF

save_path: github-repos/MATLAB/cortex-lab-allenCCF