text:
function res = isdigit(s)
% res = isdigit(s);
% res(i) = 1 is s(1) is a digit between 0 and 9 (inclusive)
%
id = find(s >= '0' & s <= '9');
res = zeros(size(s));
res(id) = 1;

return;


text:
function count = subset_sum_count ( w, t, r )

%*****************************************************************************80
%
%% SUBSET_SUM_COUNT counts the solutions to the subset sum problem in a given range.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 May 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer W(N), a set of weights.  The length of this
%    array must be no more than 31.
%
%    Input, integer T, the target value.
%
%    Input, integer R(2), the lower and upper limits to be searched.
%    If this argument is omitted, the entire range, [0, 2^N-1 ] will
%    be searched.
%
%    Output, integer COUNT, the number of solutions found in this range.
%
  c = [];
  index = [];
%
%  Using a single integer to track the subsets only works if the number
%  of objects is small enough.  MATLAB can pack no more than 31 bits of
%  information into a nonnegative integer.
%
  if ( nargin < 1 )
    w = input ( 1, '  Enter the vector of weights [ w(1), w(2), ..., w(n)]: ' );
  end

  n = length ( w );

  if ( 31 < n )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SUBSET_SUM_FIND - Fatal error!\n' );
    fprintf ( 1, '  This function is restricted to N <= 31.\n' );
    error ( 'SUBSET_SUM_FIND - Fatal error!' );
  end

  if ( nargin < 2 )
    t = input ( '  Enter the target value T: ' );
  end
%
%  Make sure the range is reasonable.
%
  if ( nargin < 3 )
    r(1) = 0;
    r(2) = 2^n - 1;
  else
    r(1) = max ( r(1), 0 );
    r(2) = min ( r(2), 2^n - 1 );
  end

  fprintf ( 1, '\n' )
  fprintf ( 1, '  Searching indices %d through %d\n', r(1), r(2) );
%
%  Run through the range.
%
  count = 0;

  for index = r(1) : r(2)
%
%  Convert INDEX into vector of indices in W.
%
    c = find ( bitget ( index, 1:n ) );
%
%  If the sum of those weights matches the target, increment the count.
%
    if ( sum ( w(c) ) == t )
      count = count + 1;
    end

  end

  index = [];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%DEMEV2	Demonstrate Bayesian classification for the MLP.
%
%	Description
%	A synthetic two class two-dimensional dataset X is sampled  from a
%	mixture of four Gaussians.  Each class is associated with two of the
%	Gaussians so that the optimal decision boundary is non-linear. A 2-
%	layer network with logistic outputs is trained by minimizing the
%	cross-entropy error function with isotroipc Gaussian regularizer (one
%	hyperparameter for each of the four standard weight groups), using
%	the scaled conjugate gradient optimizer. The hyperparameter vectors
%	ALPHA and BETA are re-estimated using the function EVIDENCE. A graph
%	is plotted of the optimal, regularised, and unregularised decision
%	boundaries.  A further plot of the moderated versus unmoderated
%	contours is generated.
%
%	See also
%	EVIDENCE, MLP, SCG, DEMARD, DEMMLP2
%

%	Copyright (c) Ian T Nabney (1996-2001)

clc;

disp('This program demonstrates the use of the evidence procedure on')
disp('a two-class problem.  It also shows the improved generalisation')
disp('performance that can be achieved with moderated outputs; that is')
disp('predictions where an approximate integration over the true')
disp('posterior distribution is carried out.')
disp(' ')
disp('First we generate a synthetic dataset with two-dimensional input')
disp('sampled from a mixture of four Gaussians.  Each class is')
disp('associated with two of the Gaussians so that the optimal decision')
disp('boundary is non-linear.')
disp(' ')
disp('Press any key to see a plot of the data.')
pause;

% Generate the matrix of inputs x and targets t.

rand('state', 423);
randn('state', 423);

ClassSymbol1 = 'r.';
ClassSymbol2 = 'y.';
PointSize = 12;
titleSize = 10;

fh1 = figure;
set(fh1, 'Name', 'True Data Distribution');
whitebg(fh1, 'k');

%
% Generate the data
%
n=200;

% Set up mixture model: 2d data with four centres
% Class 1 is first two centres, class 2 from the other two
mix = gmm(2, 4, 'full');
mix.priors = [0.25 0.25 0.25 0.25];
mix.centres = [0 -0.1; 1.5 0; 1 1; 1 -1];
mix.covars(:,:,1) = [0.625 -0.2165; -0.2165 0.875];
mix.covars(:,:,2) = [0.25 0; 0 0.25];
mix.covars(:,:,3) = [0.2241 -0.1368; -0.1368 0.9759];
mix.covars(:,:,4) = [0.2375 0.1516; 0.1516 0.4125];

[data, label] = gmmsamp(mix, n);

%
% Calculate some useful axis limits
%
x0 = min(data(:,1));
x1 = max(data(:,1));
y0 = min(data(:,2));
y1 = max(data(:,2));
dx = x1-x0;
dy = y1-y0;
expand = 5/100;			% Add on 5 percent each way
x0 = x0 - dx*expand;
x1 = x1 + dx*expand;
y0 = y0 - dy*expand;
y1 = y1 + dy*expand;
resolution = 100;
step = dx/resolution;
xrange = [x0:step:x1];
yrange = [y0:step:y1];
%
% Generate the grid
%
[X Y]=meshgrid([x0:step:x1],[y0:step:y1]);
%
% Calculate the class conditional densities, the unconditional densities and
% the posterior probabilities
%
px_j = gmmactiv(mix, [X(:) Y(:)]);
px = reshape(px_j*(mix.priors)',size(X));
post = gmmpost(mix, [X(:) Y(:)]);
p1_x = reshape(post(:, 1) + post(:, 2), size(X));
p2_x = reshape(post(:, 3) + post(:, 4), size(X));

plot(data((label<=2),1),data(label<=2,2),ClassSymbol1, 'MarkerSize', ...
PointSize)
hold on
axis([x0 x1 y0 y1])
plot(data((label>2),1),data(label>2,2),ClassSymbol2, 'MarkerSize', ...
    PointSize)

% Convert targets to 0-1 encoding
target=[label<=2];
disp(' ')
disp('Press any key to continue')
pause; clc;

disp('Next we create a two-layer MLP network with 6 hidden units and')
disp('one logistic output.  We use a separate inverse variance')
disp('hyperparameter for each group of weights (inputs, input bias,')
disp('outputs, output bias) and the weights are optimised with the')
disp('scaled conjugate gradient algorithm.  After each 100 iterations')
disp('the hyperparameters are re-estimated twice.  There are eight')
disp('cycles of the whole algorithm.')
disp(' ')
disp('Press any key to train the network and determine the hyperparameters.')
pause;

% Set up network parameters.
nin = 2;		% Number of inputs.
nhidden = 6;		% Number of hidden units.
nout = 1;		% Number of outputs.
alpha = 0.01;		% Initial prior hyperparameter.
aw1 = 0.01;
ab1 = 0.01;
aw2 = 0.01;
ab2 = 0.01;

% Create and initialize network weight vector.
prior = mlpprior(nin, nhidden, nout, aw1, ab1, aw2, ab2);
net = mlp(nin, nhidden, nout, 'logistic', prior);

% Set up vector of options for the optimiser.
nouter = 8;			% Number of outer loops.
ninner = 2;			% Number of innter loops.
options = foptions;		% Default options vector.
options(1) = 1;			% This provides display of error values.
options(2) = 1.0e-5;		% Absolute precision for weights.
options(3) = 1.0e-5;		% Precision for objective function.
options(14) = 100;		% Number of training cycles in inner loop.

% Train using scaled conjugate gradients, re-estimating alpha and beta.
for k = 1:nouter
  net = netopt(net, options, data, target, 'scg');
  [net, gamma] = evidence(net, data, target, ninner);
  fprintf(1, '\nRe-estimation cycle %d:\n', k);
  disp(['  alpha = ', num2str(net.alpha')]);
  fprintf(1, '  gamma =  %8.5f\n\n', gamma);
  disp(' ')
  disp('Press any key to continue.')
  pause;
end

disp(' ')
disp('Network training and hyperparameter re-estimation are now complete.')
disp('Notice that the final error value is close to the number of data')
disp(['points (', num2str(n), ') divided by two.'])
disp('Also, the hyperparameter values differ, which suggests that a single')
disp('hyperparameter would not be so effective.')
disp(' ')
disp('First we train an MLP without Bayesian regularisation on the')
disp('same dataset using 400 iterations of scaled conjugate gradient')
disp(' ')
disp('Press any key to train the network by maximum likelihood.')
pause;
% Train standard network
net2 = mlp(nin, nhidden, nout, 'logistic');
options(14) = 400;
net2 = netopt(net2, options, data, target, 'scg');
y2g = mlpfwd(net2, [X(:), Y(:)]);
y2g = reshape(y2g(:, 1), size(X));

disp(' ')
disp('We can now plot the function represented by the trained networks.')
disp('We show the decision boundaries (output = 0.5) and the optimal')
disp('decision boundary given by applying Bayes'' theorem to the true')
disp('data model.')
disp(' ')
disp('Press any key to add the boundaries to the plot.')
pause;

% Evaluate predictions.
[yg, ymodg] = mlpevfwd(net, data, target, [X(:) Y(:)]);
yg = reshape(yg(:,1),size(X));
ymodg = reshape(ymodg(:,1),size(X));

% Bayesian decision boundary
[cB, hB] = contour(xrange,yrange,p1_x,[0.5 0.5],'b-');
[cNb, hNb] = contour(xrange,yrange,yg,[0.5 0.5],'r-');
[cN, hN] = contour(xrange,yrange,y2g,[0.5 0.5],'g-');
set(hB, 'LineWidth', 2);
set(hNb, 'LineWidth', 2);
set(hN, 'LineWidth', 2);
Chandles = [hB(1) hNb(1) hN(1)];
legend(Chandles, 'Bayes', ...
  'Reg. Network', 'Network', 3);

disp(' ')
disp('Note how the regularised network predictions are closer to the')
disp('optimal decision boundary, while the unregularised network is')
disp('overtrained.')

disp(' ')
disp('We will now compare moderated and unmoderated outputs for the');
disp('regularised network by showing the contour plot of the posterior')
disp('probability estimates.')
disp(' ')
disp('The first plot shows the regularised (moderated) predictions')
disp('and the second shows the standard predictions from the same network.')
disp('These agree at the level 0.5.')
disp('Press any key to continue')
pause
levels = 0:0.1:1;
fh4 = figure;
set(fh4, 'Name', 'Moderated outputs');
hold on
plot(data((label<=2),1),data(label<=2,2),'r.', 'MarkerSize', PointSize)
plot(data((label>2),1),data(label>2,2),'y.', 'MarkerSize', PointSize)

[cNby, hNby] = contour(xrange, yrange, ymodg, levels, 'k-');
set(hNby, 'LineWidth', 1);

fh5 = figure;
set(fh5, 'Name', 'Unmoderated outputs');
hold on
plot(data((label<=2),1),data(label<=2,2),'r.', 'MarkerSize', PointSize)
plot(data((label>2),1),data(label>2,2),'y.', 'MarkerSize', PointSize)

[cNbm, hNbm] = contour(xrange, yrange, yg, levels, 'k-');
set(hNbm, 'LineWidth', 1);

disp(' ')
disp('Note how the moderated contours are more widely spaced.  This shows')
disp('that there is a larger region where the outputs are close to 0.5')
disp('and a smaller region where the outputs are close to 0 or 1.')
disp(' ')
disp('Press any key to exit')
pause
close(fh1);
close(fh4);
close(fh5);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%ASAGLOB_MEAN_ADJ Subtraction of the signal mean identifier
%   Declaring the variable ASAGLOB_MEAN_ADJ in the base workspace, for
%   example by typing at the command prompt,
%     ASAglob_mean_adj = 1;
%   will enable the ARMASA functions, ARMASEL, SIG2AR, SIG2MA and
%   SIG2ARMA to read out the value of ASAGLOB_MEAN_ADJ. A value 1 will
%   inform these functions, that the input data SIG have been level
%   adjusted using the mean of SIG. This will have a small effect on the
%   estimation of the Prediction Error, that is determined in these
%   functions.
%
%   By default the functions subtract the mean from the input data SIG
%   and ASAGLOB_MEAN_ADJ is automatically assigned. If one choses to
%   subtract the mean beforehand, or when it is known that the mean has
%   already been subtracted by some other proces, assign 1 to
%   ASAGLOB_MEAN_ADJ. To control the subtraction of the mean itself, use
%   ASAGLOB_SUBTR_MEAN. ASAGLOB_SUBTR_MEAN will NOT be effected by any
%   choice of ASAGLOB_MEAN_ADJ.
%
%   See also: ASAGLOB_SUBTR_MEAN

disp('  Variable ASAGLOB_MEAN_ADJ is currently undefined.')
disp('  Type ''help ASAglob_mean_adj'' for more information.')

author: Sable

repo: mcbench-benchmarks


text:
function a_cum = i4vec_cum ( n, a )

%*****************************************************************************80
%
%% I4VEC_CUM computes the cumulative sum of the entries of an I4VEC.
%
%  Discussion:
%
%    An I4VEC is a vector of I4's.
%
%  Example:
%
%    Input:
%
%      A = (/ 1, 2, 3, 4 /)
%
%    Output:
%
%      A_CUM = (/ 1, 3, 6, 10 /)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, integer A(N), the vector to be summed.
%
%    Output, integer A_CUM(1:N), the cumulative sum of the entries of A.
%
  a_cum(1) = a(1);

  for i = 2 : n
    a_cum(i) = a_cum(i-1) + a(i);
  end

  return
end


text:
function pred = ml_predictgmm(trials, model)
% Prediction function for Gaussian Mixture Models.
% Prediction = ml_predictgmm(Trials, Model)
%
% In:
%   Trials  : the data a matrix, as in ml_predict
%
%   Model   : predictive model as produced by ml_trainlogreg
%
% Out:
%   Prediction  : discrete probability distribution, formatted as
%                 {'disc' [NxC] [Cx1]}, with element #2 being the per-class probability and
%                 element #3 the original target values per class
%                 thus, the expected target values are Prediction{2}*Prediction{3}
%
% Examples:
%   targets might look like this: [-1 -1 1 -1 1 -1 -1 1 -1 -1 1 -1 -1 1 ...]'
%
%   model = ml_traingmm(data,targets)
%   p = ml_predictgmm(data, model); expectation = p{2}*p{3};
%   now expectation might look like this: [-0.6 -0.9 +0.4 -0.7 +0.8 -0.1 +0.5 +1.0 -0.9 +1.0 -1.0 -1.0 +1.0 ...]'
%
% See also:
%   ml_traingmm
%
%                           Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                           2010-04-05

% scale data
trials = hlp_applyscaling(trials,model.sc_info);

% predict probabilities
switch model.variant
    case {'avdp','vdp','bj','bjrnd','cdp','csb','vb'}
        for c=1:length(model.classes)
            model.class{c}.test_data = trials';
            % obtain log probabilities under each component
            tmp = vdpgm(model.data{c}',model.class{c});
            pdfmat(:,c) = exp(tmp.predictive_posterior');
        end
        % normalize the log probabilities, incorporate the prior, and turn them into true conditional probabilities
        probs = gmmb_normalize(pdfmat);
    case {'em','fj','gem'}
        pdfmat = gmmb_pdf(trials, model.class);
        probs = gmmb_normalize(gmmb_weightprior(pdfmat, model.class));
end
pred = {'disc', probs, model.classes};

repo: REST


text:
function jed = ymdf_to_jed_islamic_a2 ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_TO_JED_ISLAMIC_A2 converts an Islamic A YMDF date to a JED.
%
%  Discussion:
%
%    The algorithm has the beauty of being comprehensible%
%
%    Count the days up to the day before the start of the calendar,
%    plus the days in the current month, the 29 days guaranteed
%    in the previous months of this year, the (months/2) 30th days,
%    the 354 days in each of the previous years, plus the total number
%    of leap days in the preceding years.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Reingold, Nachum Dershowitz,
%    Calendrical Calculations I,
%    Software - Practice and Experience,
%    Volume 20, Number 9, September 1990, pages 899-928.
%
%  Parameters:
%
%    Input, integer Y, M, D, real F, the YMDF date.
%
%    Output, real JED, the corresponding Julian Ephemeris Date.
%

%
%  Check the date.
%
  [ y, m, d, ierror ] = ymd_check_islamic ( y, m, d );

  if ( ierror ~= 0 )
    jed = -1.0;
    return
  end

  jed_epoch = epoch_to_jed_islamic_a ( );

  jed = jed_epoch - 1 + d + 29 * ( m - 1 ) + floor ( m / 2 ) ...
    + 354 * ( y - 1 ) + floor ( ( 11 * y + 3 ) / 30 ) + f;

  return
end


text:
function g = fileKernGradient(kern, x, varargin)

% FILEKERNGRADIENT Gradient of FILE kernel's parameters.
% FORMAT
% DESC computes the gradient of functions with respect to the
% stored file
% kernel's parameters. As well as the kernel structure and the
% input positions, the user provides a matrix PARTIAL which gives
% the partial derivatives of the function with respect to the
% relevant elements of the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x : the input locations for which the gradients are being
% computed.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The argument takes
% the form of a square matrix of dimension  numData, where numData is
% the number of rows in X.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters. The ordering of the vector should match
% that provided by the function kernExtractParam.
%
% FORMAT
% DESC computes the derivatives as above, but input locations are
% now provided in two matrices associated with rows and columns of
% the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x1 : the input locations associated with the rows of the
% kernel matrix.
% ARG x2 : the input locations associated with the columns of the
% kernel matrix.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The matrix should
% have the same number of rows as X1 and the same number of columns
% as X2 has rows.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters.
%
% SEEALSO fileKernParamInit, kernGradient, fileKernDiagGradient, kernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006

% KERN

% The last argument is covGrad
if nargin < 4
  [k, sk] = fileKernCompute(kern, x);
else
  [k, sk] = fileKernCompute(kern, x, varargin{1});
end
g(1) =  sum(sum(varargin{end}.*sk));
%/~
if any(isnan(g))
  warning('g is NaN')
end
%~/

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function U = fda(X, t, q)
% Fisher (linear) discriminant analysis
% Input:
%   X: d x n data matrix
%   t: 1 x n class label
%   d: target dimension
% Output:
%   U: projection matrix y=U'*x
% Written by Mo Chen (sth4nth@gmail.com).
n = size(X,2);
k = max(t);

E = sparse(1:n,t,true,n,k,n);  % transform label into indicator matrix
nk = full(sum(E));

m = mean(X,2);
Xo = bsxfun(@minus,X,m);
St = (Xo*Xo')/n;                   % 4.43

mk = bsxfun(@times,X*E,1./nk);
mo = bsxfun(@minus,mk,m);
mo = bsxfun(@times,mo,sqrt(nk/n));
Sb = mo*mo';                       % 4.46
% Sw = St-Sb;                        % 4.45

[U,A] = eig(Sb,St,'chol');
[~,idx] = sort(diag(A),'descend');
U = U(:,idx(1:q));

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function X = mohsst5_remove_land(X)

data = mohsst5_loaddata();

land = sum(~isnan(data.observations),2) == 0;

X(land,:) = [];

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
% Driver script for solving the 1D Burgers equations
Globals1D;

% Orer of polymomials used for approximation
N = 6;

% Generate simple mesh
xL = -1.0; xR = 1.0;
[Nv, VX, K, EToV] = MeshGen1D(xL,xR,10);

% Initialize solver and construct grid and metric
StartUp1D;

% Set initial conditions
epsilon = 0.1;
%u = -tanh((x+0.5)/(2*epsilon)) + 1.0;
%u=2*sin(x);
xmid=(x(end)+x(1))/2; x1=x<xmid; x2=x>=xmid; u=1*x1+0.1*x2;

% Solve Problem
FinalTime = 0.15;
[u] = Burgers1D(u,epsilon,xL,xR,FinalTime);

repo: Aero-matlab


text:
function p02_default ( dim_num )

%*****************************************************************************80
%
%% P02_DEFAULT sets default values for problem 02.
%
%  Discussion:
%
%    If a problem uses vector parameters, then the spatial dimension
%    DIM_NUM is needed as input, so that the vector parameters can be
%    properly dimensioned.
%
%    Every problem keeps a count of the number of function calls.  Calling
%    this routine causes that count to be reset to 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the problem.
%
  i4 = 0;
  p02_i4 ( 'D', '*', i4 );

  return
end


text:
function r82row_print_part_test ( )

%*****************************************************************************80
%
%% R82ROW_PRINT_PART_TEST tests R82ROW_PRINT_PART.
%
%  Discussion:
%
%    An R82ROW is a (2,N) array of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 April 2015
%
%  Author:
%
%    John Burkardt
%
  n = 10;
%
%  Set up the array by rows.
%
  v = [ ...
    11.0,  12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0; ...
    21.0,  22.0, 23.0, 24.0, 25.0, 26.0, 27.0, 28.0, 29.0, 30.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R82ROW_PRINT_PART_TEST\n' );
  fprintf ( 1, '  R82ROW_PRINT_PART prints part of an R82ROW.\n' );

  max_print = 2;
  r82row_print_part ( n, v, max_print, '  Output with MAX_PRINT = 2' );

  max_print = 5;
  r82row_print_part ( n, v, max_print, '  Output with MAX_PRINT = 5' );

  max_print = 25;
  r82row_print_part ( n, v, max_print, '  Output with MAX_PRINT = 25' );

  return
end


text:
function [m] = cellssq(x, dim)

% [S] = CELLSSQ(X, DIM) computes the sum of squares, across all cells in x along
% the dimension dim.
%
% X should be an linear cell-array of matrices for which the size in at
% least one of the dimensions should be the same for all cells

nx = size(x);
if ~iscell(x) || length(nx)>2 || all(nx>1),
  error('incorrect input for cellssq');
end

if nargin==1,
  scx1 = cellfun('size', x, 1);
  scx2 = cellfun('size', x, 2);
  if     all(scx2==scx2(1)), dim = 2; %let second dimension prevail
  elseif all(scx1==scx1(1)), dim = 1;
  else   error('no dimension to compute sum of squares for');
  end
end

nx   = max(nx);
ssmp = cellfun(@sumsq,   x, repmat({dim},1,nx), 'UniformOutput', 0);
m    = sum(cell2mat(ssmp), dim);

function [s] = sumsq(x, dim)

s = sum(x.^2, dim);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [a,T] = compute_tf_model(xyu_trimmed, p_drone, p_physics)
% COMPUTE_TF_MODEL - Function computing the transfer function model
%
% Syntax: out = airdata(x)
%
% Inputs:
%   x_trim - desired trimmed state
%   u_trim - desired trimmed input delta
%   y_trim - desired trimmed air data
%   P      - parameters structure
%
% Outputs:
%   T_phi_da        - transfer function (tf) from delta_aileron to roll
%   T_chi_phi       - tf from chi to roll
%   T_theta_de      - tf from delta_elevator to pitch
%   T_h_theta       - tf from pitch to altitude
%   T_h_Va          - tf from Va to altitude
%   T_Va_dt         - tf from delta_thrust to Va
%   T_Va_theta      - tf from pitch to Va
%   T_vy_dr         - tf from delta_rudder to vy

x_trim = xyu_trimmed(1:12);
y_trim = xyu_trimmed(13:15);
u_trim = xyu_trimmed(16:19);

% Relabel inputs
% Position in NED
% pn          = x_trim(1);
% pe          = x_trim(2);
% pd          = x_trim(3); % -altitude

% UAV velocity wrt inertial frame in body frame
% vx          = x_trim(4);
% vy          = x_trim(5);
% vz          = x_trim(6);
% v_xyz       = x_trim(4:6);

% Euler angles
% phi         = x_trim(7);
theta       = x_trim(8);
% psi         = x_trim(9);

% Rotation rates
% p           = x_trim(10);
% q           = x_trim(11);
% r           = x_trim(12);

% Actuators
de = u_trim(1); % delta elevator
% da = u_trim(2); % delta aileron
% dr = u_trim(3); % delta rudder
dt = u_trim(4); % delta thrust

% Wind
Va =    y_trim(1);
alpha = y_trim(2);
% beta =  y_trim(3);

% Compute parametrs
p_dyn = 0.5*p_physics.rho*Va^2;
kb = p_drone.kb(Va);
if isinf(kb)
    kb =0;
end
kc = p_drone.kc(Va);
if isinf(kc)
    kc =0;
end

a.phi1 = - p_dyn*p_drone.S_wing*p_drone.b*p_drone.Cp_p*kb;
a.phi2 =   p_dyn*p_drone.S_wing*p_drone.b*p_drone.Cp_da;

a.beta1 = - 0.5 * p_physics.rho * Va * p_drone.S_wing / p_drone.mass * p_drone.CY_beta;
a.beta2 =   0.5 * p_physics.rho * Va * p_drone.S_wing / p_drone.mass * p_drone.CY_dr;

a.theta1 = - p_dyn * p_drone.c * p_drone.S_wing / p_drone.Jy * p_drone.Cm_q * kc;
a.theta2 = - p_dyn * p_drone.c * p_drone.S_wing / p_drone.Jy * p_drone.Cm_alpha;
a.theta3 =   p_dyn * p_drone.c * p_drone.S_wing / p_drone.Jy * p_drone.Cm_de;

a.va1 = p_physics.rho * Va * p_drone.S_wing / p_drone.mass * ( p_drone.CD0 + p_drone.CD_alpha*alpha + p_drone.CD_de*de ) + ...
       p_physics.rho * p_drone.S_prop / p_drone.mass * p_drone.C_prop * Va;
a.va2 = p_physics.rho * p_drone.S_prop / p_drone.mass * p_drone.C_prop * p_drone.k_motor^2 * dt;
a.va3 = p_physics.gravity * cos(theta - alpha);

% Define transfer functions
T.phi_da   = tf([a.phi2],[1,a.phi1,0]);
T.chi_phi       = tf([p_physics.gravity/Va],[1,0]);
T.theta_de = tf(a.theta3,[1,a.theta1,a.theta2]);
T.h_theta       = tf([Va],[1,0]);
T.h_Va          = tf([theta],[1,0]);
T.Va_dt    = tf([a.va2],[1,a.va1]);
T.Va_theta      = tf([-a.va3],[1,a.va1]);
T.vy_dr     = tf([Va*a.beta2],[1,a.beta1]);

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab


text:
function [nfft,tableout]=floor23(n)
%FLOOR23  Previous number with only 2,3 factors
%   Usage: nceil=floor23(n);
%
%   `floor23(n)` returns the first number less than or equal to *n*,
%   which can be written as a product of powers of *2* and *3*.
%
%   The algorithm will look up the best size in a table, which is computed
%   the first time the function is run. If the input size is larger than the
%   largest value in the table, the input size will be reduced by factors of
%   *2*, until it is in range.
%
%   `[nceil,table]=floor23(n)` additionally returns the table used for lookup.
%
%   Examples:
%   ---------
%
%   Return the first number smaller or equal to *26* that can be written
%   solely as products of powers of *2* and *3*:::
%
%     floor23(26)
%
%   This plot shows the behaviour of |floor23| and |ceil23| for numbers
%   up to 100:::
%
%     x=1:100;
%     plot(x,floor23(x),x,ceil23(x));
%     legend('floor23','ceil23','Location','Northwest');
%
%   See also: ceil23, floor235, nextfastfft

%   AUTHOR: Peter L. Søndergaard

persistent table;

maxval=2^20;

if isempty(table)
    % Compute the table for the first time, it is empty.
    l2=log(2);
    l3=log(3);
    l5=log(5);
    lmaxval=log(maxval);
    table=zeros(143,1);
    ii=1;
    prod2=1;
    for i2=0:floor(lmaxval/l2)
        prod3=prod2;
        for i3=0:floor((lmaxval-i2*l2)/l3)
            table(ii)=prod3;
            prod3=prod3*3;
            ii=ii+1;
        end;
        prod2=prod2*2;
    end;
    table=sort(table);
end;

% Copy input to output. This allows us to efficiently work in-place.
nfft=n;

% Handle input of any shape by Fortran indexing.
for ii=1:numel(n)
  n2reduce=0;

  if n(ii)>maxval
    % Reduce by factors of 2 to get below maxval
    n2reduce=ceil(log2(nfft(ii)/maxval));
    nfft(ii)=nfft(ii)/2^n2reduce;
  end;

  % Use a simple bisection method to find the answer in the table.
  from=1;
  to=numel(table);
  while from<=to
    mid = round((from + to)/2);
    diff = table(mid)-nfft(ii);
    if diff<0
      from=mid+1;
    else
      to=mid-1;
    end
  end
  if nfft(ii)~=table(from)
      nfft(ii)=table(from-1);
  end;

  % Add back the missing factors of 2 (if any)
  nfft(ii)=nfft(ii)*2^n2reduce;

end;

tableout=table;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
clear; clc; close all;
wkdir = '../'; % The root foler of FM-Bench
addpath([wkdir 'vlfeat-0.9.21/toolbox/']);
vl_setup;

% Datasets = {'TUM', 'KITTI', 'Tanks_and_Temples', 'CPC'};
Datasets = {'KITTI'};

matcher='SIFT-RT'; % SIFT with Ratio Test
estimator='RANSAC';

for s = 1 : length(Datasets)
     dataset = Datasets{s};

    % An example for DoG detector
    FeatureDetection(wkdir, dataset);

    % An example for SIFT descriptor
    FeatureExtraction(wkdir, dataset);

    % An example for exhaustive nearest neighbor matching with ratio test
    FeatureMatching(wkdir, dataset, matcher);

    % An example for RANSAC based FM estimation
    GeometryEstimation(wkdir, dataset, matcher, estimator);

end

author: JiawangBian

repo: FM-Bench


text:
function [C_df,Df] = extract_DF_F_new(A,C,b,f,P,options)

% extract DF/F signals after performing NMF
% inputs:  Y raw data (d X T matrix, d # number of pixels, T # of timesteps)
%          A matrix of spatial components (d x K matrix, K # of components)
%          C matrix of temporal components (K x T matrix)
%          P neuron structure, used to read the baseline activity for each
%                    component of C
%          options structure used for specifying method for determining DF
%           default method is the median of the trace. By changing
%           options.df_prctile an arbitray percentile can be used (between 0 and 100).
%           a moving window can also be established by specifying options.df_window

% outputs:  C_df temporal components in the DF/F domain
%           Df   background for each component to normalize the filtered raw data

% Written by:
% Eftychios A. Pnevmatikakis, Simons Foundation, 2016

%memmaped = isobject(Y);
defoptions = CNMFSetParms;
if nargin < 6 || isempty(options)
    options = defoptions;
end
if ~isfield(options,'df_prctile') || isempty(options.df_prctile)
    options.df_prctile = defoptions.df_prctile;
end
if ~isfield(options,'df_window') || isempty(options.df_window)
    options.df_window = defoptions.df_window;
end
if ~isfield(options,'full_A') || isempty(options.full_A); full_A = defoptions.full_A; else full_A = options.full_A; end

[K,T] = size(C);

Bas = zeros(K,T);

if ~(nargin < 5 || isempty(P))
    bas_val = cell2mat(P.b);
    Ntr = size(bas_val,2);
    if Ntr > 1
        ln = diff(P.cs_frtrs);
        for i = 1:Ntr
            Bas(:,:,P.cs_frtrs(i)+1:P.cs_frtrs(i+1)) = repmat(bas_val(:,i),1,ln(i));
        end
    else
        Bas = repmat(bas_val,1,T);
    end
end

nA = sqrt(sum(A.^2));

AA = A'*A;
AA(1:K+1:end) = 0;

Cf = bsxfun(@times,C - Bas,nA(:).^2);
C2 = repmat(AA*bas_val,1,T) + (A'*b)*f;

if isempty(options.df_window) || (options.df_window > size(C,2))
    if options.df_prctile == 50
        Df = median(C2,2);
    else
        Df = prctile(C2,options.df_prctile,2);
    end
    C_df = bsxfun(@times,Cf,1./Df(:));
else
    if options.df_prctile == 50
        if verLessThan('matlab','2015b')
            warning('Median filtering at the boundaries might be inaccurate due to zero padding.')
            Df = medfilt1(C2,options.df_window,[],2);
        else
            Df = medfilt1(C2,options.df_window,[],2,'truncate');
        end
    else
        Df = zeros(size(C2));
        for i = 1:size(Df,1);
            df_temp = running_percentile(C2(i,:), options.df_window, options.df_prctile);
            Df(i,:) = df_temp(:)';
        end
    end
    C_df = Cf./Df;
end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function[C] = Cp(rho)
% C = Cp(rho)
% C_p(rho) = sqrt{\frac{1}{2\pi\int_{0}^1 r rho^2 J_1(rho r)^2 dr}}
% This constant is chosen so that e_p(r) = C_p(rho) J_0(rho_p r) has a unit
% H_0^1(B) norm. (that is \int_{B}|\nabla e_p|^2 = 1).

C = sqrt(1./(pi*rho.*(besselj(0, rho).^2.*rho+besselj(1, rho).^2.*rho-2*besselj(0, rho).*besselj(1, rho))));
C(rho==0) = 1;

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%  fix_point_iterator solves the fix point of a function
%
%  ::
%
%    [T,itercode,retcode]=fix_point_iterator(iterate_func,T0,options,varargin)
%
%  Args:
%
%     - iterate_func : [func_handle]: a function handle which returns 2
%       elements: The improved value of T and F0, the value of the function
%       whose discrepancy is to be minmized
%
%     - T0 : is the initial guess for the solution
%
%     - options : [struct] with fields
%       - fix_point_explosion_limit : [positive scalar |{1e+6}] : maximum
%         divergence of F0
%       - fix_point_TolFun : [positive scalar |{sqrt(eps)}] : tolerance
%         criterion
%       - fix_point_maxiter : [positive scalar |{1000}] : maximum number of
%         iterations
%       - fix_point_verbose : [true|{false}] : show iterations or not
%       - fix_point_valid_func : [function_handle|{@utils.error.valid}] :
%         function that checks the validity (nans, complex, inf, etc.) of
%         the calculations. The default function checks for complex numbers.
%         If you want to allow for complex numbers, use instead
%         @utils.error.validComplex
%
%     - varargin
%
%  Returns:
%     :
%     - T : final solution
%     - itercode : final number of iterations
%     - retcode : return code
%       - 0 : successful
%       - 21 : maximum number of iterations reached
%       - 22 : nan or inf in F0 or T
%       - 23 : divergence
%
%  Note:
%
%  Example:
%
%     See also:
%

author: jmaih


text:
function [Y, spectrum] = slgembed(G, Gc, d, fm, varargin)
%SLGEMBED Solves the general graph-based embedding
%
% $ Syntax $
%   - Y = slgembed(G, Gc, d, fm, ...)
%   - [Y, spectrum] = slgembed(G, Gc, d, fm, ...)
%
% $ Arguments $
%   - G:        The graph to be optimized
%   - Gc:       The constraint graph
%   - d:        The dimension of the embedding
%   - fm:       The type of formulation
%   - Y:        The embedding sample coordinates
%
% $ Description $
%   - Y = slgembed(G, Gc, d, fm, ...) solves the general graph-based
%     embedding of dimension d. In mathematics, it is to solve the
%     following optimization problem:
%           min/max y^T M y,  s.t. y^T C y = I
%     Based on different fm, the formulations of M and C are different:
%     fm is a cell array of two char string elements:
%           fm = {fg, fc}
%     fg indicates the formulation of M, fc indicates the formulation of
%     the constraint matrix C.
%       fg has the following different values:
%           - 'minW':   minimization using M = W, (W is the adjmat of G)
%           - 'maxW':   maximization using M = W
%           - 'minL':   minimization using M = L = D - W
%           - 'maxL':   maximization using M = L = D - W
%       fc has the following different values:
%           - 'I':      use C = I, that is y^T * y = I
%           - 'D':      use C = D, that is y^T * D * y = I (based on G)
%           - 'WC':      use C = W, adjacency matrix (based on Gc)
%           - 'LC':      use C = L = D - W, (based on Gc)
%     For example, if you specify fm = {'maxW', 'D'}, then the function
%     will solve the following optimization problem:
%           maximize y^T W y,   s.t. y^T * D * y = I
%     You can also specify the following properties:
%       - 'inv':        The parameters do to eigenvalue inverse
%                       {method, ...}
%                       (refer to slinvevals for details)
%                       This parameter take effects only when fc = 'WC" or
%                       fc = 'LC'.
%       - 'skip':       How many eigen-components to skip. default = 0.
%                       (In some algorithms, it is necessary to skip
%                        the first or first several eigen-components).
%
%   - [Y, spectrum] = slgembed(G, Gc, d, fm, ...) additionally return
%     the spectrum of the embedding, the eigenvalues of the whitened
%     C^(-1/2) M C^(-1/2).
%
% $ Remarks $
%   - The fc can be 'WC' or 'LC' only when Gc is full matrix, in current
%     version of implementation.
%
% $ History $
%   - Created by Dahua Lin, on Sep 12, 2006
%

%% parse and verify input arguments

if nargin < 4
    raise_lackinput('slgembed', 4);
end

gi = slgraphinfo(G, {'square'});
n = gi.n;
W = sladjmat(G, 'valtype', 'numeric', 'sparse', issparse(G));

[fg, fc] = deal(fm{:});

Wc = [];
if strcmp(fc, 'WC') || strcmp(fc, 'LC')
    if ~isempty(Gc)
        if ~isnumeric(Gc) || issparse(Gc) || ~isequal(size(Gc), [n n])
            error('sltoolbox:invalidarg', ...
                'In current implementation, Gc must be an n x n full numeric matrix for fc = WC or GC');
        end
        Wc = Gc;
    else
        error('sltoolbox:invalidarg', ...
            'When fc is WC or LC, Wc should be specified');
    end
end

opts.inv = {};
opts.skip = 0;
opts = slparseprops(opts, varargin{:});

if d + opts.skip > n
    error('sltoolbox:invalidarg', ...
        'The embedding dimension d plus the skip dimension should not exceed n');
end

%% main skeleton

W = W + W';
if ~isempty(Wc)
    Wc = Wc + Wc';
end

% parse formulation and decide scheme

switch fg
    case 'minW'
        Mfunc = @make_Wmat;
        optype = 'max';
    case 'maxW'
        Mfunc = @make_Wmat;
        optype = 'max';
    case 'minL'
        Mfunc = @make_Lmat;
        optype = 'min';
    case 'maxL'
        Mfunc = @make_Lmat;
        optype = 'max';
    otherwise
        error('sltoolbox:invalidarg', 'Invalid fg type: %s', fg);
end

switch fc
    case 'I'
        whfunc = [];
    case 'D'
        whfunc = @(M) whM_by_D(M, W);
        nyfunc = @(Y, T) nY_by_S(Y, T);
    case 'WC'
        whfunc = @(M) whM_by_W(M, Wc, opts.inv);
        nyfunc = @(Y, T) nY_by_T(Y, T);
    case 'LC'
        whfunc = @(M) whM_by_L(M, Wc, opts.inv);
        nyfunc = @(Y, T) nY_by_T(Y, T);
    otherwise
        error('sltoolbox:invalidarg', 'Invalid fc type: %s', fc);
end

% make M
M = Mfunc(W);

% whiten M
if ~isempty(whfunc)
    MTcell = whfunc(M);
    M = MTcell{1};
    T = MTcell{2};
    clear MTcell;
end

% optimize to solve embedding
[Y, spectrum] = optim_embed(M, optype, d, opts.skip);

% normalize the embedding to satisfy the constraint
if ~isempty(whfunc)
    Y = nyfunc(Y, T);
end
Y = Y';

%% Core routines

function M = make_Wmat(W)

M = W;

function vD = make_Dvec(W)

vD = sum(W, 1)';
if issparse(vD)
    vD = full(vD);
end

function L = make_Lmat(W)

vD = make_Dvec(W);
n = size(W, 1);
if issparse(W)
    D = sparse((1:n)', (1:n)', vD, n, n, n);
    L = D - W;
else
    L = -W;
    dinds = (1:n)'*(n+1)-n;
    L(dinds) = L(dinds) + vD;
end

function MTcell = whM_by_D(M, W)

vD = make_Dvec(W);
vD(vD < eps) = eps;
cv = 1 ./ sqrt(vD);

n = size(M, 1);
if issparse(M)
    Mw = M;
    for i = 1 : n
        Mw(:,i) = Mw(:,i) * cv(i);
    end
    for i = 1 : n
        Mw(i,:) = Mw(i,:) * cv(i);
    end
else
    Mw = slmulrowcols(M, cv', cv);
end

MTcell = {Mw, cv};

function MTcell = whM_by_W(M, W, invparams)

T = slwhiten_from_cov(W, invparams{:});
Mw = T' * M * T;
MTcell = {Mw, T};

function MTcell = whM_by_L(M, W, invparams)

L = make_Lmat(W);
MTcell = whM_by_W(M, L, invparams);

function [Y, spectrum] = optim_embed(M, optype, d, dskip)

switch optype
    case 'min'
        ord = 'ascend';
    case 'max'
        ord = 'descend';
end

d0 = d + dskip;
[spectrum, Y] = slsymeig(M, d0, ord);
if dskip > 0
    spectrum = spectrum(dskip+1:d0);
    Y = Y(:, dskip+1:d0);
end

function Y = nY_by_S(Y, T)

Y = slmulvec(Y, T, 1);

function Y = nY_by_T(Y, T)

Y = T * Y;

author: lmthang


text:
% closest larva, based on dhead2central
function [data,units,mind] = compute_closestlarva_head2central(trx,n,dosave_d)

 if nargin < 3,
   dosave_d = true;
 end

larvae = trx.exp2flies{n};
nlarvae = numel(larvae);
closestlarva = cell(1,nlarvae);
mind = cell(1,nlarvae);

for i1 = 1:nlarvae,
  larva1 = larvae(i1);
  fprintf('fly1 = %d\n',larva1);
  larvae2 = larvae;
  d = nan(numel(larvae2),trx(larva1).nframes);

  for i2 = 1:numel(larvae2),
    larva2 = larvae2(i2);
    if larva1 == larva2,
      continue;
    end
    d(i2,:) = dhead2central_pair(trx,larva1,larva2);
  end
  [mind{i1},closesti] = min(d,[],1);
  closestlarva{i1} = larvae2(closesti);
  closestlarva{i1}(isnan(mind{i1})) = nan;
end

% so that we don't compute dcenter twice
 if dosave_d,
   data = mind; %#ok<NASGU>
   units = parseunits('mm'); %#ok<NASGU>
   filename = trx.GetPerFrameFile('dhead2central',n);

   save(filename,'data','units');
 end

data = closestlarva;
units = parseunits('unit');

author: kristinbranson

repo: JAABA


text:
function seq = ScanMatch_NumToDoubleStr(num, modulus)
%SCANMATCH_NUMTODOUBLESTR transform an array of numbers into a sequence of
% double characters (each number is represented by a residue of two letter) in
% function of the modulus. The modulus is the length of the
% alphabet used for the second letter of a residue.
% i.e.  seq = ScanMatch_NumToDoubleStr([1 2 27], 26) % modulus = 26 (all the alphabet is used)
%  returns aAaBBa
% i.e.  seq = ScanMatch_NumToDoubleStr([1 2 27], 10) % modulus = 26 (all the alphabet is used)
%  returns aAaBcG
%
% The modulus can be any number between 1 and 26.
% The Input numbers have to be between 1 and 676 (which is the maximum you
% can code with two characters: 676 is xX with a modulus of 26 )
%
%   Part of the ScanMatch toolbox
%   Written by Filipe Cristino
%   $Version: 1.00 $  $Date: 10/09/2009

% ---- Check Modulus ----
if modulus > 26 || modulus < 1
    error('ScanMatch:WrongInput', 'The modulus has to be between 1 and 26!')
end

% ---- Check input numbers ----
max_num = 26 * modulus; % Max number you can convert in function of the modulus
if sum(num <1 | num > max_num)
    error('ScanMatch:WrongInput', 'The input number can not be less than 1 or more than (26 * modulus)')
end

% ---- substract one to the array as the conversion works from zero
num = num - 1;

% ---- Start conversion ----
num_l = length(num);
ind= 1;
for i=1:num_l
    seq(ind) = lower(char((fix(num(i)/modulus))+65));
    seq(ind+1) = char(rem(num(i),modulus)+65);
    ind = ind+2;
end

author: Allopart

repo: rbpf-gmapping


text:
function weights = update_weight_fcn(step_size, err_sig, delayed_signal, filter_coefficients, reset_weights)
% This function updates the adaptive filter weights based on LMS algorithm

% NOTE: To instruct Embedded MATLAB to compile an external function,
% add the following compilation directive or pragma to the function code
%#eml

fm = fimath(step_size);

step_sig = fi(step_size .* err_sig, 1, 32, 20, fm);
correction_factor = fi(delayed_signal .* step_sig, 1, 32, 20, fm);
updated_weight = fi(correction_factor + filter_coefficients, 1, 16, 16, fm);

if reset_weights
    weights = fi(zeros(1,40), 1, 16, 16, fm);
else
    weights = updated_weight;
end


text:
%% Multiple comparison correction with Threshold-Free Cluster Enhancement
%
% This example demonstrates cosmo_cluster_neighborhood and
% cosmo_montecarlo_cluster_stat
%
% Note: this example shows multiple-comparison for a single subject, but
% the same logic can be applied to a group of subjects to do a group
% analysis.
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

%% Define data
config=cosmo_config();
data_path=fullfile(config.tutorial_data_path,'ak6','s01');

targets=repmat(1:6,1,10);
chunks=floor(((1:60)-1)/6)+1;

ds = cosmo_fmri_dataset(fullfile(data_path,'glm_T_stats_perrun.nii'),...
                        'mask',fullfile(data_path, 'brain_mask.nii'), ...
                                'targets',targets,'chunks',chunks);

% There are 10 chunks, for which the data is assumed to be independent.
% Construct a dataset with 10 samples corresponding to each chunk, with
% the average value across all six targets. Each sample is considered to be
% the same condition, namely the effect of the stimulus-versus-baseline
% effect; thus all target values must be set to 1.
%
% Use either:
%   - cosmo_split and cosmo_stack
%   - cosmo_average_samples
%   - (advanced) cosmo_fx
%
% Assign the result to a variable 'ds_stim'
% >@@>
ds_stim=cosmo_fx(ds,@(x)mean(x,1),{'chunks'});

% % alternative:
%     ds_stim=ds;
%     ds_stim.sa.targets(:)=1;
%     ds_split=cosmo_split(ds_stim,{'chunks'});
%
%     for k=1:numel(ds_split)
%         ds_avg_k=cosmo_slice(ds_split{k},1);
%         ds_avg_k.samples=mean(ds_split{k}.samples,1);
%         ds_split{k}=ds_avg_k;
%     end
%
%     ds_stim=cosmo_stack(ds_split);
%
% <@@<

%% Define a cluster neighborhood for this dataset and assign the result to
% a variable 'cl_nh'.
% hint: use cosmo_cluster_neighborhood

% >@@>
cl_nh=cosmo_cluster_neighborhood(ds_stim);
% <@@<

% Show a plot with the sorted number of neighbors
% for each voxel

% >@@>
n_neighbors_per_feature=cellfun(@numel,cl_nh.neighbors);
plot(sort(n_neighbors_per_feature))
% <@@<

%% Run cosmo_montecarlo_cluster_stat

% There is one condition per chunk; all targets are set to 1.
% Thus the subsequent anaylsis is a one-sample t-test.
% Note: if this was a group analysis, then each sample (row in ds.samples)
% would contain data from one subject; each unique value in .sa.chunks
% would correspond to one subject; and each unique value in .sa.targets
% would correspond to a condition of interest.

% Since this is a one-sample t-test against a mean of zero, we set this as
% a (required) option

opt=struct();
opt.h0_mean=0;

% set the number of iterations ('niter' option).
% At least 10000 is adviced for publication-quality analyses; because that
% takes quite a while to compute, here we use 200

% >@@>
% Note: for publication-quality analyses, niter=10000 or more is
% recommended
opt.niter=200;
% <@@<

% using cosmo_montecarlo_cluster_stat, compute a map with z-scores
% against the null hypothesis of a mean of zero, corrected for multiple
% comparisons. Store the result in a variable named 'tfce_z_ds_stim'

% >@@>
tfce_z_ds_stim=cosmo_montecarlo_cluster_stat(ds_stim,cl_nh,opt);
% <@@<
cosmo_plot_slices(tfce_z_ds_stim);

%% Using the same logic, run a two-sample t-test for primates versus bugs

% >@@>
primates_insects_mask=cosmo_match(ds.sa.targets,[1 2 5 6]);
ds_primates_insects=cosmo_slice(ds, primates_insects_mask);

% set primates=1, insects=2
ds_primates_insects.sa.targets(cosmo_match(...
                            ds_primates_insects.sa.targets,[1 2]))=1;
ds_primates_insects.sa.targets(cosmo_match(...
                            ds_primates_insects.sa.targets,[5 6]))=2;

% compute average for each unique combination of targets and chunks
ds_avg_primate_insects=cosmo_average_samples(ds_primates_insects);

cl_nh=cosmo_cluster_neighborhood(ds_avg_primate_insects);

opt=struct();

% set the number of iterations.
% At least 10000 is adviced for publication-quality analyses; because that
% takes quite a while to compute, here we use 200

% Note: for publication-quality analyses, niter=10000 or more is
% recommended
opt.niter=200;

tfce_z_ds_primate_vs_insects=cosmo_montecarlo_cluster_stat(...
                                    ds_avg_primate_insects,cl_nh,opt);
% <@@<
cosmo_plot_slices(tfce_z_ds_primate_vs_insects);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function dicomrt_plotctramp(name,materials,mat_ct_up_bound,density_lo_bound,density_up_bound)
% dicomrt_plotctramp(name,materials,mat_ct_up_bound,density_lo_bound,density_up_bound)
%
% Plot ct conversion ramp used to create ct-phantom.
%
% materials			    materials name
% mat_ct_up_bound		materials upper bound (see dosxyz user's manual)
% density_lo_bound		density lower bound (see dosxyz user's manual)
% density_up_bound		density upper bound (see dosxyz user's manual)
%
% See also dicomrt_ctcreate, dicomrt_createwphantom
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Plot
figure
set(gcf,'Name',['dicomrt_plotctramp: CT conversion ramp for ',name]);
hold
for i=1:size(materials,1)
    if i==1
        plot([0 mat_ct_up_bound(i)], [0 density_up_bound(i)],'LineWidth',2);
        text(mat_ct_up_bound(i)/2,density_up_bound(i)/2,['\leftarrow', materials(i,:)],'FontSize',12,'FontWeight','bold');
    else
        locateXtext=(mat_ct_up_bound(i)-mat_ct_up_bound(i-1))./2+mat_ct_up_bound(i-1);
        locateYtext=(density_up_bound(i)-density_lo_bound(i))./2+density_lo_bound(i);
        plot([mat_ct_up_bound(i-1) mat_ct_up_bound(i)], [density_lo_bound(i) density_up_bound(i)],...
            'LineWidth',2);
        plot(mat_ct_up_bound(i-1),density_lo_bound(i),'Linestyle','none','Marker','d','MarkerFaceColor','auto');
        text(locateXtext,locateYtext,['  \leftarrow ', materials(i,:)],'FontSize',12,'FontWeight','bold');
    end
end
grid on
title('CT conversion ramp','FontSize',18);
ylabel('Density (g/cm^3)','FontSize',14);
xlabel('CT number','FontSize',14);
set(gca,'XLim',[0 max(mat_ct_up_bound)])
set(gca,'YLim',[0 max(density_up_bound)])

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [B,twom] = multiord_f(A,gamma,omega)
%MULTIORD_F  returns multilayer Newman-Girvan modularity matrix for ordered undirected layers, function handle version
% Only works for undirected networks
%
% Version: 2.2.0
% Date: Thu 11 Jul 2019 12:25:42 CEST
%
%   Input: A: Cell array of NxN adjacency matrices for each layer of an
%          ordered undirected multilayer network
%          gamma: intralayer resolution parameter
%          omega: interlayer coupling strength
%
%   Output: B: function handle where B(i) returns the ith column of
%          [NxT]x[NxT] flattened modularity tensor for the
%           multilayer network with uniform ordinal coupling (T is
%           the number of layers of the network)
%           twom: normalisation constant
%
%   Example of usage: [B,twom]=multiord_f(A,gamma,omega);
%          [S,Q]= genlouvain(B); % see iterated_genlouvain.m and
%          postprocess_ordinal_multilayer.m for how to improve output
%          multilayer partition
%          Q=Q/twom;
%          S=reshape(S,N,T);
%
%   [B,twom] = MULTIORD_F(A,GAMMA, OMEGA) with A a cell array of square
%   symmetric matrices of equal size each representing an undirected network
%   "layer" computes the multilayer Newman-Girvan modularity matrix using
%   the quality function described in Mucha et al. 2010, with intralayer
%   resolution parameter GAMMA, and with interlayer coupling OMEGA connecting
%   nearest-neighbor ordered layers.  The null model used for the quality
%   function is the Newman-Girvan null model (see e.g. Bazzi et al. for other
%   possible null models). Once the mulilayer modularity matrix is computed,
%   optimization can be performed by the generalized Louvain code GENLOUVAIN
%   or ITERATED_GENLOUVAIN. The output B can be used with other heuristics,
%   provided the same mapping is used to go from the multilayer tensor to
%   the multilayer flattened matrix. That is, the node-layer tuple (i,s)
%   is mapped to i + (s-1)*N. [Note that we can define a mapping between a
%   multilayer partition S_m stored as an N by T matrix and the corresponding
%   flattened partition S stored as an NT by 1 vector. In particular
%   S_m = reshape(S,N,T) and S = S_m(:).]
%
%
%   See also
%       genlouvain heuristics:      GENLOUVAIN, ITERATED_GENLOUVAIN
%       multilayer wrappers:        MULTICAT, MULTICATF, MULTIORD
%       other heuristics:           SPECTRAL23
%       Kernighan-Lin improvement:  KLNB
%
%   Notes:
%     The matrices in the cell array A are assumed to be square,
%     symmetric, and of equal size.  These assumptions are not checked here.
%
%     For smaller systems, it is potentially more efficient (and easier) to
%     directly use the sparse quality/modularity matrix B in MULTIORD. For
%     large systems with directed layer networks, use MULTIORDDIR_F.
%
%     This code serves as a template and can be modified for situations
%     with other wrinkles (e.g., different intralayer null models,
%     different numbers of nodes from layer-to-layer, or systems which are
%     both multiplex and longitudinal).  That is, this code is only a
%     starting point; it is by no means exhaustive.
%
%     By using this code, the user implicitly acknowledges that the authors
%     accept no liability associated with that use.  (What are you doing
%     with it anyway that might cause there to be a potential liability?!?)
%
%   References:
%     Blondel, Vincent D., Jean-Loup Guillaume, Renaud Lambiotte, and
%     Etienne Lefebvre, "Fast unfolding of communities in large networks,"
%     Journal of Statistical Mechanics: Theory and Experiment, P10008
%     (2008).
%
%     Fortunato, Santo, "Community detection in graphs," Physics Reports
%     486, 75-174 (2010).
%
%     Good, Benjamin H., Yves-Alexandre de Montjoye, and Aaron Clauset,
%     "Performance of modularity maximization in practical contexts,"
%     Physical Review E 81, 046106 (2010).
%
%     Newman, Mark E. J. and Michelle Girvan. "Finding and Evaluating
%     Community Structure in Networks", Physical Review E 69, 026113 (2004).
%
%     Mucha, Peter J., Thomas Richardson, Kevin Macon, Mason A. Porter, and
%     Jukka-Pekka Onnela. "Community Structure in Time-Dependent,
%     Multiscale, and Multiplex Networks," Science 328, 876-878 (2010).
%
%     Bazzi, Marya, Mason A. Porter, Stacy Williams, Mark McDonald, Daniel
%     J. Fenn, and Sam D. Howison. "Community Detection in Temporal
%     Multilayer Networks, with an Application to Correlation Networks",
%     MMS: A SIAM Interdisciplinary Journal 14, 1-41 (2016).
%
%     Porter, M. A., J. P. Onnela, and P. J. Mucha, "Communities in
%     networks," Notices of the American Mathematical Society 56, 1082-1097
%     & 1164-1166 (2009).
%
%   Acknowledgments:
%     Thank you to Dani Bassett, Jesse Blocher, Bruce Rogers, and Simi Wang
%     for their collaborative help which led to significant cleaning up
%     of earlier versions of our multilayer community detection codes.

if nargin<2||isempty(gamma)
    gamma=1;
end

if nargin<3
    omega=1;
end

N=length(A{1});
T=length(A);

if length(gamma)==1
    gamma=repmat(gamma,T,1);
end

ii=[]; jj=[]; vv=[];
ki=[]; kj=[]; kv=[];
twom=0;
for s=1:T
    indx=(1:N)'+(s-1)*N;
    [i,j,v]=find(A{s});
    ii=[ii;indx(i)]; jj=[jj;indx(j)]; vv=[vv;v];
    k=sum(A{s});
    mm=sum(k);
    twom=twom+mm;
    ki=[ki;indx];
    kj=[kj;ones(N,1)*s];
    kv=[kv;k(:)./mm];
end
AA = sparse(ii,jj,vv,N*T,N*T);
K=sparse(ki,kj,kv,N*T,T);
clear ii jj vv ki kj kv
kvec = full(sum(AA));
AA = AA + omega*spdiags(ones(N*T,2),[-N,N],N*T,N*T);
B = @(i) AA(:,i) - gamma(ceil(i/(N+eps)))*K(:,ceil(i/(N+eps)))*kvec(i);
twom=twom+2*N*(T-1)*omega;

end

author: GenLouvain

repo: GenLouvain


text:
function degree = dunavant_degree ( rule )

%*****************************************************************************80
%
%% DUNAVANT_DEGREE returns the degree of a Dunavant rule for the triangle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    David Dunavant,
%    High Degree Efficient Symmetrical Gaussian Quadrature Rules
%    for the Triangle,
%    International Journal for Numerical Methods in Engineering,
%    Volume 21, 1985, pages 1129-1148.
%
%    James Lyness, Dennis Jespersen,
%    Moderate Degree Symmetric Quadrature Rules for the Triangle,
%    Journal of the Institute of Mathematics and its Applications,
%    Volume 15, Number 1, February 1975, pages 19-32.
%
%  Parameters:
%
%    Input, integer RULE, the index of the rule.
%
%    Output, integer DEGREE, the polynomial degree of exactness of
%    the rule.
%
  if ( 1 <= rule & rule <= 20 )
    degree = rule;
  else

    degree = -1;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DUNAVANT_DEGREE - Fatal error!\n' );
    fprintf ( 1, '  Illegal RULE = %d\n', rule );
    error ( 'DUNAVANT_DEGREE - Fatal error!' )

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ magic, offset, msg ] = pnmmagic( fid, fmt )
%PNMMAGIC Find and identify magic number in PBM/PGM/PPM file.
%
%   [ MAGIC, OFFSET, MSG ] = PNMMAGIC( FID ) scans the file with
%   filehandle FID for a PNM file magic number.
%
%   [ ... ] = PNMMAGIC( FID, FMT ) allows a restriction on the image
%   formats to look for. FMT must be one of 'pbm', 'pgm', 'ppm' or any
%   combination if put in a list (cell array). The special string 'pnm'
%   indicates any of 'pbm', 'pgm' and 'ppm'. If FMT is omitted, 'pnm' is
%   used.
%
%   MAGIC is the magic number found and OFFSET is the offset of the
%   magic number within the file. MSG contains an error message if the
%   search for the magic number failed, otherwise MSG is empty.
%
%   If PNMMAGIC succeeds, the file position indicator will be right
%   after the magic number so there is no need to call FSEEK when the
%   rest of the image is to be read. If PNMMAGIC fails, the file
%   position indicator will be at the end of the file.
%
%   PNM is not an image format by itself but means any of PBM, PGM, PPM.

%   Don't call error() if there is something wrong with the image file,
%   only if there is something wrong with the way this program is
%   called.

%   Author:      Peter J. Acklam
%   Time-stamp:  1998-05-10 19:03:42
%   E-mail:      jacklam@math.uio.no (Internet)
%   URL:         http://www.math.uio.no/~jacklam

%
% Checks input arguments and assign default value to omitted argument.
%
error( nargchk( 1, 2, nargin ) );
if nargin < 2
   fmt = 'pnm';
end
if all( fid ~= fopen( 'all' ) )
   error( 'Invalid file identifier.' );
end

% Make sure the format variable is a list (cell array).
if ischar( fmt )
   fmt = { fmt };
end

% From the list of formats create a list of magic numbers to look for.
magics = {};
for i = 1:length(fmt)
   switch lower( fmt{i} )
      case 'pbm'                                % Portable bitmap.
         magics = { magics{:}, 'P1', 'P4' };
      case 'pgm'                                % Portable graymap.
         magics = { magics{:}, 'P2', 'P5' };
      case 'ppm'                                % Portable pixmap.
         magics = { magics{:}, 'P3', 'P6' };
      case 'pnm'                                % Portable anymap.
         magics = { 'P1', 'P2', 'P3', 'P4', 'P5', 'P6' };
      otherwise
         error( [ 'Invalid format string: ' fmt ] );
   end
end

% Initialize output arguments.
magic  = '';                    % Magic number string.
offset = [];                    % Offset of magic number within file.
msg    = '';                    % Error message string.

while 1

   % The magic number should be followed by whitespace, so we can use
   % FSCANF to get the magic number.
   [ str, count ] = fscanf( fid, '%s', 1 );

   % Look for magic number.
   if any( strcmp( str, magics ) )
      magic = str;
      offset = ftell( fid ) - 2;   % Subtract length of magic number.
      return
   end

   %
   % If a magic number was not found, start scanning from the beginning
   % of the next line of the file. We don't use FGETL to get the next
   % line since the file might be enormous and contain only one line.
   %
   blocksize = 8192;            % Number of bytes read at a time.
   while 1

      % Read a chunk of data from the file.
      [ data, count ] = fread( fid, blocksize );
      if count == 0
         msg = 'File ended while scanning for magic number.';
         return
      end

      % Look for newline characters.
      nlchars = ( data == 10 ) | ( data == 13 );
      if any( nlchars )
         nlpos = find( nlchars );       % Position of newline chars.
         nlpos = nlpos(1);              % In case there are several.
         fseek( fid, nlpos - count, 0 );        % Put FPS at newline.
         break
      end

   end

end

author: Sable

repo: mcbench-benchmarks


text:
function [pyr,pind] = buildWUpyr(im, Nsc, daub_order);

% [PYR, INDICES] = buildWUpyr(IM, HEIGHT, DAUB_ORDER)
%
% Construct a separable undecimated orthonormal QMF/wavelet pyramid
% on matrix (or vector) IM.
%
% HEIGHT specifies the number of pyramid levels to build. Default
% is maxPyrHt(IM,FILT).  You can also specify 'auto' to use this value.
%
% DAUB_ORDER: specifies the order of the daubechies wavelet filter used
%
% PYR is a vector containing the N pyramid subbands, ordered from fine
% to coarse.  INDICES is an Nx2 matrix containing the sizes of
% each subband.

% JPM, Univ. de Granada, 03/2003, based on Rice Wavelet Toolbox
% function "mrdwt" and on Matlab Pyrtools from Eero Simoncelli.

if Nsc < 1,
    display('Error: Number of scales must be >=1.');
else,

Nor = 3; % fixed number of orientations;
h = daubcqf(daub_order);
[lpr,yh] = mrdwt(im, h, Nsc+1); % performs the decomposition

[Ny,Nx] = size(im);

% Reorganize the output, forcing the same format as with buildFullSFpyr2

pyr = [];
pind = zeros((Nsc+1)*Nor+2,2);    % Room for a "virtual" high pass residual, for compatibility
nband = 1;
for nsc = 1:Nsc+1,
    for nor = 1:Nor,
        nband = nband + 1;
        band = yh(:,(nband-2)*Nx+1:(nband-1)*Nx);
        sh = (daub_order/2 - 1)*2^nsc;  % approximate phase compensation
        if nor == 1,        % horizontal
            band = shift(band, [sh 2^(nsc-1)]);
        elseif nor == 2,    % vertical
            band = shift(band, [2^(nsc-1) sh]);
        else
            band = shift(band, [sh sh]);    % diagonal
        end
        if nsc>2,
            band = real(shrink(band,2^(nsc-2)));  % The low freq. bands are shrunk in the freq. domain
        end
        pyr = [pyr; vector(band)];
        pind(nband,:) = size(band);
    end
end

band = lpr;
band = shrink(band,2^Nsc);
pyr = [pyr; vector(band)];
pind(nband+1,:) = size(band);

end

author: ricedsp


text:
% This is like robot1, except we only use a Kalman filter.
% The goal is to study how the precision matrix changes.

seed = 1;
rand('state', seed);
randn('state', seed);

if 0
  T = 20;
  ctrl_signal = [repmat([1 0]', 1, T/4) repmat([0 1]', 1, T/4) ...
		 repmat([-1 0]', 1, T/4) repmat([0 -1]', 1, T/4)];
else
  T = 60;
  ctrl_signal = repmat([1 0]', 1, T);
end

nlandmarks = 6;
if 0
  true_landmark_pos = [1 1;
		    4 1;
		    4 4;
		    1 4]';
else
  true_landmark_pos = 10*rand(2,nlandmarks);
end
if 0
figure(1); clf
hold on
for i=1:nlandmarks
  %text(true_landmark_pos(1,i), true_landmark_pos(2,i), sprintf('L%d',i));
  plot(true_landmark_pos(1,i), true_landmark_pos(2,i), '*')
end
hold off
end

init_robot_pos = [0 0]';

true_robot_pos = zeros(2, T);
true_data_assoc = zeros(1, T);
true_rel_dist = zeros(2, T);
for t=1:T
  if t>1
    true_robot_pos(:,t) = true_robot_pos(:,t-1) + ctrl_signal(:,t);
  else
    true_robot_pos(:,t) = init_robot_pos + ctrl_signal(:,t);
  end
  nn = argmin(dist2(true_robot_pos(:,t)', true_landmark_pos'));
  %true_data_assoc(t) = nn;
  %true_data_assoc = wrap(t, nlandmarks); % observe 1, 2, 3, 4, 1, 2, ...
  true_data_assoc  = sample_discrete(normalise(ones(1,nlandmarks)),1,T);
  true_rel_dist(:,t) = true_landmark_pos(:, nn) - true_robot_pos(:,t);
end

R = 1e-3*eye(2); % noise added to observation
Q = 1e-3*eye(2); % noise added to robot motion

% Create data set
obs_noise_seq = sample_gaussian([0 0]', R, T)';
obs_rel_pos = true_rel_dist + obs_noise_seq;
%obs_rel_pos = true_rel_dist;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create params for inference

% X(t) = A X(t-1) + B U(t) + noise(Q)

% [L1]  = [1     ]  * [L1]       + [0]  * Ut  + [0   ]
% [L2]    [  1   ]    [L2]         [0]          [ 0  ]
% [R ]t   [     1]    [R ]t-1      [1]          [   Q]

% Y(t)|S(t)=s  = C(s) X(t) + noise(R)
% Yt|St=1 = [1 0 -1]  * [L1]  + R
%                       [L2]
%                       [R ]

% Create indices into block structure
bs = 2*ones(1, nlandmarks+1); % sizes of blocks in state space
robot_block =  block(nlandmarks+1, bs);
for i=1:nlandmarks
  landmark_block(:,i) = block(i, bs)';
end
Xsz = 2*(nlandmarks+1); % 2 values for each landmark plus robot
Ysz = 2; % observe relative location
Usz = 2; % input is (dx, dy)

% create block-diagonal trans matrix for each switch
A = zeros(Xsz, Xsz);
for i=1:nlandmarks
  bi = landmark_block(:,i);
  A(bi, bi) = eye(2);
end
bi = robot_block;
A(bi, bi) = eye(2);
A = repmat(A, [1 1 nlandmarks]); % same for all switch values

% create block-diagonal system cov

Qbig = zeros(Xsz, Xsz);
bi = robot_block;
Qbig(bi,bi) = Q; % only add noise to robot motion
Qbig = repmat(Qbig, [1 1 nlandmarks]);

% create input matrix
B = zeros(Xsz, Usz);
B(robot_block,:) = eye(2); % only add input to robot position
B = repmat(B, [1 1 nlandmarks]);

% create observation matrix for each value of the switch node
% C(:,:,i) = (0 ... I ... -I) where the I is in the i'th posn.
% This computes L(i) - R
C = zeros(Ysz, Xsz, nlandmarks);
for i=1:nlandmarks
  C(:, landmark_block(:,i), i) = eye(2);
  C(:, robot_block, i) = -eye(2);
end

% create observation cov for each value of the switch node
Rbig = repmat(R, [1 1 nlandmarks]);

% initial conditions
init_x = zeros(Xsz, 1);
init_v = zeros(Xsz, Xsz);
bi = robot_block;
init_x(bi) = init_robot_pos;
%init_V(bi, bi) = 1e-5*eye(2); % very sure of robot posn
init_V(bi, bi) = Q; % simualate uncertainty due to 1 motion step
for i=1:nlandmarks
  bi = landmark_block(:,i);
  init_V(bi,bi)= 1e5*eye(2); % very uncertain of landmark psosns
  %init_x(bi) = true_landmark_pos(:,i);
  %init_V(bi,bi)= 1e-5*eye(2); % very sure of landmark psosns
end

%k = nlandmarks-1; % exact
k = 3;
ndx = {};
for t=1:T
  landmarks = unique(true_data_assoc(t:-1:max(t-k,1)));
  tmp = [landmark_block(:, landmarks) robot_block'];
  ndx{t} = tmp(:);
end

[xa, Va] = kalman_filter(obs_rel_pos, A, C, Qbig, Rbig, init_x, init_V, ...
				     'model', true_data_assoc, 'u', ctrl_signal, 'B', B, ...
		       'ndx', ndx);

[xe, Ve] = kalman_filter(obs_rel_pos, A, C, Qbig, Rbig, init_x, init_V, ...
				     'model', true_data_assoc, 'u', ctrl_signal, 'B', B);

if 0
est_robot_pos = x(robot_block, :);
est_robot_pos_cov = V(robot_block, robot_block, :);

for i=1:nlandmarks
  bi = landmark_block(:,i);
  est_landmark_pos(:,i) = x(bi, T);
  est_landmark_pos_cov(:,:,i) = V(bi, bi, T);
end
end

nrows = 10;
stepsize = T/(2*nrows);
ts = 1:stepsize:T;

if 1 % plot

clim = [0 max(max(Va(:,:,end)))];

figure(2)
if 0
  imagesc(Ve(1:2:end,1:2:end, T))
  clim = get(gca,'clim');
else
  i = 1;
  for t=ts(:)'
    subplot(nrows,2,i)
    i = i + 1;
    imagesc(Ve(1:2:end,1:2:end, t))
    set(gca, 'clim', clim)
    colorbar
  end
end
suptitle('exact')

figure(3)
if 0
  imagesc(Va(1:2:end,1:2:end, T))
  set(gca,'clim', clim)
else
  i = 1;
  for t=ts(:)'
    subplot(nrows,2,i)
    i = i+1;
    imagesc(Va(1:2:end,1:2:end, t))
    set(gca, 'clim', clim)
    colorbar
  end
end
suptitle('approx')

figure(4)
i = 1;
for t=ts(:)'
  subplot(nrows,2,i)
  i = i+1;
  Vd = Va(1:2:end,1:2:end, t) - Ve(1:2:end,1:2:end,t);
  imagesc(Vd)
  set(gca, 'clim', clim)
  colorbar
end
suptitle('diff')

end % all plot

for t=1:T
  i = 1:2*nlandmarks;
  denom = Ve(i,i,t) + (Ve(i,i,t)==0);
  Vd =(Va(i,i,t)-Ve(i,i,t)) ./ denom;
  Verr(t) = max(Vd(:));
end
figure(6); plot(Verr)
title('max relative Verr')

for t=1:T
  %err(t)=rms(xa(:,t), xe(:,t));
  err(t)=rms(xa(1:end-2,t), xe(1:end-2,t)); % exclude robot
end
figure(5);plot(err)
title('rms mean pos')

author: bayesnet

repo: bnt


text:
function [TW, TPW, TR, TPR] = ft_realtime_benchmark(target)

% FT_REALTIME_BENCHMARK times the reading and writing of data
%
% Use as
%   ft_realtime_benchmark(target)
% where target is the location of the buffer or the file, for
% example 'buffer://localhost:1972'
%
% Please note that any data in the target will be overwritten.

% Copyright (C) 2010, Stefan Klanke
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

N = 100;

hdr        = [];
hdr.Fs     = 1000; % this is not really used here
hdr.nChans = 64;

fprintf(1,'Writing header and flushing data in %s...\n', target);
ft_write_data(target, single([]), 'header', hdr, 'append', false);

blocklen = [8, 16, 32, 64, 128, 256, 512, 1024, 2048];
blocksize = hdr.nChans * blocklen;

TW  = zeros(N, length(blocklen));	% Times for writing
TPW = zeros(1, length(blocklen));	% Throughput writing
TR  = zeros(N, length(blocklen));	% Times for reading
TPR = zeros(1, length(blocklen));	% Throughput reading

for k=1:length(blocklen);
  fprintf(1,'\nDetermining WRITE throughput at %ix%i samples per block...\n', hdr.nChans, blocklen(k));
  % generate random data for writing
  X = single(randn(hdr.nChans, blocklen(k)));
  for n=1:N;
    tic;
    ft_write_data(target, X, 'header', hdr, 'append', true);
    TW(n,k) = toc;
  end
  mt = mean(TW(:,k));
  st = std(TW(:,k));
  ss = blocksize(k) * N / sum(TW(:,k));
  fprintf(1,'Time per block %f +/- %f  => %d samples/sec\n', mt, st, ss);
  TPW(k) = ss;

  hdr = ft_read_header(target);

  fprintf(1,'\nDetermining READ throughput at %ix%i samples per block...\n', hdr.nChans, blocklen(k));
  begS = 1;
  for n=1:N;
    endS = begS+blocklen(k)-1;
    tic;
    X = ft_read_data(target, 'header', hdr, 'begsample', begS, 'endsample', endS);
    TR(n,k) = toc;
    begS = endS+1;
  end
  mt = mean(TR(:,k));
  st = std(TR(:,k));
  ss = blocksize(k) * N / sum(TR(:,k));
  fprintf(1,'Time per block %f +/- %f  => %d samples/sec\n', mt, st, ss);

  TPR(k) = ss;
end

loglog(blocksize, TPW,'r+-', blocksize, TPR, 'b+-');
legend('Writing','Reading','Location','NorthWest');
xlabel('block size');
ylabel('samples per sec');

author: fieldtrip


text:
function a = line_adj_eigen_right ( n )

%*****************************************************************************80
%
%% LINE_ADJ_EIGEN_RIGHT returns the right eigenvectors of the LINE_ADJ matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the right eigenvector matrix.
%
  a = zeros ( n, n );

  for i = 1: n
    for j = 1 : n
      angle = ( i * j ) * pi / ( n + 1 );
      a(i,j) = sqrt ( 2.0 / ( n + 1 ) ) * sin ( angle );
    end
  end

  return
end


text:
function [ bootql ] = offpolicy_eval_tdlearning( qldata3, physpol, gamma, num_iter )
% V value averaged over state population
% hence the difference with mean(V) stored in recqvi(:,3)

ncl=size(physpol,1)-2;
bootql=cell(num_iter,1);
p=unique(qldata3(:,8));
prop=5000/numel(p); %5000 patients of the samples are used
prop=min([prop 0.75]);  %max possible value is 0.75 (75% of the samples are used)

ii=qldata3(:,1)==1;
a=qldata3(ii,2);
d=zeros(ncl,1);
 for i=1:ncl
  d(i)=sum(a==i);    % intitial state disctribution
 end

fprintf('Progress of Q-Learning:\n');
fprintf(['\n' repmat('.',1,num_iter) '\n\n']);

parfor i=1:num_iter
fprintf('\b|\n');

ii=floor(rand(size(p,1),1)+prop);     % select a random sample of trajectories
j=ismember(qldata3(:,8),p(ii==1));
q=qldata3(j==1,1:4);

[Qoff, ~]=OffpolicyQlearning150816( q , gamma, 0.1, 300000);

V=zeros(750,25);
for k=1:750
    for j=1:25
        V(k,j)=physpol(k,j)*Qoff(k,j);
    end
end

Vs =sum(V')';
bootql(i)={nansum(Vs(1:750).*d)/sum(d)};

% Vs=nansum((physpol.*Qoff)')';
% bootql(i)={sum(Vs(1:ncl).*d)/sum(d)};
end

bootql=cell2mat(bootql);

end

author: matthieukomorowski


text:
 % Right Stance Domain
 %
 % Contact: Right Toe
function domain = DoubleStance(model, load_path)
    % construct the right stance domain of Cassie
    %
    % Parameters:
    % model: the right body model of Cassie robot

    %% first make a copy of the robot model
    %| @note Do not directly assign with the model variable, since it is a
    %handle object.
    domain = copy(model);
    % set the name of the new copy
    domain.setName('DoubleStance');

    if nargin < 2
        load_path = [];
    end

    %% Add contact
    % left foot point contact
    [right_foot, fric_coef,geometry] = sys.frames.RightSole(model);
    right_foot_contact = ToContactFrame(right_foot,...
        'PlanarContactWithFriction');
    domain = addContact(domain,right_foot_contact,fric_coef,geometry,load_path);

    [left_foot, fric_coef,geometry] = sys.frames.LeftSole(model);
    left_foot_contact = ToContactFrame(left_foot,...
        'PlanarContactWithFriction');
    domain = addContact(domain,left_foot_contact,fric_coef,geometry,load_path);

    %% Add event
    % the force on the left sole
    fLeftSole = domain.Inputs.ConstraintWrench.fLeftSole;
    left_lift = UnilateralConstraint(domain,fLeftSole(3),'leftFootLift','fLeftSole');
    domain = addEvent(domain, left_lift);

    fRightSole = domain.Inputs.ConstraintWrench.fRightSole;
    right_lift = UnilateralConstraint(domain,fRightSole(3),'rightFootLift','fRightSole');
    domain = addEvent(domain, right_lift);

    %% Add Virtual Constraints
    %     x = domain.States.x;
    %     dx = domain.States.dx;

    % phase variable: linearized hip position
    %     r_hip_link = domain.Links(getLinkIndices(domain, 'l_uleg'));
    %     r_hip_frame = r_hip_link.Reference;
    %     r_hip_frame = domain.Joints(getJointIndices(domain,'r_leg_hpy'));
    %     r_ankle_frame = domain.Joints(getJointIndices(domain,'r_leg_aky'));
    %     p_rhp = getCartesianPosition(domain, r_hip_frame);
    %     p_rak = getCartesianPosition(domain, r_ankle_frame);
    %     % p_sf = getCartesianPosition(domain, domain.ContactPoints.RightSole);
    %
    %
    %     p_hip = p_rhp(1) - p_rak(1);
    %     deltaPhip = linearize(p_hip,x);

    % @state based phase variable
    %     p = SymVariable('p',[2,1]);
    %     tau = (deltaPhip-p(2))/(p(1)-p(2));

    % @time based phase variable
    %     t = SymVariable('t');
    %     p = SymVariable('p',[2,1]);
    %     tau = (t-p(1))/(p(2)-p(1));
    %
    %
    %
    %     % relative degree two outputs:
    %     y_rap = x('r_leg_aky');
    %     y_rar = x('r_leg_akx');
    %     y_rkp = x('r_leg_kny');
    %     y_rtp = -x('r_leg_aky') - x('r_leg_kny') - x('r_leg_hpy');
    %     y_rtr = -x('r_leg_akx') - x('r_leg_hpx');
    %     y_rhy = x('r_leg_hpz');
    %     y_lkp = x('l_leg_kny');
    %     y_llr = x('r_leg_hpx') - x('l_leg_hpx');
    %
    %     l_hip_frame = domain.Joints(getJointIndices(domain,'l_leg_hpy'));
    %     l_ankle_frame = domain.Joints(getJointIndices(domain,'l_leg_aky'));
    %     p_lhp = getCartesianPosition(domain, l_hip_frame);
    %     p_lak = getCartesianPosition(domain, l_ankle_frame);
    %     % p_lak = getCartesianPosition(domain, l_ankle_frame);
    %     nsl = (p_lak(1) - p_lhp(1))/(p_lak(3) - p_lhp(3));
    %     y_nsl = linearize(nsl,x);
    %
    %
    %     left_sole_inside = sys.frames.LeftSoleInside(model);
    %     left_sole_outside = sys.frames.LeftSoleOutside(model);
    %     left_toe = sys.frames.LeftToe(model);
    %     left_heel = sys.frames.LeftHeel(model);
    %     p_lsi = getCartesianPosition(domain, left_sole_inside);
    %     p_lso = getCartesianPosition(domain, left_sole_outside);
    %     p_ltoe = getCartesianPosition(domain, left_toe);
    %     p_lheel = getCartesianPosition(domain, left_heel);
    %
    %     y_lfr = p_lsi(3) - p_lso(3);
    %     y_lfp = p_lheel(3) - p_ltoe(3);
    %     y_lfy = p_lheel(2) - p_ltoe(2);
    %
    %     ya_2 = [y_rap;
    %         y_rar;
    %         y_rkp;
    %         y_rtp;
    %         y_rtr;
    %         y_rhy;
    %         y_lkp;
    %         y_nsl;
    %         y_llr;
    %         y_lfr;
    %         y_lfp;
    %         y_lfy];
    %
    %     % optional
    %     y2_label = {'RightAnklePitch',...
    %         'RightAnkleRoll',...
    %         'RightKneePitch',...
    %         'RightTorsoPitch',...
    %         'RightTorsoRoll',...
    %         'RightHipYaw',...
    %         'LeftKneePitch',...
    %         'LeftLinNSlopeWithBase',...
    %         'LeftLegRoll',...
    %         'LeftFootCartX',...
    %         'LeftFootCartY',...
    %         'LeftFootCartZ'};
    %
    %     y2 = VirtualConstraint(domain,ya_2,'position','DesiredType','Bezier','PolyDegree',5,...
    %         'RelativeDegree',2,'OutputLabel',{y2_label},'PhaseType','TimeBased',...
    %         'PhaseVariable',tau,'PhaseParams',p,'Holonomic',true,'LoadPath',load_path);
    %
    %
    %     domain = addVirtualConstraint(domain,y2);
end

author: ayonga

repo: frost-dev

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


text:
function [h, compUp] =  lfmComputeH3(gamma1_p, gamma1_m, sigma2, t1,t2,preFactor,...
    mode, term)
% LFMCOMPUTEH3 Helper function for computing part of the LFM kernel.
% FORMAT
% DESC computes a portion of the LFM kernel.
% ARG gamma1 : Gamma value for first system.
% ARG gamma2 : Gamma value for second system.
% ARG sigma2 : length scale of latent process.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% ARG mode: indicates in which way the vectors t1 and t2 must be transposed
% RETURN h : result of this subcomponent of the kernel for the given values.
%
% COPYRIGHT : Mauricio Alvarez, 2008
%
% MODIFICATIONS : Neil D. Lawrence, 2007
%
% SEEALSO : lfmKernParamInit, lfmXlfmKernCompute

% KERN

% Evaluation of h

if nargin<8
    term =[];
end

if ~mode
    if ~term
        if nargout >1
            compUp = lfmComputeUpsilonMatrix(gamma1_p,sigma2, t1,t2);
            h = preFactor*compUp;
        else
            h = preFactor*lfmComputeUpsilonMatrix(gamma1_p,sigma2, t1,t2);
        end
    else
        if nargout > 1
            compUp = lfmComputeUpsilonMatrix(gamma1_p,sigma2, t1,t2);
            h = -preFactor(1)*compUp + preFactor(2)*conj(compUp);
        else
            upsilon = lfmComputeUpsilonMatrix(gamma1_p,sigma2, t1,t2);
            h = -preFactor(1)*upsilon + preFactor(2)*conj(upsilon);
        end
    end
else
    if nargout>1
        compUp = cell(2,1);
        compUp{1} = lfmComputeUpsilonMatrix(gamma1_p,sigma2, t1,t2);
        compUp{2} = lfmComputeUpsilonMatrix(gamma1_m,sigma2, t1,t2);
        h = preFactor(1)*compUp{1} + preFactor(2)*compUp{2};
    else
        h = preFactor(1)*lfmComputeUpsilonMatrix(gamma1_p,sigma2, t1,t2) ...
            + preFactor(2)*lfmComputeUpsilonMatrix(gamma1_m,sigma2, t1,t2);
    end
end

author: SheffieldML

repo: GPmat


text:
function [outputMovie] = createSideBySide(primaryMovie,secondaryMovie,varargin)
	% Auto-create side-by-side, either save output as hdf5 or return a matrix
	% Biafra Ahanonu
	% started: 2014.01.04 (code taken from controllerAnalysis)
	% inputs
		% primaryMovie - string pointing to the video file (.avi, .tif, or .hdf5 supported, auto-detects based on extension) OR a matrix
		% secondaryMovie - string pointing to the video file (.avi, .tif, or .hdf5 supported, auto-detects based on extension) OR a matrix
	% outputs
		% outputMovie - horizontally concatenated movie
	% changelog
		% 2014.01.27 [22:57:19] - changed to allow input of either a path or a matrix, more generalized
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
	% TODO
		% allow option for vertical concat?
		% spatio-temporal should be one abstracted for-loop, no? - DONE 2014.01.27 [21:52:37]

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	% ========================
	% number of frames in each movie to load, [] = all, 1:500 would be 1st to 500th frame.
	options.frameList = [];
	% if want a different set of frames loaded for the second movie
	options.frameListSecondary = [];
	% whether to convert movie to double on load, not recommended
	options.convertToDouble = 0;
	% name of HDF5 dataset name to load
	options.inputDatasetName = '/1';
	% string to a movie, preferably AVI
	options.recordMovie = 0;
	% amount of pixels around the border to crop in primary movie
	options.pxToCrop = [];
	% downsample combined movie
	options.downsampleFactorFinal = 1;
	% make movies equal by using NaN instead of temporal downsampling.
	options.makeTimeEqualUsingNans = 0;
	% 1 = yes, 0 = no. whether to normalize each movies pixel distribution between [0 1]
	options.normalizeMovies = 1;
	% rotate xy dims of second movie
	options.rotateSecondMovie = 0;
	% rotate xy dims of primary movie
	options.rotatePrimaryMovie = 0;
	% Binary: 1 = increase to largest movie instead of downsampling to the smallest movie.
	options.increaseToLargestMovie = 0;
	% Binary: 1 = display info
	options.displayInfo = 1;
	% get options
	options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	try
		% ========================
		% get the movie
		if strcmp(class(primaryMovie),'char')|strcmp(class(primaryMovie),'cell')
			primaryMovie = loadMovieList(primaryMovie,'convertToDouble',options.convertToDouble,'frameList',options.frameList,'inputDatasetName',options.inputDatasetName);
		end

		% load secondary movie
		if isempty(options.frameListSecondary)
			if strcmp(class(secondaryMovie),'char')|strcmp(class(secondaryMovie),'cell')
				secondaryMovie = loadMovieList(secondaryMovie,'convertToDouble',options.convertToDouble,'frameList',options.frameList);
			end
		else
			if strcmp(class(secondaryMovie),'char')|strcmp(class(secondaryMovie),'cell')
				secondaryMovie = loadMovieList(secondaryMovie,'convertToDouble',options.convertToDouble,'frameList',options.frameListSecondary);
			end
		end
		% ========================
		% If movies are 4D rgb ([x y t c] where c = 3 for RGB colors) then call function again again for each dimension then concat the resulting individual RGB movies into a single output movie
		primaryDimNum = length(size(primaryMovie));
		secondaryDimNum = length(size(secondaryMovie));
		if primaryDimNum==4|secondaryDimNum==4
			subfxnDisp('4D movie detected...')

			% get size of 4th dim depending on which input is 4D
			if secondaryDimNum==4&primaryDimNum~=4
				nColorDims = size(secondaryMovie,4);
			else
				nColorDims = size(primaryMovie,4);
			end
			sideBySideMatrix = {};
			for colorNo = 1:nColorDims
				subfxnDisp(repmat('=',1,21))
				subfxnDisp(['color: ' num2str(colorNo) '/' num2str(nColorDims)])
				if secondaryDimNum==4
					if primaryDimNum==4
						sideBySideMatrix{colorNo} = createSideBySide(squeeze(primaryMovie(:,:,:,colorNo)),squeeze(secondaryMovie(:,:,:,colorNo)),'options',options);
					else
						sideBySideMatrix{colorNo} = createSideBySide(primaryMovie,squeeze(secondaryMovie(:,:,:,colorNo)),'options',options);
					end
				else
					sideBySideMatrix{colorNo} = createSideBySide(squeeze(primaryMovie(:,:,:,colorNo)),secondaryMovie,'options',options);
				end
			end
			outputMovie = cat(4,sideBySideMatrix{:});
			return;
		end

		if options.rotatePrimaryMovie==1
			subfxnDisp('rotating...')
			subfxnDisp(['pre-rotation dims: ' num2str(size(primaryMovie))])
			primaryMovie = permute(primaryMovie,[2 1 3]);
			subfxnDisp(['post-rotation dims: ' num2str(size(primaryMovie))])
		end
		if options.rotateSecondMovie==1
			subfxnDisp('rotating...')
			subfxnDisp(['pre-rotation dims: ' num2str(size(secondaryMovie))])
			secondaryMovie = permute(secondaryMovie,[2 1 3]);
			subfxnDisp(['post-rotation dims: ' num2str(size(secondaryMovie))])
		end

		% ========================
		% make movies single for calculations sake
		% primaryMovie = single(primaryMovie);
		% secondaryMovie = single(secondaryMovie);

		% ========================
		subfxnDisp('cropping primary movie...')
		% Get the x and y corner coordinates as integers
		if ~isempty(options.pxToCrop)
			if size(primaryMovie,2)>=size(primaryMovie,1)
				coords(1) = options.pxToCrop; %xmin
				coords(2) = options.pxToCrop; %ymin
				coords(3) = size(primaryMovie,1)-options.pxToCrop;   %xmax
				coords(4) = size(primaryMovie,2)-options.pxToCrop;   %ymax
			else
				coords(1) = options.pxToCrop; %xmin
				coords(2) = options.pxToCrop; %ymin
				coords(4) = size(primaryMovie,1)-options.pxToCrop;   %xmax
				coords(3) = size(primaryMovie,2)-options.pxToCrop;   %ymax
			end
			rowLen = size(primaryMovie,1);
			colLen = size(primaryMovie,2);
			% a,b are left/right column values
			a = coords(1);
			b = coords(3);
			% c,d are top/bottom row values
			c = coords(2);
			d = coords(4);
			cropChoice = 2;
			switch cropChoice
				case 1
					primaryMovie(1:rowLen,1:a,:) = NaN;
					primaryMovie(1:rowLen,b:colLen,:) = NaN;
					primaryMovie(1:c,1:colLen,:) = NaN;
					primaryMovie(d:rowLen,1:colLen,:) = NaN;
				case 2
					primaryMovie = primaryMovie(coords(2):coords(4), coords(1): coords(3),:);
				otherwise

			end
		end
		% ========================
		subfxnDisp('making movies spatially and temporally identical...')
		% to generalize out downsampling, create cell arrays to call that contain the dimension information
		dimensionList = {3, 2, 1};
		dimensionNameList = {'time','space','space'};
		% loop over each of the dimensions to resize to
		loopList = [1 3]; %1:length(dimensionList)
		for i=loopList
			thisDim = dimensionList{i};
			thisDimName = dimensionNameList{i};
			lengthPrimary = size(primaryMovie,thisDim);
			lengthSecond = size(secondaryMovie,thisDim);
			subfxnDisp(['length primary: ' num2str(lengthPrimary)])
			subfxnDisp(['length secondary: ' num2str(lengthSecond)])
			if options.makeTimeEqualUsingNans==1&strcmp(thisDimName,'time')
				if lengthPrimary>lengthSecond
					subfxnDisp('adding NaNs to end of second movie...')
					movieDiff = lengthPrimary-lengthSecond;
					secondaryMovie(:,:,end+movieDiff) = NaN;
				elseif lengthSecond>lengthPrimary
					subfxnDisp('adding NaNs to end of first movie...')
					movieDiff = lengthSecond-lengthPrimary;
					primaryMovie(:,:,end+movieDiff) = NaN;
				end
				continue
			end
			if lengthPrimary>lengthSecond
				% display(['downsampling: ' movieList{1}]);
				if options.increaseToLargestMovie==1
					downsampleFactor = lengthSecond/lengthPrimary;
					secondaryMovie = downsampleMovie(secondaryMovie,'downsampleDimension',thisDimName,'downsampleFactor',downsampleFactor,'downsampleZ',lengthPrimary,'downsampleX',lengthPrimary);
				else
					downsampleFactor = lengthPrimary/lengthSecond;
					primaryMovie = downsampleMovie(primaryMovie,'downsampleDimension',thisDimName,'downsampleFactor',downsampleFactor,'downsampleZ',lengthSecond,'downsampleX',lengthSecond);
				end
			elseif lengthSecond>lengthPrimary
				% display(['downsampling: ' vidList{1}]);
				if options.increaseToLargestMovie==1
					downsampleFactor = lengthPrimary/lengthSecond;
					primaryMovie = downsampleMovie(primaryMovie,'downsampleDimension',thisDimName,'downsampleFactor',downsampleFactor,'downsampleZ',lengthSecond,'downsampleX',lengthSecond);
				else
					downsampleFactor = lengthSecond/lengthPrimary;
					secondaryMovie = downsampleMovie(secondaryMovie,'downsampleDimension',thisDimName,'downsampleFactor',downsampleFactor,'downsampleZ',lengthPrimary,'downsampleX',lengthPrimary);
				end
			end
		end
		%========================
		% normalize movies between 0 and 1 so they display correctly together, better if their distributions are the same
		if options.normalizeMovies==1
			subfxnDisp('normalizing movies...')
			[primaryMovie] = normalizeVector(single(primaryMovie),'normRange','zeroToOne');
			[primaryMovie] = normalizeMovie(primaryMovie,'normalizationType','meanSubtraction');
			[secondaryMovie] = normalizeVector(single(secondaryMovie),'normRange','zeroToOne');
			[secondaryMovie] = normalizeMovie(secondaryMovie,'normalizationType','meanSubtraction');
		end
		% ========================
		% horizontally concat the movies
		subfxnDisp('concatenating movies...')
		subfxnDisp(['size primary: ' num2str(size(primaryMovie))])
		subfxnDisp(['size secondary: ' num2str(size(secondaryMovie))])
		outputMovie = horzcat(primaryMovie,secondaryMovie);
		clear primaryMovie secondaryMovie
		if options.downsampleFactorFinal>1
			subfxnDisp('downsampling final movie...')
			outputMovie = downsampleMovie(outputMovie,'downsampleDimension','space','downsampleFactor',options.downsampleFactorFinal);
		end

		% if user ask to save movie, do so.
		if options.recordMovie~=0
			writerObj = VideoWriter(options.recordMovie);
			open(writerObj);
			nFrames = size(outputMovie,3);
			reverseStr = '';
			for frame=1:nFrames
				writeVideo(writerObj,squeeze(outputMovie(:,:,frame)));
				if frame==1||mod(frame,5)==0||frame==nFrames
					reverseStr = cmdWaitbar(frame,nFrames,reverseStr,'inputStr','writing movie');drawnow;
				end
			end
			close(writerObj);
		end
	catch err
		display(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		display(repmat('@',1,7))
	end
	% subfxnDisp: function description
	function subfxnDisp(txt)
		if options.displayInfo==1
			display(txt)
		end
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah