text:
function texturedVRML2(outputFname,vertex, face, textureUV, textureImageFile)
%{
  outputFname is the file name for VRML2 file. usualy it is *.wrl
      when it is done, open the file in meshlab
  vertex is 3xN points
  face is KxM faces (index for faces with K vertices and M face in total.
        K=3 for triangulated mesh. K=4 for rectangular mesh.)
  textureUV is KxMx2 uv coordinates for each vertice of the faces
        (all values should be from 0 to 1)
  texture coordinate system

  V (0,1)      (1,1)
  ^
  |
  |
  |
  +--------> U (1,0)
 (0,0)

  more document: http://www.c3.hu/cryptogram/vrmltut/part6.html
%}

fp = fopen(outputFname,'w');
fprintf(fp,'#VRML V2.0 utf8\n');
fprintf(fp,'Shape{ geometry IndexedFaceSet{\n');
fprintf(fp,'coord Coordinate{point[\n');
for i=1:size(vertex,2)
    fprintf(fp,'%f %f %f,\n',vertex(1,i),vertex(2,i),vertex(3,i));
end
fprintf(fp,']}\n');
fprintf(fp,'coordIndex[\n');
for i=1:size(face,2)
    for j=1:size(face,1)
        fprintf(fp,'%d,',face(j,i));
    end
    fprintf(fp,'-1,\n');
end
fprintf(fp,']\n');
fprintf(fp,'texCoord TextureCoordinate{point[\n');
for i=1:size(textureUV,2)
    for j=1:size(textureUV,1)
        fprintf(fp,'%f %f\n',textureUV(j,i,1),textureUV(j,i,2));
    end
end
fprintf(fp,']}\n');
fprintf(fp,'texCoordIndex[\n');
cnt =0;
for i=1:size(textureUV,2)
    for j=1:size(textureUV,1)
        fprintf(fp,'%d ',cnt);
        cnt=cnt+1;
    end
    fprintf(fp,'-1\n');
end
fprintf(fp,']\n');
fprintf(fp,'}\n');
fprintf(fp,'appearance Appearance{material Material{diffuseColor 1 1 1 }  texture ImageTexture { url "%s" } }\n',textureImageFile);
fprintf(fp,'}\n');
fclose(fp);
end

author: jianxiongxiao

repo: ProfXkit


text:
function [ n_data, x, p, cdf ] = geometric_cdf_values ( n_data )

%*****************************************************************************80
%
%% GEOMETRIC_CDF_VALUES returns values of the geometric CDF.
%
%  Discussion:
%
%    The geometric or Pascal probability density function gives the
%    probability that the first success will happen on the X-th Bernoulli
%    trial, given that the probability of a success on a single trial is P.
%
%    The value of CDF ( X, P ) is the probability that the first success
%    will happen on or before the X-th trial.
%
%    In Mathematica, the function can be evaluated by:
%
%      Needs["Statistics`DiscreteDistributions`]
%      dist = GeometricDistribution [ p ]
%      CDF [ dist, x ]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%    Daniel Zwillinger and Stephen Kokoska,
%    CRC Standard Probability and Statistics Tables and Formulae,
%    Chapman and Hall / CRC Press, 2000.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer X, the number of trials.
%
%    Output, real P, the probability of success
%    on one trial.
%
%    Output, real CDF, the cumulative density function.
%
  n_max = 14;

  cdf_vec = [ ...
     0.1900000000000000E+00, ...
     0.2710000000000000E+00, ...
     0.3439000000000000E+00, ...
     0.6861894039100000E+00, ...
     0.3600000000000000E+00, ...
     0.4880000000000000E+00, ...
     0.5904000000000000E+00, ...
     0.9141006540800000E+00, ...
     0.7599000000000000E+00, ...
     0.8704000000000000E+00, ...
     0.9375000000000000E+00, ...
     0.9843750000000000E+00, ...
     0.9995117187500000E+00, ...
     0.9999000000000000E+00 ];

  p_vec = [ ...
     0.1E+00, ...
     0.1E+00, ...
     0.1E+00, ...
     0.1E+00, ...
     0.2E+00, ...
     0.2E+00, ...
     0.2E+00, ...
     0.2E+00, ...
     0.3E+00, ...
     0.4E+00, ...
     0.5E+00, ...
     0.5E+00, ...
     0.5E+00, ...
     0.9E+00 ];

  x_vec = [ ...
    1,  2,  3, 10, 1, ...
    2,  3, 10,  3, 3, ...
    3,  5, 10,  3 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0;
    p = 0.0;
    cdf = 0.0;
  else
    x = x_vec(n_data);
    p = p_vec(n_data);
    cdf = cdf_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function	X = make_nn_diag(N,K)
% Make N x N matrix with K-principal dimension
% X = make_nn_diag(N,K)
Adiag = 0.2;

Z = rand(N,K);
Z = repmultiply( Z, 1./sqrt(sum(Z.^2,1)) );
X = Adiag * eye(N) + Z * Z';

return
%%%--- END ---%%%

% Make N x N matrix whose n-th diag elements have correlation with neighbor
% X(n,n+k) = 1, k=-K:K

X  = speye(N,N);

%ix = randperm(N);
%ix = ix(1:K);

Nd = max(fix(N/K),1);
ix = 1:Nd:N;

for n = ix
	%a = exp(-(n/NN).^2);
%	X1= spdiags(a*ones(N,1),n,N,N);
	X1= spdiags(rand(N,1),n,N,N);
%	X = X + X1 ;
	X = X + X1 + X1';
end

X = X ./ repmat( 1./sqrt(sum(X.^2,2)) , [1, N]);

author: KamitaniLab

repo: GenericObjectDecoding


text:
function varargout = area(f, varargin)
%AREA   Filled CHEBFUN area plot.
%   AREA(X, F) or AREA(F) is the same as PLOT(X, F) or PLOT(F) except that the
%   area between 0 and F is filled. When F is array-valued, AREA(F) plots the
%   columns of Y as filled areas.
%
%   AREA(F, LEVEL) specifies the base level for the area plot to be at the value
%   y = LEVEL. The default value is LEVEL = 0.
%
%   AREA(..., 'Prop1', VALUE1, 'Prop2', VALUE2,...) sets the specified
%   properties of the underlying areaseries objects.
%
%   AREA(AX, ...) plots into axes with handle AX. Use GCA to get the handle to
%   the current axes or to create one if none exist.
%
%   H = AREA(...) returns a vector of handles to areaseries objects.
%
% See also PLOT.

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

% TODO: Quasimatrix support?

% Deal with empty case:
if ( isempty(f) )
    varargout{1} = plot(); %#ok<LTARG>
    return
end

if ( numel(f) > 1 )
    try
        f = quasi2cheb(f);
    catch ME
        error('CHEBFUN:CHEBFUN:area:quasi', ...
            'AREA() does not support quasimatrices.');
    end
end

% Get the (x, y) data from PLOTDATA():
if ( (nargin > 1) && isa(varargin{1}, 'chebfun') )
    g = varargin{1};
    if ( numel(g) > 1 )
        try
            g = quasi2cheb(g);
        catch ME
            error('CHEBFUN:CHEBFUN:area:quasi', ...
                'AREA() does not support quasimatrices.');
        end
    end
    varargin(1) = [];
    data = plotData(f, g);
else
    data = plotData(f);
end

x = data.xLine;
y = data.yLine;

% Remove NaNs from jumps:
if ( isnan(x(1)) )
    x(1) = x(2);
    y(1,:) = y(2,:);
end
if ( isnan(x(end)) )
    x(end) = x(end-1);
    y(end,:) = y(end-1,:);
end
mask = find(isnan(x));
% Take average of left/right for y:
y(mask,:) = (y(mask+1,:) + y(mask-1,:))/2;
% Take subsequent point for x:
x(mask) = x(mask+1);

% Call built-in AREA():
h = area(x, y, 'LineStyle', 'None', varargin{:});

% Output handle to plot window:
if ( nargout == 1 )
    varargout{1} = h;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function test_pull1271

% MEM 5gb
% WALLTIME 00:20:00
% DEPENDENCY ft_read_data ft_read_header ft_read_event

% the DCCNPATH function will take the file from the present working directory if possible
% otherwise it will search on /home/common/matlab (a DCCN specific linux network drive)
% and if that fails it will search on H:\common\matlab (the windows version of that network drive)
%
% cd '/Users/alan/Dropbox/Lab-BML/Code/[2020.04]Fieldtrip_Neuroomega_nsx_PR';

fname_nsx         = dccnpath('/home/common/matlab/fieldtrip/data/test/original/eeg/blackrock/test_data_blackrock.ns2');
fname_neuroomega  = dccnpath('/home/common/matlab/fieldtrip/data/test/original/eeg/neuroomega/LT1D0.000F0000.mat');

%% loading blackrock header and data
hdr_nsx=ft_read_header(fname_nsx);
ft_read_data(fname_nsx,'header',hdr_nsx);

%% loading neuroomega header, events and data
hdr_neuroomega=ft_read_header(fname_neuroomega,'chantype','micro');
ft_read_data(fname_neuroomega,'header',hdr_neuroomega);
ft_read_event(fname_neuroomega);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% fft_compress.m

% 作用：使用离散傅里叶变换对输入的 JPG 图片 image 按照指定的滤波比例 ratio 进行压缩
% 返回值：返回低通滤波之后的频率分布和压缩之后的图片
function [z, k] = fft_compress(image, ratio)

    % 对原 JPG 图片在三个颜色上分别做二维离散傅里叶变换，得到频率分布
    z(:,:,1) = fft2(image(:,:,1));
    z(:,:,2) = fft2(image(:,:,2));
    z(:,:,3) = fft2(image(:,:,3));

    % 获取图片的尺寸大小
    [a, b, ~] = size(image);

    % 低通滤波
    for i = 1 : a
        for j = 1 : b
            if (i + j > (a+b) * ratio)
                z(i, j, 1) = 0;
                z(i, j, 2) = 0;
                z(i, j, 3) = 0;
            end
        end
    end

    % 对过滤之后的结果在三个颜色上分别做进行二维反离散傅里叶变换
    k(:,:,1) = ifft2(z(:,:,1));
    k(:,:,2) = ifft2(z(:,:,2));
    k(:,:,3) = ifft2(z(:,:,3));

    % 类型转换，转换为 0-255 范围内的颜色值
    k = uint8(k);
end

author: yandld


text:
function xdot = MathieuEquation(t,x)
% Mathieu Equation is y''(z)+eta.y'(z)+(a-2qcos(2z))sin(y) = 0
% Written into two first order differential equations
% y'(z) = x
% x'(z) = -eta.y'(z)-(a-2qcos(2z))sin(y)
%
n = length(x) ;
xdot=zeros(n,1);
theta = x(1) ;
Dtheta = x(2) ;
q = x(3) ;
a = x(4) ;
eta = x(5) ;
xdot(1) = Dtheta;
xdot(2) = -eta*Dtheta-(a-2*q*cos(2*t))*sin(theta);

author: Sable

repo: mcbench-benchmarks


text:
function [cc] = chaincode(b,unwrap)
%   Freeman Chain Code
%
%   Description: Give Freeman chain code 8-connected representation of a
%                boundary
%   Author.....: Alessandro Mannini <mannini@esod.it>
%   Date.......: 2010, september
%
%   usage:
%   --------------------------------------------------------
%   [cc] = chaincode(b,u)
%
%   INPUT:
%   --------------------------------------------------------
%   b       - boundary as np-by-2 array;
%             np is the number of pixels and each element is a pair (y,x) of
%             pixel coordinates
%   unwrap  - (optional, default=false) unwrap code;
%             if enable phase inversions are eliminated
%
%
%   OUTPUT:
%   --------------------------------------------------------
%   cc is structure with the following fields:
%
%   cc.code     - 8-connected Freeman chain code as 1-by-np array (or
%                 1-by-(np-1) if the boundary isn't close)
%   cc.x0,cc.y0 - respectively the abscissa and ordinate of start point
%   cc.ucode    - unwrapped 8-connected Freeman chain code (if required)
%

%
%
%   used direction-to-code convention is:       3  2  1
%                                                \ | /
%                                             4 -- P -- 0
%                                                / | \
%                                               5  6  7
%
%   and in terms of deltax,deltay if next pixel compared to the current:
%   --------------------------
%   | deltax | deltay | code |
%   |------------------------|
%   |    0   |   +1   |   2  |
%   |    0   |   -1   |   6  |
%   |   -1   |   +1   |   3  |
%   |   -1   |   -1   |   5  |
%   |   +1   |   +1   |   1  |
%   |   +1   |   -1   |   7  |
%   |   -1   |    0   |   4  |
%   |   +1   |    0   |   0  |
%   --------------------------
%

% check input arguments
if nargin>2
    error('Too many arguments');
elseif nargin==0
    error('Too few arguments');
elseif nargin==1
    unwrap=false;
end
% compute dx,dy  by a circular shift on coords arrays by 1 element
sb=circshift(b,[-1 0]);
delta=sb-b;
% check if boundary is close, if not cut last element
if abs(delta(end,1))>1 || abs(delta(end,2))>1
    delta=delta(1:(end-1),:);
end
% check if boundary is 8-connected
n8c=find(abs(delta(:,1))>1 | abs(delta(:,2))>1);
if size(n8c,1)>0
    s='';
    for i=1:size(n8c,1)
        s=[s sprintf(' idx -> %d \n',n8c(i))];
    end
    error('Curve isn''t 8-connected in elements: \n%s',s);
end

% convert dy,dx pairs to scalar indexes thinking to them (+1) as base-3 numbers
% according to: idx=3*(dy+1)+(dx+1)=3dy+dx+4 (adding 1 to have idx starting
% from 1)
% Then use a mapping array cm
%   --------------------------------------
%   | deltax | deltay | code | (base-3)+1 |
%   |-------------------------------------|
%   |    0   |   +1   |   2  |      8     |
%   |    0   |   -1   |   6  |      2     |
%   |   -1   |   +1   |   3  |      7     |
%   |   -1   |   -1   |   5  |      1     |
%   |   +1   |   +1   |   1  |      9     |
%   |   +1   |   -1   |   7  |      3     |
%   |   -1   |    0   |   4  |      4     |
%   |   +1   |    0   |   0  |      6     |
%   ---------------------------------------

idx=3*delta(:,1)+delta(:,2)+5;
cm([1 2 3 4 6 7 8 9])=[5 6 7 4 0 3 2 1];

% finally the chain code array and the starting point
cc.x0=b(1,2);
cc.y0=b(1,1);
cc.code=(cm(idx))';

% If unwrapping is required, use the following algorithm
%
% if a(k), k=1..n is the original code and u(k) the unwrapped:
%
% - u(1)=a(1)
% - u(k)=g(k),
%        g(k) in Z | (g(k)-a(k)) mod 8=0 and |g(k)-u(k-1)| is minimized
%
if (unwrap)
    a=cc.code;
    u(1)=a(1);
    la=size(a,1);
    for i=2:la
        n=round((u(i-1)-a(i))/8);
        u(i)=a(i)+n*8;
    end
    cc.ucode=u';
end

end


text:
function idx = getDealiasingIndexes(~, N, nVars)
%GETGDEALIASINGINDEXES  % Returns indexes for dealiasing procedure (2/3-rule).
%
% See also SPINOPSPHERE.

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

toOne = floor(N/2) + 1 - ceil(N/6):floor(N/2) + ceil(N/6);
idx = false(N, N);
idx(toOne, toOne) = 1;
idx = repmat(idx, nVars, 1);

end


text:
function e = maxdiff(a,b,rel)
% MAXDIFF(A,B) returns the maximum difference in any field or element.
% Matching infinities or NaNs do not count.
%
% MAXDIFF(A,B,REL) measures the per-element relative difference (A-B)/(REL + A)
%
% Examples:
%   maxdiff([1 2 3 nan inf -inf],[1 2 4 nan inf -inf]) % = 1

% Written by Tom Minka
% (c) Microsoft Corporation. All rights reserved.

if nargin < 3
  rel = [];
end

e = 0;
if ~isequal(class(a), class(b))
  fprintf('maxdiff: incompatible types\n');
  e = Inf;
  return
end
if isa(a,'struct')
  for f = fieldnames(a)'
    field = char(f);
    if ~isfield(b,field)
      fprintf('maxdiff: second argument lacks field %s\n', field);
      e = Inf;
      return
    end
    e = max(e,maxdiff(a.(field), b.(field), rel));
  end
  return
end
if ~isequal(size(a),size(b))
  fprintf('maxdiff: size mismatch\n');
  e = Inf;
  return
end
a = a(:);
b = b(:);
if iscell(a)
  for i = 1:numel(a)
    e = max(e,maxdiff(a{i},b{i}));
  end
  return
end
i = isnan(a);
if any(i ~= isnan(b))
  % mismatched NaNs
  e = Inf;
  return
elseif sum(i) > 0
  a = a(~i);
  b = b(~i);
end
i = ~isfinite(a);
if any(i ~= ~isfinite(b))
  % mismatched infs
  e = Inf;
  return
elseif ~isequal(a(i),b(i))
  e = Inf;
  return
else
  a = a(~i);
  b = b(~i);
end
if isempty(a)
  e = 0;
  return
end
e = abs(a(:) - b(:));
if ~isempty(rel)
  e = e ./ (rel + abs(a(:)));
end
e = max(e);

author: andrewssobral

repo: lrslibrary


text:
%% Mean-Shift Video Tracking
% by Sylvain Bernhardt
% July 2008
%% Description
% Measures the similarity between two
% density estimations q and p done with
% a kernel which profile is k.
% q is the estimation of a reference patch
% and p the estimation of a candidate one 'T2'
% which size is H,W.
% The outputs are the similarity value f
% and the weight mask w for the gradient ascent
% in the extended Mean-Shift algorithm.
%
% [f,w] = Simil_func(q,p,T2,k,H,W)

function [f,w] = Simil_func(q,p,T2,k,H,W)

w = zeros(H,W);
f = 0;
for i=1:H
    for j=1:W
        w(i,j) = sqrt(q(T2(i,j)+1)/p(T2(i,j)+1));
        f = f+w(i,j)*k(i,j);
    end
end
% Normalization of f
f = f/(H*W);


text:
%% clear the workspace and select data
clear; clc; close all;

%% choose multiple datasets or just one
neuron = Sources2D();
nams = {'./data_1p.tif'};          % you can put all file names into a cell array; when it's empty, manually select files
nams = neuron.select_multiple_files(nams);  %if nam is [], then select data interactively

%% parameters
% -------------------------    COMPUTATION    -------------------------  %
pars_envs = struct('memory_size_to_use', 8, ...   % GB, memory space you allow to use in MATLAB
    'memory_size_per_patch', 0.5, ...   % GB, space for loading data within one patch
    'patch_dims', [64, 64],...  %GB, patch size
    'batch_frames', 1000);           % number of frames per batch
  % -------------------------      SPATIAL      -------------------------  %
gSig = 3;           % pixel, gaussian width of a gaussian kernel for filtering the data. 0 means no filtering
gSiz = 13;          % pixel, neuron diameter
ssub = 1;           % spatial downsampling factor
with_dendrites = true;   % with dendrites or not
if with_dendrites
    % determine the search locations by dilating the current neuron shapes
    updateA_search_method = 'dilate';  %#ok<UNRCH>
    updateA_bSiz = 5;
    updateA_dist = neuron.options.dist;
else
    % determine the search locations by selecting a round area
    updateA_search_method = 'ellipse'; %#ok<UNRCH>
    updateA_dist = 5;
    updateA_bSiz = neuron.options.dist;
end
spatial_constraints = struct('connected', true, 'circular', false);  % you can include following constraints: 'circular'
spatial_algorithm = 'hals';

% -------------------------      TEMPORAL     -------------------------  %
Fs = 10;             % frame rate
tsub = 1;           % temporal downsampling factor
deconv_options = struct('type', 'ar1', ... % model of the calcium traces. {'ar1', 'ar2'}
    'method', 'foopsi', ... % method for running deconvolution {'foopsi', 'constrained', 'thresholded'}
    'smin', -5, ...         % minimum spike size. When the value is negative, the actual threshold is abs(smin)*noise level
    'optimize_pars', true, ...  % optimize AR coefficients
    'optimize_b', true, ...% optimize the baseline);
    'max_tau', 100);    % maximum decay time (unit: frame);

nk = 3;             % detrending the slow fluctuation. usually 1 is fine (no detrending)
% when changed, try some integers smaller than total_frame/(Fs*30)
detrend_method = 'spline';  % compute the local minimum as an estimation of trend.

% -------------------------     BACKGROUND    -------------------------  %
bg_model = 'ring';  % model of the background {'ring', 'svd'(default), 'nmf'}
nb = 1;             % number of background sources for each patch (only be used in SVD and NMF model)
bg_neuron_factor = 1.4;
ring_radius = round(bg_neuron_factor * gSiz);  % when the ring model used, it is the radius of the ring used in the background model.
%otherwise, it's just the width of the overlapping area
num_neighbors = 50; % number of neighbors for each neuron

% -------------------------      MERGING      -------------------------  %
show_merge = false;  % if true, manually verify the merging step
merge_thr = 0.65;     % thresholds for merging neurons; [spatial overlap ratio, temporal correlation of calcium traces, spike correlation]
method_dist = 'max';   % method for computing neuron distances {'mean', 'max'}
dmin = 5;       % minimum distances between two neurons. it is used together with merge_thr
dmin_only = 2;  % merge neurons if their distances are smaller than dmin_only.
merge_thr_spatial = [0.8, 0.4, -inf];  % merge components with highly correlated spatial shapes (corr=0.8) and small temporal correlations (corr=0.1)

% -------------------------  INITIALIZATION   -------------------------  %
K = [];             % maximum number of neurons per patch. when K=[], take as many as possible.
min_corr = 0.8;     % minimum local correlation for a seeding pixel
min_pnr = 8;       % minimum peak-to-noise ratio for a seeding pixel
min_pixel = gSig^2;      % minimum number of nonzero pixels for each neuron
bd = 0;             % number of rows/columns to be ignored in the boundary (mainly for motion corrected data)
frame_range = [];   % when [], uses all frames
save_initialization = false;    % save the initialization procedure as a video.
use_parallel = true;    % use parallel computation for parallel computing
show_init = true;   % show initialization results
choose_params = true; % manually choose parameters
center_psf = true;  % set the value as true when the background fluctuation is large (usually 1p data)
% set the value as false when the background fluctuation is small (2p)

% -------------------------  Residual   -------------------------  %
min_corr_res = 0.7;
min_pnr_res = 6;
seed_method_res = 'auto';  % method for initializing neurons from the residual
update_sn = true;

% ----------------------  WITH MANUAL INTERVENTION  --------------------  %
with_manual_intervention = false;

% -------------------------  FINAL RESULTS   -------------------------  %
save_demixed = true;    % save the demixed file or not
kt = 3;                 % frame intervals

% -------------------------    UPDATE ALL    -------------------------  %
neuron.updateParams('gSig', gSig, ...       % -------- spatial --------
    'gSiz', gSiz, ...
    'ring_radius', ring_radius, ...
    'ssub', ssub, ...
    'search_method', updateA_search_method, ...
    'bSiz', updateA_bSiz, ...
    'dist', updateA_bSiz, ...
    'spatial_constraints', spatial_constraints, ...
    'spatial_algorithm', spatial_algorithm, ...
    'tsub', tsub, ...                       % -------- temporal --------
    'deconv_options', deconv_options, ...
    'nk', nk, ...
    'detrend_method', detrend_method, ...
    'background_model', bg_model, ...       % -------- background --------
    'nb', nb, ...
    'ring_radius', ring_radius, ...
    'num_neighbors', num_neighbors, ...
    'merge_thr', merge_thr, ...             % -------- merging ---------
    'dmin', dmin, ...
    'method_dist', method_dist, ...
    'min_corr', min_corr, ...               % ----- initialization -----
    'min_pnr', min_pnr, ...
    'min_pixel', min_pixel, ...
    'bd', bd, ...
    'center_psf', center_psf);
neuron.Fs = Fs;

%% distribute data and be ready to run source extraction
neuron.getReady_batch(pars_envs);

%% initialize neurons in batch mode
neuron.initComponents_batch(K, save_initialization, use_parallel);

%% udpate spatial components for all batches
neuron.update_spatial_batch(use_parallel);

%% udpate temporal components for all bataches
neuron.update_temporal_batch(use_parallel);

%% update background
neuron.update_background_batch(use_parallel);

%% delete neurons

%% merge neurons

%% get the correlation image and PNR image for all neurons
neuron.correlation_pnr_batch();

%% concatenate temporal components
neuron.concatenate_temporal_batch();
neuron.viewNeurons([],neuron.C_raw);

%% save workspace
neuron.save_workspace_batch();

author: zhoupc


text:
% See also: hsvargplvm_init

function [options, optionsDyn] = hsvargplvmOptions(globalOpt, timeStampsTraining, labelsTrain)

if nargin < 2
    timeStampsTraining = [];
end

if nargin < 3
    labelsTrain = [];
end

%-- One options structure where there are some parts shared for all
% models/layers and some parts specific for a few layers / submodels.

options = vargplvmOptions('dtcvar');

% Taken from globalOpt
options.H = globalOpt.H;
options.baseKern = globalOpt.baseKern;
options.Q = globalOpt.Q;
options.K = globalOpt.K;
options.enableDgtN = globalOpt.DgtN;
options.initial_X = globalOpt.initial_X;
options.initX = globalOpt.initX;
options.multOutput = globalOpt.multOutput;
if options.multOutput > 2
    warning('Multoutput > 2 has wrong derivatives!!')
end
%
options.optimiser = 'scg2';

% !!!!! Be careful to use the same type of scaling and bias for all models!!!
% scale = std(Ytr);
% scale(find(scale==0)) = 1;
%options.scaleVal = mean(std(Ytr));
% options.scaleVal = sqrt(var(Ytr{i}(:))); %%% ??
options.scale2var1 = globalOpt.scale2var1;

options.fixInducing = globalOpt.fixInducing;

%----- Parent prior (we call priors "dynamics", but it can actually be some
% other type of prior, eg labels etc.).
% The relevant fields of globalOpt are coming from the svargplvm_init
% (called within hsvargplvm_init).
if isempty(globalOpt.dynamicsConstrainType) || nargout < 2
    optionsDyn = [];
else
    if ~isempty(labelsTrain)
        optionsDyn.labelsTrain = labelsTrain;
    end
    % This does not mean it needs time inputs, it's just saying that it'll
    % use the kernel types and reparametrization used in VGPDS (regressive
    % "dynamics" etc).
    optionsDyn.type = 'vargpTime';
    optionsDyn.inverseWidth=30;
    optionsDyn.vardistCovars = globalOpt.vardistCovarsMult;
    if iscell(globalOpt.initX)
        optionsDyn.initX = globalOpt.initX{end};
    else
        optionsDyn.initX = globalOpt.initX;
    end
    optionsDyn.constrainType = globalOpt.dynamicsConstrainType;
    if ~isempty(timeStampsTraining)
        optionsDyn.t = timeStampsTraining;
    end
    optionsDyn.kern = globalOpt.dynamicKern;
    optionsDyn.initCovarMedian = globalOpt.dynamicsInitCovarMedian;
    optionsDyn.initCovarMedianLowest = globalOpt.initCovarMedianLowest;
    optionsDyn.initCovarMedianHighest = globalOpt.initCovarMedianHighest;
end

author: SheffieldML

repo: deepGP

save_path: github-repos/MATLAB/SheffieldML-deepGP