function [sigma,C] = calcBSImpVol(cp,P,S,K,T,r,q)
%
% [sigma,C] = calcBSImpVol(cp,P,S,K,T,r,q)
%
% Calculates Black-Scholes Implied Volatility Surface for an Option Price Matrix.
% Uses Li's Rational Function Approximator for the Initial Estimate, followed by
% 3rd-Order Householder's Root Finder (i.e. using vega,vomma & ultima) for greater
% convergence rate and wider domain-of-convergence relative to Newton-Raphson. Both
% Li's Approximator and the Root Finder are calculated matrix-wise (i.e.
% fully vectorized) for increased efficiency.
%
%
%   Input Parameters
%       cp      Call[+1],Put[-1]  [m x n],[1 x 1]
%       P       Option Price Matrix [m x n]
%       S       Underlying Price [1 x 1]
%       K       Strike Price [m x n]
%       T       Time to Expiry [m x n]
%       r       Continuous Risk-Free Rate [m x n],[1 x 1]
%       q       Continuos Div Yield [m x n],[1 x 1]
%
%   Output Parameters
%       sigma   Implied Volatility  [m x n]
%       C       Convergence Flag  [m x n]
%
%
%
%   Example 1
%     S = 100; K = (40:25:160)'; T = (0.25:0.25:1)'; % Define Key Variables
%     r = 0.01; q = 0.03;
%     cp = 1; % i.e. call
%     P = ...
%         [[59.3526805861312,34.4154741312210,10.3406451776045,0.501199199160055,0.0101623685145268;];
%         [58.7107005379958,33.8563481863964,10.9917759513981,1.36915029885860,0.143324063090580;];
%         [58.0742593358310,33.3567195106962,11.5012247391034,2.12859686975881,0.400045353619436;];
%         [57.4444414750070,32.9126689586500,11.9027988146544,2.77274776123341,0.708059729236718;];];
%     [mK,mT] = meshgrid(K,T);
%     [sigma,C] = calcBSImpVol(cp,P,S,mK,mT,r,q);
%     mesh(mK,mT,sigma);
%
%     Example 2
%       S = 100; K = (40:25:160)'; T = (0.25:0.25:1)'; % Define Key Variables
%       cp = [ones(4,3),-ones(4,2)]; % [Calls[4,3],Puts[4,2]]
%       R = 0.01*repmat([1.15;1.10;1.05;1],1,5); %
%       Q = 0.03*repmat([1.3;1.2;1.1;1],1,5);
%       P = ...
%           [[59.1445725607811,34.2167401269277,10.1798771553458,16.1224863211251,40.5779719086946];
%           [58.4355054500906,33.5945826994415,10.7977275764632,17.4776751735401,41.1533978186314];
%           [57.8694061672804,33.1636044111551,11.3636963648521,18.6294544130139,41.7369813312724];
%           [57.4444414750070,32.9126689586500,11.9027988146694,19.5839252875422,42.2704830992694]];
%       [mK,mT] = meshgrid(K,T);
%       [sigma,C] = calcBSImpVol(cp,P,S,mK,mT,R,Q);
%       mesh(mK,mT,sigma);
%       hold on; scatter3(mK(:),mT(:),sigma(:),60,[0,0,0],'filled'); hold off
%       xlabel('Strike'); ylabel('Expiry'); zlabel('Volatility');
%
% References:
%   1)  Li, 2006, "You Don't Have to Bother Newton for Implied Volatility"
%       http://papers.ssrn.com/sol3/papers.cfm?abstract_id=952727
%   2)  http://en.wikipedia.org/wiki/Householder's_method
%   3)  http://en.wikipedia.org/wiki/Greeks_(finance)
%
%% APPLY LI's RATIONAL-FUNCTION APPROXIMATOR
[g,h] = size(P);
if isscalar(r); r = r*ones(g,h); end
if isscalar(q); q = q*ones(g,h); end
if isscalar(cp); cp = cp*ones(g,h); end
p = [-0.969271876255; 0.097428338274; 1.750081126685];
m = ones(1,1,14);
m(:) = [...
    6.268456292246;
    -6.284840445036;
    30.068281276567;
    -11.780036995036;
    -2.310966989723;
    -11.473184324152;
    -230.101682610568;
    86.127219899668;
    3.730181294225;
    -13.954993561151;
    261.950288864225;
    20.090690444187;
    -50.117067019539;
    13.723711519422];
m = m(ones(g,1),ones(1,h),:); % Repmat to size [g,h]
n = ones(1,1,14);
n(:) = [...
    -0.068098378725;
    0.440639436211;
    -0.263473754689;
    -5.792537721792;
    -5.267481008429;
    4.714393825758;
    3.529944137559;
    -23.636495876611;
    -9.020361771283;
    14.749084301452;
    -32.570660102526;
    76.398155779133;
    41.855161781749;
    -12.150611865704];
n = n(ones(g,1),ones(1,h),:); % Repmat to size [g,h]
i = ones(1,1,14);
i(:) = [0,1,0,1,2,0,1,2,3,0,1,2,3,4];
i = i(ones(g,1),ones(1,h),:); % Repmat to size [g,h]
j = ones(1,1,14);
j(:) = [1,0,2,1,0,3,2,1,0,4,3,2,1,0];
j = j(ones(g,1),ones(1,h),:); % Repmat to size [g,h]
x = log(S.*exp((r-q).*(T))./K); % Calculate Normalized Moneyness Measure
P(cp==-1) = P(cp==-1) + S.*exp(-q(cp==-1).*T(cp==-1)) - K(cp==-1).*exp(-r(cp==-1).*T(cp==-1)); % Convert Put to Call by Parity Relation
c = P./(S.*exp(-q.*(T))); % Normalized Call Price
x = x(:,:,ones(1,14)); % Repmat to 3d size 14
c = c(:,:,ones(1,14)); % Repmat to 3d size 14
% Rational Function -  Eqn(19) of Li 2006
fcnv = @(p,m,n,i,j,x,c)(p(1).*x(:,:,1) + p(2).*sqrt(c(:,:,1)) + p(3).*c(:,:,1) + (sum(n.*((x.^i).*(sqrt(c).^j)),3))./(1 + sum(m.*((x.^i).*(sqrt(c).^j)),3)));
v1 = fcnv(p,m,n,i,j,x,c); % D- Domain (x<=-1)
v2 = fcnv(p,m,n,i,j,-x,exp(x).*c + 1 -exp(x)); % Reflection for D+ Domain (x>1)
v = zeros(g,h); v(x(:,:,1)<=0)=v1(x(:,:,1)<=0); v(x(:,:,1)>0)=v2(x(:,:,1)>0);
% Domain-of-Approximation is x={-0.5,+0.5},v={0,1},x/v={-2,2}
domainFilter = x(:,:,1)>=-0.5 & x(:,:,1)<=0.5 & v > 0 & v <1 & (x(:,:,1)./v)<=2 & (x(:,:,1)./v)>=-2;
sigma = v./sqrt(T); % v = sigma.*(sqrt(T));
sigma(~domainFilter) = 0.8; % use 0.8 arbtrarily as best vol-guess for out-of-domain points
%% HOUSEHOLDER ROOT-FINDER FOR INCREASED CONVERGENCE
d1fcn = @(sig,C)((log(S./K(C)) + (r(C)-q(C)+sig(C).^2*0.5).*(T(C)))./(sig(C).*sqrt(T(C))));
d2fcn = @(sig,C)((log(S./K(C)) + (r(C)-q(C)-sig(C).^2*0.5).*(T(C)))./(sig(C).*sqrt(T(C))));
callfcn = @(sig,C)( +exp(-q(C).*T(C)).*S.*fcnN(d1fcn(sig,C)) - exp(-r(C).*T(C)).*K(C).*fcnN(d2fcn(sig,C)) );
vegafcn = @(sig,C)(S.*exp(-q(C).*(T(C))).*fcnn(d1fcn(sig,C)).*(sqrt(T(C))));
vommafcn = @(sig,C)(S.*exp(-q(C).*(T(C))).*fcnn(d1fcn(sig,C)).*(sqrt(T(C))).*d1fcn(sig,C).*d2fcn(sig,C)./sig(C));
ultimafcn = @(sig,C)(-S.*exp(-q(C).*(T(C))).*fcnn(d1fcn(sig,C)).*(sqrt(T(C))).*(d1fcn(sig,C).*d2fcn(sig,C).*(1-d1fcn(sig,C).*d2fcn(sig,C))+d1fcn(sig,C).^2+d2fcn(sig,C).^2)./(sig(C).^2));
tolMat=1e-12;
k_max = 10; % 10 Householder Iterations
objfcn = @(sig,C)(P(C) - callfcn(sig,C));
C = true(size(P(:))); err = objfcn(sigma,C); % calculate initial error
C = abs(err)>tolMat; % Convergence Matrix
k = 1; % Initialize Count
while any(C) && k<=k_max % Iterate until sooner of Convergence or Count-limit

    % Calculate Derivatives (Greeks)
    vega = vegafcn(sigma,C); %f'(x_n)
    vomma = vommafcn(sigma,C); %f''(x_n)
    ultima = ultimafcn(sigma,C); %f'''(x_n)
%     % Newton Raphson Method x_n+1 = x_n + f(x_n)/f'(x_n)
%     sigma = sigma  + (err(C)./vega) ;
%     % Halley Method x_n+1 = x_n - f(x_n)/( f'(x_n) - f(x_n)*f''(x_n)/2*f'(x_n))
%     sigma = sigma  - err(C)./(-vega-(-err(C).*vomma./(-2.*vega)));
    % Householder Method x_n+1 = x_n - f(x_n)/( f'(x_n) - f(x_n)*f''(x_n)/2*f'(x_n))
    sigma(C) = sigma(C) - (6.*err(C).*vega.^2 + 3.*err(C).^2.*vomma)./(-6.*vega.^3 - 6.*err(C).*vega.*vomma - err(C).^2.*ultima);

    % Update Error
    err(C) = objfcn(sigma,C); %

    % Ascertain Convergence to Tolerance
    C = abs(err)>tolMat; % Convergence Matrix

    % Increment Count
    k = k + 1;
end
sigma(C) = NaN; % any remaining sigma are not worth calculating
end
%% Gaussian Subfunctions
function p=fcnN(x)
p=0.5*(1.+erf(x./sqrt(2)));
end
%
function p=fcnn(x)
p=exp(-0.5*x.^2)./sqrt(2*pi);
end


text:
 function pgd_step_test
%function pgd_step_test
% test the PGD algorithm using 2D LS cost function

% cost function terms
kap = 4;
A = [1 0; 0 sqrt(kap)];
W = eye(2);
M = eye(2);
yy = 0;

f.niter = 10;
x = [-kap; 1];

% run PSD
xpsd = qpwls_psd(x, A, W, yy, 0, 'precon', 1, 'niter', f.niter, 'isave', 'all');

% run PGD
data = {yy, A, W};
f.step = 3.3; % best step is 0.4, so illustrate suboptimal step
[xpgd1 steps] = pgd_step(x, data, @costgrad, 'precon', M, ...
		'step_method', 'user', 'step', f.step, ...
		'niter', f.niter, 'isave', 'all', 'chat', 0);

%pr steps
[xpgd2 steps] = pgd_step(x, data, @costgrad, 'precon', M, ...
		'step_method', 'der2', ...
		'niter', f.niter, 'isave', 'all', 'chat', 0);
%pr steps

if im
	clf, subplot(211)
	plot(	xpsd(1,:), xpsd(2,:), 'y-o', ...
		xpgd1(1,:), xpgd1(2,:), 'b--+', ...
		xpgd2(1,:), xpgd2(2,:), 'm:.')
	xlabel x1, ylabel x2
	title 'QPWLS example'
	leg = {'PSD', sprintf('PGD:%g', f.step), 'PGD:Newton'};
	legend(leg{:}, 3)
end

% cost function map
if ~isvar('qq')
	x1 = linspace(-4.5,0.5,41)';
	x2 = linspace(-1,1.5,43)';
	[xx1 xx2] = ndgrid(x1,x2);
	qq = 0 * xx1;
	for i1=1:length(x1)
		for i2=1:length(x2)
			x = [x1(i1) x2(i2)]';
			qq(i1,i2) = norm(sqrtm(W) * (yy - A * x)).^2;
		end
	end
end

if im
	hold on
	contour(x1, x2, qq', 8)
	plot(0,0, 'rx')
	hold off
	axis equal
	axis([-4.5 0.5 -1 1.5])
	xtick(-4:0)
	ytick(-1:1)
	colormap(0.3+0.7*gray)
end

if im
	subplot(212)
	plot(	0:f.niter, sum(abs(xpsd), 1), 'y-o', ...
		0:f.niter, sum(abs(xpgd1)), 'b--x', ...
		0:f.niter, sum(abs(xpgd2)), 'm:.')
	xlabel 'iteration', ylabel '||x||_1'
	legend(leg{:}, 1)
end

% cost function and gradient for basic WLS
function [cost, grad] = costgrad(x, data)
y = data{1};
A = data{2};
W = data{3};

p = A * x;
cost = 0.5 * norm(sqrtm(W) * (y - p)).^2;
if nargout > 1
	grad = -A' * (W * (y - p));
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function out = find_closest_region(atl, input_xyz_coord)
% Find the closest atlas region to an input [x, y, z] mm coordinate
%
% First line: One-line summary description of function
%
% :Usage:
% ::
%
%     out = find_closest_region(atlas_obj, input_xyz_coord)
%
%
% :Inputs:
%
%   **atlas_obj:**
%        An atlas-class object
%
%   **input_xyz_coord:**
%        A row vector with x, y, z coordinates in mm
%
% :Outputs:
%
%   **out:**
%        A structure containing:
%
%
% Example:
% Find the Glasser 2016 parcel containing a coordinate and display it
% --------------------------------------------------------------------
% atl = load_atlas('glasser');
% input_xyz_coord = [-18 12 -22];
% out = find_closest_region(atl, input_xyz_coord);
% disp(out)
% orthviews(out.closest_region_obj);
% spm_orthviews('reposition', input_xyz_coord)
%
% Note: requires dist.m from Matlab neural net toolbox to work correctly.
% It is very fast.
% There are many replacements that could be added, but this is the current
% implementation.

% ..
%     Author and copyright information:
%
%     Copyright (C) 2021 Tor Wager
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
% ..

out = struct();

if ~isrow(input_xyz_coord)
    input_xyz_coord = input_xyz_coord';
end

if ~isrow(input_xyz_coord)
    error('Input a row vector with [x y z] coordinates in mm');
end

% transform all coordinates in object into mm
mmcoords = voxel2mm(atl.volInfo.xyzlist', atl.volInfo.mat);

% find closest
out.d = dist(input_xyz_coord , mmcoords(:,~atl.removed_voxels));

[out.mind, out.wh] = min(out.d);

out.region_number = atl.dat(out.wh);

if length(unique(out.region_number)) > 1
    warning('More than one region is equally close to input coordinate.');
end

out.closest_label = atl.labels{out.region_number};

out.closest_region_obj = select_atlas_subset(atl, out.region_number);

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%ESTIMATEAFFINE3D  Computes an optimal affine transformation between two 3D point sets
%
%     M = cv.estimateAffine3D(src, dst)
%     [M, inliers] = cv.estimateAffine3D(...)
%     [...] = cv.estimateAffine3D(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __src__ First input 3D point set containing `(X,Y,Z)`. Cell array of
%   3-element vectors `{[x,y,z],...}` or Nx3/Nx1x3/1xNx3 numeric array.
% * __dst__ Second input 3D point set containing `(x,y,z)`. Same size and type
%   as `src`.
%
% ## Output
% * __M__ Output 3D affine transformation matrix 3x4 of the form
%   `[a11 a12 a13 b1; a21 a22 a23 b2; a31 a32 a33 b3]`
% * __inliers__ Output vector of same length as number of points, indicating
%   which points are inliers (1-inlier, 0-outlier).
% * __result__ success flag.
%
% ## Options
% * __RansacThreshold__ Maximum reprojection error in the RANSAC algorithm to
%   consider a point as an inlier. default 3.0.
% * __Confidence__ Confidence level, between 0 and 1, for the estimated
%   transformation. Anything between 0.95 and 0.99 is usually good enough.
%   Values too close to 1 can slow down the estimation significantly. Values
%   lower than 0.8-0.9 can result in an incorrectly estimated transformation.
%   default 0.99.
%
% It computes:
%
%     [x; y; z] = [a11 a12 a13; a21 a22 a23; a31 a32 a33] * [X; Y; Z] + [b1; b2; b3]
%
% The function estimates an optimal 3D affine transformation between two 3D
% point sets using the RANSAC algorithm.
%
% See also: cv.estimateAffine2D
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function dose = matRad_calcPhotonDoseBixel(SAD,m,betas,Interp_kernel1,...
                  Interp_kernel2,Interp_kernel3,radDepths,geoDists,...
                  isoLatDistsX,isoLatDistsZ)
% matRad photon dose calculation for an individual bixel
%
% call
%   dose = matRad_calcPhotonDoseBixel(SAD,m,betas,Interp_kernel1,...
%                  Interp_kernel2,Interp_kernel3,radDepths,geoDists,...
%                  isoLatDistsX,isoLatDistsZ)
%
% input
%   SAD:                source to axis distance
%   m:                  absorption in water (part of the dose calc base
%                       data)
%   betas:              beta parameters for the parameterization of the
%                       three depth dose components
%   Interp_kernel1/2/3: kernels for dose calculation
%   radDepths:          radiological depths
%   geoDists:           geometrical distance from virtual photon source
%   isoLatDistsX:       lateral distance in X direction in BEV from central
%                       ray at iso center plane
%   isoLatDistsZ:       lateral distance in Z direction in BEV from central
%                       ray at iso center plane
%
% output
%   dose:               photon dose at specified locations as linear vector
%
% References
%   [1] http://www.ncbi.nlm.nih.gov/pubmed/8497215
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Define function_Di
func_Di = @(beta,x) beta/(beta-m) * (exp(-m*x) - exp(-beta*x));

% Calulate lateral distances using grid interpolation.
lat1 = Interp_kernel1(isoLatDistsX,isoLatDistsZ);
lat2 = Interp_kernel2(isoLatDistsX,isoLatDistsZ);
lat3 = Interp_kernel3(isoLatDistsX,isoLatDistsZ);

% now add everything together (eq 19 w/o inv sq corr -> see below)
dose = lat1 .* func_Di(betas(1),radDepths) + ...
       lat2 .* func_Di(betas(2),radDepths) + ...
       lat3 .* func_Di(betas(3),radDepths);

% inverse square correction
dose = dose .* (SAD./geoDists(:)).^2;

% check if we have valid dose values and adjust numerical instabilities
% from fft convolution
dose(dose < 0 & dose > -1e-14) = 0;
if any(isnan(dose)) || any(dose<0)
   error('Error in photon dose calculation.');
end

repo: matRad


text:
function C = multiTimes( A, B, type )
% Matrix multiply each submatrix of two 3D arrays without looping.
%
% type controls the matrix multiplication to perform (for each i):
%  1    - C(:,:,i) = A(:,:,i)*B(:,:)
%  1.1  - C(i,:,:) = A(:,:,i)*B(:,:)
%  1.2  - C(:,:,i) = A(:,:)*B(:,:,i)
%  2    - C(:,:,i) = A(:,:,i)*B(:,:,i)
%  2.1  - C(:,:,i) = A(:,:,i)'*B(:,:,i)
%  2.2  - C(:,:,i) = A(:,:,i)*B(:,:,i)'
%  3    - C(:,i) = A(:,:,i)*B(:,i)
%  3.1  - C(:,i) = A(:,:,i)'*B(:,i)
%  3.2  - C(:,i) = A(:,i)'*B(:,:,i)
%  4.1  - C(i) = trace(A(:,:,i)'*B(:,:,i))
% Corresponding dimensions of A and B must match appropriately.
%
% USAGE
%  C = multiTimes( A, B, type )
%
% INPUTS
%  A         - [ma x na x oa] matrix
%  B         - [mb x nb x ob] matrix
%  type      - multiplication type (see above)
%
% OUTPUTS
%  C         - result of the multiplication
%
% EXAMPLE
%  n=10000; A=randn(2,2,n); B=randn(2);
%  tic, C1=multiTimes(A,B,1); toc
%  tic, C2=zeros(size(A)); for i=1:n, C2(:,:,i)=A(:,:,i)*B; end; toc
%
% See also BSXFUN, MTIMES
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.52
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

ma = size(A,1); na = size(A,2); oa = size(A,3);
mb = size(B,1); nb = size(B,2); ob = size(B,3);

% just to simplify the reading
%if( ma==mb ); m = ma; end
%if( na==nb ); n = na; end
if( oa==ob ); o = oa; end

switch type
  case 1    % C(:,:,i) = A(:,:,i)*B(:,:)
    C = permute(reshape(reshape(permute(A,[1 3 2]),ma*oa,na)*B,...
      ma,oa,nb),[1 3 2]);
  case 1.1  % C(i,:,:) = A(:,:,i)*B(:,:)
    C = reshape(reshape(permute(A,[3 1 2]),ma*oa,na)*B,oa,ma,nb);
  case 1.2  % C(:,:,i) = A(:,:)*B(:,:,i)
    C = reshape(A*reshape(B,mb,nb*ob),ma,nb,ob);
  case 2    % C(:,:,i) = A(:,:,i)*B(:,:,i)
    C = reshape(sum(bsxfun(@times,reshape(A,...
      [ma na 1 o]),reshape(B,[1 mb nb o])),2),[ma nb o]);
  case 2.1  % C(:,:,i) = A(:,:,i)'*B(:,:,i)
    C = reshape(sum(bsxfun(@times,reshape(permute(A,[2,1,3]),...
      [na ma 1 o]),reshape(B,[1 mb nb o])),2),[na nb o]);
  case 2.2  % C(:,:,i) = A(:,:,i)*B(:,:,i)'
    C = reshape(sum(bsxfun(@times,reshape(A,...
      [ma na 1 o]),reshape(permute(B,[2,1,3]),[1 nb mb o])),2),[ma mb o]);
  case 3    % C(:,i) = A(:,:,i)*B(:,i)
    C = reshape(sum(bsxfun(@times, A, reshape(B,[1 mb nb ])),2),ma,nb);
  case 3.1  % C(:,i) = A(:,:,i)'*B(:,i)
    C = reshape(sum(bsxfun(@times, permute(A,[2,1,3]), ...
      reshape(B,[1 mb nb ])),2),na,nb);
  case 3.2  % C(:,i) = A(:,i)'*B(:,:,i)
    C = reshape(sum(bsxfun(@times, reshape(A,ma,1,na), B),1), nb,na);
  case 4.1  % C(i) = tr(A(:,:,i)'*B(:,:,i))
    C = reshape(sum(sum(A.*B,1),2),1,o);
  otherwise
    error('unknown type: %f',type);
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function ekg_filt(up)
%EKG_FILT extracts respiratory signals using various filtering techniques
% from the ECG signal as specified in PC's literature review.
%	            ekg_filt(option, up)
%
%	Inputs:
%		option          the option which has led to this function being used
%       up              universal parameters structure
%
%	Outputs:
%       ...
%

fprintf('\n--- Extracting Respiratory Signals from ECG using Filtering Techniques ');
log_int_respSig = 0;             % Has value 1 unless this is a final respiratory signal

for subj = up.paramSet.subj_list

    %% Cycle through each ecg signal
    for sig_no = 1 : length(up.paramSet.ekg_sigs)

        %% Cycle through each method
        for filt_no = 1 : length(up.al.options.ekg_filt)

            %% Skip if this processing has been done previously
            eval(['save_name = ''' up.paramSet.ekg_sigs{sig_no}, up.paths.filenames.filt '_' up.al.options.ekg_filt{filt_no} ''';']);
            iden_resp_sig_file_ending
            savepath = [up.paths.data_save_folder, num2str(subj), ending];
            exist_log = check_exists(savepath, save_name);
            if exist_log
                continue
            end

            %% Load relevant data
            if ~exist('data', 'var')
                load([up.paths.data_load_folder, up.paths.data_load_filename]);
            end
            % Extract EKG data
            eval(['rel_data = data(subj).' up.paramSet.ekg_sigs{sig_no} ';']);

            %% Filter the raw signal using this method
            respWave = feval(up.al.options.ekg_filt{filt_no}, rel_data, up);

            %% Band-pass filter
            filtered_data = bpf_signal_to_remove_non_resp_freqs(respWave, respWave.fs, up);
            eval([save_name ' = filtered_data;']);

            %% Save processed data
            save_or_append_data
        end

    end

end

end

author: peterhcharlton

repo: RRest


text:
function YESNO = isnan(F)
YESNO=0;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [candidates, score] = run_edge_boxes90(im, num_candidates)

  if ~isdeployed
    old_path = path;

    addpath('../edge_box');
    addpath('../edge_box/toolbox/matlab/');
    addpath('../edge_box/toolbox/channels/');
  end

  % Demo for Edge Boxes (please see readme.txt first).

  %% load pre-trained edge detection model (see edgesDemo.m)
  model=load('/path/to/edge_box/models/forest/modelBsds'); model=model.model;

  %% set up opts for edgeBoxes (see edgeBoxes.m)
  opts = edgeBoxes;
  opts.alpha = .85;     % step size of sliding window search
  opts.beta  = .95;     % nms threshold for object proposals
  opts.minScore = .01;  % min score of boxes to detect
  opts.maxBoxes = 1e4;  % max number of boxes to detect

  %% detect bbs (no visualization code for now)
  tic, bbs=edgeBoxes(im,model,opts); toc
  candidates = double(bbs(:,1:4));
  candidates(:,3:4) = candidates(:,3:4) + candidates(:,1:2);
  score = double(bbs(:,end));

  if ~isdeployed
    path(old_path);
  end
end

author: hosang

repo: detection-proposals


text:
function E = ECGModelError(X,ECGmn,Phasemn,flag)
% E = ECGModelError(X,ECGmn,Phasemn,flag)
% Overview: Generate Synthetic ECG beat with parameters given in the variable X. Generate the artificial ECG
% over phase given in Phasemn. Compute the ECGModelError as the difference between the synthetic
% ECG beat and ECG beat stored in ECGmn.
%
% Inputs:
%           X - contains the amplitude (alphai), standard deviation (bi) and phase (tetai) for
%           each Gaussian function used for generating the synthetic ECG beat.
%
%           ECGmn - Average ECG beat estimated from the ECG analysis
%           window.
%
%           Phasemn - array contains the phase points over which the synthetic
%           ECG beat is generated.
%
%           flag - set to 1 to return generated sythetic ECG beat. Set to 0 to
%           return the difference between generated synthetic ECG beat and average ECG beat for analysis window.
%
% Outputs:
%           E - synthetic ECG beat generated using Gaussian parameter estimates if flag set to 1.
%               difference between synthetic ECG beat and average ECG beat for analysis window if flag set to 0.
%
% Open Source ECG Toolbox, version 1.0, November 2006
% Released under the GNU General Public License
% Copyright (C) 2006  Reza Sameni
% Sharif University of Technology, Tehran, Iran -- LIS-INPG, Grenoble, France
% reza.sameni@gmail.com
% last modified:
% on 11/30/2020 by Ismail Sadiq
%
% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2 of the License, or (at your
% option) any later version.
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
% Public License for more details. You should have received a copy of the
% GNU General Public License along with this program; if not, write to the
% Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
% MA  02110-1301, USA.

L = (length(X)/3);

alphai = X(1:L);
bi = X(L+1:2*L);
tetai = X(2*L+1:3*L);

Z = zeros(size(ECGmn));
for j = 1:length(alphai),
    dtetai = rem(Phasemn - tetai(j) + pi,2*pi)-pi;
    Z = Z + alphai(j) .* exp(-dtetai .^2 ./ (2*bi(j) .^ 2));
end

if(flag==0)
    E = (Z-ECGmn);
elseif(flag==1)
    E = Z;
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
% pop_plotdata() - Plot average of EEG channels or independent components in
%                  a rectangular array. Else, (over)plot single trials.
% Usage:
%   >> avg = pop_plotdata(EEG, typeplot, indices, trials, title, singletrials, ydir, ylimits);
%
% Inputs:
%   EEG        - Input dataset
%   typeplot   - Type data to plot (1=channels, 0=components) {Default:1}
%   indices    - Array of channels (or component) indices to plot
%                     {Default: all}
%   trials     - Array of trial indices. sum specific trials in the average
%                     {Default: all}
%   title      - Plot title. {Default: []}.
%   singletrials - [0|1], Plot average or overplot single trials
%                      0 plot average, 1 plot single trials {Default: 0}
%   ydir       - [1|-1] y-axis polarity (pos-up = 1; neg-up = -1)
%                {def command line-> pos-up; def GUI-> neg-up}
%   ylimits    - [ymin ymax] plotting limits {default [0 0] -> data range}
%
% Outputs:
%   avg        - [matrix] Data average
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: plotdata(), eeglab()

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% 01-25-02 reformated help & license -ad
% 03-08-02 add eeglab options -ad
% 03-18-02 added title -ad & sm
% 03-30-02 added single trial capacities -ad

function [sigtmp, com] = pop_plotdata(EEG, typeplot, indices, trials, plottitle, singletrials, ydir,ylimits);

ylimits = [0 0];% default ylimits = data range

% warning signal must be in a 3D form before averaging

sigtmp = [];

com = '';
if nargin < 1
   help pop_plotdata;
   return;
end;
if nargin < 2
	typeplot = 1; % 1=signal; 0=component
end;
if exist('plottitle') ~= 1
	plottitle = '';
end;

if nargin <3
	if typeplot % plot signal channels
		result = inputdlg2({  'Channel number(s):'  ...
                                      'Plot title:'  ...
                                      'Vertical limits ([0 0]-> data range):' ...
                                   }, ...
                                  'Channel ERPs in rect. array -- pop_plotdata()', 1, ...
                                   {['1:' int2str(EEG.nbchan)] [fastif(isempty(EEG.setname), '',[EEG.setname ' ERP'])] ['0 0'] }, ...
                                   'pop_plotdata' );
	else % plot components
		result = inputdlg2({  'Component number(s):' ...
                                      'Plot title:' ...
                                      'Vertical limits ([0 0]-> data range):' ...
                                   }, ...
                                   'Component ERPs in rect. array -- pop_plotdata()', 1, ...
                                   {   ['1:' int2str(size(EEG.icawinv,2))] [fastif(isempty(EEG.setname), '',[EEG.setname ' ERP'])] ['0 0'] }, ...
                                   'pop_plotdata' );
	end;
	if length(result) == 0 return; end;

	indices   	 = eval( [ '[' result{1} ']' ] );

	plottitle    = result{2};

	singletrials = 0;

	ylimits   	 = eval( [ '[' result{3} ']' ] );
        if length(ylimits) ~= 2
            ylimits = [0 0]; % use default if 2 values not given
        end
end;
if ~(exist('trials') == 1)
	trials = 1:EEG.trials;
end;
if exist('plottitle') ~= 1
    plottitle = '';
end;
if exist('singletrials') ~= 1
    singletrials = 0;
end;

if EEG.trials > 1 & singletrials == 0
    fprintf('Selecting trials and components...\n');
	if typeplot == 1
	   sigtmp = nan_mean(EEG.data(indices,:,trials),3);
 	else
	   if isempty(EEG.icasphere)
	      error('no ICA data for this set, first run ICA');
	   end;
       tmpdata = eeg_getdatact(EEG, 'component', indices, 'trialindices', trials);
	   fprintf('Averaging...\n');
	   sigtmp = nan_mean(tmpdata,3);
	end;
else
	if typeplot == 1
	   sigtmp = EEG.data(indices,:,trials);
	else
	   if isempty(EEG.icasphere)
	      error('no ICA data for this set, first run ICA');
	   end;
       sigtmp = eeg_getdatact(EEG, 'component', indices, 'trialindices', trials);
	end;
end;
figure;
try, icadefs; set(gcf, 'color', BACKCOLOR); catch, end;
if exist('YDIR') ~= 1
     ydir = 1;
else
     ydir = YDIR;
end

%
%%%%%%%%%%%%%%%%%%%%%%%%%% make the plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
sigtmp = reshape( sigtmp, size(sigtmp,1),  size(sigtmp,2)*size(sigtmp,3));
ymin = ylimits(1);
ymax = ylimits(2);

if ischar(ymin)
   ymin = str2num(ymin);
   ymax = str2num(ymax);
end

% com = sprintf('pop_plotdata(%s, %d, %s, [1:%d], ''%s'', %d, %d, [%f %f]);', ...
% inputname(1), typeplot, vararg2str(indices), EEG.trials, plottitle, singletrials,ydir,ymin,ymax);
% fprintf([com '\n']);
if ~isempty(EEG.chanlocs) && typeplot == 1
    tmpchanlocs = EEG.chanlocs;
    chanlabels = strvcat({ tmpchanlocs(indices).labels });
else
    chanlabels = num2str(indices(:));
end;
plottopo( sigtmp, 'frames', EEG.pnts, 'limits', [EEG.xmin*1000 EEG.xmax*1000 ymin ymax], 'title', plottitle, 'chans', 1:size(sigtmp,1), 'ydir', ydir, 'channames', chanlabels);
%plotdata(sigtmp, EEG.pnts, [EEG.xmin*1000 EEG.xmax*1000 ymin ymax], plottitle, indices,0,0,ydir);

%
%%%%%%%%%%%%%%%%%%%%%%%%%% add figure title %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if typeplot == 1
	set(gcf, 'name', 'Plot > Channel ERPs > In rect. array -- plotdata()');
else
	set(gcf, 'name', 'Plot > Component ERPs > In rect. array -- plotdata()');
end;
%
%%%%%%%%%%%%%%%%%%%%%%%%%% set y-axis direction %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if exist('ydir') ~= 1
  if exist('YDIR') ~= 1
    ydir = 1;  % default positive up
  else
    ydir = YDIR; % icadefs.m system-wide default
  end
end

if ydir==1
   set(gca,'ydir','normal');
else % ydir == -1
   set(gca,'ydir','reverse');
end

switch nargin
	case {0, 1, 2, 3}, com = sprintf('pop_plotdata(%s, %d, %s, [1:%d], ''%s'', %d, %d, [%g %g]);', inputname(1), typeplot, vararg2str(indices), EEG.trials, plottitle, singletrials,ydir,ymin,ymax);
	case 4, com = sprintf('pop_plotdata(%s, %d, %s, %s, ''%s'', %d, %d, [%g %g]);', inputname(1), typeplot, vararg2str(indices), vararg2str(trials), plottitle, singletrials,ydir,ymin,ymax);
end;

fprintf([com '\n']);
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = nan_mean(in, dim)
    tmpin = in;
    tmpin(find(isnan(in(:)))) = 0;
    out = sum(tmpin, dim) ./ sum(~isnan(in),dim);

author: PatternRecognition

repo: OpenBMI


text:
clear; close all
        % Impunerea primei solutii initiale
sol_init = bvpinit(linspace(0,4,5),[1 0]);
        % Gasirea primei solutii
sol = bvp4c(@fbound,@ffront,sol_init);
        % Evaluarea primei solutii
x = linspace(0,4);
y1 = bvpval(sol,x);
        % Impunerea solutiei initiale a doua
sol_init = bvpinit(linspace(0,4,5),[-1 0]);
        % Gasirea solutiei a doua
sol = bvp4c(@fbound,@ffront,sol_init);
        % Evaluarea solutiei a doua
y2 = bvpval(sol,x);
        % Reprezentarea grafica a celor doua solutii
plot(x,y1(1,:),x,y2(1,:),'Linewidth',1.5);
xlabel('x');
ylabel('y');
grid
legend('Prima solutie','A doua solutie')

author: Sable

repo: mcbench-benchmarks


text:
% sim_spinecho_xN.m
% Robin Simpson and Jamie Near, 2014.
%
% USAGE:
% out = sim_spinecho(n,sw,Bfield,linewidth,sys,tau,Nechoes)
%
% DESCRIPTION:
% This function simulates a multi-echo spin-echo experiment with 'Nechoes' instantaneous RF pulses.
%
% INPUTS:
% n         = number of points in fid/spectrum
% sw        = desired spectral width in [Hz]
% Bfield    = main magnetic field strength in [T]
% linewidth = linewidth in [Hz]
% sys       = spin system definition structure
% tau       = echo time in [ms]
% Nechoes   = number of spin echoes (optional.  Default = 10);
%
% OUTPUTS:
% out       = simulated spectrum, in FID-A structure format, using multi-echo
%             spin-echo sequence.

function out = sim_spinecho_xN(n,sw,Bfield,linewidth,sys,tau,Nechoes)

if nargin<7
    Nechoes=10;
end

%Set water to centre
centreFreq=4.65;
for k=1:length(sys)
    sys(k).shifts=sys(k).shifts-centreFreq;
end

%Calculate Hamiltonian matrices and starting density matrix.
[H,d]=sim_Hamiltonian(sys,Bfield);

delay=tau/(2*Nechoes);

%BEGIN PULSE SEQUENCE************
d=sim_excite(d,H,'x');                            %EXCITE
for k=1:Nechoes
    d=sim_evolve(d,H,delay/1000);                        %Evolve by tau/2
    d=sim_rotate(d,H,180,'y');                       %180 degree refocusing pulse about y' axis.
    d=sim_evolve(d,H,delay/1000);                        %Evolve by tau/2
end
[out,dout]=sim_readout(d,H,n,sw,linewidth,90);  %Readout along y (90 degree phase);
%END PULSE SEQUENCE**************

%Correct the ppm scale:
out.ppm=out.ppm-(4.65-centreFreq);

%Fill in structure header fields:
out.seq=['spinecho_x' num2str(Nechoes)];
out.te=tau;
out.sim='ideal';

%Additional fields for compatibility with FID-A processing tools.
out.sz=size(out.specs);
out.date=date;
out.dims.t=1;
out.dims.coils=0;
out.dims.averages=0;
out.dims.subSpecs=0;
out.dims.extras=0;
out.averages=1;
out.rawAverages=1;
out.subspecs=1;
out.rawSubspecs=1;
out.flags.writtentostruct=1;
out.flags.gotparams=1;
out.flags.leftshifted=0;
out.flags.filtered=0;
out.flags.zeropadded=0;
out.flags.freqcorrected=0;
out.flags.phasecorrected=0;
out.flags.averaged=1;
out.flags.addedrcvrs=1;
out.flags.subtracted=1;
out.flags.writtentotext=0;
out.flags.downsampled=0;
out.flags.isFourSteps=0;

author: CIC-methods

repo: FID-A

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


text:
function [tests,pass,perf]=test_woody(varargin)

inputs={'verbose'};
verbose=0;
for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};']);
    end
end

%Test the examples
% Note: This requires the Signal Processing toolbox
str=['t=[0:1/1000:1];N=1001;x=sin(2*pi*t)+sin(4*pi*t)+sin(8*pi*t);' ...
    'y=exp(0.01*[-1*[500:-1:1] 0 -1*[1:500]]);'...
  's=x.*y;sig1=0;sig2=0.1;M=100;S=zeros(N,M);center=501;TAU=round((rand(1,M)-0.5)*160);'...
  'for i=1:M;tau=TAU(i);if(tau<0)' ...
  'S(:,i)=[s(-1*tau:end)''; zeros(-1*(tau+1),1)];'...
  'else; S(:,i)=[zeros(tau,1);s(1:N-tau)''; ];end;'...
  'if(i<50) S(:,i)=S(:,i) + randn(N,1).*sig1; else; S(:,i)=S(:,i) + randn(N,1).*sig2; end; end; '...
  '[wood]=woody(S,[],[],''woody'',''biased'');[thor]=woody(S,[],[],''thornton'',''biased'');'...
  'figure;subplot(211);' ...
  'plot(s,''b'',''LineWidth'',2);hold on;plot(S,''r'');plot(s,''b'',''LineWidth'',2);legend(''Signal'',''Measurements'');'...
  'subplot(212);plot(s);hold on;plot(mean(S,2),''r'');plot(wood,''g'');plot(thor,''k'');'...
  'legend(''Signal'',''Normal Ave'',''Woody Ave'',''Thornton Ave'');close all;'];

test_string={str};

clean_up={};
[tests,pass,perf]=test_wrapper(test_string,clean_up,verbose);

author: ikarosilva

repo: wfdb-app-toolbox


text:
function fem1d_bvp_linear_test03 ( )

%*****************************************************************************80
%
%% FEM1D_BVP_LINEAR_TEST03 carries out test case #3.
%
%  Location:
%
%    http://people.sc.fsu.edu/~jburkardt/m_src/fem1d_bvp_linear/fem1d_bvp_linear_test03.m
%
%  Discussion:
%
%    Use A3, C3, F3, EXACT3, EXACT_UX3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dianne O'Leary,
%    Scientific Computing with Case Studies,
%    SIAM, 2008,
%    ISBN13: 978-0-898716-66-5,
%    LC: QA401.O44.
%
  n = 11;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_BVP_LINEAR_TEST03\n' );
  fprintf ( 1, '  Solve -( A(x) U''(x) )'' + C(x) U(x) = F(x)\n' );
  fprintf ( 1, '  for 0 < x < 1, with U(0) = U(1) = 0.\n' );
  fprintf ( 1, '  A3(X)  = 1.0\n' );
  fprintf ( 1, '  C3(X)  = 2.0 * X\n' );
  fprintf ( 1, '  F3(X)  = - X * ( 2 * X * X - 3 * X - 3 ) * exp ( X )\n' );
  fprintf ( 1, '  U3(X)  = X * ( 1 - X ) * exp ( X )\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nodes = %d\n', n );
%
%  Geometry definitions.
%
  x_first = 0.0;
  x_last = 1.0;
  x = linspace ( x_first, x_last, n );
  x = x(:);

  u = fem1d_bvp_linear ( n, @a3, @c3, @f3, x );

  uexact = exact3 ( x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     I    X         U         Uexact    Error\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %4d  %8f  %8f  %8f  %8e\n', ...
      i, x(i), u(i), uexact(i), abs ( u(i) - uexact(i) ) );
  end

  e1 = l1_error ( n, x, u, @exact3 );
  e2 = l2_error_linear ( n, x, u, @exact3 );
  h1s = h1s_error_linear ( n, x, u, @exact_ux3 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  l1 norm of error  = %g\n', e1 );
  fprintf ( 1, '  L2 norm of error  = %g\n', e2 );
  fprintf ( 1, '  Seminorm of error = %g\n', h1s  );

  return
end
function value = a3 ( x )

%*****************************************************************************80
%
%% A3 evaluates A function #3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of A(X).
%
  value = 1.0;

  return
end
function value = c3 ( x )

%*****************************************************************************80
%
%% C3 evaluates C function #3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of C(X).
%
  value = 2.0 * x;

  return
end
function value = exact3 ( x )

%*****************************************************************************80
%
%% EXACT3 evaluates exact solution #3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of U(X).
%
  value = x .* ( 1.0 - x ) .* exp ( x );

  return
end
function value = exact_ux3 ( x )

%*****************************************************************************80
%
%% EXACT_UX3 evaluates the derivative of exact solution #3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of dUdX(X).
%
  value = ( 1.0 - x - x .* x ) .* exp ( x );

  return
end
function value = f3 ( x )

%*****************************************************************************80
%
%% F3 evaluates right hand side function #3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 March 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of F(X).
%
  value = - x .* ( 2.0 * x .* x - 3.0 * x - 3.0 ) .* exp ( x );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
maxmag = ceil(10*max(newt2.Magnitude))/10;
mima = min(newt2.Magnitude);
if mima > 0 ; mima = 0 ; end

[bval,xt2] = hist(newt2.Magnitude,(mima:0.1:maxmag));
% normalise to annula rates
bval = bval/(max(newt2.Date)-min(newt2.Date));
bvalsum = cumsum(bval); % N for M <=
bval2 = bval(length(bval):-1:1);
bvalsum3 = cumsum(bval(length(bval):-1:1));    % N for M >= (counted backwards)
xt3 = (maxmag:-0.1:mima);

backg_ab = log10(bvalsum3);

figure_w_normalized_uicontrolunits(bfig);delete(gca);delete(gca); delete(gca); delete(gca)
rect = [0.22,  0.3, 0.65, 0.6];           % plot Freq-Mag curves
axes('position',rect);

%%
% plot the cum. sum in each bin  %%
%%

pl =semilogy(xt3,bvalsum3,'sb');
set(pl,'LineWidth',1.0,'MarkerSize',6,...
    'MarkerFaceColor','w','MarkerEdgeColor','k');
hold on
%pl1 =semilogy(xt3,bval2,'^b');
%set(pl1,'LineWidth',1.0,'MarkerSize',4,...
%    'MarkerFaceColor',[0.7 0.7 .7],'MarkerEdgeColor','k');

bv2 = [];bv3 = [] ; me = [];BV = [];
ni2 = 300;
BB = [];
for i = 1:ni2/1:length(newt2)-ni2
    [bv magco stan ] =  bvalca2(newt2(i:i+ni2,:));
    nn2 = newt2(i:i+ni2,:);
    l = nn2(:,6) >= magco+0.0;
    nn2 = nn2(l,:);

    [bval,xt2] = hist(nn2(:,6),(mima:0.1:maxmag));
    % normalise to annual rates
    bval = bval/(max(nn2(:,3))-min(nn2(:,3)));
    k = mima:0.1:magco;
    bval(1:length(k)) = nan;
    BB = [BB ; bval];
    bvalsum = cumsum(bval); % N for M <=
    bval2 = bval(length(bval):-1:1);
    bvalsum3 = cumsum(bval(length(bval):-1:1));    % N for M >= (counted backwards)

    hold on
    pl =semilogy(xt3,bvalsum3,'sb');
    set(pl,'LineWidth',1.0,'MarkerSize',4,...
        'MarkerFaceColor',[rand(1,1) rand(1,1) rand(1,1)],'MarkerEdgeColor',[rand(1,1) rand(1,1) rand(1,1)]);
    hold on

end

allsum = mean(BB, 'omitnan');
bvalsum = cumsum(allsum); % N for M <=
bvalsum3 = cumsum(allsum(length(allsum):-1:1));    % N for M >= (counted backwards)
pl =semilogy(xt3,bvalsum3,'hb');
set(pl,'LineWidth',1.0,'MarkerSize',10,...
    'MarkerFaceColor','y','MarkerEdgeColor','k');
hold on

author: CelsoReyes


text:
function [H,mu,SR]=standardEpanechnikovKernelBW(param1,N)
%%STANDARDEPANECHNIKOVKERNELBW Estimate the (univariate or multivariate)
%           kernel bandwidth when approximating a set of samples whose
%           underlying density is (or is approximated as) Gaussian with an
%           Epanechnikov kernel. The Epanechnikov kernel has certain
%           desirable optimality properties.
%
%INPUTS: param1 If inputType=0, then this is:
%              xi A numDimXN set of samples of the distribution from which
%                 a Gaussian kernel estimator should interpolate a PDF. For
%                 multidimensional estimation, there must be >=numDim
%                 samples. Specifically, the sample covariance matrix must
%                 be positive definite.
%              If inputType=1, then this is:
%              SR A root covariance matrix such that SR*SR' equals the
%                 (exact or approximated) covariance matrix of the PDF
%                 being approximated.
%            N If this parameter is provided, then param1 is assumed to be
%              SR and N is the number of samples that are being smoothed.
%              Otherwise, if this is omitted or an empty matrix is passed,
%              param1 is xi.
%
%OUTPUTS: H The numDimXnumDim kernel bandwidth.
%         mu If xi is given, then this is the sample mean. Otherwise, this
%            is an empty matrix.
%         SR The root covariance matrix from the input or computed from xi.
%
%The formulae for the optimal bandwidth of an Epanechnikov kernel given
%samples from an underlying Gaussian distribution is given in Equation 2.3
%of [1]. The function EpanechnikovKernel implements the Epanechnikov
%kernel.
%
%REFERENCES:
%[1] C. Musso, N. Oudjane, and F. LeGland, "Improving regularised particle
%    filters," in Sequential Monte Carlo Methods in Practice, A. Doucet, J.
%    F. G. de Freitas, and N. J. Gordon, Eds., Ch. 12, New York:
%    Springer-Verlag, 2001.
%
%March 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2)
    N=[];
end

nx=size(param1,1);%Dimensionality.

if(isempty(N))%xi is given
    xi=param1;
    N=size(xi,2);
    [mu,S]=calcMixtureMoments(xi);
    SR=cholSemiDef(S,'lower',1);%Not fully triangular. but S=SR*SR'
else%SR is given and N must be given.
    mu=[];
    SR=param1;
end

cn=hypersphereVolume(1,nx);
A=((8/cn)*(nx+4)*2*sqrt(pi)^nx)^(1/(nx+4));
hOpt=A*N^(-1/(nx+4));

H=hOpt*SR*eye(nx);

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function flattenFromROI(view, ROI)
% Initiates the mrFlatMesh GUI, filling in some fields
%
%  flattenFromROI(view, ROI)
%
% Wrapper for mrFlatMesh GUI. Opens the GUI, filling in some
% default starting values. Uses the center of mass of the specified
% ROI as the start point.
%
% GLOBALS: uses mrSESSION fields to initialize the mrFlatMesh GUI.
%
% Example:
%   flattenFromROI(VOLUME{1});    % Uses selected (current) ROI
%
% HISTORY:
% 2002.03.04 RFD (bob@white.stanford.edu): wrote it
% 7/16/02 djh, replaced mrSESSION.vAnatomyPath with global vANATOMYPATH
% 2002.07.19 RFD: cleaned up default filenames a bit
% 2006.02.05 MMS: made it work with ROIs with only one coordinate

mrGlobals

switch view.viewType
    case 'Gray'
        leftGrayPath  = view.leftPath;
        rightGrayPath = view.rightPath;
    otherwise
        error([mfilename,': unsupported view type "',view.viewType,'".']);
end

if(~exist('ROI','var'))
    if(view.selectedROI>0)
        ROI = view.ROIs(view.selectedROI);
    else
        myErrorDlg('No ROI selected');
        return;
    end
end

mmPerPixXYZ = readVolAnatHeader(vANATOMYPATH);

startXYZ = round(mean(ROI.coords',1)');
startXYZ = [startXYZ(2);startXYZ(1);startXYZ(3)];

% Now how do we find which hemisphere this coord belongs to?
hemi = viewGet(view, 'roihemi', ROI.coords);
if isempty('hemi'), hemi = 'left'; end;

if strcmpi('hemi', 'left')
    grayPath = fullfile(leftGrayPath,'');
else
    grayPath = fullfile(rightGrayPath,'');
end

[grayPathDir,grayPathFileName] = fileparts(grayPath);
if(exist(grayPath)~=2)
    if(exist(grayPathDir,'dir'))
        % at least get them close
        grayPath = grayPathDir;
    else
        grayPath = '';
        grayPathDir = '';
    end
end

% If the Gray coords are saved in the project as 'coords.mat', then use
% this to get the gray graph
coordsFile = fullfile(HOMEDIR, 'Gray', 'coords.mat');
if exist(coordsFile, 'file'),
    grayPath = coordsFile;
end

% Try to guess the mesh file path
[p f e] = fileparts(leftGrayPath);

% if the gray path is a nifti class file, then we use this and not a mesh
if strcmpi(e, '.gz')
    if strcmpi('hemi', 'left')
        meshPath = fullfile(leftGrayPath,'');
    else
        meshPath = fullfile(rightGrayPath,'');
    end
else
    % otherwise we look for a mesh
    subDirs = {'','3dMeshes','3DMeshes','3Dmeshes'};
    extensions = {'.MrM','.mrm','.Mrm','.MRm','.MRM', 'nii.gz'};
    meshPath = grayPathDir;
    for ii=1:length(subDirs)
        for jj=1:length(extensions)
            if(exist(fullfile(grayPathDir,subDirs{ii},[grayPathFileName,extensions{jj}]))==2)
                meshPath = fullfile(grayPathDir,subDirs{ii},[grayPathFileName,extensions{jj}]);
                break;
            end
        end
    end
end
ts = now;
saveFileName = ['unfold_',datestr(ts,'mm'),datestr(ts, 'dd'),datestr(ts,'yy')];
if(exist(fullfile(grayPathDir,'Unfolds'), 'dir'))
    savePath = fullfile(grayPathDir, 'Unfolds', saveFileName);
else
    savePath = fullfile(grayPathDir, saveFileName);
end

unfoldRadiusMM = 60;

disp('calling mrFlatMesh...');

mrFlatMesh({'grayPath', grayPath, ...
    'meshPath', meshPath, ...
    'savePath', savePath, ...
    'scaleXYZ', mmPerPixXYZ , ...
    'startXYZ', startXYZ, ...
    'unfoldRadiusMM', unfoldRadiusMM, ...
    'hemi', hemi});

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function w = uminus(w)
   %-  Unary minus for waveforms.
   %   -A negates the elements of A.
   %
   %   B = UMINUS(W) is called for the syntax '-W' when W is an waveform.

   % AUTHOR: Celso Reyes, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   for n=1:numel(w)
      w(n).data = -w(n).data;
   end
end

author: geoscience-community-codes

repo: GISMO


text:
function out = cfg_example_run_add2(job)
% Example function that returns the sum of two numbers given in job.a in
% out. The output is referenced as out(1), this is defined in
% cfg_example_vout_add2.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_example_run_add2.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

out = sum(job.a);

author: spm


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function y = sabrtmax(asabr,b,r,nsabr,f,fbar,eps)
% calclates the maximum time for the sabr approximation formula
% to be valid (see Risk paper by Doust)

    a = asabr / eps;
    n = nsabr / eps;

    zF = f^(1-b)/asabr/(1-b);
    zFbar = fbar.^(1-b)/asabr/(1-b);

    fav = sqrt(f*fbar);
    gamma1 = b./fav;
    gamma2 = b*(b-1)./fav.^2;

    kH = 0.125 * (2*gamma2-gamma1.^2)*a^2.*fav.^(2*b) ...
        + 0.75 * r*n*a*gamma1.*fav.^b ...
        + 0.125*(2-3*r^2)*n^2;

    z = zF - zFbar;
    integral = (f^(1-b)-fbar.^(1-b)) / (1-b);

    xz = 1/nsabr *log((sqrt(1-2*nsabr*r*z+nsabr^2*z.^2)-r+nsabr*z)/(1-r));
    if f==fbar
       y = 12 ./(eps^2*8 * kH);
    else
       y = 12 ./(eps^2*(8 * kH + a^2*(log(f./fbar)./integral .* (z./xz)).^2));
    end
end


text:
function [num, den, gamma] = designIIR(phi, D, m0, h, M)

% designIIR: design filters F_i(z) (as in Section III.B.)
%
% Usage:    [num, den, gamma] = designIIR(phi, D, m0, h, M)
%
% INPUTS:
%   phi: input rational transfer functions
%   D:   vector of fractional delays
%   m0:  system delay tolerance
%   h:   fast sampling interval
%   M:   superresolution factor (integer)
%
% OUTPUTS:
%   gamma:  the H infinity norm of the induced error system K
%   num: numerator vectors
%   den: denominator vectors
%
% filter F_i(z) will be an IIR filter with num{i} is the coefficients of
% the numerator and den{i} is the coefficients of the denominator.
%
% See also: getF, demo

% Get the ingeter and residual of the delays
m = floor(D(:)/h);
d = D(:) - m * h;

% Get the digital system as in Prop. 1 and 2
Ad = getAd(phi, h, d);
Bd = getBd(phi, h, d);
Cd = getCd(phi, h, d);
Dd = zeros(size(Cd,1), size(Bd,2));

% The Integer Delay Operator
sysd = IntDelayOp([m0; m]);

% Get the digital system as in Prop. 3 by taking into account the integer
% delay operators
sys = sminreal( sysd * ss(Ad, Bd, Cd, Dd, -1) );    % sminreal to reduce the system's dimension on the fly
% sys = sysd * ss(Ad, Bd, Cd, Dd, -1);

% Get the system P (see Fig. 6)
[Ap, Bp, Cp, Dp, p, q] = getP(sys, M);

% Convert P into analog, since hinfsyn works in analog domain
[Ac, Bc, Cc, Dc] = Digital2Analog(Ap, Bp, Cp, Dp);

% Design in continuous time
P = pck(Ac, Bc, Cc, Dc);

ubd = 1;
lbd = 0;
inc = .0001;

[F, g, gamma] = hinfsyn(P,p,q,lbd,ubd,inc);

% Get the synthesis system F
[Af, Bf, Cf, Df] = unpck(F);

% Convert F back to digital domain
[Af, Bf, Cf, Df] = Analog2Digital(Af, Bf, Cf, Df);

% Convert to filter coefficients
[num, den] = getF(Af, Bf, Cf, Df, M);


text:
% Example 7.2: Maximum entropy distribution
% Section 7.2, Figures 7.2-7.3
% Boyd & Vandenberghe, "Convex Optimization"
% Originally by Lieven Vandenberghe
% Adapted for CVX by Michael Grant 4/11/06
%
% We consider a probability distribution on 100 equidistant points in the
% interval [-1,1]. We impose the following prior assumptions:
%
%    -0.1 <= E(X) <= +0.1
%    +0.5 <= E(X^2) <= +0.6
%    -0.3 <= E(3*X^3-2*X) <= -0.2
%    +0.3 <= Pr(X<0) <= 0.4
%
% Along with the constraints sum(p) == 1, p >= 0, these constraints
% describe a polyhedron of probability distrubtions. In the first figure,
% the distribution that maximizes entropy is computed. In the second
% figure, we compute upper and lower bounds for Prob(X<=a_i) for each
% point -1 <= a_i <= +1 in the distribution, as well as the maximum
% entropy CDF.

%
% Represent the polyhedron as follows:
%     A * p <= b
%     sum( p ) == 1
%     p >= 0
%

n  = 100;
a  = linspace(-1,1,n);
a2 = a .^ 2;
a3 = 3 * ( a.^ 3 ) - 2 * a;
ap = +( a < 0 );
A  = [ a   ; -a  ; a2 ; -a2  ; a3 ; -a3 ; ap ; -ap ];
b  = [ 0.1 ; 0.1 ;0.5 ; -0.5 ; -0.2 ; 0.3 ; 0.4 ; -0.3 ];

%
% Compute the maximum entropy distribution
%

cvx_expert true
cvx_begin
    variables pent(n)
    maximize( sum(entr(pent)) )
    sum(pent) == 1;
    A * pent <= b;
cvx_end

%
% Compute the bounds on Prob(X<=a_i), i=1,...,n
%

Ubnds = zeros(1,n);
Lbnds = zeros(1,n);
for t = 1 : n,
    cvx_begin quiet
        variable p( n )
        minimize sum( p(1:t) )
        p >= 0;
        sum( p ) == 1;
        A * p <= b;
    cvx_end
    Lbnds(t) = cvx_optval;
    cvx_begin quiet
        variable p( n )
        maximize sum( p(1:t) )
        p >= 0;
        sum( p ) == 1;
        A * p <= b;
    cvx_end
    Ubnds(t) = cvx_optval;
    disp( sprintf( '%g <= Prob(x<=%g) <= %g', Lbnds(t), a(t), Ubnds(t) ) );
end

%
% Generate the figures
%

figure( 1 )
stairs( a, pent );
xlabel( 'x' );
ylabel( 'PDF( x )' );

figure( 2 )
stairs( a, cumsum( pent ) );
grid on;
hold on
d = stairs(a, Lbnds,'r-');  set(d,'Color',[0 0.5 0]);
d = stairs(a, Ubnds,'r-');  set(d,'Color',[0 0.5 0]);
d = plot([-1,-1], [Lbnds(1), Ubnds(1)],'r-');
set(d,'Color',[0 0.5 0]);
axis([-1.1 1.1 -0.1 1.1]);
xlabel( 'x' );
ylabel( 'CDF( x )' );
hold off

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
% Calculate statistics of a rate map that depend on probability distribution function (PDF)
%
% Calculates information, sparsity and selectivity of a rate map. Calculations are done
% according to 1993 Skaggs et al. "An Information-Theoretic Approach to Deciphering the Hippocampal Code"
% paper. Another source of information is 1996 Skaggs et al. paper called
% "Theta phase precession in hippocampal neuronal populations and the compression of temporal sequences".
%
%  USAGE
%   [information, sparsity, selectivity] = analyses.mapStatsPDF(map)
%   map             Structure with rate map, output of analyses.map
%   information     Structure with fields:
%       rate        information rate [bits/sec]
%       content     Spatial information content [bits/spike]
%   sparsity        Sparsity value
%   selectivity     Selectivity value
%
%  SEE
%   See also analyses.map
%
function [information, sparsity, selectivity] = mapStatsPDF(map)
    if ~isstruct(map)
        error('BNT:arg', 'Incorrect argument. Map should be a structure. You are probably relying on old code when mapStatsPDF accepted map as matrix.');
    end
    T = nansum(nansum(map.time)); % overall trial duration
    posPDF = map.time / (T+eps); % probability of animal being in bin x:
                                 % duration of time animal spent in the bin divided by overall trial duration.

    meanrate = nansum(nansum(map.z .* posPDF));
    meansquarerate = nansum(nansum( (map.z .^ 2) .* posPDF ));
    if meansquarerate == 0
       sparsity = NaN;
    else
        sparsity = meanrate^2 / meansquarerate;
    end

    maxrate = nanmax(nanmax(map.z));
    if meanrate == 0;
       selectivity = NaN;
       information.content = nan;
       information.rate = nan;
    else
       selectivity = maxrate / meanrate;
       logArg = map.z / meanrate;
       logArg(logArg < 1) = 1;

       information.rate = nansum(nansum(posPDF .* map.z .* log2(logArg)));
       information.content = information.rate / meanrate;
    end
end

author: kavli-ntnu


text:
function cnn_imagenet_minimal()
% CNN_IMAGENET_MINIMAL   Minimalistic demonstration of how to run an ImageNet CNN model

% setup toolbox
run(fullfile(fileparts(mfilename('fullpath')), ...
  '..', 'matlab', 'vl_setupnn.m')) ;

% download a pre-trained CNN from the web
if ~exist('imagenet-vgg-f.mat', 'file')
  urlwrite('http://www.vlfeat.org/matconvnet/models/imagenet-vgg-f.mat', ...
    'imagenet-vgg-f.mat') ;
end
net = load('imagenet-vgg-f.mat') ;

% obtain and preprocess an image
im = imread('peppers.png') ;
im_ = single(im) ; % note: 255 range
im_ = imresize(im_, net.normalization.imageSize(1:2)) ;
im_ = im_ - net.normalization.averageImage ;

% run the CNN
res = vl_simplenn(net, im_) ;

% show the classification result
scores = squeeze(gather(res(end).x)) ;
[bestScore, best] = max(scores) ;
figure(1) ; clf ; imagesc(im) ;
title(sprintf('%s (%d), score %.3f',...
   net.classes.description{best}, best, bestScore)) ;

author: huangzehao

repo: caffe-vdsr


text:
function [blockSizes, blockStarts, blockSymbols] = getContigBlocks( ids )
%  Given a sequence "ids" of discrete positive integers
%     find all contiguous blocks within it
%USAGE:
%  [S, L, A] = getContigBlocks( [1 1 1 2 2 2 2 1 3 3] );
%  will return [B,L,S], where
%      S : [3 4 1 2] is vector of block sizes
%      L : [1 4 8 9] is vector of starting locations
%      A : [1 2 1 3] is vector of block labels (original vals in "ids")

blockSizes = [];
blockStarts = [];
blockSymbols = [];
aa = 1;
while aa <= length( ids )

    blockStarts(end+1) = aa;
    blockSymbols(end+1) = ids(aa);

    if aa == length( ids )
        % Reached the end of the input vector
        blockSizes(end+1) = 1;
        break;
    else
        % Search forwards
        didFindMismatch = 0;
        for bb = aa+1:length(ids)
            if ids(bb) ~= ids(aa);
                didFindMismatch = 1;
                bb = bb - 1;
                break;
            end
        end
        if didFindMismatch
            blockSizes(end+1) = bb - aa + 1;
        else
            blockSizes(end+1) = length(ids) - aa + 1;
            break;
        end

    end

    aa = bb+1;

end

author: michaelchughes

repo: NPBayesHMM


text:
function g = p03_g ( n, x )

%*****************************************************************************80
%
%% P03_G evaluates the gradient for problem 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real G(N), the gradient of the objective function.
%
  y = p03_yvec ( );

  g = zeros ( n, 1 );

  for i = 1 : 15

    d1 = 0.5 * ( i - 1 );
    d2 = 3.5 - d1 - x(3);
    arg = - 0.5 * x(2) * d2 * d2;
    t = x(1) * exp ( arg ) - y(i);

    g(1) = g(1) + 2.0 * exp ( arg ) * t;
    g(2) = g(2) - x(1) * exp ( arg ) * t * d2 * d2;
    g(3) = g(3) + 2.0 * x(1) * x(2) * exp ( arg ) * t * d2;

  end

  return
end


text:
function engine = belprop_inf_engine(fg, max_iter, momentum, tol, maximize)

if nargin < 2, max_iter = length(fg.G); end
if nargin < 3, momentum = 0; end
if nargin < 4, tol = 1e-3; end
if nargin < 5, maximize = 0; end

engine.fgraph = fg;
engine.max_iter = max_iter;
engine.momentum = momentum;
engine.tol = tol;
engine.maximize = maximize;

% store results computed by enter_evidence here
ndoms = length(fg.doms);
nvars = length(fg.vars);
engine.marginal_domains = cell(1, ndoms);

% to compute the marginal on each variable, we need to know which domain to marginalize
% so we represent each domain as a bit vector, and compute its (pre-evidence) weight
engine.dom_weight = [];

% engine.dom_bitv = sparse(ndoms, nvars);
% ns = fg.node_sizes;
% for i=1:ndoms
%   engine.dom_bitv(i, fg.doms{i}) = 1;
%   engine.dom_weight(i) = prod(ns(fg.doms{i}));
% end

engine = class(engine, 'belprop_inf_engine');

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
classdef Triplets < handle
    %TRIPLET Summary of this class goes here
    %   Detailed explanation goes here

    properties (Access = public)
    nmax,m,n,nnz;
    nzval,ri,ci;
    end

    methods (Access = public)
        function obj=Triplets()
            obj.nmax = 0;
            obj.m = 0;
            obj.n = 0;
            obj.nnz = 0;
            obj.ri = 0;
            obj.ci = 0;
        end
        function clear(obj)
            obj.nmax = 0;
            obj.m = 0;
            obj.n = 0;
            obj.nnz = 0;
            obj.ri = 0;
            obj.ci = 0;
        end
        function reserve(obj,rows,cols,num)
           clear(obj);
           obj.m = rows;
           obj.n = cols;
           obj.nmax = num;
           obj.nzval = zeros(num);
           obj.ri = zeros(num);
           obj.ci = zeros(num);
        end

        function add_row(obj,r,c,val)
           obj.ri(obj.nnz+1) = r;
           obj.ci(obj.nnz+1) = c;
           obj.nzval(obj.nnz+1) = val;
           obj.nnz = obj.nnz + 1;
           obj.nnz

        end
    end
end

author: SuTanTank

repo: VideoStitchingViaShakinessRemoving


text:
classdef ciatah < dynamicprops
	% Performs analysis on calcium imaging data. Also performs analysis on behavior (response signals) compared to stimulus or other continuous signals during an imaging session.
	% Biafra Ahanonu
	% started: 2014.07.31 (2020.10.25 [22:30:29])
	% This is a re-write of my old code from controllerAnalysis. Encapsulating the functions and variables as methods and properties in a class should allow easier maintenance/flexibility.
	% inputs
		%
	% outputs
		%

	% changelog
		% updated: 2017.01.15 [01:31:54]
		% 2019.05.08 [10:59:59] - Added check for required toolboxes used in class.
		% 2019.07.25 [09:39:16] - Updated loading so that users only need to be in the root calciumImagingAnalysis path but required folders do not need to be loaded.
		% 2019.08.20 [09:33:32] - Improved loading of folders and Miji to save time.
		% 2019.08.30 [12:58:37] - Added Java heap space memory check on initialization.
		% 2019.09/10 - Added GUI font option and signalExtractionTraceOutputType.
		% 2019.10.15 [21:57:45] - Improved checking for directories that should not be loaded, remove need for verLessThan('matlab','9.0') check.
		% 2019.12.22 [09:01:38] - Re-vamped selection list dialog boxes for methods, cell-extraction, folders, etc. Now loads as a figure and for some (like methods list) has tooltips that describe each list item.
		% 2020.03.23 [22:36:36] - Re-vamped selection dialog so everything is on a single panel and users can interactively alter how they filter folders.
		% 2020.05.07 [15:13:53] - Refactored a bit to remove methods that do not need to be in the main class definition to allow easier updating later.
		% 2021.03.25 [20:15:31] - Use ciapkg.getDir() to ensure using root directory consistently rather than ciapkgRoot.
		% 2021.03.25 [22:46:38] - Moved many smaller helper functions in ciatah.m to their own M-files to make maintenance easier.
		% 2021.06.18 [20:28:28] - Added modelVarsFromFilesCheck support.
		% 2021.06.30 [11:14:12] - Added FRAMES_PER_SECOND_PLAYBACK.
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
		% 2021.10.20 [21:46:52] - Added fontSizeGuiScale.
		% 2022.02.09 [23:51:18] - Misc code fixes to conform to better Matlab language standards.
	% TODO
		%

	% dynamicprops is a subclass of handle, allowing addition of properties

	properties(GetAccess = 'public', SetAccess = 'public')
		% public read and write access.

		% FPS of movie(s) being analyzed
		FRAMES_PER_SECOND =  5;
		% FPS of movie playback
		FRAMES_PER_SECOND_PLAYBACK =  20;
		% Int: what factor temporally are analyzed movie from raw data
		DOWNSAMPLE_FACTOR =  4;
		% Float: estimated um per pixel
		MICRON_PER_PIXEL =  2.51; % 2.37;

		% Int: set the default UI font size
		fontSizeGui = 10;

		% Float: allows font size in GUI to be scaled
		fontSizeGuiScale = 1;

		defaultObjDir = pwd;
		classVersion = ciapkg.version();
		serverPath = '';
		privateSettingsPath = [ciapkg.getDir() filesep 'private' filesep 'settings' filesep 'privateLoadBatchFxns.m'];
		% place where functions can temporarily story user settings
		functionSettings = struct(...
			'null', NaN...
		);

		% user information
		userName = 'USA';
		% Taxonomic species name, defaults to mouse.
		species = 'mus musculus';

		% counters and stores
		% index of current folder
		fileNum = 1;
		% same as fileNum, will transfer to this since more clear
		folderNum = 1;
		% number to current stimulus index
		stimNum = 1;

		% Github repo path, for updating the repository code
		githubUrl = 'https://github.com/bahanonu/calciumImagingAnalysis';
		% Folder to put downloaded external programs
		externalProgramsDir = ciapkg.getDirExternalPrograms();

		% String: name of the analysis file to put in a folder to indicate to other computers the current computer is analyzing the folder and they should skip
		concurrentAnalysisFilename = '_currentlyAnalyzingFolderCheck.mat';

		% Cell array strings: List of methods to fast track folder to analyze dialog or skip altogether
		methodExcludeList = {'showVars','showFolders','setMainSettings','loadDependencies','saveObj','setStimulusSettings','modelDownsampleRawMovies','setMovieInfo','setup','update'};
		methodExcludeListVer2 = {'modelEditStimTable','behaviorProtocolLoad','modelPreprocessMovie','modelModifyMovies','removeConcurrentAnalysisFiles'};
		methodExcludeListStimuli = {'modelVarsFromFiles'};

		% 0 = load variables from disk, reduce RAM usage. 1 = load from disk to ram, faster for analysis.
		loadVarsToRam = 0;
		% show GUI for view functions?
		guiEnabled = 1;
		% indices for folders to analyze, [] = all
		foldersToAnalyze = [];
		% indices for stimuli to analyze, [] = all
		discreteStimuliToAnalyze = [];
		% io settings
		fileFilterRegexp = 'crop';
		% Regular expression for alternative file
		fileFilterRegexpAlt = 'crop';
		% Regular expression for alternative file during cell extraction
		fileFilterRegexpAltCellExtraction = '';
		% raw movie
		fileFilterRegexpRaw = 'concat';
		% behavior video regexp
		behaviorVideoRegexp = '';
		% loop over all files during analysis? 'individual' or 'group'
		analysisType  = 'group';
		% 1 = perform certain analysis on dF/F instead of peaks
		dfofAnalysis = 0;
		% 'filtered' returns auto/manually filtered signals/images, 'raw' returns raw
		modelGetSignalsImagesReturnType = 'filtered'
		% name of input dataset name for preprocessing
		inputDatasetName = '/1';
		% name of output dataset name for preprocessing
		outputDatasetName = '/1';
		%
		stimTriggerOnset = 0;
		% paths for specific types of files
		currentDateTimeSessionStr = datestr(now,'yyyymmdd_HHMMSS_FFF','local');
		currentDateTimeStr = datestr(now,'yyyymmdd','local');
		picsSavePath = [ciapkg.getDir() filesep 'private' filesep 'pics' filesep datestr(now,'yyyymmdd','local') filesep];
		dataSavePath = [ciapkg.getDir() filesep 'private' filesep 'data' filesep datestr(now,'yyyymmdd','local') filesep];
		dataSavePathFixed = [ciapkg.getDir() filesep 'private' filesep 'data' filesep];
		logSavePath = [ciapkg.getDir() filesep 'private' filesep 'logs' filesep datestr(now,'yyyymmdd','local') filesep];
		settingsSavePath = [ciapkg.getDir() filesep 'private' filesep 'settings'];
		%
		dataSaveFilenameModifier = '';
		% table save
		delimiter = ',';
		% name of i/o HDF5 dataset names
		hdf5Datasetname = '/1';
		% type of images to save analysis as '-dpng','-dmeta','-depsc2'
		imgSaveTypes = {'-dpng'};
		% Custom subject regexp find string to pass to getFileInfo
		subjectRegexp = '(m|M|f|F|Mouse|mouse|Mouse_|mouse_)\d+';
		% colormap to be used
		% colormap = customColormap([]);
		% colormap = customColormap({[0 0 0.7],[1 1 1],[0.7 0 0]});
		% colormap = customColormap({[0 0 1],[1 1 1],[0.5 0 0],[1 0 0]});
		% colormapAlt = customColormap({[0 0 0.7],[1 1 1],[0.7 0 0]});
		% colormapAlt2 = diverging_map(linspace(0,1,100),[0 0 0.7],[0.7 0 0]);
		colormap = parula(100);
		colormapAlt = parula(100);
		colormapAlt2 = parula(100);
		% colormap = customColormap({[27 52 93]/256,[1 1 1],[106 41 50]/256},'nPoints',50);
		% use for stimulus related viewing functions
		% frames before/after stimulus to look
		timeSequence = -25:25;
		postStimulusTimeSeq = 0:10;
		% bin analysis, integer only
		binDownsampleAmount = 1;
		% number of standard deviation for threshold crossing
		numStdsForThresh = 2.5;
		%
		stimulusTableValueName = 'frameSessionDownsampled';
		stimulusTableFrameName = 'frameSessionDownsampled';
		stimulusTableResponseName = 'response';
		stimulusTableTimeName = 'time';
		stimulusTableSessionName = 'trial';

		% methods
		currentMethod = 'modelAddNewFolders';

		% Region analysis
		regionModSaveStr = '_regionModSelectUser.mat'

		usrIdxChoiceStr = {...
			'PCAICA',...
			'CNMF',...
			'CNMFE',...
			'ROI',...
			'CELLMax',...
			'EXTRACT',...
			'EM'};
		usrIdxChoiceDisplay = {...
			'PCAICA (Mukamel, 2009)',...
			'CNMF (Giovannucci, 2019)',...
			'CNMF-E (Zhou, 2018)',...
			'ROI',...
			'CELLMax (Kitch/Ahanonu)',...
			'EXTRACT (Inan, 2021)',...
			'CELLMax [EM] (Kitch/Ahanonu)'...
			};

		% PCAICA, EM, EXTRACT, CNMF, CNMFE
		signalExtractionMethod = 'PCAICA';
		% signalExtractionMethod = 'EM';
		% Int: indicates which trace output to use from 1 = primary trace, 2 = secondary trace, etc.
		signalExtractionTraceOutputType = 1;

		extractionMethodSaveStr = struct(...
			'PCAICA', '_pcaicaAnalysis',...
			'EM', '_emAnalysis',...
			'CELLMax', '_cellmaxAnalysis',...
			'EXTRACT', '_extractAnalysis',...
			'CNMF', '_cnmfAnalysis',...
			'CNMFE', '_cnmfeAnalysis',...
			'ROI', '_roiAnalysis'...
		);
		extractionMethodStructSaveStr = struct(...
			'PCAICA', '_pcaicaAnalysis.mat',...
			'EM', '_emAnalysis.mat',...
			'CELLMax', '_cellmaxAnalysis.mat',...
			'EXTRACT', '_extractAnalysis.mat',...
			'CNMF', '_cnmfAnalysis.mat',...
			'CNMFE', '_cnmfeAnalysis.mat',...
			'ROI', '_roiAnalysis.mat'...
		);
		extractionMethodSortedSaveStr = struct(...
			'PCAICA', '_ICdecisions.mat',...
			'EM', '_emAnalysisSorted.mat',...
			'CELLMax', '_cellmaxAnalysisSorted.mat',...
			'EXTRACT', '_extractAnalysisSorted.mat',...
			'CNMF', '_cnmfAnalysisSorted.mat',...
			'CNMFE', '_cnmfeAnalysisSorted.mat',...
			'ROI', '_roiAnalysisSorted.mat'...
		);
		extractionMethodClassifierSaveStr = struct(...
			'PCAICA', '_ICclassifierDecisions.mat',...
			'EM', '_emAnalysisClassifierDecisions.mat',...
			'CELLMax', '_cellmaxAnalysisClassifierDecisions.mat',...
			'EXTRACT', '_extractAnalysisClassifierDecisions.mat',...
			'CNMF', '_cnmfAnalysisClassifierDecisions.mat',...
			'CNMFE', '_cnmfeAnalysisClassifierDecisions.mat',...
			'ROI', '_roiAnalysisClassifierDecisions.mat'...
		);
		extractionMethodValidVarname = struct(...
			'PCAICA', 'valid',...
			'EM', 'validCellMax',...
			'CELLMax', 'validCellMax',...
			'EXTRACT', 'validEXTRACT',...
			'CNMF', 'validCNMF',...
			'CNMFE', 'validCNMFE',...
			'ROI', 'validROI'...
		);
		extractionMethodStructVarname = struct(...
			'PCAICA', 'pcaicaAnalysisOutput',...
			'EM', 'emAnalysisOutput',...
			'CELLMax', 'cellmaxAnalysisOutput',...
			'EXTRACT', 'extractAnalysisOutput',...
			'CNMF', 'cnmfAnalysisOutput',...
			'CNMFE', 'cnmfeAnalysisOutput',...
			'ROI', 'roiAnalysisOutput'...
		);

		% PCAICA names
		rawPCAICAStructSaveStr = '_pcaicaAnalysis.mat';
		rawICfiltersSaveStr = '_ICfilters.mat';
		rawICtracesSaveStr = '_ICtraces.mat';
		sortedICfiltersSaveStr = '_ICfilters_sorted.mat';
		sortedICtracesSaveStr = '_ICtraces_sorted.mat';
		sortedICdecisionsSaveStr = '_ICdecisions.mat';
		classifierICdecisionsSaveStr = '_ICclassifierDecisions.mat';
		structPCAICAVarname = 'pcaicaAnalysisOutput';
		validPCAICAStructVarname = 'valid';
		% ROI names
		rawROIStructSaveStr = '_roiAnalysis.mat';
		sortedRoiStructSaveStr = '_roiAnalysisSorted.mat';
		structROIVarname = 'roiAnalysisOutput';
		classifierRoiStructSaveStr = '_roiAnalysisClassifierDecisions.mat';
		rawROItracesSaveStr = '_ROItraces.mat';
		validRoiStructVarname = 'validROI';
		% EM names
		rawEMStructSaveStr = '_emAnalysis.mat';
		sortedEMStructSaveStr = '_emAnalysisSorted.mat';
		classifierEMStructSaveStr = '_emAnalysisClassifierDecisions.mat';
		structEMVarname = 'emAnalysisOutput';
		validEMStructVarname = 'validCellMax';
		% EXTRACT names
		rawEXTRACTStructSaveStr = '_extractAnalysis.mat';
		sortedEXTRACTStructSaveStr = '_extractAnalysisSorted.mat';
		classifierEXTRACTStructSaveStr = '_extractAnalysisClassifierDecisions.mat';
		validEXTRACTStructVarname = 'validEXTRACT';
		structEXTRACTVarname = 'extractAnalysisOutput';
		% CNMF names
		rawCNMFStructSaveStr = '_cnmfAnalysis.mat';
		sortedCNMFStructSaveStr = '_cnmfAnalysisSorted.mat';
		classifierCNMFStructSaveStr = '_cnmfAnalysisClassifierDecisions.mat';
		validCNMFStructVarname = 'validCNMF';
		structCNMRVarname = 'cnmfAnalysisOutput';

		% NWB-related properties
		% Whether to write or load NWB files, e.g. from modelExtractSignalsFromMovie
		nwbLoadFiles = 0;
		% Str: sub-folder where NWB files are stored. Leave blank to load from current folder.
		nwbFileFolder = 'nwbFiles';
		% Str: blank, use calciumImagingAnalysis regexp, else force use of this regexp for NWB files
		nwbFileRegexp = '';
		% Name of H5 group for movies, images, and signal series in NWB files
		nwbGroupImages = '/processing/ophys/ImageSegmentation/PlaneSegmentation';
		nwbGroupSignalSeries = '/processing/ophys/Fluorescence/RoiResponseSeries';
		nwbGroupMovie = '/acquisition/TwoPhotonSeries/data';

		settingOptions = struct(...
			'analysisType',  {{'group','individual'}},...
			'loadVarsToRam', {{0,1}},...
			'guiEnabled', {{0,1}},...
			'dfofAnalysis', {{0,1}},...
			'picsSavePath', {{['private' filesep 'pics' filesep datestr(now,'yyyymmdd','local') filesep]}},...
			'delimiter', {{',','tab'}},...
			'imgSaveTypes', {{'-dpng','-dmeta','-depsc2'}}...
		);

		filterImageOptions = struct(...
			'minNumPixels', 10,...
			'maxNumPixels', 100,...
			'SNRthreshold', 1.45,...
			'minPerimeter', 5,...
			'maxPerimeter', 50,...
			'minSolidity', 0.8,...
			'minEquivDiameter', 3,...
			'maxEquivDiameter', 30,...
			'slopeRatioThreshold', 0.04...
		);

		downsampleRawOptions = struct(...
			'folderListInfo','USER_PATH',...
			'downsampleSaveFolder','USER_PATH',...
			'downsampleSrcFolder','USER_PATH',...
			'downsampleFactor','4',...
			'fileFilterRegexp','(recording.*.hdf5|recording.*.tif|.*.isxd)',...
			'datasetName','/images',...
			'maxChunkSize','25000',...
			'srcFolderFilterRegexp','201\d',...
			'srcSubfolderFileFilterRegexp','(recording.*.(txt|xml)|*.gpio)',...
			'srcSubfolderFileFilterRegexpExt','(.txt|.xml|.gpio)',...
			'downsampleSaveFolderTwo','',...
			'downsampleFactorTwo','2',...
			'outputDatasetName','/1'...
		);

		% Pre-processing options
		% Empty by default, if user opts to keep pre-process settings, this is saved as a struct.
		preprocessSettings = [];
		saveLoadPreprocessingSettings = 1;
		motionCorrectionRefFrame = 100;

		% io folders
		inputFolders = {};
		videoDir = '';
		videoSaveDir = '';
		trackingDir = '';
		stimulusDir = '';
		% if want to automatically save object to a specific location.
		objSaveLocation = [];

		% signal related
		% either the raw signals (traces) or
		rawSignals = {};
		% secondary either the raw signals (traces) or
		rawSignals2 = {};
		%
		rawImages = {};
		% computed signal peaks/locations, to reduce computation in functions
		signalPeaks = {};
		%
		signalPeaksArray = {};
		% computed centroid locations {[x y],...}
		objLocations = {};
		% mean correlation coefficient between image and movie
		imageMovieCorr = {};
		% cellmaps indicating which cells were filtered
		rawImagesFiltered = {};
		% structure of classifier structures, each field in the property should be named after the folder's subject, e.g. classifierStructs.m667 is the classification structure for subject m667
		classifierStructs = {};
		% structure of classifier structures for each folder, e.g. after running classification on each
		classifierFolderStructs = {};
		% cell array with {signalNo}.signalFeatures, {signalNo}.imageFeatures
		classifierFeatures = {};
		% cell array with {signalNo}.signalFeatures, {signalNo}.imageFeatures
		classifierImageFeaturesNames = {'Eccentricity','EquivDiameter','Area','Orientation','Perimeter','Solidity'};
		% structure for all valid classifications to go
		valid = {};
		% Automated or manual classification
		validManual = {};
		% from automated classification
		validAuto = {};
		% valid cells based on a regional modification
		validRegionMod = {};
		% polygon vertices from previously selected regions
		validRegionModPoly = {};
		% whether or not rawSignals/rawImages have been replaced by only valid signals, hence, ignore validManual/validAuto
		validPurge = 0;
		% ROI to use for exclusion analysis
		analysisROIArray = {};
		% number of expected [PCs ICs] for PCA-ICA, alter for other procedures
		numExpectedSignals = {};

		% subject info
		% all are cell array of strings or numbers as specified in the name
		dataPath = {};
		subjectNum = {};
		subjectStr = {};
		subjectProtocolStr = {};
		assay = {};
		protocol = {};
		assayType = {};
		assayNum = {};
		imagingPlane = {};
		imagingPlaneNum = {};
		date = {};
		fileIDArray = {};
		fileIDNameArray = {};
		folderBaseSaveStr = {};
		folderBasePlaneSaveStr = {};
		folderBaseDisplayStr = {};
		folderBaseSaveStrUnique = {};

		% path to CSV/TAB file or matlab table containing trial information and frames when stimuli occur
		discreteStimulusTable = {};
		% cell array of strings
		stimulusNameArray = {};
		% cell array of strings, used for saving pictures/etc.
		stimulusSaveNameArray = {};
		% cell array of numbered values for stimulus, e.g. {65,10}
		stimulusIdArray = {};
		% [1 numTrialFrames] vectors with 1 for when stimulus occurs
		stimulusVectorArray = {};
		% vector sequence before/after to analyze stimulus
		stimulusTimeSeq = {};

		% path to a CSV/TAB file
		continuousStimulusTable = {};
		% cell array of strings
		continuousStimulusNameArray = {};
		% cell array of strings, used for saving pictures/etc.
		continuousStimulusSaveNameArray = {};
		% cell array of numbered values for stimulus, e.g. {65,10}
		continuousStimulusIdArray = {};
		% [1 numTrialFrames] vectors with 1 for when stimulus occurs
		continuousStimulusVectorArray = {};
		% vector sequence before/after to analyze stimulus
		continuousStimulusTimeSeq = {};

		% session trial names
		sessionTrialNames = '';

		% behavior metrics
		behaviorMetricTable = {};
		behaviorMetricNameArray = {};
		behaviorMetricIdArray = {};

		% List of available calciumImagingAnalysis methods
		methodsList = {...
		'------- SETUP -------';
		'modelAddNewFolders';
		'help';
		'setup';
		'update';
		'loadDependencies';
		'setMovieInfo';
		'resetMijiClass';
		'';
		'------- CLASS/BEHAVIOR -------';
		'showVars';
		'showFolders';
		'saveObj';
		'initializeObj';
		'setMainSettings';
		'';
		'------- DATA CHECK/LOAD/EXPORT -------';
		'modelGetFileInfo';
		'modelVerifyDataIntegrity';
		'modelBatchCopyFiles';
		'modelLoadSaveData';
		'modelExportData';
		'';
		'------- VISUALIZATION/PREPROCESS CHECK -------';
		'viewMovieFiltering';
		'viewMovieRegistrationTest';
		'viewMovie';
		'';
		'------- PREPROCESS -------';
		'modelDownsampleRawMovies';
		'modelPreprocessMovie';
		'modelModifyMovies';
		'removeConcurrentAnalysisFiles';
		'';
		'------- CELL/SIGNAL EXTRACTION -------';
		'modelExtractSignalsFromMovie';
		'viewCellExtractionOnMovie';
		'removeConcurrentAnalysisFiles';
		'';
		'------- LOAD CELL-EXTRACTION/SIGNAL DATA -------';
		'modelVarsFromFiles';
		'';
		'------- SIGNAL SORTING -------';
		'computeManualSortSignals';
		'modelModifyRegionAnalysis';
		'';
		'------- PREPROCESS VERIFICATION -------';
		'viewMovie';
		'viewObjmaps';
		'viewCreateObjmaps';
		'viewMovieCreateSideBySide';
		'';
		'------- ACROSS SESSION ANALYSIS: COMPUTE/VIEW -------';
		'viewSubjectMovieFrames';
		'computeMatchObjBtwnTrials';
		'';
		'viewMatchObjBtwnSessions';
		'modelSaveMatchObjBtwnTrials';
		'computeCellDistances';
		'computeCrossDayDistancesAlignment';
		'';
		'------- TRACKING -------';
		'modelTrackingData';
		'viewOverlayTrackingToVideo'
		};

		% Tooltips for calciumImagingAnalysis methods
		tts = struct(...
			'SETUP', 'Methods for setting up the class.',...
			'modelAddNewFolders', 'Add new folders to calciumImagingAnalysis.',...
			'help', 'Some tips if running into trouble using calciumImagingAnalysis.',...
			'setup', 'Setup calciumImagingAnalysis if running for the 1st time or need to re-install dependencies (e.g. Fiji).',...
			'update', 'Update calciumImagingAnalysis, direct user toward relevant web address.',...
			'loadDependencies', 'Download and setup calciumImagingAnalysis dependencies.',...
			'setMovieInfo', 'Set information',...
			'resetMijiClass', 'Reset Fiji (e.g. Miji) if having problems loading.',...
			'CLASS_BEHAVIOR', 'Methods for modifying calciumImagingAnalysis behavior or saving.',...
			'showVars', 'Show all properties for this instance of calciumImagingAnalysis.',...
			'showFolders', 'Show all folders users have loaded into the class.',...
			'saveObj', 'Save this instance of the object for later loading.',...
			'initializeObj', 'IGNORE.',...
			'setMainSettings', 'IGNORE.',...
			'DATA_CHECK_LOAD_EXPORT', 'Methods for verify data and updating calciumImagingAnalysis meta-data.',...
			'modelGetFileInfo', 'Update information about each folder loaded into the class. Information derived from folder name.',...
			'modelVerifyDataIntegrity', 'Various sub-methods to verify data is valid and that certain files are present in each folder.',...
			'modelBatchCopyFiles', 'Batch copy files from one location to another or to move files to sub-folders within each folder (e.g. if users want to archive particular files).',...
			'modelLoadSaveData', 'IGNORE.',...
			'modelExportData', 'IGNORE.',...
			'PREPROCESS', 'Methods for pre-processing imaging movies.',...
			'modelDownsampleRawMovies', 'Spatially downsample large raw movies. Done in chunks so can process movies larger than available RAM.',...
			'viewMovieFiltering', 'IGNORE.',...
			'viewMovieRegistrationTest', 'Allows users to test several pre-processing setting to choose one that works best for their dataset.<br>Mainly focused on motion correction and spatial filtering.',...
			'viewMovie', 'Allows users to view their movie in several ways for each folder.',...
			'modelPreprocessMovie', 'Main pre-processing method to do motion correction, spatial filtering, calculate relative fluorescence, etc.',...
			'modelModifyMovies', 'Used to add user-defined borders to the movie.<br>e.g. if want to blank out areas outside GRIN lens so they do not interfere with cell extraction, etc.',...
			'removeConcurrentAnalysisFiles', 'Removes temporary files created when running analysis across multiple workstations in "modelPreprocessMovie" or "modelExtractSignalsFromMovie".',...
			'CELL_SIGNAL_EXTRACTION', 'Methods related to cell extraction.',...
			'modelExtractSignalsFromMovie', 'Main method to allow pre-processing',...
			'viewCellExtractionOnMovie', 'Overlays cell-extraction outputs on the imaging movie to allow users to verify cell shape, that cells were not missed, etc.',...
			'LOAD_CELL_EXTRACTION_SIGNAL_DATA', 'Methods relating to loading cell-extraction data into calciumImagingAnalysis.',...
			'modelVarsFromFiles', 'Load cell-extraction data into calciumImagingAnalysis.<br>RUN EACH TIME re-loading calciumImagingAnalysis and have already run cell extraction on the folders added to calciumImagingAnalysis.',...
			'SIGNAL_SORTING', 'Methods relating to sorting/classification of cell-extraction outputs.',...
			'computeManualSortSignals', 'A GUI that allows ',...
			'modelModifyRegionAnalysis', '',...
			'PREPROCESS_VERIFICATION', '',...
			'viewObjmaps', 'Creates several figures displaying cell-extraction outputs from each selected folder.',...
			'viewCreateObjmaps', 'Creates several figures displaying cell-extraction outputs from each selected folder.',...
			'viewSubjectMovieFrames', 'Takes frames or cell-extraction cell maps from each folder associated with an animal and displays them as a movie to help aid in visual inspection of cross-session alignment possibilities.',...
			'viewMovieCreateSideBySide', 'Combines synchronized imaging and behavior or other experimental videos to allow side-by-side comparisons.',...
			'TRACKING', 'Methods related to tracking animals.',...
			'modelTrackingData', 'After running ImageJ-based tracking, use this module to clean-up the data and overlay it on the behavior movie.',...
			'viewOverlayTrackingToVideo', 'Creates a video where the animal''s computed centroid and velocity is overlaid on the video.',...
			'ACROSS_SESSION_ANALYSIS__COMPUTE_VIEW', 'Methods related to aligning cells across imaging sessions.',...
			'computeMatchObjBtwnTrials', 'Main method to match cell-extraction outputs (e.g. cells) across imaging sessions.',...
			'viewMatchObjBtwnSessions', 'Allows users to visualize cross session matches and outputs videos to help with assessment as well.',...
			'modelSaveMatchObjBtwnTrials', 'Will allow users to save output of cross-session matching.',...
			'computeCellDistances', 'Computes the cell-cell distance for all cells in each imaging session and outputs as a CSV table.<br>Should be used to assess the largest distance to still count cells as matching during cross-session matching.',...
			'computeCrossDayDistancesAlignment', 'Computes the cell-cell distance of all cross-session matched cells.<br>e.g. ideally all will be below the cross-session cell distance cutoff.'...
		);
	end
	properties(GetAccess = 'public', SetAccess = 'private')
		% public read access, but private write access.

		% summary statistics data save stores
		sumStats = {};
		detailStats = {};
		tmpStats = {};
		saveDataType = {};

		% counters and stores
		figNames = {};
		figNo = {};
		figNoAll = 777;

		% signal related
		nSignals = {};
		nFrames = {};
		signalPeaksCopy = {};
		alignedSignalArray = {};
		alignedSignalShuffledMeanArray = {};
		alignedSignalShuffledStdArray = {};

		% stimulus
		% reorganize discreteStimulusTable into stimulus structures to reduce memory footprint
		discreteStimulusArray = {};
		discreteStimMetrics = {};

		% stimulus
		% reorganize discreteStimulusTable into stimulus structures to reduce memory footprint
		continuousStimulusArray = {};
		continuousStimMetrics = {};

		% behavior metric
		% reorganize behaviorMetricTable into stimulus structures to reduce memory footprint/provide common IO
		behaviorMetricArray = {};

		% for each folder have a structure containing different metrics.
		neuroBehaviorMetric = {};

		% distance metrics
		distanceMetric = {};
		distanceMetricShuffleMean = {};
		distanceMetricShuffleStd = {};

		% correlation metrics
		corrMatrix = {};

		% significant signals, different variables for controlling which signals are statistically significant, given some test
		currentSignificantArray = [];
		significantArray = {};
		significantStatsArray = {};
		sigModSignals = {};
		sigModSignalsAll = {};
		ttestSignSignals = {};

		% cross session alignment
		globalIDs = [];
		globalIDCoords = {};
		globalIDFolders = {};
		globalIDImages = {};
		globalRegistrationCoords = {};
		globalObjectMapTurboreg = [];
		globalStimMetric = [];
		globalIDStruct = {};
	end
	properties(GetAccess = 'private', SetAccess = 'private')
		% private read and write access
	end
	properties(Constant = true)
		% cannot be changed after object is created
		% MICRON_PER_PIXEL =  2.37;
	end

	methods
		% methods, including the constructor are defined in this block
		function obj = ciatah(varargin)
			% CLASS CONSTRUCTOR
			warning on;
			clc
			% ' Calcium Imaging Analysis Class
			disp([...
			ciapkg.pkgName ' (calciumImagingAnalysis)' 10 ...
			'A software package for analyzing one- and two-photon calcium imaging datasets.' 10 10 ...
			'Biafra Ahanonu <<a href="emailto:bahanonu@alum.mit.edu">bahanonu@alum.mit.edu</a>>' 10 ...
			'Version ' obj.classVersion 10 10 ...
			'Made in USA']);

			obj.usaflag;

			disp(repmat('*',1,42))
			disp(['Constructing ' ciapkg.pkgName ' imaging analysis object...'])

			% Ensure that default directory is the CIAtah repository root
			functionLocation = dbstack('-completenames');
			functionLocation = functionLocation(1).file;
			[functionDir,~,~] = fileparts(functionLocation);
			[functionDir,~,~] = fileparts(functionDir);
			obj.defaultObjDir = functionDir;
			clear functionDir functionLocation;

			% Make sure all functions are loaded
			obj.loadBatchFunctionFolders();
			disp(repmat('*',1,42))

			% Make sure defaults set
			try
				setFigureDefaults();
				obj.classVersion = ciapkg.version();
			catch
			end

			% Because the obj
			%========================
			% obj.exampleOption = '';
			% get options
			obj = ciapkg.io.getOptions(obj,varargin);
			% disp(options)
			% unpack options into current workspace
			% fn=fieldnames(options);
			% for i=1:length(fn)
			%	 eval([fn{i} '=options.' fn{i} ';']);
			% end
			%========================

			obj = initializeObj(obj);

			ciapkg.io.updatePkg('updatePackage',0);

			disp(repmat('*',1,42))
			disp(['Done initializing ' ciapkg.pkgName '!'])
			disp(repmat('*',1,42))

			disp([...
			'Run processing pipeline by typing below (or clicking link) into command window (no semi-colon!):' 10 ...
			'<a href="matlab: obj">obj</a>'])

			% disp([...
			% 'Run processing pipeline by typing into command window:' 10 ...
			% '<a href="">obj.runPipelineProcessing</a>' 10 ...
			% 'or for advanced features: ' 10 ...
			% '<a href="">obj.runPipeline</a>' 10])
		end
		% getter and setter functions
		function dataPath = get.dataPath(obj)
			dataPath = obj.dataPath;
		end
	end
	methods(Static = true)
		% functions that are related but not dependent on instances of the class
		% function obj = loadObj(oldObj)
		%	  [filePath,folderPath,~] = uigetfile('*.*','select text file that points to analysis folders','example.txt');
		%	  % exit if user picks nothing
		%	  % if folderListInfo==0; return; end
		%	  load([folderPath filesep filePath]);
		%	  oldObj = obj;
		%	  obj = calciumImagingAnalysis;
		%	  obj = getOptions(obj,oldObj);
		% end
	end
	methods(Access = private)
		% methods only executed by other class methods

		% model methods, usually for input-output like saving information to files
		% for obtaining the current stim from tables
		[behaviorMetric] = modelGetBehaviorMetric(obj,inputID)
	end
	methods(Access = protected)
		% methods only executed by other class methods, also available to subclasses
	end
	methods(Access = public)
		% these are in separate M-files

		% calciumImagingAnalysis main class methods
		[idNumIdxArray, validFoldersIdx, ok] = ciatahMainGui(obj,fxnsToRun,inputTxt,currentIdx)
		[idNumIdxArray, ok] = pipelineListBox(obj,fxnsToRun,inputTxt,currentIdx)
		obj = computeCrossDayDistancesAlignment(obj)
		obj = computeCellDistances(obj)
		obj = resetMijiClass(obj)
		obj = runPipeline(obj,varargin)
		obj = initializeObj(obj)
		obj = loadBatchFunctionFolders(obj)
		obj = loadDependencies(obj)
		obj = removeConcurrentAnalysisFiles(obj)
		obj = setMovieInfo(obj)
		[validFoldersIdx] = pipelineFolderFilter(obj,useAltValid,validFoldersIdx)

		[output, output2] = modelGetStim(obj,idNum,varargin)

		% view help about the object
		obj = help(obj)

		% specific experiments

		% compute methods, performs some computation and returns calculation to class property
		obj = computeMatchObjBtwnTrials(obj,varargin)
		obj = computeManualSortSignals(obj)
		obj = computeClassifyTrainSignals(obj)

		% just need stimulus files

		% view methods, for displaying charts
		% no prior computation
		obj = viewCreateObjmaps(obj,varargin)

		% pre-processing checking
		obj = viewMovie(obj)
		obj = viewMovieFiltering(obj)
		obj = viewMovieRegistrationTest(obj)
		obj = viewMovieCreateSideBySide(obj)
		obj = modelModifyMovies(obj)
		obj = viewCellExtractionOnMovie(obj,varargin)

		% require pre-computation, individual

		% require pre-computation, group, global alignment
		obj = viewMatchObjBtwnSessions(obj,varargin)
		% movies

		% require pre-computation and behavior metrics, individual

		% requires manual sorting

		% need tracking
		obj = viewOverlayTrackingToVideo(obj)

		% model methods, usually for input-output like saving information to files
		obj = modelAddNewFolders(obj,varargin)
		obj = modelExportData(obj,varargin)
		obj = modelReadTable(obj,varargin)
		obj = modelTableToStimArray(obj,varargin)
		obj = modelGetFileInfo(obj)
		obj = modelVerifyDataIntegrity(obj)
		obj = modelSaveImgToFile(obj,saveFile,thisFigName,thisFigNo,thisFileID,varargin)
		obj = modelSaveSummaryStats(obj,varargin)
		obj = modelSaveDetailedStats(obj)
		obj = modelVarsFromFiles(obj)
		obj = modelModifyRegionAnalysis(obj,varargin)
		obj = modelExtractSignalsFromMovie(obj,varargin)
		obj = modelPreprocessMovie(obj)
		obj = modelDownsampleRawMovies(obj)
		obj = modelBatchCopyFiles(obj)
		obj = modelLoadSaveData(obj)
		obj = modelSaveMatchObjBtwnTrials(obj,varargin)

		% helps clean and load tracking data
		obj = modelTrackingData(obj)

		% helper
		[inputSignals, inputImages, signalPeaks, signalPeaksArray, valid, validType, inputSignals2] = modelGetSignalsImages(obj,varargin)
		[fileIdxArray, idNumIdxArray, nFilesToAnalyze, nFiles] = getAnalysisSubsetsToAnalyze(obj)
		[preprocessSettingStruct, preprocessingSettingsAll] = getRegistrationSettings(obj,inputTitleStr,varargin)

		% set methods, for IO to specific variables in a controlled manner
		obj = setMainSettings(obj)

		% Helper functions
		% obj = delete(obj)
		obj = setup(obj)
		obj = update(obj)
		obj = display(obj)
		obj = showVars(obj)
		obj = downloadLatestGithubVersion(obj)
		obj = showProtocolSubjectsSessions(obj)
		obj = showFolders(obj)
		[valid] = getValid(obj,validType)
		obj = changeCaxis(obj)
		obj = changeFont(obj)
		obj = checkToolboxes(obj)
		obj = GetSize(obj)
		obj = makeFolderDirs(obj)
		[flagOut] = usaflag(obj)
		obj = saveObj(obj)
		obj = modelVarsFromFilesCheck(obj,folderNo,varargin)
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function h_out= axis_getQuietly(varargin)
%AXIS_GETQUIETLY - Make an axis the current axis without interfering the
% 'Visible' property

vis= get(gcf, 'Visible');
if nargin==1,
  axes(varargin{:});
else
  out= axes(varargin{:});
end
set(gcf, 'Visible',vis);

if nargout>0,
  h_out= out;
end

author: bbci