text:
function [M,scal] = spm_affreg(VG,VF,flags,M,scal)
% Affine registration using least squares.
% FORMAT [M,scal] = spm_affreg(VG,VF,flags,M0,scal0)
%
% VG        - Vector of template volumes.
% VF        - Source volume.
% flags     - a structure containing various options.  The fields are:
%             WG       - Weighting volume for template image(s).
%             WF       - Weighting volume for source image
%                        Default to [].
%             sep      - Approximate spacing between sampled points (mm).
%                        Defaults to 5.
%             regtype  - regularisation type.  Options are:
%                        'none'  - no regularisation
%                        'rigid' - almost rigid body
%                        'subj'  - inter-subject registration (default).
%                        'mni'   - registration to ICBM templates
%             globnorm - Global normalisation flag (1)
% M0        - (optional) starting estimate. Defaults to eye(4).
% scal0     - (optional) starting estimate.
%
% M         - affine transform, such that voxels in VF map to those in
%             VG by   VG.mat\M*VF.mat
% scal      - scaling factors for VG
%
% When only one template is used, then the cost function is approximately
% symmetric, although a linear combination of templates can be used.
% Regularisation is based on assuming a multi-normal distribution for the
% elements of the Henckey Tensor. See:
% "Non-linear Elastic Deformations". R. W. Ogden (Dover), 1984.
% Weighting for the regularisation is determined approximately according
% to:
% "Incorporating Prior Knowledge into Image Registration"
% J. Ashburner, P. Neelin, D. L. Collins, A. C. Evans & K. J. Friston.
% NeuroImage 6:344-352 (1997).
%
%_______________________________________________________________________
% @(#)spm_affreg.m	2.3 John Ashburner 03/02/18

if nargin<5, scal = ones(length(VG),1); end;
if nargin<4, M    = eye(4);             end;

def_flags = struct('sep',5, 'regtype','subj','WG',[],'WF',[],'globnorm',1,'debug',0);
if nargin < 2 | ~isstruct(flags),
	flags = def_flags;
else,
	fnms = fieldnames(def_flags);
	for i=1:length(fnms),
		if ~isfield(flags,fnms{i}),
			flags = setfield(flags,fnms{i},getfield(def_flags,fnms{i}));
		end;
	end;
end;

% Check to ensure inputs are valid...
% ---------------------------------------------------------------
if length(VF)>1, error('Can not use more than one source image'); end;
if ~isempty(flags.WF),
	if length(flags.WF)>1,
		error('Can only use one source weighting image');
	end;
	if any(any(VF.mat-flags.WF.mat)),
		error('Source and its weighting image must have same orientation');
	end;
	if any(any(VF.dim(1:3)-flags.WF.dim(1:3))),
		error('Source and its weighting image must have same dimensions');
	end;
end;
if ~isempty(flags.WG),
	if length(flags.WG)>1,
		error('Can only use one template weighting image');
	end;
	tmp = reshape(cat(3,VG(:).mat,flags.WG.mat),16,length(VG)+length(flags.WG));
else,
	tmp = reshape(cat(3,VG(:).mat),16,length(VG));
end;
if any(any(diff(tmp,1,2))),
	error('Reference images must all have the same orientation');
end;
if ~isempty(flags.WG),
	tmp = cat(1,VG(:).dim,flags.WG.dim);
else,
	tmp = cat(1,VG(:).dim);
end;
if any(any(diff(tmp(:,1:3),1,1))),
	error('Reference images must all have the same dimensions');
end;
% ---------------------------------------------------------------

% Generate points to sample from, adding some jitter in order to
% make the cost function smoother.
% ---------------------------------------------------------------
rand('state',0); % want the results to be consistant.
dg   = VG(1).dim(1:3);
df   = VF(1).dim(1:3);

if length(VG)==1,
	skip = sqrt(sum(VG(1).mat(1:3,1:3).^2)).^(-1)*flags.sep;
	[x1,x2,x3]=ndgrid(1:skip(1):dg(1)-.5, 1:skip(2):dg(2)-.5, 1:skip(3):dg(3)-.5);
	x1   = x1 + rand(size(x1))*0.5; x1 = x1(:);
	x2   = x2 + rand(size(x2))*0.5; x2 = x2(:);
	x3   = x3 + rand(size(x3))*0.5; x3 = x3(:);
end;

skip = sqrt(sum(VF(1).mat(1:3,1:3).^2)).^(-1)*flags.sep;
[y1,y2,y3]=ndgrid(1:skip(1):df(1)-.5, 1:skip(2):df(2)-.5, 1:skip(3):df(3)-.5);
y1   = y1 + rand(size(y1))*0.5; y1 = y1(:);
y2   = y2 + rand(size(y2))*0.5; y2 = y2(:);
y3   = y3 + rand(size(y3))*0.5; y3 = y3(:);
% ---------------------------------------------------------------

if flags.globnorm,
	% Scale all images approximately equally
	% ---------------------------------------------------------------
	for i=1:length(VG),
		VG(i).pinfo(1:2,:) = VG(i).pinfo(1:2,:)/spm_global(VG(i));
	end;
	VF(1).pinfo(1:2,:) = VF(1).pinfo(1:2,:)/spm_global(VF(1));
end;
% ---------------------------------------------------------------

if length(VG)==1,
	[G,dG1,dG2,dG3]  = spm_sample_vol(VG(1),x1,x2,x3,1);
	if ~isempty(flags.WG),
		WG = abs(spm_sample_vol(flags.WG,x1,x2,x3,1))+eps;
	end;
end;

[F,dF1,dF2,dF3]  = spm_sample_vol(VF(1),y1,y2,y3,1);
if ~isempty(flags.WF),
	WF = abs(spm_sample_vol(flags.WF,y1,y2,y3,1))+eps;
end;
% ---------------------------------------------------------------
n_main_its = 0;
ss         = Inf;
W          = [Inf Inf Inf];
est_smo    = 1;
% ---------------------------------------------------------------

for iter=1:256,
	pss   = ss;
	p0    = [0 0 0  0 0 0  1 1 1  0 0 0];

	% Initialise the cost function and its 1st and second derivatives
	% ---------------------------------------------------------------
	n     = 0;
	ss    = 0;
	Beta  = zeros(12+length(VG),1);
	Alpha = zeros(12+length(VG));

	if length(VG)==1,
		% Make the cost function symmetric
		% ---------------------------------------------------------------

		% Build a matrix to rotate the derivatives by, converting from
		% derivatives w.r.t. changes in the overall affine transformation
		% matrix, to derivatives w.r.t. the parameters p.
		% ---------------------------------------------------------------
		dt  = 0.0001;
		R   = eye(13);
		MM0 = inv(VG.mat)*inv(spm_matrix(p0))*VG.mat;
		for i1=1:12,
			p1          = p0;
			p1(i1)      = p1(i1)+dt;
			MM1         = (inv(VG.mat)*inv(spm_matrix(p1))*(VG.mat));
			R(1:12,i1)  = reshape((MM1(1:3,:)-MM0(1:3,:))/dt,12,1);
		end;
		% ---------------------------------------------------------------
		[t1,t2,t3] = coords((M*VF(1).mat)\VG(1).mat,x1,x2,x3);
		msk        = find((t1>=1 & t1<=df(1) & t2>=1 & t2<=df(2) & t3>=1 & t3<=df(3)));
		if length(msk)<32, error_message; end;
		t1         = t1(msk);
		t2         = t2(msk);
		t3         = t3(msk);
		t          = spm_sample_vol(VF(1), t1,t2,t3,1);

		% Get weights
		% ---------------------------------------------------------------
		if ~isempty(flags.WF) | ~isempty(flags.WG),
			if isempty(flags.WF),
				wt = WG(msk);
			else,
				wt = spm_sample_vol(flags.WF(1), t1,t2,t3,1)+eps;
				if ~isempty(flags.WG), wt = 1./(1./wt + 1./WG(msk)); end;
			end;
			wt = sparse(1:length(wt),1:length(wt),wt);
		else,
			wt = speye(length(msk));
			wt = [];
		end;
		% ---------------------------------------------------------------
		clear t1 t2 t3

		% Update the cost function and its 1st and second derivatives.
		% ---------------------------------------------------------------
		[AA,Ab,ss1,n1] = costfun(x1,x2,x3,dG1,dG2,dG3,msk,scal^(-2)*t,G(msk)-(1/scal)*t,wt);
		Alpha = Alpha + R'*AA*R;
		Beta  = Beta  + R'*Ab;
		ss    = ss    + ss1;
		n     = n     + n1;
		t     = G(msk) - (1/scal)*t;
	end;

	if 1,
		% Build a matrix to rotate the derivatives by, converting from
		% derivatives w.r.t. changes in the overall affine transformation
		% matrix, to derivatives w.r.t. the parameters p.
		% ---------------------------------------------------------------
		dt = 0.0001;
		R  = eye(12+length(VG));
		MM0 = inv(M*VF.mat)*spm_matrix(p0)*M*VF.mat;
		for i1=1:12,
			p1          = p0;
			p1(i1)      = p1(i1)+dt;
			MM1         = (inv(M*VF.mat)*spm_matrix(p1)*M*VF.mat);
			R(1:12,i1)  = reshape((MM1(1:3,:)-MM0(1:3,:))/dt,12,1);
		end;
		% ---------------------------------------------------------------
		[t1,t2,t3] = coords(VG(1).mat\M*VF(1).mat,y1,y2,y3);
		msk        = find((t1>=1 & t1<=dg(1) & t2>=1 & t2<=dg(2) & t3>=1 & t3<=dg(3)));
		if length(msk)<32, error_message; end;

		if length(msk)<32, error_message; end;
		t1 = t1(msk);
		t2 = t2(msk);
		t3 = t3(msk);
		t  = zeros(length(t1),length(VG));

		% Get weights
		% ---------------------------------------------------------------
		if ~isempty(flags.WF) | ~isempty(flags.WG),
			if isempty(flags.WG),
				wt = WF(msk);
			else,
				wt = spm_sample_vol(flags.WG(1), t1,t2,t3,1)+eps;
				if ~isempty(flags.WF), wt = 1./(1./wt + 1./WF(msk)); end;
			end;
			wt = sparse(1:length(wt),1:length(wt),wt);
		else,
			wt = speye(length(msk));
		end;
		% ---------------------------------------------------------------

		if est_smo,
			% Compute derivatives of residuals in the space of F
			% ---------------------------------------------------------------
			[ds1,ds2,ds3] = transform_derivs(VG(1).mat\M*VF(1).mat,dF1(msk),dF2(msk),dF3(msk));
			for i=1:length(VG),
				[t(:,i),dt1,dt2,dt3] = spm_sample_vol(VG(i), t1,t2,t3,1);
				ds1   = ds1 - dt1*scal(i); clear dt1
				ds2   = ds2 - dt2*scal(i); clear dt2
				ds3   = ds3 - dt3*scal(i); clear dt3
			end;
			dss   = [ds1'*wt*ds1 ds2'*wt*ds2 ds3'*wt*ds3];
			clear ds1 ds2 ds3
		else,
			for i=1:length(VG),
				t(:,i)= spm_sample_vol(VG(i), t1,t2,t3,1);
			end;
		end;

		clear t1 t2 t3

		% Update the cost function and its 1st and second derivatives.
		% ---------------------------------------------------------------
		[AA,Ab,ss2,n2] = costfun(y1,y2,y3,dF1,dF2,dF3,msk,-t,F(msk)-t*scal,wt);
		Alpha = Alpha  + R'*AA*R;
		Beta  = Beta   + R'*Ab;
		ss    = ss     + ss2;
		n     = n      + n2;
	end;

	if est_smo,
		% Compute a smoothness correction from the residuals and their
		% derivatives.  This is analagous to the one used in:
		% 	"Analysis of fMRI Time Series Revisited"
		% 	Friston KJ, Holmes AP, Poline JB, Grasby PJ, Williams SCR,
		% 	Frackowiak RSJ, Turner R.  Neuroimage 2:45-53 (1995).
		% ---------------------------------------------------------------
		vx     = sqrt(sum(VG(1).mat(1:3,1:3).^2));
		pW     = W;
		W      = (2*dss/ss2).^(-.5).*vx;
		W      = min(pW,W);
		if flags.debug, fprintf('\nSmoothness FWHM:  %.3g x %.3g x %.3g mm\n', W*sqrt(8*log(2))); end;
		if length(VG)==1, dens=2; else, dens=1; end;
		smo    = prod(min(dens*flags.sep/sqrt(2*pi)./W,[1 1 1]));
		est_smo=0;
		n_main_its = n_main_its + 1;
	end;

	% Update the parameter estimates
	% ---------------------------------------------------------------
	nu      = n*smo;
	sig2    = ss/nu;
	[d1,d2] = reg(M,12+length(VG),flags.regtype);

	soln    = (Alpha/sig2+d2)\(Beta/sig2-d1);
	scal    = scal - soln(13:end);
	M       = spm_matrix(p0 + soln(1:12)')*M;

	if flags.debug,
		fprintf('%d\t%g\n', iter, ss/n);
		piccies(VF,VG,M,scal,b)
	end;

	% If cost function stops decreasing, then re-estimate smoothness
	% and try again.  Repeat a few times.
	% ---------------------------------------------------------------
	ss = ss/n;
	if iter>1, spm_chi2_plot('Set',ss); end;
	if (pss-ss)/pss < 1e-6,
		est_smo = 1;
	end;
	if n_main_its>3, break; end;

end;
return;
%_______________________________________________________________________

%_______________________________________________________________________
function [X1,Y1,Z1] = transform_derivs(Mat,X,Y,Z)
% Given the derivatives of a scalar function, return those of the
% affine transformed function
%_______________________________________________________________________

t1 = Mat(1:3,1:3);
t2 = eye(3);
if sum((t1(:)-t2(:)).^2) < 1e-12,
        X1 = X;Y1 = Y; Z1 = Z;
else,
        X1    = Mat(1,1)*X + Mat(1,2)*Y + Mat(1,3)*Z;
        Y1    = Mat(2,1)*X + Mat(2,2)*Y + Mat(2,3)*Z;
        Z1    = Mat(3,1)*X + Mat(3,2)*Y + Mat(3,3)*Z;
end;
return;
%_______________________________________________________________________

%_______________________________________________________________________
function [d1,d2] = reg(M,n,typ)
% Analytically compute the first and second derivatives of a penalty
% function w.r.t. changes in parameters.

if nargin<3, typ = 'subj'; end;
if nargin<2, n   = 13;     end;

[mu,isig] = priors(typ);
ds  = 0.000001;
d1  = zeros(n,1);
d2  = zeros(n);
p0  = [0 0 0  0 0 0  1 1 1  0 0 0];
h0  = penalty(p0,M,mu,isig);
for i=7:12, % derivatives are zero w.r.t. rotations and translations
	p1    = p0;
	p1(i) = p1(i)+ds;
	h1    = penalty(p1,M,mu,isig);
	d1(i) = (h1-h0)/ds; % First derivative
	for j=7:12,
		p2    = p0;
		p2(j) = p2(j)+ds;
		h2    = penalty(p2,M,mu,isig);
		p3    = p1;
		p3(j) = p3(j)+ds;
		h3    = penalty(p3,M,mu,isig);
		d2(i,j) = ((h3-h2)/ds-(h1-h0)/ds)/ds; % Second derivative
	end;
end;
return;
%_______________________________________________________________________

%_______________________________________________________________________
function h = penalty(p,M,mu,isig)
% Return a penalty based on the elements of an affine transformation,
% which is given by:
% 	spm_matrix(p)*M
%
% The penalty is based on the 6 unique elements of the Hencky tensor
% elements being multinormally distributed.
%_______________________________________________________________________

% Unique elements of symmetric 3x3 matrix.
els = [1 2 3 5 6 9];

T = spm_matrix(p)*M;
T = T(1:3,1:3);
T = 0.5*logm(T'*T);
T = T(els)' - mu;
h = T'*isig*T;
return;
%_______________________________________________________________________

%_______________________________________________________________________
function [mu,isig] = priors(typ)
% The parameters for this distribution were derived empirically from 227
% scans, that were matched to the ICBM space.
%_______________________________________________________________________

mu   = zeros(6,1);
isig = zeros(6);
switch deblank(lower(typ)),

case 'mni', % For registering with MNI templates...
	mu   = [0.0667 0.0039 0.0008 0.0333 0.0071 0.1071]';
	isig = 1e4 * [
	    0.0902   -0.0345   -0.0106   -0.0025   -0.0005   -0.0163
	   -0.0345    0.7901    0.3883    0.0041   -0.0103   -0.0116
	   -0.0106    0.3883    2.2599    0.0113    0.0396   -0.0060
	   -0.0025    0.0041    0.0113    0.0925    0.0471   -0.0440
	   -0.0005   -0.0103    0.0396    0.0471    0.2964   -0.0062
	   -0.0163   -0.0116   -0.0060   -0.0440   -0.0062    0.1144];

case 'rigid', % Constrained to be almost rigid...
	mu   = zeros(6,1);
	isig = eye(6)*1e9;

case 'isochoric', % Volume preserving...
	error('Not implemented');

case 'isotropic', % Isotropic zoom in all directions...
	error('Not implemented');

case 'subj', % For inter-subject registration...
	mu   = zeros(6,1);
	isig = 1e3 * [
	    0.8876    0.0784    0.0784   -0.1749    0.0784   -0.1749
	    0.0784    5.3894    0.2655    0.0784    0.2655    0.0784
	    0.0784    0.2655    5.3894    0.0784    0.2655    0.0784
	   -0.1749    0.0784    0.0784    0.8876    0.0784   -0.1749
	    0.0784    0.2655    0.2655    0.0784    5.3894    0.0784
	   -0.1749    0.0784    0.0784   -0.1749    0.0784    0.8876];

case 'none', % No regularisation...
	mu   = zeros(6,1);
	isig = zeros(6);

otherwise,
	error(['"' typ '" not recognised as type of regularisation.']);
end;
return;

%_______________________________________________________________________
function [y1,y2,y3]=coords(M,x1,x2,x3)
% Affine transformation of a set of coordinates.
%_______________________________________________________________________

y1 = M(1,1)*x1 + M(1,2)*x2 + M(1,3)*x3 + M(1,4);
y2 = M(2,1)*x1 + M(2,2)*x2 + M(2,3)*x3 + M(2,4);
y3 = M(3,1)*x1 + M(3,2)*x2 + M(3,3)*x3 + M(3,4);
return;
%_______________________________________________________________________

%_______________________________________________________________________
function A = make_A(x1,x2,x3,dG1,dG2,dG3,t)
% Generate part of a design matrix using the chain rule...
% df/dm = df/dy * dy/dm
% where
% 	df/dm is the rate of change of intensity w.r.t. affine parameters
% 	df/dy is the gradient of the image f
% 	dy/dm crange of position w.r.t. change of parameters
%_______________________________________________________________________

A  = [x1.*dG1 x1.*dG2 x1.*dG3 ...
	  x2.*dG1 x2.*dG2 x2.*dG3 ...
	  x3.*dG1 x3.*dG2 x3.*dG3 ...
	      dG1     dG2     dG3    t];
return;
%_______________________________________________________________________

%_______________________________________________________________________
function [AA,Ab,ss,n] = costfun(x1,x2,x3,dG1,dG2,dG3,msk,lastcols,b,wt)
chunk = 10240;
lm    = length(msk);
AA    = zeros(12+size(lastcols,2));
Ab    = zeros(12+size(lastcols,2),1);
ss    = 0;
n     = 0;

for i=1:ceil(lm/chunk),
	ind  = (((i-1)*chunk+1):min(i*chunk,lm))';
	msk1 = msk(ind);

	A1   = make_A(x1(msk1),x2(msk1),x3(msk1),dG1(msk1),dG2(msk1),dG3(msk1),lastcols(ind,:));
	b1   = b(ind);
	if ~isempty(wt),
		wt1   = wt(ind,ind);
		AA    = AA  + A1'*wt1*A1;
		%Ab   = Ab  + A1'*wt1*b1;
		Ab    = Ab  + (b1'*wt1*A1)';
		ss    = ss  + b1'*wt1*b1;
		n     = n   + trace(wt1);
		clear wt1
	else,
		AA    = AA  + spm_atranspa(A1);
		%Ab   = Ab  + A1'*b1;
		Ab    = Ab  + (b1'*A1)';
		ss    = ss  + b1'*b1;
		n     = n   + length(msk1);
	end;
	clear A1 b1 msk1 ind
end;
return;
%_______________________________________________________________________

%_______________________________________________________________________
function error_message
% Display an error message for when things go wrong.
str = {	'There is not enough overlap in the images',...
	'to obtain a solution.',...
	' ',...
	'Please check that your header information is OK.'};
spm('alert*',str,mfilename,sqrt(-1));
error('insufficient image overlap')
return
%_______________________________________________________________________

%_______________________________________________________________________
function piccies(VF,VG,M,scal,b)
% This is for debugging purposes.
% It shows the linear combination of template images, the affine
% transformed source image, the residual image and a histogram of the
% residuals.
%_______________________________________________________________________

figure(2);
Mt = spm_matrix([0 0 (VG(1).dim(3)+1)/2]);
M  = (M*VF(1).mat)\VG(1).mat;
t  = zeros(VG(1).dim(1:2));
for i=1:length(VG);
	t  = t + spm_slice_vol(VG(i),  Mt,VG(1).dim(1:2),1)*scal(i);
end;
u  = spm_slice_vol(VF(1),M*Mt,VG(1).dim(1:2),1);
subplot(2,2,1);imagesc(t');axis image xy off
subplot(2,2,2);imagesc(u');axis image xy off
subplot(2,2,3);imagesc(u'-t');axis image xy off
subplot(2,2,4);hist(b,50); % Entropy of residuals may be a nice cost function?
drawnow;
return;
%_______________________________________________________________________

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function suite = test_assertVectorsAlmostEqual
initTestSuite;

%===============================================================================
function test_happyCase

A = [1 1e10];
B = [2 1e10];
% All code here should execute with no error.
assertVectorsAlmostEqual(A, B);
assertVectorsAlmostEqual(A, B, 'custom message');

%===============================================================================
function test_failedAssertion

A = [1 1e6];
B = [2 1e6];

f = @() assertVectorsAlmostEqual(A, B);
assertExceptionThrown(f, 'assertVectorsAlmostEqual:tolExceeded');

%===============================================================================
function test_failedAssertionWithCustomMessage

A = [1 1e6];
B = [2 1e6];
f = @() assertVectorsAlmostEqual(A, B, 'my message');
assertExceptionThrown(f, 'assertVectorsAlmostEqual:tolExceeded');

%===============================================================================
function test_nonFloatInputs()
assertExceptionThrown(@() assertVectorsAlmostEqual('hello', 'world'), ...
    'assertVectorsAlmostEqual:notFloat');

%===============================================================================
function test_sizeMismatch()
assertExceptionThrown(@() assertVectorsAlmostEqual(1, [1 2]), ...
    'assertVectorsAlmostEqual:sizeMismatch');

%===============================================================================
function test_finiteAndInfinite()
assertExceptionThrown(@() assertVectorsAlmostEqual([1 2], [1 Inf]), ...
    'assertVectorsAlmostEqual:tolExceeded');

%===============================================================================
function test_infiniteAndInfinite
assertExceptionThrown(@() assertVectorsAlmostEqual([1 Inf], [1 Inf]), ...
    'assertVectorsAlmostEqual:tolExceeded');

%===============================================================================
function test_finiteAndNaN
assertExceptionThrown(@() assertVectorsAlmostEqual([1 2], [1 NaN]), ...
    'assertVectorsAlmostEqual:tolExceeded');

%===============================================================================
function test_NanAndNan
assertExceptionThrown(@() assertVectorsAlmostEqual([1 NaN], [1 NaN]), ...
    'assertVectorsAlmostEqual:tolExceeded');

%===============================================================================
function test_oppositeSignInfs
assertExceptionThrown(@() assertVectorsAlmostEqual([1 Inf], [1 -Inf]), ...
    'assertVectorsAlmostEqual:tolExceeded');

author: ricedsp


text:
function tapas_BLT_runITI(params, iTrial)
% tapas_BLT_runITI shows the Fixation Cross
%
% About:
%   Coded by F.Petzschner 19. April 2017
%	Last change: F.Petzschner 19. April 2017

% Show Cue
Screen('DrawTexture', params.screen.window, params.text.imgiti, [], params.screen.fit, 0);
Screen('Flip', params.screen.window);
tapas_BLT_logEvent('iti_on',params);

% Open CO2 valve if specified (and using valves)
if params.valves == 1 && (sum(ismember(params.co2, iTrial)) == 1)
    outp(params.port_address, 2);
    params.CO2valve_open = GetSecs;
end

% Wait
params = wait2_escapeOption(params.dur.iti(iTrial), params);
tapas_BLT_logEvent('iti_off',params);

% Ensure valves are closed
if params.valves == 1
    outp(params.port_address, 0);
end

end

author: translationalneuromodeling

repo: tapas


text:
%% Copyright (C) 2014-2016, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod  @@sym limit (@var{expr}, @var{x}, @var{a}, @var{dir})
%% @defmethodx @@sym limit (@var{expr}, @var{x}, @var{a})
%% @defmethodx @@sym limit (@var{expr}, @var{a})
%% @defmethodx @@sym limit (@var{expr})
%% Evaluate symbolic limits.
%%
%% The limit of @var{expr} as @var{x} tends to @var{a} from
%% @var{dir}.  @var{dir} can be @code{left} or @code{right}.
%%
%% Examples:
%% @example
%% @group
%% syms x
%% L = limit(sin(x)/x, x, 0)
%%   @result{} L = (sym) 1
%% L = limit(1/x, x, sym(inf))
%%   @result{} L = (sym) 0
%% L = limit(1/x, x, 0, 'left')
%%   @result{} L = (sym) -∞
%% L = limit(1/x, x, 0, 'right')
%%   @result{} L = (sym) ∞
%% @end group
%% @end example
%%
%% If @var{x} is omitted, @code{symvar} is used to determine the
%% variable.  If @var{a} is omitted, it defaults to 0.
%%
%% @var{dir} defaults to @code{right}.  Note this is different from
%% Matlab's Symbolic Math Toolbox which returns @code{NaN} for
%% @code{limit(1/x, x, 0)}
%% (and @code{+/-inf} if you specify @code{left/right}).  I'm not
%% sure how to get this nicer behaviour from SymPy.
%% FIXME: this is https://github.com/cbm755/octsympy/issues/74
%%
%% @seealso{@@sym/diff}
%% @end defmethod

function L = limit(f, x, a, dir)

  if (nargin > 4 || nargin < 1)
    print_usage ();
  end

  f = sym(f);
  if (nargin < 4)
    dir= 'right';
  end
  if (nargin == 2)
    a = x;
    x = symvar(f, 1);
  end
  if (nargin == 1)
    x = symvar(f, 1);
    a = 0;
  end

  switch (lower (dir))
    case {'left' '-'}
      pdir = '-';
    case {'right' '+'}
      pdir = '+';
    otherwise
      print_usage ();
  end

  if (isempty (x))
    L = f;
    return
  end

  L = elementwise_op ('lambda f, x, a, dir: f.limit(x, a, dir=dir)', ...
                      sym(f), sym(x), sym(a), pdir);
end

%!error limit (sym(1), 2, 3, 4, 5)

%!shared x, oo
%! syms x
%! oo = sym(inf);

%!assert (isa (limit(x, x, pi), 'sym'))

%!assert (isequal (limit(x, x, pi), sym(pi)))

%!assert (isequal (limit(sin(x)/x, x, 0), 1))

%!test
%! % left/right-hand limit
%! assert (isequal (limit(1/x, x, 0, 'right'), oo))
%! assert (isequal (limit(1/x, x, 0), oo))
%! assert (isequal (limit(1/x, x, 0, 'left'), -oo))
%! assert (isequal (limit(1/x, x, oo), 0))
%! assert (isequal (limit(sign(x), x, 0, 'left'), -1))
%! assert (isequal (limit(sign(x), x, 0, 'right'), 1))
%! assert (isequal (limit(sign(x), x, 0, '-'), -1))
%! assert (isequal (limit(sign(x), x, 0, '+'), 1))

%!test
%! % matrix
%! syms y
%! A = [x 1/x x*y];
%! B = sym([3 sym(1)/3 3*y]);
%! assert (isequal (limit(A, x, 3), B))

%!test
%! % omitting arguments
%! syms a
%! assert (isequal (limit(a), 0))
%! assert (isequal (limit(a*x+a+2), a+2))
%! assert (isequal (limit(a*x+a+2, 6), 7*a+2))

%!test
%! % constants
%! assert (isequal (limit(sym(6)), 6))
%! assert (isequal (limit(sym(6), 7), 6))
%! assert (isequal (limit([sym(6) sym(2)], 7), [6 2]))

%!test
%! % double constant, with sym limit
%! a = limit (6, sym(0));
%! assert (isa (a, 'sym'))
%! assert (isequal (a, sym(6)))

repo: octsympy


text:
function [dm] = m2dm(m)
% Convert length from meters to decimeters.
% Chad A. Greene 2012
dm = m*10;

author: Sable

repo: mcbench-benchmarks


text:
function [error_train, error_val] = ...
    learningCurve(X, y, Xval, yval, lambda)
%LEARNINGCURVE Generates the train and cross validation set errors needed
%to plot a learning curve
%   [error_train, error_val] = ...
%       LEARNINGCURVE(X, y, Xval, yval, lambda) returns the train and
%       cross validation set errors for a learning curve. In particular,
%       it returns two vectors of the same length - error_train and
%       error_val. Then, error_train(i) contains the training error for
%       i examples (and similarly for error_val(i)).
%
%   In this function, you will compute the train and test errors for
%   dataset sizes from 1 up to m. In practice, when working with larger
%   datasets, you might want to do this in larger intervals.
%

% Number of training examples
m = size(X, 1);

% You need to return these values correctly
error_train = zeros(m, 1);
error_val   = zeros(m, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return training errors in
%               error_train and the cross validation errors in error_val.
%               The vector numex_vec contains the number of training
%               examples to use for each calculation of training error and
%               cross validation error, i.e, error_train(i) and
%               error_val(i) should give you the errors
%               obtained after training on i examples.
%
% Note: You should evaluate the training error on the first i training
%       examples (i.e., X(1:i, :) and y(1:i)).
%
%       For the cross-validation error, you should instead evaluate on
%       the _entire_ cross validation set (Xval and yval).
%
% Note: If you are using your cost function (linearRegCostFunction)
%       to compute the training and cross validation error, you should
%       call the function with the lambda argument set to 0.
%       Do note that you will still need to use lambda when running
%       the training to obtain the theta parameters.
%
% Hint: You can loop over the examples with the following:
%
%       for i = 1:m
%           % Compute train/cross validation errors using training examples
%           % X(1:i, :) and y(1:i), storing the result in
%           % error_train(i) and error_val(i)
%           ....
%
%       end
%

% ---------------------- Sample Solution ----------------------

% linearRegCostFunction(X, y, theta, lambda)

for i = 1:m,
  X_train = X(1:i, :);
  y_train = y(1:i);
  theta = trainLinearReg(X_train, y_train, lambda);
  error_train(i)  = linearRegCostFunction(X_train, y_train, theta, 0);
  error_val(i)    = linearRegCostFunction(Xval, yval, theta, 0);
end

% -------------------------------------------------------------

% =========================================================================

end

author: schneems

repo: Octave

save_path: github-repos/MATLAB/schneems-Octave


text:
clear; clc;
% P2.4
x = [1,-2,4,6,-5,8,10];
n = [-4:2];

% P2.4a
[x11,n11] = sigshift(x,n,-2);
[x12,n12] = sigshift(x,n,4);
[x1,n1] = sigadd(3*x11,n11,x12,n12);
[x1,n1] = sigadd(x1,n1,-2*x,n);
[xe,xo,m] = evenodd(x1,n1);
figure(1);
subplot(2,2,1); stem(n1,x1,'ko'); title('Original Sequence')
xlabel('n'); ylabel('x(n)');
subplot(2,2,2); stem(m,xe,'ko'); title('Even Part')
xlabel('n'); ylabel('xe(n)');
subplot(2,2,4); stem(m,xo,'ko'); title('Odd Part')
xlabel('n'); ylabel('xo(n)');

% P2.4b
[x21,n21] = sigshift(x,n,-5);
[x22,n22] = sigshift(x,n,-4);
[x2,n2] = sigadd(5*x21,n21,4*x22,n22);
[x2,n2] = sigadd(x2,n2,3*x,n);
[xe,xo,m] = evenodd(x2,n2);
figure(2);
subplot(2,2,1); stem(n2,x2,'ko'); title('Original Sequence')
xlabel('n'); ylabel('x(n)');
subplot(2,2,2); stem(m,xe,'ko'); title('Even Part')
xlabel('n'); ylabel('xe(n)');
subplot(2,2,4); stem(m,xo,'ko'); title('Odd Part')
xlabel('n'); ylabel('xo(n)');

% P2.4c
[x31,n31] = sigshift(x,n,-4);
[x32,n32] = sigshift(x,n,1);
[x33,n33] = sigmult(x31,n31,x32,n32);
[x31,n31] = sigfold(x,n);    [x31,n31] = sigshift(x31,n31,2);
[x32,n32] = sigmult(x31,n31,x,n);
[x3,n3] = sigadd(x33,n33,x32,n32);
[xe,xo,m] = evenodd(x3,n3);
figure(3);
subplot(2,2,1); stem(n3,x3,'ko'); title('Original Sequence')
xlabel('n'); ylabel('x(n)');
subplot(2,2,2); stem(m,xe,'ko'); title('Even Part')
xlabel('n'); ylabel('xe(n)');
subplot(2,2,4); stem(m,xo,'ko'); title('Odd Part')
xlabel('n'); ylabel('xo(n)');

% P2.4d
[x41,n41] = sigshift(x,n,-2);
x42 = cos(0.1*pi*n);
[x42,n42] = sigmult(x42,n,x41,n41);
x43 = 2*exp(0.5*n);
[x41,n41] = sigmult(x43,n,x,n);
[x4,n4] = sigadd(x42,n42,x41,n41);
[xe,xo,m] = evenodd(x4,n4);
figure(4);
subplot(2,2,1); stem(n4,x4,'ko'); title('Original Sequence')
xlabel('n'); ylabel('x(n)');
subplot(2,2,2); stem(m,xe,'ko'); title('Even Part')
xlabel('n'); ylabel('xe(n)');
subplot(2,2,4); stem(m,xo,'ko'); title('Odd Part')
xlabel('n'); ylabel('xo(n)');

% P2.4e
x5 = zeros(size(n));
for k = 1:5,
    [x51,n51] = sigshift(x,n,k);
    x52 = n.*x51;
    x5 = sigadd(x5,n51,x52,n51);
end
[xe,xo,m] = evenodd(x5,n51);
figure(5);
subplot(2,2,1); stem(n51,x5,'ko'); title('Original Sequence')
xlabel('n'); ylabel('x(n)');
subplot(2,2,2); stem(m,xe,'ko'); title('Even Part')
xlabel('n'); ylabel('xe(n)');
subplot(2,2,4); stem(m,xo,'ko'); title('Odd Part')
xlabel('n'); ylabel('xo(n)');


text:
function [ar,ASAsellog,ASAcontrol] = sig2ar(sig,cand_order,last)
%SIG2AR AR model identification
%   [AR,SELLOG] = SIG2AR(SIG) estimates autoregressive models from the
%   data vector SIG and selects a model with optimal predictive
%   qualities. The selected model is returned in the parameter vector AR.
%   The structure SELLOG provides additional information on the selection
%   process.
%
%   SIG2AR(SIG,CAND_ORDER) selects only from candidate models whose
%   orders are entered in CAND_ORDER. CAND_ORDER must either be a row of
%   ascending orders, or a single order (in which case no true order
%   selection is performed).
%
%   Without user intervention, the mean of SIG is subtracted from the
%   data. To control the subtraction of the mean, see the help topics on
%   ASAGLOB_SUBTR_MEAN and ASAGLOB_MEAN_ADJ.
%
%   SIG2AR is an ARMASA main function.
%
%   See also: SIG2MA, SIG2ARMA, ARMASEL.

%   References: P. M. T. Broersen, Facts and Fiction in Spectral
%               Analysis, IEEE Transactions on Instrumentation and
%               Measurement, Vol. 49, No. 4, August 2000, pp. 766-772.

%Header
%=========================================================================

%Declaration of variables
%------------------------

%Declare and assign values to local variables
%according to the input argument pattern
switch nargin
case 1
   if isa(sig,'struct'), ASAcontrol=sig; sig=[];
   else, ASAcontrol=[];
   end
   cand_order=[];
case 2
   if isa(cand_order,'struct'), ASAcontrol=cand_order; cand_order=[];
   else, ASAcontrol=[];
   end
case 3
   if isa(last,'struct'), ASAcontrol=last;
   else, error(ASAerr(39))
   end
otherwise
   error(ASAerr(1,mfilename))
end

if isequal(nargin,1) & ~isempty(ASAcontrol)
      %ASAcontrol is the only input argument
   ASAcontrol.error_chk = 0;
   ASAcontrol.run = 0;
end

%Declare ASAglob variables
ASAglob = {'ASAglob_subtr_mean';'ASAglob_mean_adj'; ...
      'ASAglob_rc';'ASAglob_ar'};

%Assign values to ASAglob variables by screening the
%caller workspace
for ASAcounter = 1:length(ASAglob)
   ASAvar = ASAglob{ASAcounter};
   eval(['global ' ASAvar]);
   if evalin('caller',['exist(''' ASAvar ''',''var'')'])
      eval([ASAvar '=evalin(''caller'',ASAvar);']);
   else
      eval([ASAvar '=[];']);
   end
end

%ARMASA-function version information
%-----------------------------------

%This ARMASA-function is characterized by
%its current version,
ASAcontrol.is_version = [2000 12 30 20 0 0];
%and its compatability with versions down to,
ASAcontrol.comp_version = [2000 12 30 20 0 0];

%This function calls other functions of the ARMASA
%toolbox. The versions of these other functions must
%be greater than or equal to:
ASAcontrol.req_version.burg = [2000 12 30 20 0 0];
ASAcontrol.req_version.cic = [2000 12 30 20 0 0];
ASAcontrol.req_version.rc2arset = [2000 12 30 20 0 0];

%Checks
%------

if ~any(strcmp(fieldnames(ASAcontrol),'error_chk')) | ASAcontrol.error_chk
      %Perform standard error checks
   %Input argument format checks
   ASAcontrol.error_chk = 1;
   if ~isnum(sig)
      error(ASAerr(11,'sig'))
   elseif ~isavector(sig)
      error([ASAerr(14) ASAerr(15,'sig')])
   elseif size(sig,2)>1
      sig = sig(:);
      warning(ASAwarn(25,{'row';'sig';'column'},ASAcontrol))
   end
   if ~isempty(cand_order)
      if ~isnum(cand_order) | ~isintvector(cand_order) |...
            cand_order(1)<0 | ~isascending(cand_order)
         error(ASAerr(12,{'candidate';'cand_order'}))
      elseif size(cand_order,1)>1
         cand_order = cand_order';
         warning(ASAwarn(25,{'column';'cand_order';'row'},ASAcontrol))
      end
   end

   %Input argument value checks
   if ~isreal(sig)
      error(ASAerr(13))
   end
   if max(cand_order) > length(sig)-1
      error(ASAerr(21))
   end
end

if ~any(strcmp(fieldnames(ASAcontrol),'version_chk')) | ASAcontrol.version_chk
      %Perform version check
   ASAcontrol.version_chk = 1;

   %Make sure the requested version of this function
   %complies with its actual version
   ASAversionchk(ASAcontrol);

   %Make sure the requested versions of the called
   %functions comply with their actual versions
   burg(ASAcontrol);
   cic(ASAcontrol);
   rc2arset(ASAcontrol);
end

if ~any(strcmp(fieldnames(ASAcontrol),'run')) | ASAcontrol.run
      %Run the computational kernel
   ASAcontrol.run = 1;
   ASAcontrol.version_chk = 0;
   ASAcontrol.error_chk = 0;
   ASAtime = clock;
   ASAdate = now;

%Main
%=====================================================

%Initialization of variables
%---------------------------

if isempty(ASAglob_subtr_mean) | ASAglob_subtr_mean
   sig = sig-mean(sig);
   if isempty(ASAglob_mean_adj)
      ASAglob_mean_adj = 1;
   end
elseif isempty(ASAglob_mean_adj)
   ASAglob_mean_adj = 0;
end

n_obs = size(sig,1);

%Determination of the maximum candidate AR order
%-----------------------------------------------

if ~isempty(cand_order)
   max_order = cand_order(end);
else
   max_order = ...
      min(fix(n_obs/2),fix(200*log10(n_obs)));
   if max_order > 1000;
      max_order = 1000;
   end
end

%Estimation procedure
%--------------------

[rc,var] = burg(sig,max_order,ASAcontrol);

%AR model order selection
%------------------------

rc(1) = 0;
res = var*cumprod(1-rc.^2);
rc(1) = 1;
[cicar,pe_est] = cic(res,n_obs,cand_order,ASAcontrol);
[min_value,sel_location] = min(cicar);
if isempty(cand_order)
   cand_order = 0:max_order;
end
sel_order = cand_order(sel_location);

%Arranging output arguments
%--------------------------

ar = rc2arset(rc(1:sel_order+1),ASAcontrol);

ASAglob_rc = rc;
ASAglob_ar = ar;

if nargout>1
   ASAsellog.funct_name = mfilename;
   ASAsellog.funct_version = ASAcontrol.is_version;
   ASAsellog.date_time = ...
      [datestr(ASAdate,8) 32 datestr(ASAdate,0)];
   ASAsellog.comp_time = etime(clock,ASAtime);
   ASAsellog.ar = ar;
   ASAsellog.rcarlong = rc;
   ASAsellog.mean_adj = ASAglob_mean_adj;
   ASAsellog.cand_order = cand_order;
   ASAsellog.cic = cicar;
   ASAsellog.pe_est = pe_est;
end

%Footer
%=====================================================

else %Skip the computational kernel
   %Return ASAcontrol as the first output argument
   if nargout>1
      warning(ASAwarn(9,mfilename,ASAcontrol))
   end
   ar = ASAcontrol;
   ASAcontrol = [];
end

%Program history
%======================================================================
%
% Version                Programmer(s)          E-mail address
% -------                -------------          --------------
% former versions        P.M.T. Broersen        p.m.t.broersen@tudelft.nl
% [2000 12 30 20 0 0]    W. Wunderink           wwunderink01@freeler.nl


text:
%%% SetupBipedRobot.m
%%% Set biped robot structure of Figure 2.19, 2.20
%%% Field definition: Table 2.1 Link Parameters

global uLINK

ToDeg = 180/pi;
ToRad = pi/180;
UX = [1 0 0]';
UY = [0 1 0]';
UZ = [0 0 1]';

uLINK    = struct('name','BODY'    , 'm', 10, 'sister', 0, 'child', 2, 'b',[0  0    0.7]','a',UZ,'q',0);

uLINK(2) = struct('name','RLEG_J0' , 'm',  5, 'sister', 8, 'child', 3, 'b',[0 -0.1 0]'   ,'a',UZ,'q',0);
uLINK(3) = struct('name','RLEG_J1' , 'm',  1, 'sister', 0, 'child', 4, 'b',[0  0   0]'   ,'a',UX,'q',0);
uLINK(4) = struct('name','RLEG_J2' , 'm',  5, 'sister', 0, 'child', 5, 'b',[0  0   0]'   ,'a',UY,'q',0);
uLINK(5) = struct('name','RLEG_J3' , 'm',  1, 'sister', 0, 'child', 6, 'b',[0  0  -0.3]' ,'a',UY,'q',0);
uLINK(6) = struct('name','RLEG_J4' , 'm',  6, 'sister', 0, 'child', 7, 'b',[0  0  -0.3]' ,'a',UY,'q',0);
uLINK(7) = struct('name','RLEG_J5' , 'm',  2, 'sister', 0, 'child', 0, 'b',[0  0   0  ]' ,'a',UX,'q',0);

uLINK(8) = struct('name','LLEG_J0' , 'm',  5, 'sister', 0, 'child', 9, 'b',[0  0.1 0]'   ,'a',UZ,'q',0);
uLINK(9) = struct('name','LLEG_J1' , 'm',  1, 'sister', 0, 'child',10, 'b',[0  0   0]'   ,'a',UX,'q',0);
uLINK(10)= struct('name','LLEG_J2' , 'm',  5, 'sister', 0, 'child',11, 'b',[0  0   0]'   ,'a',UY,'q',0);
uLINK(11)= struct('name','LLEG_J3' , 'm',  1, 'sister', 0, 'child',12, 'b',[0  0  -0.3]' ,'a',UY,'q',0);
uLINK(12)= struct('name','LLEG_J4' , 'm',  6, 'sister', 0, 'child',13, 'b',[0  0  -0.3]' ,'a',UY,'q',0);
uLINK(13)= struct('name','LLEG_J5' , 'm',  2, 'sister', 0, 'child', 0, 'b',[0  0   0  ]' ,'a',UX,'q',0);

[uLINK(1).vertex,uLINK(1).face]   = MakeBox([0.1 0.3 0.5]  ,[0.05 0.15 -0.05] );    % BODY
[uLINK(7).vertex,uLINK(7).face]   = MakeBox([0.2 0.1 0.02] ,[0.05  0.05 0.05]);     % Foot
[uLINK(13).vertex,uLINK(13).face] = MakeBox([0.2 0.1 0.02] ,[0.05  0.05 0.05]);     % Foot

FindMother(1);   % Find mother link from sister and child data

%%% Substitute the ID to the link name variables. For example, BODY=1.
for n=1:length(uLINK)
    eval([uLINK(n).name,'=',num2str(n),';']);
end

uLINK(BODY).p = [0.0, 0.0, 0.65]';
uLINK(BODY).R = eye(3);
ForwardKinematics(1);

uLINK(BODY).v = [0 0 0]';
uLINK(BODY).w = [0 0 0]';
for n=1:length(uLINK)
    uLINK(n).dq     = 0;            % joitn speed   [rad/s]
end

author: s-kajita

repo: IntroductionToHumanoidRobotics


text:
function h = genHazFun(t,b)

% Sample general hazard function; depends on three parameters

h = b(1) + b(2)*(t./b(3)).*exp(-(t./b(3)));


text:
 % Copyright 2001, Brown University, Providence, Rhode Island.
 %
 % All Rights Reserved
 %
 % Permission to use this software for noncommercial research and
 % educational purposes is hereby granted without fee.
 % Redistribution, sale, or incorporation of this software into a
 % commercial product is prohibited.
 %
 % BROWN UNIVERSITY DISCLAIMS ANY AND ALL WARRANTIES WITH REGARD TO
 % THIS SOFTWARE,INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY
 % AND FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL BROWN
 % UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
 % DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE,
 % DATA OR PROFITS.

% ----------------------------------------------------------------
%  jacobd() - derivative of jacobi polynomials - vector version
%
%  Get a vector 'poly' of values of the derivative of the n_th order
%  Jacobi polynomial P^(alpha,beta)_N(z) at the np points z.
%
%  To do this we have used the relation
%
%  d   alpha,beta   1                  alpha+1,beta+1
%  -- P (z)       = -(alpha+beta+n+1) P  (z)
%  dz  n            2                  n-1
%  ----------------------------------------------------------------*/

function jd = jacobd(z, n, alpha, beta)

 jd = zeros(size(z));

  one = 1.0;
  if(n == 0)
    jd(:,:) = 0.0;
  else
    jd = am282jacobi1d(z,n-1,alpha+one,beta+one);
    jd = jd*0.5*(alpha + beta + n + one);
  end

repo: Aero-matlab


text:
function [ y, ierror ] = y_check_common ( y )

%*****************************************************************************80
%
%% Y_CHECK_COMMON checks a Common year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, the year, which must not be 0.
%
%    Output, integer IERROR, is 0 if Y is legal, and 1 otherwise.
%
  if ( y ~= 0 )
    ierror = 0;
  else
    ierror = 1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [gt,dt] = evalRes( gt0, dt0, thr, mul )
% Evaluates detections against ground truth data.
%
% Uses modified Pascal criteria that allows for "ignore" regions. The
% Pascal criteria states that a ground truth bounding box (gtBb) and a
% detected bounding box (dtBb) match if their overlap area (oa):
%  oa(gtBb,dtBb) = area(intersect(gtBb,dtBb)) / area(union(gtBb,dtBb))
% is over a sufficient threshold (typically .5). In the modified criteria,
% the dtBb can match any subregion of a gtBb set to "ignore". Choosing
% gtBb' in gtBb that most closely matches dtBb can be done by using
% gtBb'=intersect(dtBb,gtBb). Computing oa(gtBb',dtBb) is equivalent to
%  oa'(gtBb,dtBb) = area(intersect(gtBb,dtBb)) / area(dtBb)
% For gtBb set to ignore the above formula for oa is used.
%
% Highest scoring detections are matched first. Matches to standard,
% (non-ignore) gtBb are preferred. Each dtBb and gtBb may be matched at
% most once, except for ignore-gtBb which can be matched multiple times.
% Unmatched dtBb are false-positives, unmatched gtBb are false-negatives.
% Each match between a dtBb and gtBb is a true-positive, except matches
% between dtBb and ignore-gtBb which do not affect the evaluation criteria.
%
% In addition to taking gt/dt results on a single image, evalRes() can take
% cell arrays of gt/dt bbs, in which case evaluation proceeds on each
% element. Use bbGt>loadAll() to load gt/dt for multiple images.
%
% Each gt/dt output row has a flag match that is either -1/0/1:
%  for gt: -1=ignore,  0=fn [unmatched],  1=tp [matched]
%  for dt: -1=ignore,  0=fp [unmatched],  1=tp [matched]
%
% USAGE
%  [gt, dt] = bbGt( 'evalRes', gt0, dt0, [thr], [mul] )
%
% INPUTS
%  gt0  - [mx5] ground truth array with rows [x y w h ignore]
%  dt0  - [nx5] detection results array with rows [x y w h score]
%  thr  - [.5] the threshold on oa for comparing two bbs
%  mul  - [0] if true allow multiple matches to each gt
%
% OUTPUTS
%  gt   - [mx5] ground truth results [x y w h match]
%  dt   - [nx6] detection results [x y w h score match]
%
% EXAMPLE
%
% See also bbGt, bbGt>compOas, bbGt>loadAll

% get parameters
if(nargin<3 || isempty(thr)), thr=.7; end
if(nargin<4 || isempty(mul)), mul=0; end

% if gt0 and dt0 are cell arrays run on each element in turn
if( iscell(gt0) && iscell(dt0) ), n=length(gt0);
  assert(length(dt0)==n); gt=cell(1,n); dt=gt;
  for i=1:n, [gt{i},dt{i}] = evalRes(gt0{i},dt0{i},thr,mul); end; return;
end

% check inputs
if(isempty(gt0)), gt0=zeros(0,5); end
if(isempty(dt0)), dt0=zeros(0,5); end
assert( size(dt0,2)==5 ); nd=size(dt0,1);
assert( size(gt0,2)==5 ); ng=size(gt0,1);

% sort dt highest score first, sort gt ignore last
[~,ord]=sort(dt0(:,5),'descend'); dt0=dt0(ord,:);
[~,ord]=sort(gt0(:,5),'ascend'); gt0=gt0(ord,:);
gt=gt0;  dt=dt0; dt=[dt zeros(nd,1)];
gt(:,5)=-gt(:,5);
% Attempt to match each (sorted) dt to each (sorted) gt
oa = compOas( dt(:,1:4), gt(:,1:4), gt(:,5)==-1 );
for d = 1:nd
  bstOa = thr;
  bstg = 0;
  bstm = 0; % info about best match so far
  for g = 1:ng
    % if this gt already matched, continue to next gt
    m = gt(g,5);
    if(m==1 && ~mul)
        continue;
    end
    % if dt already matched, and on ignore gt, nothing more to do
    if(bstm~=0 && m==-1)
        break;
    end
    % compute overlap area, continue to next gt unless better match made
    if(oa(d,g)<bstOa)
        continue;
    end
    % match successful and best so far, store appropriately
    bstOa = oa(d,g);
    bstg = g;
    if(m==0)
        bstm = 1;
    else
        bstm = -1;
    end
  end;
  g = bstg;
  m = bstm;
  % store type of match for both dt and gt
  if(m==-1)
      dt(d,6) = m;
  elseif(m==1)
      gt(g,5) = m;
      dt(d,6) = m;
  end
end

end

author: VisDrone


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

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

%Test the examples
test1_str=['wfdb2mat(''mitdb/200'');' ...
           '[tm,signal,Fs,siginfo]=rdmat(''200m'');'...
           '[signal2,Fs2,tm2]=rdsamp(''200m'');' ...
           'if(sum(abs(signal-signal2)) ~=0);error(''data reading failed'');end' ...
           ];

test_string={test1_str};
clean_up={['delete([pwd filesep ''mitdb'' filesep ''*'']);' ...
          'if(exist([pwd filesep ''mitdb''],''dir''));rmdir([pwd filesep ''mitdb''],''s'');end;' ...
          'delete([pwd filesep ''200m*'']);' ...
          ]};
[tests,pass,perf]=test_wrapper(test_string,clean_up,verbose);

author: ikarosilva

repo: wfdb-app-toolbox


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   PARAMETERS Returns a data structure containing the parameters of the
%   ABB IRB840A.
%
%   Author: Arturo Gil. Universidad Miguel Hernandez de Elche.
%   email: arturo.gil@umh.es date:   07/11/2020
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function robot = parameters()

robot.name= 'ABB_IRB840_A';

robot.DH.theta= '[pi/2 pi/2 0 q(4)]';
robot.DH.d='[q(1) q(2) q(3) 0]';
robot.DH.a='[2.0 0 0 0]';
robot.DH.alpha= '[pi/2 -pi/2 0 0]';
robot.J=[];

robot.inversekinematic_fn = 'inversekinematic_irb840(robot, T)';
robot.directkinematic_fn = 'directkinematic(robot, q)';

%number of degrees of freedom
robot.DOF = 4;

%rotational: 0, translational: 1
robot.kind=['T' 'T' 'T' 'R'];

%minimum and maximum rotation angle in rad
robot.maxangle =[-1 1; %Axis 1, minimum, maximum
                -1 1; %Axis 2, minimum, maximum
                -1 1; %Axis 3
                deg2rad(-400) deg2rad(400)]; %Axis 4: Unlimited (400� default)

%maximum absolute speed of each joint rad/s or m/s
robot.velmax = [1; %Axis 1, m/s
                1; %Axis 2, m/s
                1; %Axis 3, m/s
                deg2rad(360)]; %Axis 4, m/s

robot.accelmax=robot.velmax/0.1; % 0.1 is here an acceleration time

% end effectors maximum velocity
robot.linear_velmax = 2.5; %m/s

%base reference system
robot.T0 = eye(4);

robot.T0 = [0 0 1 0;
            1 0 0 0;
            0 1 0 0;
            0 0 0 1];

%INITIALIZATION OF VARIABLES REQUIRED FOR THE SIMULATION
%position, velocity and acceleration
robot=init_sim_variables(robot);
robot.path = pwd;

% GRAPHICS
robot.graphical.has_graphics=0;
robot.graphical.color = [255 102 51]./255;
%for transparency
robot.graphical.draw_transparent=0;
%draw DH systems
robot.graphical.draw_axes=1;
%DH system length and Font size, standard is 1/10. Select 2/20, 3/30 for
%bigger robots
robot.graphical.axes_scale=1;
%adjust for a default view of the robot
robot.axis=[-0.5 0.75 -0.75 0.75 0 1.1];
%read graphics files
robot = read_graphics(robot);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DYNAMIC PARAMETERS
%   WARNING! These parameters do not correspond to the actual IRB 140
%   robot. They have been introduced to demonstrate the necessity of
%   simulating the robot and should be used only for educational purposes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
robot.has_dynamics=0;

%consider friction in the computations
robot.dynamics.friction=0;

%link masses (kg)
robot.dynamics.masses=[25 27 15 10 2.5 1.5];

%COM of each link with respect to own reference system
robot.dynamics.r_com=[0       0          0; %(rx, ry, rz) link 1
                     -0.05	 0.006	 0.1; %(rx, ry, rz) link 2
                    -0.0203	-0.0141	 0.070;  %(rx, ry, rz) link 3
                     0       0.019       0;%(rx, ry, rz) link 4
                     0       0           0;%(rx, ry, rz) link 5
                     0       0         0.032];%(rx, ry, rz) link 6

%Inertia matrices of each link with respect to its D-H reference system.
% Ixx	Iyy	Izz	Ixy	Iyz	Ixz, for each row
robot.dynamics.Inertia=[0      0.35	0   	0	0	0;
    .13     .524	.539	0	0	0;
    .066	.086	.0125	0	0	0;
    1.8e-3	1.3e-3	1.8e-3	0	0	0;
    .3e-3	.4e-3	.3e-3	0	0	0;
    .15e-3	.15e-3	.04e-3	0	0	0];

robot.motors=load_motors([5 5 5 4 4 4]);
%Speed reductor at each joint
robot.motors.G=[300 300 300 300 300 300];

repo: ARTE


text:
function value = halton_base_check ( dim_num, base )

%*****************************************************************************80
%
%% HALTON_BASE_CHECK is TRUE if BASE is legal.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer BASE(1:DIM_NUM), the Halton bases.
%    Each base should be greater than 1.  Only the integer
%    part of a base is used.
%
%    Output, logical VALUE.
%
  if ( any ( base(1:dim_num) <= 1 ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HALTON_BASE_CHECK - Fatal error!\n' );
    fprintf ( 1, '  At least one of the input bases is <= 1!\n' );
    value = 0;
  else
    value = 1;
  end

  return
end


text:
function R = getHPfilter(freq,numsamps)

% freq is number of frequencies to cut off.
base = 1:2*pi/numsamps:2*pi

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef IncreaseLast_LineSearchInitiator < LineSearchInitiator

    properties (Access = private)
        incrementFactor
        maxValue
        minValue
    end

    methods (Access = public)

        function obj = IncreaseLast_LineSearchInitiator(cParams)
            obj.init(cParams);
            obj.incrementFactor = cParams.incrementFactor;
            obj.maxValue        = cParams.maxValue;
            obj.minValue        = cParams.minValue;
        end

        function initStep = compute(obj,lastStep)
            if lastStep <= obj.minValue
                initStep = 5*lastStep;
            else
                k = obj.incrementFactor;
                initStep = k*lastStep;
                initStep = min(initStep,obj.maxValue);
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function n = ndims(t)
%NDIMS Number of dimensions of a sparse tensor.
%
%   NDIMS(T) returns the number of dimensions of sparse tensor T.
%
%   Examples:
%   T = sptenrand([3 2 2],5);
%   ndims(T) %<-- should return 3
%
%   See also SPTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

n = size(t.size,2);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
clear all; close all;
I=imread('pout.tif');
M=stretchlim(I);
J=imadjust(I, M, [ ]);
figure;
subplot(121);
imshow(uint8(I));
subplot(122);
imshow(uint8(J));

author: HuangCongQing


text:
function item = cfg_intree(varargin)

% This is currently only a "marker" class that should be inherited by all
% within-tree classes. It does not add fields and does not have methods.
%
% 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_intree.m 1716 2008-05-23 08:18:45Z volkmar $

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

item = class(struct('unused',[]), mfilename);


text:
function [Wrot, CC] = get_whitening_matrix_faster(rez)
% based on a subset of the data, compute a channel whitening matrix
% this requires temporal filtering first (gpufilter)

ops = rez.ops;
Nbatch = ops.Nbatch;
twind = ops.twind;
NchanTOT = ops.NchanTOT;
NT = ops.NT;
NTbuff = ops.NTbuff;
chanMap = ops.chanMap;
Nchan = rez.ops.Nchan;
xc = rez.xc;
yc = rez.yc;

% load data into patches, filter, compute covariance
fprintf('Getting channel whitening matrix... \n');

ibatch = 1;
ibatch = ibatch:ops.nSkipCov:Nbatch;

% scrappy progress bar in command window
% allBatches = 1:Nbatch;
pb = progBar(ibatch, 20);

if getOr(ops, 'useMemMapping',0)
    %% Parallelize raw data loading (needs memmapfile implement for parallel read access)
        % slice ops for parallel
        ntbuff = ops.ntbuff; % NOTE: this is just the buffer length, not ".NTbuff" (which == the actual batch + buffer length)
        NT = ops.NT;
        bufferedNT = NT + 2*ntbuff;
        tend = ops.tend;
        twind = ops.twind;
        doCAR = logical(ops.CAR);
        CC = cell(1, numel(ibatch));
        [CC{:}] = deal(zeros(Nchan,  Nchan, 'single')); % we'll estimate the covariance from data batches, then add to this variable

        % set up the parameters of the filter % just one instance of filters
        if isfield(ops,'fslow')&&ops.fslow<ops.fs/2
            [b1, a1] = butter(3, [ops.fshigh/ops.fs,ops.fslow/ops.fs]*2, 'bandpass'); % butterworth filter with only 3 nodes (otherwise it's unstable for float32)
        else
            [b1, a1] = butter(3, ops.fshigh/ops.fs*2, 'high'); % the default is to only do high-pass filtering at 150Hz
        end

        % Memory map the raw data file to m.Data.x of size [nChan, nSamples)
        mmfRaw = dir(ops.fbinary);
        nsamp = mmfRaw.bytes/2/ops.NchanTOT;
        mmfRaw = memmapfile(ops.fbinary, 'Format',{'int16', [ops.NchanTOT, nsamp], 'x'}); % ignore tstart here b/c we'll account for it on each mapped read (using .twind)
        %parmmf = parallel.pool.Constant(mmf);
        tic
        parfor i = 1:length(ibatch)
        % for i = 1:length(ibatch)
            ii = ibatch(i);
            %offset = max(0, twind + 2*NchanTOT*((NT - ntbuff) * (ii-1) - 2*ntbuff));
            % index batch timepoints to read:  [ntbuff + NT + ntbuff]
            twindow = [max(1, twind + NT*(ii-1) - ntbuff): min(tend, twind + NT*(ii) + ntbuff)];

            % read from memory mapped file
            buff = mmfRaw.Data.x(:,twindow);

            % % --- step inside gpufilter.m operations ---%

            % subsample only good channels & transpose for filtering
            buff = single(buff(chanMap,:))';    % buff dims now: [time, channel]

            % --- Demean before padding ---
            % subtract within-channel means from each channel
            buff = buff - mean(buff, 1); % subtract mean of each channel

            % CAR, common average referencing by median
            if doCAR %getOr(ops, 'CAR', 1)
                buff = buff - median(buff, 2); % subtract median across channels
            end

            % Now can pad first & last batches with zeros
            nsampcurr = size(buff,1);
            if nsampcurr<bufferedNT
                if ii~=1
                    % append zeros if end batch
                    buff = [buff; zeros(bufferedNT-nsampcurr, size(buff,2))];
                else
                    % prepend zeros if first batch
                    buff = [zeros(bufferedNT-nsampcurr, size(buff,2)); buff];
                end
            end

            % apply high/low pass filtering
            % next four lines should be equivalent to filtfilt (which cannot be used because it requires float64)
            buff = filter(b1, a1, buff); % causal forward filter
            buff = flipud(buff); % reverse time
            buff = filter(b1, a1, buff); % causal forward filter again
            buff = flipud(buff); % reverse time back

            % % --- end of gpufilter.m operations ---%

            buff = buff(ntbuff+1:end-ntbuff,:);   % remove timepoints used as buffers

            CC{i}        = (buff' * buff)/NT; % sample covariance
            %pb.check(i); % progress bar updates not meaningful when running on parallel pool
        end
        toc
        % scale & combine parallel outputs
        CC = cell2mat(cellfun(@(x) permute(x,[3,1,2]), CC', 'uni',0)); % shift dims of cell so that cell2mat maintains array dimensions
        CC = squeeze(mean(CC));
    %%

else
    % old single-thread method
    fid = fopen(ops.fbinary, 'r');
    tic
    CC = gpuArray.zeros( Nchan,  Nchan, 'single'); % we'll estimate the covariance from data batches, then add to this variable

    for i = 1:length(ibatch) %allBatches
        ii = ibatch(i);
        offset = max(0, twind + 2*NchanTOT*((NT - ops.ntbuff) * (ii-1) - 2*ops.ntbuff));
        fseek(fid, offset, 'bof');
        buff = fread(fid, [NchanTOT NTbuff], '*int16');

        if isempty(buff)
            break;
        end
        nsampcurr = size(buff,2);
        if nsampcurr<NTbuff
            buff(:, nsampcurr+1:NTbuff) = repmat(buff(:,nsampcurr), 1, NTbuff-nsampcurr);
        end

        datr    = gpufilter(buff, ops, chanMap); % apply filters and median subtraction

        CC        = CC + (datr' * datr)/NT; % sample covariance

        pb.check(ii) % update progress bar in command window
    end
    CC = CC / ceil((Nbatch-1)/ops.nSkipCov); % normalize by number of batches
    toc
    fclose(fid);
end

% CC = diag(diag(CC));
% mu = mean(diag(CC));

% CC = gpuArray(eye(size(CC), 'single'));
% CC = CC * mu;

% xp = cat(2, rez.xc, rez.yc);
% CC = mu * (.05 + .95 * kernel2D(xp, xp, 30));

% plot(diag(CC))

if ops.whiteningRange<Inf
    % if there are too many channels, a finite whiteningRange is more robust to noise in the estimation of the covariance
    ops.whiteningRange = min(ops.whiteningRange, Nchan);
    Wrot = whiteningLocal(gather(CC), yc, xc, ops.whiteningRange); % this function performs the same matrix inversions as below, just on subsets of channels around each channel
else
    Wrot = whiteningFromCovariance(CC);
end
Wrot    = ops.scaleproc * Wrot; % scale this from unit variance to int 16 range. The default value of 200 should be fine in most (all?) situations.

fprintf('Channel-whitening matrix computed. \n');

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
% Test file for bndfun/poly.m

function pass = test_poly(pref)

% Get preferences.
if ( nargin < 1 )
    pref = chebfunpref();
end

% Set the domain
dom = [-2 7];

%%
% Check a few simple examples.

f = bndfun(@(x) zeros(size(x)), struct('domain', dom), pref);
p = poly(f);
pass(1) = (norm(p, inf) <= get(f, 'vscale')*eps);

f = bndfun(@(x) 3*ones(size(x)), struct('domain', dom), pref);
p = poly(f);
pass(2) = (norm(p - 3, inf) < get(f, 'vscale')*eps);

f = bndfun(@(x) 6.4*x - 3i, struct('domain', dom), pref);
p = poly(f);
pass(3) = (norm(p - [6.4 (-3i)], inf) < get(f, 'vscale')*eps);

f = bndfun(@(x) 2i*x.^5 - 3.2*x.^4 + 2*x.^2 - (1.2 + 3i), ...
    struct('domain', dom), pref);
p = poly(f);
pass(4) = (norm(p - [2i (-3.2) 0 2 0 -(1.2 + 3i)], inf) ...
    < get(f, 'vscale')*eps);

%%
% Verify operation for array-valued bndfun objects.

f = bndfun(@(x) [3*ones(size(x)), (6.4*x - 3i), (4*x.^2 - 2i*x + 3.7)], ...
    struct('domain', dom), pref);
p = poly(f);
p_exact = [0 0     3;
           0 6.4   (-3i);
           4 (-2i) 3.7];
pass(5) = (norm(p(:) - p_exact(:), inf) < ...
    max(get(f, 'vscale')*eps));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function outConMatStack=legaliseConMatRow(inConMatStack, badRowIndex)
% Makes sure that there are only two entries on any row
% Input is a cell array of connection matrices with a row 'badRow' in
% common
% badRow has more than 2 entries (say, t)
% there are x= (nchoosek([...],2) new rows that can be generated from
% badRow - each one having 2 entries.
% For each conMat in the input stack, generate x output conMats with each
% of the possible new rows.
% Last edited $Date: 2007/07/05 19:50:13 $
defaultConMat=sparse(inConMatStack{1});
badRow=defaultConMat(badRowIndex,:);
badRowVals=find(badRow);
badRowCombinations=nchoosek(badRowVals,2);
nCombinations=length(badRowCombinations);
nConMats=length(inConMatStack);
counter=1;
fprintf('Generating %d conmats',nConMats*nCombinations);
for thisConMat=1:nConMats
	for thisRowComb=1:nCombinations;
		outConMatTemp=sparse(inConMatStack{thisConMat});
		outConMatTemp(badRowIndex,:)=0;
		outConMatTemp(badRowIndex,badRowCombinations(thisRowComb,:))=1;
		outConMatStack{counter}=sparse(outConMatTemp);
		counter=counter+1;
	end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ an, bn ] = lcrg_anbn ( a, b, c, n )

%*****************************************************************************80
%
%% LCRG_ANBN computes the "N-th power" of a linear congruential generator.
%
%  Discussion:
%
%    We are considering a linear congruential random number generator.
%    The LCRG takes as input an integer value called SEED, and returns
%    an updated value of SEED,
%
%      SEED(out) = ( a * SEED(in) + b ) mod c.
%
%    and an associated pseudorandom real value
%
%      U = SEED(out) / c.
%
%    In most cases, a user is content to call the LCRG repeatedly, with
%    the updating of SEED being taken care of automatically.
%
%    The purpose of this routine is to determine the values of AN and BN
%    that describe the LCRG that is equivalent to N applications of the
%    original LCRG.
%
%    One use for such a facility would be to do random number computations
%    in parallel.  If each of N processors is to compute many random values,
%    you can guarantee that they work with distinct random values
%    by starting with a single value of SEED, using the original LCRG to generate
%    the first N-1 "iterates" of SEED, so that you now have N "seed" values,
%    and from now on, applying the N-th power of the LCRG to the seeds.
%
%    If the K-th processor starts from the K-th seed, it will essentially
%    be computing every N-th entry of the original random number sequence,
%    offset by K.  Thus the individual processors will be using a random
%    number stream as good as the original one, and without repeating, and
%    without having to communicate.
%
%    To evaluate the N-th value of SEED directly, we start by ignoring
%    the modular arithmetic, and working out the sequence of calculations
%    as follows:
%
%      SEED(0)   =     SEED.
%      SEED(1)   = a * SEED      + b
%      SEED(2)   = a * SEED(1)   + b = a^2 * SEED           + a * b + b
%      SEED(3)   = a * SEED(2)   + b = a^3 * SEED + a^2 * b + a * b + b
%      ...
%      SEED(N-1) = a * SEED(N-2) + b
%
%      SEED(N) = a * SEED(N-1) + b = a^N * SEED
%                                    + ( a^(n-1) + a^(n-2) + ... + a + 1 ) * b
%
%    or, using the geometric series,
%
%      SEED(N) = a^N * SEED + ( a^N - 1) / ( a - 1 ) * b
%              = AN * SEED + BN
%
%    Thus, from any SEED, we can determine the result of N applications of the
%    original LCRG directly if we can solve
%
%      ( a - 1 ) * BN = ( a^N - 1 ) * b in modular arithmetic,
%
%    and evaluate:
%
%      AN = a^N
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 May 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Barry Wilkinson, Michael Allen,
%    Parallel Programming:
%    Techniques and Applications Using Networked Workstations and Parallel Computers,
%    Prentice Hall,
%    ISBN: 0-13-140563-2,
%    LC: QA76.642.W54.
%
%  Parameters:
%
%    Input, integer A, the multiplier for the LCRG.
%
%    Input, integer B, the added value for the LCRG.
%
%    Input, integer C, the base for the modular arithmetic.
%    For 32 bit arithmetic, this is often 2^31 - 1, or 2147483647.  It is
%    required that 0 < C.
%
%    Input, integer N, the "index", or number of times that the
%    LCRG is to be applied.  It is required that 0 <= N.
%
%    Output, integer AN, BN, the multiplier and added value for
%    the LCRG that represent N applications of the original LCRG.
%
  if ( n < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LCRG_ANBN - Fatal error!\n' );
    fprintf ( 1, '  Illegal input value of N = %d\n', n );
    error ( 'LCRG_ANBN - Fatal error!' );
  end

  if ( c <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LCRG_ANBN - Fatal error!\n' );
    fprintf ( 1, '  Illegal input value of C = %d\n', c );
    error ( 'LCRG_ANBN - Fatal error!' );
  end

  if ( n == 0 )
    an = 1;
    bn = 0;
  elseif ( n == 1 )
    an = a;
    bn = b;
  else
%
%  Compute A^N.
%
    an = power_mod ( a, n, c );
%
%  Solve
%    ( a - 1 ) * BN = ( a^N - 1 ) mod B
%  for BN.
%
    am1 = a - 1;
    anm1tb = ( an - 1 ) * b;

    [ bn, ierror ] = congruence ( am1, c, anm1tb );

    if ( ierror ~= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'LCRG_ANBN - Fatal error!\n' );
      fprintf ( 1, '  An error occurred in the CONGRUENCE routine.\n' );
      fprintf ( 1, '  The error code was IERROR = %d\n', ierror );
      error ( 'LCRG_ANGN - Fatal error!' );
    end

  end

  return
end


text:
function [ output_args ] = test_image_registration( input_args )
%TEST_IMAGE_REGISTRATION Summary of this function goes here
%   Detailed explanation goes here
Imoving=imread('images/lenag1.png');
Istatic=imread('images/lenag2.png');

% Register the images
[Ireg,O_trans,Spacing,M,B,F] = image_registration(Imoving,Istatic);

% Show the registration result
figure,
subplot(2,2,1), imshow(Imoving); title('moving image');
subplot(2,2,2), imshow(Istatic); title('static image');
subplot(2,2,3), imshow(Ireg); title('registerd moving image');
% Show also the static image transformed to the moving image
Ireg2=movepixels(Istatic,F);
subplot(2,2,4), imshow(Ireg2); title('registerd static image');

% Show the transformation fields
figure,
subplot(2,2,1), imshow(B(:,:,1),[]); title('Backward Transf. in x direction');
subplot(2,2,2), imshow(F(:,:,2),[]); title('Forward Transf. in x direction');
subplot(2,2,3), imshow(B(:,:,1),[]); title('Backward Transf. in y direction');
subplot(2,2,4), imshow(F(:,:,2),[]); title('Forward Transf. in y direction');

% Calculate strain tensors
E = strain(F(:,:,1),F(:,:,2));
% Show the strain tensors
figure,
subplot(2,2,1), imshow(E(:,:,1,1),[-1 1]); title('Strain Tensors Exx');
subplot(2,2,2), imshow(E(:,:,1,2),[-1 1]); title('Strain Tensors Exy');
subplot(2,2,3), imshow(E(:,:,2,1),[-1 1]); title('Strain Tensors Eyx');
subplot(2,2,4), imshow(E(:,:,2,2),[-1 1]); title('Strain Tensors Eyy');

end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function res = nconditions(this)
% Method for getting the number of unique conditions in the file
% FORMAT res = nconditions(obj)
% _______________________________________________________________________
% Copyright (C) 2008-2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: nconditions.m 5025 2012-10-31 14:44:13Z vladimir $

res = numel(condlist(this));

author: spm


text:
classdef TonemapDrago < handle
    %TONEMAPDRAGO  Tonemapping algorithm used to map HDR image to 8-bit range
    %
    % Adaptive logarithmic mapping is a fast global tonemapping algorithm that
    % scales the image in logarithmic domain.
    %
    % Since it's a global operator the same function is applied to all the
    % pixels, it is controlled by the bias parameter.
    %
    % Optional saturation enhancement is possible as described in [FL02].
    %
    % For more information see [DM03].
    %
    % ## References
    % [FL02]:
    % > Raanan Fattal, Dani Lischinski, and Michael Werman.
    % > "Gradient domain high dynamic range compression". In ACM Transactions
    % > on Graphics (TOG), volume 21, pages 249-256. ACM, 2002.
    %
    % [DM03]:
    % > Frederic Drago, Karol Myszkowski, Thomas Annen, and Norishige Chiba.
    % > "Adaptive logarithmic mapping for displaying high contrast scenes".
    % > In Computer Graphics Forum, volume 22, pages 419-426. Wiley Online
    % > Library, 2003.
    %
    % See also: cv.Tonemap, cv.TonemapDurand, cv.TonemapReinhard,
    %  cv.TonemapMantiuk, tonemap
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    properties (Dependent)
        % positive value for gamma correction.
        %
        % Gamma value of 1.0 implies no correction, gamma equal to 2.2 is
        % suitable for most displays. Generally gamma > 1 brightens the image
        % and gamma < 1 darkens it.
        Gamma
        % positive saturation enhancement value.
        %
        % 1.0 preserves saturation, values greater than 1 increase saturation
        % and values less than 1 decrease it.
        Saturation
        % value for bias function in [0, 1] range.
        %
        % Values from 0.7 to 0.9  usually give best results, default value
        % is 0.85
        Bias
    end

    %% TonemapDrago
    methods
        function this = TonemapDrago(varargin)
            %TONEMAPDRAGO  Creates TonemapDrago object
            %
            %     obj = cv.TonemapDrago()
            %     obj = cv.TonemapDrago('OptionName',optionValue, ...)
            %
            % ## Options
            % * __Gamma__ default 1.0
            % * __Saturation__ default 1.0
            % * __Bias__ default 0.85
            %
            % See also: cv.TonemapDrago.process
            %
            this.id = TonemapDrago_(0, 'new', varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.TonemapDrago
            %
            if isempty(this.id), return; end
            TonemapDrago_(this.id, 'delete');
        end
    end

    %% Tonemap
    methods
        function dst = process(this, src)
            %PROCESS  Tonemaps image
            %
            %     dst = obj.process(src)
            %
            % ## Input
            % * __src__ source RGB image, 32-bit `single` 3-channel array.
            %
            % ## Output
            % * __dst__ destination image of same size as `src`, 32-bit
            %   `single` 3-channel array with values in [0,1] range.
            %
            % See also: cv.TonemapDrago.TonemapDrago
            %
            dst = TonemapDrago_(this.id, 'process', src);
        end
    end

    %% Algorithm
    methods (Hidden)
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.TonemapDrago.empty, cv.TonemapDrago.load
            %
            TonemapDrago_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Returns true if the algorithm is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.TonemapDrago.clear, cv.TonemapDrago.load
            %
            b = TonemapDrago_(this.id, 'empty');
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.TonemapDrago.save, cv.TonemapDrago.load
            %
            name = TonemapDrago_(this.id, 'getDefaultName');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.TonemapDrago.load
            %
            TonemapDrago_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.TonemapDrago.save
            %
            TonemapDrago_(this.id, 'load', fname_or_str, varargin{:});
        end
    end

    %% Getters/Setters
    methods
        function value = get.Gamma(this)
            value = TonemapDrago_(this.id, 'get', 'Gamma');
        end
        function set.Gamma(this, value)
            TonemapDrago_(this.id, 'set', 'Gamma', value);
        end

        function value = get.Saturation(this)
            value = TonemapDrago_(this.id, 'get', 'Saturation');
        end
        function set.Saturation(this, value)
            TonemapDrago_(this.id, 'set', 'Saturation', value);
        end

        function value = get.Bias(this)
            value = TonemapDrago_(this.id, 'get', 'Bias');
        end
        function set.Bias(this, value)
            TonemapDrago_(this.id, 'set', 'Bias', value);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [clout]=cl_cl_intersect(clin1,clin2,varargin)

% Usage:
% clout=cl_cl_intersect(clin1,clin2)
% clout=cl_cl_intersect(clin1,clin2,'switch')
%
% clout is a cl structure containing clusters of only those voxels that
% exist in both cl structures clin1 and clin2.
%
% If the clin1 and clin2 are not in the same image space, clin2 will be
% converted to the space of clin1 by finding the nearest mm point in the
% space of clin1 to each mm point in clin2. NOTE: ***** This may be
% undiserable if the voxel sizes in clin1 are substantially smaller than
% those in clin2--eg., a 4x4x4 voxel will be translated into a single 2x2x2
% voxel, resulting in 3 'empty' 2x2x2 voxels in the space covered by the
% 4x4x4 voxel. For this reason, it is recommended that you place the lowest
% resolution cl structure in the first input.
%
% The cl.Z values and the threshold field are taken from clin1, and the
% values from clin2 are dropped. If the string 'switch' is included as in
% input, the Z values from clin2 will be used instead.
%

sw=0;
if ~isempty(varargin)
    for k=1:length(varargin)
        if strcmp(varargin{k},'switch')
            sw=1;
        end
    end
end

XYZmm1=[];
XYZmm2=[];
Z=[];
for k=1:length(clin1)
    if size(clin1(k).XYZmm,1)==3
        XYZmm1(:,end+1:end+clin1(k).numVox)=clin1(k).XYZmm;
    else
        XYZmm1(:,end+1:end+clin1(k).numVox)=clin1(k).XYZmm';
    end
    if ~sw
        Z(end+1:end+clin1(k).numVox)=clin1(k).Z;
    end
end

for k=1:length(clin2)
    if size(clin2(k).XYZmm,1)==3
        XYZmm2(:,end+1:end+clin2(k).numVox)=clin2(k).XYZmm;
    else
        XYZmm2(:,end+1:end+clin2(k).numVox)=clin2(k).XYZmm';
    end
    if sw
        Z(end+1:end+clin2(k).numVox)=clin2(k).Z;
    end
end

if clin1(1).voxSize(:)~=clin2(1).voxSize(:)
    [XYZmm2,p_ind]=mmToVoxel(XYZmm2,clin1(1).M,'valid');
    if sw
        for m=1:length(p_ind)
            ind=find(p_ind(:)==p_ind(m));
            Z(ind)=mean(Z(ind));
        end
        Z=Z(unique(p_ind));
    end
    XYZmm2=voxelToMm(XYZmm2,clin1(1).M);
end

XYZmm_out=[];
Zout=[];

for k=1:size(XYZmm1,2)
    a=find(XYZmm2(1,:)==XYZmm1(1,k));
    if ~isempty(a)
        b=find(XYZmm2(2,a)==XYZmm1(2,k));
        if ~isempty(b)
            c=find(XYZmm2(3,a(b))==XYZmm1(3,k));
            if ~isempty(c)
                XYZmm_out(:,end+1)=XYZmm2(:,a(b(c)));
                Zout(:,end+1)=Z(a(b(c)));
            end
        end
    end
end
if isempty(XYZmm_out)
	clout=struct('XYZ',{},'XYZmm',{},'name',{},'numVox',{},'M',{},'voxSize',{},'dim',{},'title',{},'threshold',{},'Z',{},'mm_center',{});
else
    XYZout=mmToVoxel(XYZmm_out,clin1(1).M,'valid');
    cl_ind=spm_clusters(XYZout);
    for k=1:max(cl_ind)
        clout(k).XYZ=XYZout(:,cl_ind==k);
        clout(k).XYZmm=voxelToMm(clout(k).XYZ,clin1(1).M);
        clout(k).name='';
        clout(k).numVox=length(find(cl_ind==k));
        clout(k).M=clin1(1).M;
        clout(k).voxSize=clin1(1).voxSize;
        clout(k).dim=clin1(1).dim;
        clout(k).title=['Cluster of ' num2str(clout(k).numVox) ' voxels from the intersection of two clusters, created by cl_cl_intersect.m.'];
        clout(k).threshold=clin1(1).threshold;
        clout(k).Z=Zout(cl_ind==k);
        clout(k).mm_center=center_of_mass(clout(k).XYZmm,clout(k).Z);
    end
end


text:
function fh = decomp_reconst_W(im,Nsc,filter,block,noise,parent,covariance,optim,sig);

% Decompose image into subbands, denoise using BLS-GSM method, and recompose again.
%		fh = decomp_reconst(im,Nsc,filter,block,noise,parent,covariance,optim,sig);
%       im:         image
%       Nsc:        number of scales
%       filter:     type of filter used (see namedFilters)
%       block:      2x1 vector indicating the dimensions (rows and columns) of the spatial neighborhood
%       noise:      signal with the same autocorrelation as the noise
%       parent:     include (1) or not (0) a coefficient from the immediately coarser scale in the neighborhood
%       covariance:	 are we considering covariance or just variance?
%       optim:		 for choosing between BLS-GSM (optim = 1) and MAP-GSM (optim = 0)
%       sig:        standard deviation (scalar for uniform noise or matrix for spatially varying noise)
% Version using a critically sampled pyramid (orthogonal wavelet), as implemented in MatlabPyrTools (Eero).

% JPM, Univ. de Granada, 3/03

if (block(1)/2==floor(block(1)/2))|(block(2)/2==floor(block(2)/2)),
   error('Spatial dimensions of neighborhood must be odd!');
end

if ~exist('parent'),
        parent = 1;
end

if ~exist('covariance'),
        covariance = 1;
end

if ~exist('optim'),
        optim = 1;
end

if ~exist('sig'),
        sig = sqrt(mean(noise.^2));
end

Nor = 3;    % number of orientations: vertical, horizontal and mixed diagonals (for compatibility)

[pyr,pind] = buildWpyr(im,Nsc,filter,'circular');
[pyrN,pind] = buildWpyr(noise,Nsc,filter,'circular');
pyrh = pyr;
Nband = size(pind,1);
for nband = 1:Nband-1, % everything except the low-pass residual
  fprintf('%d % ',round(100*(nband-1)/(Nband-1)))
  aux = pyrBand(pyr, pind, nband);
  auxn = pyrBand(pyrN, pind, nband);
  prnt = parent & (nband < Nband-Nor);   % has the subband a parent?
  BL = zeros(size(aux,1),size(aux,2),1 + prnt);
  BLn = zeros(size(aux,1),size(aux,2),1 + prnt);
  BL(:,:,1) = aux;
  BLn(:,:,1) = auxn;
  if prnt,
  	aux = pyrBand(pyr, pind, nband+Nor);
    auxn = pyrBand(pyrN, pind, nband+Nor);
    aux = real(expand(aux,2));
    auxn = real(expand(auxn,2));
    BL(:,:,2) = aux;
  	BLn(:,:,2) = auxn;
  end

  sy2 = mean2(BL(:,:,1).^2);
  sn2 = mean2(BLn(:,:,1).^2);
  if sy2>sn2,
     SNRin = 10*log10((sy2-sn2)/sn2);
  else
     disp('Signal is not detectable in noisy subband');
  end

  % main
  BL = denoi_BLS_GSM_band(BL,block,BLn,prnt,covariance,optim,sig);
  pyrh(pyrBandIndices(pind,nband)) = BL(:)';
end
fh = reconWpyr(pyrh,pind,filter,'circular');

author: ricedsp


text:
function [ x, w ] = rule06 ( n )

%*****************************************************************************80
%
%% RULE06 returns the rule of degree 6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 July 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, integer N, the number of nodes.
%
%    Output, real X(2,N), the coordinates of the nodes.
%
%    Output, real W(N), the weights.
%
  xs = [ ...
    -.7347550761673839,0.8662152988634034, ...
    0.1596873653424614,-.8905137714296896, ...
    -.1469707846748791,0.9240009259977663, ...
    -.8463324986375500,-.4086308482879689, ...
    0.5175294652720337,0.4801002492857063 ];
  ys = [ ...
    0.8933891941643415,-.7037359670513631, ...
    -.9085856749287847,0.1644347368502312, ...
    0.5352177835541986,0.4879643750888035, ...
    -.8394767448218339,-.4262330870004397, ...
    0.9176357850707917,-.1009764561823168 ];
  ws = [ ...
    0.1541850714382379,0.1900556513689156, ...
    0.2246942645703077,0.2465847648329768, ...
    0.5062382287542438,0.1829226437278864, ...
    0.1373586623279704,0.4754388545735908, ...
    0.1856913242244974,0.5252576589275637 ];

  x(1,1:n) = xs(1:n);
  x(2,1:n) = ys(1:n);
  w(1:n) = ws(1:n);

  return
end


text:
function overlap_ratio = get_overlap_1toN(rect1, rect2, cover)
% rect1: 1*4
% rect2: N*4 (left top right bottom)
% cover: 0: o/(s1+s2-o) 1: o/s1 2:o/min(s1,s2) 3: 0/s2
    if nargin < 3
        cover = 0;
    end
    rect1 = repmat(rect1, [size(rect2, 1) 1]);

    area1 = (rect1(:, 3) - rect1(:, 1)) .* (rect1(:, 4) - rect1(:, 2));
    area2 = (rect2(:, 3) - rect2(:, 1)) .* (rect2(:, 4) - rect2(:, 2));

    l = max(rect1(:, 1), rect2(:, 1));
    r = min(rect1(:, 3), rect2(:, 3));
    t = max(rect1(:, 2), rect2(:, 2));
    b = min(rect1(:, 4), rect2(:, 4));

    w = r - l;
    h = b - t;
    overlap = w .* h;
    overlap(w < 0 | h < 0) = 0;
    if cover == 1
        overlap_ratio = overlap ./ area1;
    elseif cover == 2
        overlap_ratio = overlap ./ min(area1, area2);
    elseif cover == 3
        overlap_ratio = overlap ./ area2;
    else
        overlap_ratio = overlap ./ (area1 + area2 - overlap);
    end
end

author: liuyuisanai

repo: RSA-for-object-detection


text:
list_M_species={};

for i=1:2; % reaction
    for t=1:size(list_M(i,:),2) % nubmer of metabolites
        if ~isempty(list_M{i,t})
            met(i,t)=retrieveMet(parsed_fatty_acid_new,list_M(i,t))
            list_M_species(i,2*t-1)={met(i,t).speciesAlliens}

        end
    end
end

for i=1:length(num)

    for r=2:2:length(baseR.(list_nodes{i}))

        list_R(i,r/2)=baseR.(list_nodes{i})(r);
    end

    for d=1:2:length(baseR.(list_nodes{i}))-1
        list_R_M(i,(d+1)/2)=baseR.(list_nodes{i})(d)
    end

    for r=2:2:length(baseP.(list_nodes{i}))

        list_P(i,r/2)=baseP.(list_nodes{i})(r);
    end

    for d=1:2:length(baseP.(list_nodes{i}))-1
        list_P_M(i,(d+1)/2)=baseP.(list_nodes{i})(d)
    end

    for r=2:2:length(baseC_R.(list_nodes{i}))
        list_C_R(i,r/2)=baseC_R.(list_nodes{i})(r); % each reaction has the same color
    end

    for r=2:2:length(baseC_P.(list_nodes{i}))
        list_C_P(i,r/2)=baseC_P.(list_nodes{i})(r); % each reaction has the same color
    end

end

author: opencobra

repo: cobratoolbox


text:
function [y,x,x0,eta,e,u] = simulateNLSS(n_t,f_fname,g_fname,theta,phi,u,alpha,sigma,options,x0)
% legacy code
s = warning ('on');
warning ('*** The function `simulateNLSS` is now deprecated. Please use `VBA_simulate` instead (same syntax).')
warning (s);

% fallback
try
    [y,x,x0,eta,e,u] = VBA_simulate (n_t,f_fname,g_fname,theta,phi,u,alpha,sigma,options,x0);
catch
    [y,x,x0,eta,e,u] = VBA_simulate (n_t,f_fname,g_fname,theta,phi,u,alpha,sigma,options);
end

author: MBB-team

repo: VBA-toolbox


text:
%[2015]-"TSA: Tree-seed algorithm for continuous optimization"

% (9/12/2020)

function TSA = jTreeSeedAlgorithm(feat,label,opts)
% Parameters
lb    = 0;
ub    = 1;
thres = 0.5;
ST    = 0.1;    % switch probability

if isfield(opts,'T'), max_Iter = opts.T; end
if isfield(opts,'N'), N = opts.N; end
if isfield(opts,'ST'), ST = opts.ST; end
if isfield(opts,'thres'), thres = opts.thres; end

% Objective function
fun = @jFitnessFunction;
% Number of dimensions
dim = size(feat,2);
% Initial (5)
X   = zeros(N,dim);
for i = 1:N
  for d = 1:dim
    X(i,d) = lb + (ub - lb) * rand();
  end
end
% Fitness
fit = zeros(1,N);
for i = 1:N
  fit(i) = fun(feat,label,(X(i,:) > thres),opts);
end
% Best solution (6)
[fitG, idx] = min(fit);
Xgb         = X(idx,:);
% Maximum & minimum number of seed
Smax = round(0.25 * N);
Smin = round(0.1 * N);
% Pre
curve = zeros(1,max_Iter);
curve(1) = fitG;
t = 2;
% Iteration
while t <= max_Iter
  for i = 1:N
    % Random number of seed
    num_seed = round(Smin + rand()* (Smax - Smin));
    Xnew     = zeros(num_seed, dim);
    for j = 1:num_seed
      % Random select a tree, but not i
      RN = randperm(N);
      RN(RN == i) = [];
      r  = RN(1);
      for d = 1:dim
        % Alpha in [-1,1]
        alpha = -1 + 2 * rand();
        if rand() < ST
          % Generate seed (3)
          Xnew(j,d) = X(i,d) + alpha * (Xgb(d) - X(r,d));
        else
          % Generate seed (4)
          Xnew(j,d) = X(i,d) + alpha * (X(i,d) - X(r,d));
        end
      end
      % Boundary
      XB = Xnew(j,:); XB(XB > ub) = ub; XB(XB < lb) = lb;
      Xnew(j,:) = XB;
    end
    % Fitness
    for j = 1:num_seed
      % Fitness
      Fnew = fun(feat,label,(Xnew(j,:) > thres),opts);
      % Greedy selection
      if Fnew < fit(i)
        fit(i) = Fnew;
        X(i,:) = Xnew(j,:);
      end
    end
  end
  % Best solution (6)
  [fitG_new, idx] = min(fit);
  Xgb_new         = X(idx,:);
  % Best update
  if fitG_new < fitG
    fitG = fitG_new;
    Xgb  = Xgb_new;
  end
  % Store
  curve(t) = fitG;
  fprintf('\nIteration %d Best (TSA)= %f',t,curve(t))
  t = t + 1;
end
% Select features
Pos   = 1:dim;
Sf    = Pos((Xgb > thres) == 1);
sFeat = feat(:,Sf);
% Store results
TSA.sf = Sf;
TSA.ff = sFeat;
TSA.nf = length(Sf);
TSA.c  = curve;
TSA.f  = feat;
TSA.l  = label;
end

author: JingweiToo

repo: Wrapper-Feature-Selection-Toolbox


text:
function [ y, next ] = p39_equil ( neqn, next )

%*****************************************************************************80
%
%% P39_EQUIL returns equilibrium solutions of problem p39.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, integer NEXT, the index of the previous
%    equilibrium, which should be 0 on first call.
%
%    Output, real Y(NEQN), the "next" equilibrium solution, if any.
%
%    Output, integer NEXT, the index of the current equilibrium,
%    or 0 if there are no more.
%
  next = 0;
  y = [];

  return
end

author: johannesgerer


text:
function opt = hmc2_opt(opt)
%HMC2_OPT  Default options for Hybrid Monte Carlo sampling.
%
%  Description
%    OPT = HMC2_OPT
%      return default options
%    OPT = HMC2_OPT(OPT)
%      fill empty options with default values
%
%  The options and defaults are
%    display (0)
%      1 to display the energy values and rejection threshold at
%        each step of the Markov chain
%      2 to display also position vectors at each step
%    checkgrad (0)
%      1 to check the user defined gradient function
%    steps (10)
%      Defines the trajectory length (i.e. the number of leapfrog
%      steps at each iteration)
%    nsamples (1)
%      the number of samples retained from the Markov chain
%    nomit (0)
%      the number of samples omitted from the start of the chain
%    persistence (0)
%      0 for complete replacement of momentum variables
%      1 if momentum persistence is used
%    decay (0.9)
%      defines the decay used when a persistent update of
%      (leap-frog) momentum is used. Bounded to the interval [0, 1.)
%    stepadj (0.1)
%      the step adjustment used in leap-frogs
%    stepsf ([])
%      the step size function
%    window (1)
%      the size of the acceptance window
%
%  See also
%    HMC2

%	Copyright (c) Aki Vehtari (1998)

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if nargin < 1
  opt=[];
end

if ~isfield(opt,'display')
  opt.display=0;
end
if ~isfield(opt,'checkgrad')
  opt.checkgrad=0;
end
if ~isfield(opt,'steps') | opt.steps < 1
  opt.steps=10;
end
if ~isfield(opt,'nsamples') | opt.nsamples < 1
  opt.nsamples=1;
end
if ~isfield(opt,'nomit') | opt.nomit < 0
  opt.nomit=0;
end
if ~isfield(opt,'persistence')
  opt.persistence=0;
end
if ~isfield(opt,'decay') | opt.decay < 0 | opt.decay > 1
  opt.decay=0.9;
end
if ~isfield(opt,'stepadj')
  opt.stepadj=0.1;
end
if ~isfield(opt,'stepsf')
  opt.stepsf=[];
end
if ~isfield(opt,'window') | opt.window < 0
  opt.window=1;
end
if opt.window > opt.steps
  opt.window=opt.steps;
end

author: fieldtrip

repo: fieldtrip