text:
function x = spm_normrnd(m, C, N)
% Random samples from Gaussian distribution
% FORMAT x = spm_normrnd(m, C, N)
% m        - [d x 1] mean
% C        - [d x d] covariance or cell array {dC, vC} so that
%            [vC, diag(dC)] = eig(C)
% N        - number of samples
%
% x        - [d x N] matrix of samples
%__________________________________________________________________________
% Copyright (C) 2009 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_normrnd.m 3603 2009-11-30 18:56:50Z guillaume $

if iscell(C)
    deig = C{1};
    evec = C{2};
else
    [evec, eval] = eig(C);
    deig         = diag(eval);
end

i = (abs(imag(deig))>0) | (deig<0);
if any(i)
  %warning('Covariance matrix is not positive semi-definite: redefined');
  deig(i)   = [];
  evec(:,i) = [];
end

proj = randn(N,length(deig)) * diag(sqrt(deig));
x    = repmat(m(:),1,N) + evec*proj';

author: spm


text:
%
%      This HDR Toolbox demo creates an LDR stack from an HDR image:
%	   1) Read a stack of LDR images
%	   2) Read exposure values from the EXIF
%	   3) Estimate the Camera Response Function (CRF)
%	   4) Build the radiance map using the stack and stack_exposure
%	   5) Save the radiance map in .hdr format
%	   6) Show the tone mapped version of the radiance map
%
%       Author: Francesco Banterle
%       Copyright 2016 (c)
%

clear all;

img = hdrimread('Bottles_Small.hdr');

name_folder_h = 'demos/output/hdr_to_stack_histogram';
name_folder_u = 'demos/output/hdr_to_stack_uniform';
name_folder_s = 'demos/output/hdr_to_stack_selected';

format = 'jpg';

%
% Histogram
%

mkdir(name_folder_h);

[stack, stack_exposure] = CreateLDRStackFromHDR( img, 2, 'histogram', 'gamma', 2.2);

for i=1:size(stack, 4)
    name_out = [name_folder_h, '/exp_', num2str(stack_exposure(i)), '.', format];
    imwrite(stack(:,:,:,i), name_out);
end

%
% Uniform
%

mkdir(name_folder_u);

[stack, stack_exposure] = CreateLDRStackFromHDR( img, 2, 'uniform', 'gamma', 2.2);

for i=1:size(stack, 4)
    name_out = [name_folder_u, '/exp_', num2str(stack_exposure(i)), '.', format];
    imwrite(stack(:,:,:,i), name_out);
end

%
% Selected
%

mkdir(name_folder_s);

selected_exps = [-1.0 0.0 1.0];

[stack, stack_exposure] = CreateLDRStackFromHDR( img, selected_exps, 'selected', 'gamma', 2.2);

for i=1:size(stack, 4)
    name_out = [name_folder_s, '/exp_', num2str(stack_exposure(i)), '.', format];
    imwrite(stack(:,:,:,i), name_out);
end

author: banterle


text:
% Reg1_test.m
% test of Reg1 regularization object

if 1 % test single-pixel case
	R0 = Reg1(true(1));
	x = 7;
	jf_equal(0, R0.penal(R0, x))
	jf_equal(0, R0.cgrad(R0, x))
	jf_equal(0, R0.denom(R0, x))
	jf_equal(full(R0.C), nan(0,1))
end

if 1 % test tiny image
	R0 = Reg1(true(3,2), 'offsets', '2d:hv');
	x = reshape(1:3*2, 3, 2);
	R0.penal(R0, x);
	R0.cgrad(R0, x);
	R0.denom(R0, x);
	full(R0.C);

	tmp = {true(3,2), 'offsets', '2d:hv'};
	if 0
		rng(7)
		tmp = {tmp{:}, 'user_wt', ones(3,2,2)};
	end
	R0t = Reg1(tmp{:}, 'type_penal', 'mat', 'type_diff', 'ind');
	R0x = Reg1(tmp{:}, 'type_penal', 'mex');
	jf_equal( full(R0t.C), full(R0x.C) )
	jf_equal( full(R0t.C1), full(R0x.C1) )
%	jf_equal( R0t.cgrad(R0t, x), R0x.cgrad(R0x, x) ) % not ?
%	full(R0t.C1)
%	full(R0x.C1)
end

if 1 || ~isvar('Rt'), printm 'Rt'
	if 1 % 3d
		ig = image_geom('nx', 512, 'ny', 480, 'nz', 2^5, 'fov', 500, ...
			'down', 2^2);
	else % 2d
		ig = image_geom('nx', 512, 'ny', 448, 'fov', 500, 'down', 2^3);
	end
	tmp = ig.circ(ig.fov/2*1.1) > 0;
	tmp([1 end],:,:) = 0; tmp(:, [1 end],:) = 0; % zero border for 3d
	ig.mask = tmp; clear tmp
	kappa = 2 + 1 * cos(2*pi*ig.x/ig.fov) * sin(2*pi*ig.y/ig.fov)';
%	kappa = ig.mask;
	f.offsets = ''; f.n_offset = 4;
	if ig.is3
		kappa = repmat(ig.mask_or, [1 1 ig.nz]) .* ig.mask;
		f.offsets = '3d:26';
		f.n_offset = 13;
	else
		kappa = single(kappa) .* ig.mask;
	end

	if 0
		profile on
		cpu etic
		pn = jf_protected_names;
		for qq=1:29*2
			ii = 1:prod(ig.dim);
%			[jx jy jz] = ind2sub(ig.dim, ii);
			jj = pn.ind2sub(ig.dim, ii);
		end
		cpu etoc
		profile report
	return
	end

	f.l2b = 3;
%	f.pot = 'hyper3'; f.delta = 1; f.pot_arg = {f.pot, f.delta};
%	f.pot = 'quad'; f.delta = inf; f.pot_arg = {f.pot};
%	f.pot = 'qgg2'; f.delta = 2; f.pot_arg = {f.pot, f.delta, 1.2};
	f.pot = 'gf1'; f.delta = 2; ...
		tmp = potential_fun('gf1-fit', nan, {'qgg2', [1 10], 1.2});
		f.pot_arg = {f.pot, f.delta, tmp};

%	x = ig.unitv(ig.nx, ig.ny) + ig.unitv(ig.nx/4, 5+0*ig.ny/2);
	x = single(ig.mask + ig.circ(ig.fov/4));
	rng(0)
	x = single(rand(ig.dim)) .* ig.mask;
	if ig.is3
		x = x + ig.unitv ...
			+ ig.unitv(3*ig.nx/8, ig.ny/2, 1) ...
			+ ig.unitv(ig.nx/2, 5*ig.ny/8, ig.nz);
	end
	xm = x(ig.mask);

	for order = 1:2
		pr order
		f.arg1 = {kappa, 'offsets', f.offsets, 'beta', 2^f.l2b, ...
			'edge_type', 'tight', 'order', order};
		if 0 && order == 1 % test user_wt
			rng(7)
%			tmp = ones([ig.dim f.n_offset]);
%			tmp(1:end/4,:,:) = 4;
			tmp = rand([ig.dim f.n_offset]);
			f.arg1 = {f.arg1{:}, 'user_wt', tmp};
		end
		f.arg = {f.arg1{:}, 'pot_arg', f.pot_arg};
		Rt = Reg1(f.arg{:}, 'type_penal', 'mat');
		Rx = Reg1(f.arg{:}, 'type_penal', 'mex'); % default control
%		R1 = Reg1(f.arg{:}, 'type_penal', 'mex', 'control', 1);

		if 0 % test timing of loop control 1 and 2
			% in 3d for order=2 with no user_wt,
			% control=2 is about 5% faster
			g1 = ig.embed(R1.cgrad(R1, xm)); % warm up
			g2 = ig.embed(Rx.cgrad(Rx, xm)); % warm up

			cpu etic
			g1 = ig.embed(R1.cgrad(R1, xm));
			cpu etoc 'Rx cgrad control=1 time'

			cpu etic
			g2 = ig.embed(Rx.cgrad(Rx, xm));
			cpu etoc 'Rx cgrad control=2 time'
			equivs(g1, g2)
		continue
		end

		if has_mex_jf
			switch f.pot
			case {'qgg2', 'gf1'}
				tmp = {f.delta, f.pot_arg{3}};
			otherwise
				tmp = f.delta;
			end
			Ro = Robject(f.arg1{:}, 'type_denom', 'matlab', ...
				'potential', f.pot, 'delta', tmp, ...
				'distance_power', Rt.distance_power);
		end

		if 1 && ig.is3 && has_mex_jf
			zxy = @(x) permute(x, [3 1 2]);
			xyz = @(x) permute(x, [2 3 1]);
			% convert offsets to zxy!
			f.offsets_zxy = reg_offset_xyz_to_zxy(f.offsets, size(kappa));
			Rz = Reg1(zxy(kappa), f.arg{2:end}, ...
				'user_wt', [], ...
				'offsets', f.offsets_zxy, ...
				'offsets_is_zxy', true, ...
				'type_penal', 'zxy');
		end

%		Rm = Rmem(kappa, 'pot_arg', f.pot_arg, 'beta', 2^f.l2b, ...
%			'edge_type', 'tight', ...
%			'distance_power', Rg.distance_power, ...
%			'offsets', f.offsets);

		if 1 % check C1 *
			jf_equal(Rt.C1 * x, Rx.C1 * x)
			jf_equal(Rt.C1 * xm, Rx.C1 * xm)
			jf_equal(ig.shape(Rt.C1 * xm), Rt.C1 * x)
			if has_mex_jf
				jf_equal(Ro.C1 * x, Rx.C1 * x)
			end
		end

		if 0 % check C1 (only for tiny cases)
			c1 = Rt.C1; c1 = c1(:,:);
			c2 = Rx.C1; c2 = c1(:,:);
			jf_equal(c1, c2)
			if has_mex_jf
				co = Ro.C1; co = co(:,:);
				jf_equal(c1, co)
			end
%			clear c1 c2
		end

		if 1 % check dercurv
			[t.dt t.ct] = feval(Rt.dercurv, Rt, Rt.C1 * x);
			[t.dx t.cx] = feval(Rx.dercurv, Rx, Rx.C1 * x);
%			[t.d1 t.c1] = feval(R1.dercurv, R1, R1.C1 * x);
			jf_equal(t.dt, t.dx)
			jf_equal(t.ct, t.cx)
%			jf_equal(t.dt, t.d1)
%			jf_equal(t.ct, t.c1)

			[t.dx t.cx] = feval(Rx.dercurv, Rx, Rx.C1 * xm);
			jf_equal(ig.shape(t.dx), t.dt)
			jf_equal(ig.shape(t.cx), t.ct)

			if has_mex_jf
				[t.do t.co] = feval(Ro.dercurv, Ro, Ro.C1 * xm);
				t.do = ig.shape(t.do);
				t.co = ig.shape(t.co);
				equivs(t.dt, t.do)
				% trick: curv matches *except* at outer edge(s)
				equivs(t.ct(2:end-1,2:end-1,:), t.co(2:end-1,2:end-1,:))
			end
		end

		if 1 && streq(f.pot, 'quad') % check C *
%			tmp = Rt.C * x;
			jf_equal(ig.shape(Rt.C * xm), Rt.C * x)
		end

		if 1 % check penalty value
			tmp = Rt.penal(Rt, x);
			jf_equal(tmp, Rt.penal(Rt, xm))
			equivs(tmp, Rx.penal(Rx, x))
%			equivs(tmp, R1.penal(R1, x))
			jf_equal(Rx.penal(Rx, x), Rx.penal(Rx, xm))
			if has_mex_jf
%				jf_equal(Rt.penal(Rt, xm), Ro.penal(Ro, xm))
				equivs(tmp, Ro.penal(Ro, xm), 'thresh', 18e-6)
			end
		end

		if 1 % check cgrad
			g1 = ig.embed(Rt.cgrad(Rt, xm));
			w1 = ig.embed(Rt.denom(Rt, xm));
			w1s = ig.embed(Rt.denom_sqs1(Rt, xm));
			if 1
				cpu etic
				g2 = Rt.cgrad(Rt, x);
				cpu etoc 'Rt cgrad time'
				jf_equal(g1, g2)
			end

			if has_mex_jf
				go = ig.embed(Ro.cgrad(Ro, xm));
				equivs(g1, go)

				g3 = ig.embed(Rx.cgrad(Rx, xm));
				equivs(g1, g3)

				if 1 % check mex version w/o mask
					cpu etic
					g4 = Rx.cgrad(Rx, x);
					cpu etoc 'Rx cgrad time'
					jf_equal(g3, g4)
				end

				if ig.is3 && order == 1
					tmp = zxy(x);
					% check zxy version w/o mask
					cpu etic
					[g5 w5] = feval(Rz.cgrad_denom, Rz, tmp);
					cpu etoc 'Rz cgrad/denom time'
					g5 = xyz(g5);
					w5 = xyz(w5);
					if isempty(Rx.user_wt) % trick
						equivs(g1, g5)
						equivs(w1, w5)
					end
%					im(w1-w5), return

					% zxy version with mask
					mask_zxy = zxy(ig.mask);
					[g5m w5m] = feval(Rz.cgrad_denom, Rz, tmp(mask_zxy));
					g5m = xyz(embed(g5m, Rz.mask));
					w5m = xyz(embed(w5m, Rz.mask));
					jf_equal(g5, g5m)
					jf_equal(w5, w5m)
				end
			end

			if 0 % check empirical gradient (slow)
				ge = ig.embed(Rx.egrad(Rx, xm, 0.001));
				max_percent_diff(ge, g1)
				im plc 2 3
				im(1, x, 'x'), cbar
				im(2, g1, 'g mat'), cbar
				im(3, g3, 'g mex'), cbar
				im(4, ge, 'g emp'), cbar
				im(5, g1-ge, 'mat-emp'), cbar
				im(6, g3-ge, 'mex-emp'), cbar
			prompt
			end

%			im(g3(:,:,[1 end])-g1(:,:,[1 end]))
%			im clf, im([g1; g3; g1-g3]), cbar

%			clear g1 g2 g3 g4
		end

		if 1 % check denom
			dentf = Rt.denom_sqs1(Rt, x);
			dentm = Rt.denom_sqs1(Rt, xm);
			dentm = ig.embed(dentm);
			equivs(dentf, dentm)

			if has_mex_jf
				denxf = Rx.denom_sqs1(Rx, x);
				denxm = Rx.denom_sqs1(Rx, xm);
				denxm = ig.embed(denxm);
				equivs(dentf, denxf)
				equivs(denxm, denxf)
%				im clf, im([denxf; denxm; denxf-denxm]), cbar

%				denx1 = R1.denom_sqs1(R1, x);
%				equivs(dentf, denx1)

				dentf = Rt.denom(Rt, x);
				denxf = Rx.denom(Rx, x);
%				denx1 = R1.denom(R1, x);
				equivs(dentf, denxf)
%				equivs(dentf, denx1)

				deno = ig.embed(Ro.denom(Ro, xm));
				% trick: denom matches except along outer edges
				ix = (1+order):(ig.nx-order);
				iy = (1+order):(ig.ny-order);
				equivs(dentf(ix,iy), deno(ix,iy))

				if order == 1 && ig.is3 && isempty(Rx.user_wt)
					equivs(dentf, w5) % check zxy
				end
			end
		end

		if 0 % check diag ?
			%d = Rm.diag(R);
		end

		% test threads
	end
prompt
end

return %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Robject designed to match Rmem'
if ~isvar('Ro'), printm 'Ro'
	tmp = repmat(kappa.^2, [1 1 1 length(Rm.offsets)]);
	Ro = Robject(ig.mask, 'edge_type', 'leak', 'beta', 2^f.l2b, ...
		'offsets', Rm.offsets, ...
		'potential', f.pot, 'delta', f.delta, ...
		'type_denom', 'matlab', ...
		'distance_power', 0, 'user_wt', tmp);
	clear tmp
end

cpu etic
g1 = Rm.cgrad(Rm, x);
cpu etoc 'Rm cgrad'
cpu etic
g2 = Ro.cgrad(Ro, x(ig.mask));
g2 = embed(g2, ig.mask);
cpu etoc 'Ro cgrad'
max_percent_diff(g1, g2)
if ~isequal(g1, g2), error 'cgrad mismatch', end, clear g2

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%ISSTACKED Test on stacked mapping
%
%   N = ISSTACKED(W)
%   ISSTACKED(W)
%
% INPUT
%   W  Mapping
%
% OUTPUT
%   N  Scalar, 1 if W is a stacked mapping, 0 otherwise
%
% DESCRIPTION
% Returns 1 for stacked mappings. If no output is requested, false outputs
% are turned into errors. This may be used for assertion.
%
% SEE ALSO (PRTools Guide)
% ISMAPPING, ISPARALLEL

% $Id: isstacked.m,v 1.2 2006/03/08 22:06:58 duin Exp $

function n = isstacked(w)

		if (isa(w,'prmapping')) & (strcmp(w.mapping_file,'stacked'))
%			strcmp(w.mapping_file,'dyadicm')
		n = 1;
	else
		n = 0;
	end

	% Generate an error if the input is not a stacked mapping and no output
	% is requested (assertion).

	if (nargout == 0) & (n == 0)
		error([newline '---- Stacked mapping expected -----'])
	end

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function squareStokesP1P0
%STOKE Summary of this function goes here
%   Detailed explanation goes here

close all; clear all; clc;
%figure(1); set(gcf,'Units','normal'); set(gcf,'Position',[0,0,0.7,0.75]);
%---------------------- Parameters ----------------------------------------
maxIt = 3; N = zeros(maxIt,1);
erru = zeros(maxIt,1);  errp = zeros(maxIt,1);
%---------------------- Generate initial mesh -----------------------------
node = [0,0; 1,0; 1,1; 0,1];    % nodes
elem = [2,3,1; 4,1,3];          % elements
for k = 1:3
    [node,elem] = uniformbisect(node,elem);
end
%% pde
pde.f = @f; pde.g_D = @g_D; pde.p=@exactp;
%---------------------- Finite Element Method -----------------------------
%
for k = 1:maxIt
    [node,elem] = uniformbisect(node,elem);
    [u,p,A,nodeRefine,nodeFine,elemFine] = mgStokesP1P0TwolevelNew(node,elem,pde);
    N(k) = length(u) + length(p);
    uI = exactu(nodeRefine);
    erru(k) = sqrt((u-uI(:))'*A{2}*(u-uI(:)));
    center = (nodeFine(elemFine(:,1),:)+nodeFine(elemFine(:,2),:)+nodeFine(elemFine(:,3),:))/3;
    pI = exactp(center);
    pI = pI-pI(end);
    errp(k)=max(abs(p-pI));
%     errp(k) = norm(p - pI)/sqrt(N);
    fprintf('erru = %6.3g, errp = %6.3g\n', erru(k), errp(k));
end
%% Plot convergence rates
figure(2);
showrate2(N,erru,1,'-*','||Du-Du_h||',N,errp,1,'m-+','||p-p_h||');
end
%---------------------- End of function SQUARE -----------------------------
function z = f(p)
x = p(:,1); y = p(:,2);
z(:,1) = -4*pi^2*(2*cos(2*pi*x)-1).*sin(2*pi*y)+x.^2;
z(:,2) = 4*pi^2*(2*cos(2*pi*y)-1).*sin(2*pi*x);
end
function z = exactu(p)
x = p(:,1); y = p(:,2);
z(:,1) = (1-cos(2*pi*x)).*sin(2*pi*y);
z(:,2) = -(1-cos(2*pi*y)).*sin(2*pi*x);
end
function z=exactp(p)
x = p(:,1);
z = 1/3*x.^3;
end

% %---------------------- Sub functions called by CUBE ----------------------
% function z = f(p) % load data (right hand side function)
% x = p(:,1); y = p(:,2);
% z(:,1) = 2^10*( (1-6*x+6*x.^2).*(y-3*y.^2+2*y.^3) ...
%        + (x.^2-2*x.^3+x.^4).*(-3+6*y)...
%        - (-3+6*x).*(y.^2-2*y.^3+y.^4) );
% z(:,2) = -(2^10)*( (-3+6*x).*(y.^2-2*y.^3 ...
%        + y.^4)+(x-3*x.^2+2*x.^3).*(1-6*y+6*y.^2)...
%        + (1-6*x+6*x.^2).*(y-3*y.^2+2*y.^3) );
% end
% %--------------------------------------------------------------------------
% function z = exactu(p)
% x = p(:,1); y = p(:,2);
% z(:,1) = -(2^8)*(x.^2-2*x.^3+x.^4).*(2*y-6*y.^2+4*y.^3);
% z(:,2) = 2^8*(2*x-6*x.^2+4*x.^3).*(y.^2-2*y.^3+y.^4);
% end
% %--------------------------------------------------------------------------
% function z = exactp(p)
% x = p(:,1); y = p(:,2);
% z = -(2^8)*(2-12*x+12*x.^2).*(y.^2-2*y.^3+y.^4);
% end
%--------------------------------------------------------------------------
function z = g_D(p) % Dirichlet boundary condition
z = exactu(p);
end

repo: Aero-matlab


text:
function a = abs(a)
%ABS          Taylor absolute value abs(a)
%
%   c = abs(a)
%
%Result is convex hull of left and right derivative. This allows
%also input containing zero. The expansion
%  f(x) = f(xs) + f'(zeta)*(x-xs)
%for some zeta in ch(x,xs) is still valid
%

% written  05/21/09     S.M. Rump
%

  if ~isreal(a.t)
    a.t = abs(a.t);
    a.t(2:end,:) = repmat(NaN,size(a.t(2:end,:)));
    return
  end

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  at = a.t(1,:);
  index = ( at<0 );
  a.t(:,index) = -a.t(:,index);
  if isa(at,'intval')
    index = in0(0,at(1,:));
    a.t(:,index) = infsup(-1,1)*a.t(:,index);
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function [A,S]=jade(X,m,Wprev)
% ICA implementation using JADE (used by the ICA function).
%
%   Inputs:
%       X               = Each column of X is a sample from the n sensors.
%       m               = Number of source signals (optional).
%       Wprev           = Timepoint at which to start process (default = 0 seconds).
%
%   Outputs:
%       A               = An n x m estimate of an unknown matrix.
%       S               = An m x T naive estimate of the source signals.
%
% Implementation by: JF Cardoso. (Version 1.5)
%
% Please Cite:
% Reference:
%  @article{CS_iee_94,
%   author = "Jean-Fran\c{c}ois Cardoso and Antoine Souloumiac",
%   journal = "IEE Proceedings-F",
%   title = "Blind beamforming for non {G}aussian signals",
%   number = "6",
%   volume = "140",
%   month = dec,
%   pages = {362-370},
%   year = "1993"}

% Determine the no of rows and colums of the observation matrix.
[n,T]	= size(X);
if nargin==1, m=n ; end;

nem     = m;                    % The nuumber of eigen-matrices to be diagonalized.
seuil	= 1/sqrt(T)/100;    % A statistical threshold for stopping the joint diag.

%% Whiten the Matrix:
% Assuming white noise:
if m<n,
 	[U,D]       = eig((X*X')/T);
	[pu,k]   = sort(diag(D));
 	ibl         = sqrt(pu(n-m+1:n)-mean(pu(1:n-m)));
 	bl          = ones(m,1) ./ ibl ;
 	W           = diag(bl)*U(1:n,k(n-m+1:n))';
 	IW          = U(1:n,k(n-m+1:n))*diag(ibl);
% Assuming no noise:
else
 	IW          = sqrtm((X*X')/T);
 	W           = inv(IW);
end;
Y	= W*X;

%% Cumulant Estimation:
R       = (Y*Y' )/T ;
C       = (Y*Y.')/T ;
Q       = zeros(m*m*m*m,1) ;
index	= 1;

for lx = 1:m ; Yl 	= Y(lx,:);
    for kx = 1:m ; Ykl 	= Yl.*conj(Y(kx,:));
        for jx = 1:m ; Yjkl	= Ykl.*conj(Y(jx,:));
            for ix = 1:m ;
                Q(index) = ...
                (Yjkl * Y(ix,:).')/T -  R(ix,jx)*R(lx,kx) -  R(ix,kx)*R(lx,jx) -  C(ix,lx)*conj(C(jx,kx))  ;
                index	= index + 1;
            end
        end
    end
end

%% Compute and Reshape the significant Eigen Matrices:
[U,D]	= eig(reshape(Q,m*m,m*m));
[la,K]	= sort(abs(diag(D)));

M	= zeros(m,nem*m);
Z	= zeros(m)	;
h	= m*m;
for u=1:m:nem*m,
	Z(:)            = U(:,K(h));
	M(:,u:u+m-1)	= la(h)*Z;
	h               = h-1;
end;

%% Approximate the Diagonalization of the Eigen Matrices:
B       = [ 1 0 0 ; 0 1 1 ; 0 -i i ] ;
Bt      = B' ;

encore	= 1;
if exist('Wprev','var'), V=inv(Wprev); else, V	= eye(m); end

%% Main Loop:
while encore, encore=0;
    for p = 1:m-1,
        for q = p+1:m,

            Ip = p:m:nem*m ;
            Iq = q:m:nem*m ;

            % Computing the Givens Angles:
            g	= [ M(p,Ip)-M(q,Iq)  ; M(p,Iq) ; M(q,Ip) ] ;
            [vcp,D] = eig(real(B*(g*g')*Bt));
            [la, K]	= sort(diag(D));
            angles	= vcp(:,K(3));
            if angles(1)<0 , angles= -angles ; end ;
            c	= sqrt(0.5+angles(1)/2);
            s	= 0.5*(angles(2)-j*angles(3))/c;

            % Update the matrices M and V by a Givens Rotation:
            if abs(s)>seuil,
                encore          = 1 ;
                pair            = [p;q] ;
                G               = [ c -conj(s) ; s c ] ;
                V(:,pair)       = V(:,pair)*G ;
                M(pair,:)       = G' * M(pair,:) ;
                M(:,[Ip Iq]) 	= [ c*M(:,Ip)+s*M(:,Iq) -conj(s)*M(:,Ip)+c*M(:,Iq) ] ;
            end
        end
    end
end

%% Estimation of the mixing matrix and signal separation
A	= IW*V;
S	= V'*Y ;

return ;

author: danmcduff

repo: iphys-toolbox


text:
function varargout = ensureManifoldMesh(varargin)
%ENSUREMANIFOLDMESH Apply several simplification to obtain a manifold mesh.
%
%   Try to transform an input mesh into a manifold mesh.
%
%   Not all cases of "non-manifoldity" are checked, so please use with
%   care.
%
%   [V2, F2] = ensureManifoldMesh(V, F);
%   [V2, F2] = ensureManifoldMesh(MESH);
%   MESH2 = ensureManifoldMesh(...);
%
%   Example
%   ensureManifoldMesh
%
%   See also
%    meshes3d, isManifoldMesh

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2019-02-01,    using Matlab 9.5.0.944444 (R2018b)
% Copyright 2019 INRA - Cepia Software Platform.

%% Parse input arguments

[vertices, faces] = parseMeshData(varargin{:});
verbose = true;

%% Pre-processing

% remove duplicate faces if any
faces = removeDuplicateFaces(faces);

%% Iterative processing of multiple edges
% Reduces all edges connected to more than two faces, by collapsing second
% vertex onto the first one.

% iter = 0;
% while ~isManifoldMesh(vertices, faces) && iter < 10
%     iter = iter + 1;
    [vertices, faces] = collapseEdgesWithManyFaces(vertices, faces);
% end

%% Format output
varargout = formatMeshOutput(nargout, vertices, faces);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef FontColorable < handle
    % Mixin for component with Font color (but no other editable font properties)

    % Copyright 2020-2021 The MathWorks Inc.

    %% Properties
    properties (AbortSet)

        % Font color
        FontColor (1,3) double {wt.validators.mustBeBetweenZeroAndOne} = [0 0 0]

    end %properties

    %% Internal properties
    properties (AbortSet, Transient, NonCopyable, ...
            Access = {?wt.abstract.BaseWidget, ?wt.test.BaseWidgetTest} )

        % List of graphics controls to apply to
        FontColorableComponents (:,1) matlab.graphics.Graphics

    end %properties

    %% Accessors
    methods

        function set.FontColor(obj,value)
            obj.FontColor = value;
            obj.updateFontColorableComponents()
        end

        function set.FontColorableComponents(obj,value)
            obj.FontColorableComponents = value;
            obj.updateFontColorableComponents()
        end

    end %methods

    %% Methods
    methods (Access = protected)

        function updateFontColorableComponents(obj)

                hasFontColorProp = isprop(obj.FontColorableComponents,"FontColor");

                wt.utility.fastSet(obj.FontColorableComponents(hasFontColorProp),...
                    "FontColor",obj.FontColor);

                hasForegroundColorProp = ~hasFontColorProp & ...
                    isprop(obj.FontColorableComponents,"ForegroundColor");

                wt.utility.fastSet(obj.FontColorableComponents(hasForegroundColorProp),...
                    "ForegroundColor",obj.FontColor);

        end %function

    end %methods

end %classdef


text:
function S3G = copy(S3G,ind)
% copy nodes by indece
%
% Input
%  S3G    - @SO3Grid
%  indece - int32
%
% See also
%  EBSD/delete EBSD/subGrid

 S3G.a = S3G.a(ind);
 S3G.b = S3G.b(ind);
 S3G.c = S3G.c(ind);
 S3G.d = S3G.d(ind);
 S3G.i = S3G.i(ind);

 S3G.gamma = subGrid(S3G.gamma,ind);
 S3G.alphabeta  = subGrid(S3G.alphabeta,GridLength(S3G.gamma)>0);
 S3G.gamma(GridLength(S3G.gamma)==0) = [];

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function f = plotkf(kf, plotDir);
%-------------------------------------------------------------------------------
% f = plotkf(kf, isPrintToFile);
% generate plots for Kalman filter estimates
% input:
% 1) kf is the data loaded into memory
% => use kf = readdata(solFileName, 1+9+9) to load binary file format
% => use kf = load(solFileName) to load text file format
% 2) plotDir => to control where to print a tiff file for each plot or not,
% use empty [] to turn off print
% output:
% f => figure handles
% by Dr. Yudan Yi
% modified date on: 04/08/2013
%-------------------------------------------------------------------------------
if (nargin<1) return; end;
[nn,mm] = size(kf);
if (nn<1) return; end;
isPrintToFile = false;
if nargin<2 plotDir = []; end;
if ~isempty(plotDir)
    isPrintToFile = true;
    n = length(plotDir);
    if (plotDir(n)~='\')
        plotDir = [plotDir '\'];
    end
    if ~isdir(plotDir)
        mkdir(plotDir);
    end
end
%-------------------------------------------------------------------------------
max_x = max(abs(kf(:,3)));
max_y = max(abs(kf(:,2)));
isNED = false;
if (max_y>pi||max_x>(2.0*pi))
    isNED = true;
end
f(1) = figure;
plot(kf(:,3),kf(:,2),'marker','.')
grid
axis equal
if (isNED)
    xlabel('East [m]')
    ylabel('North[m]')
else
    xlabel('Lontitude [radian]')
    ylabel('Latitude [radian]')
end
title('Ground Track');
if isPrintToFile
    print(f(1), '-dtiff', [plotDir 'ground_track']);
end
%close all
f(2) = figure;
plot(kf(:,1)-kf(1,1),kf(:,4),'marker','.')
grid
xlabel('Time [s]')
ylabel('m')
title('Height Profile');
if isPrintToFile
    print(f(2), '-dtiff', [plotDir 'ht']);
end
%close all
f(3) = figure;
plot(kf(:,1)-kf(1,1),kf(:,5:7),'marker','.')
grid
xlabel('Time [s]')
ylabel('m/s')
legend('Vx','Vy','Vz')
title('Velocity');
if isPrintToFile
    print(f(3), '-dtiff', [plotDir 'vel']);
end
%close all
% % note: roll is around pi, make to around 0
% kf(:,8) = kf(:,8);
% loc = find(kf(:,8)>pi);
% kf(loc,8) = kf(loc,8);
f(4) = figure;
plot(kf(:,1)-kf(1,1),kf(:,8:10),'marker','.')
grid
xlabel('Time [s]')
ylabel('degree')
legend('Roll','Pitch','Yaw/Heading')
title('Orientation');
if isPrintToFile
    print(f(4), '-dtiff', [plotDir 'att']);
end
%close all
if mm==19
f(5) = figure;
plot(kf(:,1)-kf(1,1),kf(:,11:13),'marker','.')
grid
xlabel('Time [s]')
ylabel('m')
legend('X','Y','Z')
title('Position RMS in ECEF XYZ');
if isPrintToFile
    print(f(5), '-dtiff', [plotDir 'rms_pos']);
end
%close all
f(6) = figure;
plot(kf(:,1)-kf(1,1),kf(:,14:16),'marker','.')
grid
xlabel('Time [s]')
ylabel('m/s')
legend('Vn','Ve','Vd')
title('Velocity RMS in NED');
if isPrintToFile
    print(f(6), '-dtiff', [plotDir 'rms_vel']);
end
%close all
f(7) = figure;
plot(kf(:,1)-kf(1,1),kf(:,17:19)*180/pi,'marker','.')
grid
xlabel('Time [s]')
ylabel('radian')
legend('Roll','Pitch','Yaw/Heading')
title('Orientation RMS in RPY');
if isPrintToFile
    print(f(7), '-dtiff', [plotDir 'rms_att']);
end
end
%close all

repo: ekfmonoslam


text:
function obj = addNodeCost(obj, func, deps, nodes, auxdata, load_path)
    % Add a cost function that only depends on variables at a particular
    % node. The input argument ''nodes'' will specify at which nodes the
    % function is defined.
    %
    % @note This function provides a simple type of cost function which
    % dependends only on variables at some particular node. If multiple
    % nodes are specified, the final cost function will be the sum of the
    % function computated at each node.
    %
    % @attention In the case the cost function dependes on variables at
    % different nodes, please use the basic ''addCost'' method directly.
    %
    % @attention This method can be used to add any terminal cost function
    % by specify ''nodes'' as either ''first'' or ''last''.
    %
    % Parameters:
    % func: a symbolic function to be integrated @type SymFunction
    % deps: a list of dependent variables @type cellstr
    % node: an indicator of 'first' or 'last' node @type char
    % auxdata: auxilary constant data to be feed in the function
    % @type double

    % basic information of NLP decision variables
    nNode  = obj.NumNode;
    vars   = obj.OptVarTable;
    if ~iscell(deps), deps = {deps}; end

    assert(isa(func,'SymFunction'),...
        'The second argument must be a SymFunction object.'); %#ok<PSIZE>

    if nargin < 5
        auxdata = [];
    else
        if ~iscell(auxdata), auxdata = {auxdata}; end
    end

    if nargin < 6
        load_path = [];
    end

    if ~isempty(load_path)
        load(func,load_path);
    end
    siz = size(func);
    assert(prod(siz)==1,...
        'The cost function must be a scalar function.'); %#ok<PSIZE>

    if ischar(nodes)
        switch nodes
            case 'first'
                node_list = 1;
            case 'last'
                node_list = nNode;
            case 'except-first'
                node_list = 2:nNode;
            case 'except-last'
                node_list = 1:nNode-1;
            case 'except-terminal'
                node_list = 2:nNode-1;
            case 'all'
                node_list = 1:nNode;
            case 'cardinal'
                node_list = 1:2:nNode;
            case 'interior'
                node_list = 2:2:nNode-1;
            otherwise
                error('Unknown node type.');
        end
    else
        if ~isnumeric(nodes)
            error(['The node must be specified as a list or following supported characters:\n',...
                '%s'],implode({'first','last','all','except-first','except-last','except-terminal', 'cardinal', 'interior'},','));
        else
            node_list = nodes;
        end
    end

    n_node = numel(node_list);
    cost(n_node) = struct();
    [cost.Name] = deal(func.Name);
    [cost.Dimension] = deal(1);
    [cost.SymFun] = deal(func);
    [cost.AuxData] = deal(auxdata);
    for i=1:n_node
        idx = node_list(i);
        dep_vars = cellfun(@(x)vars.(x)(idx),deps,'UniformOutput',false);
        cost(i).DepVariables = vertcat(dep_vars{:});
    end

    obj = addCost(obj,func.Name,nodes,cost);
end

author: ayonga

repo: frost-dev

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


text:
function [P,F] = spm_MH(L,B,y,M)
% The Rejection-Metropolis-Hastings Algorithm
% FORMAT [P,F] = spm_MH(L,P,y)
%
% L   - likelihood function: inline(P,y,M)
% B   - free parameter [structure]
% Y   - response  [stucture]
% M   - model [structure]
%
% P   - Sample from posterior p(P|y,M)
% F   - marginal likelihood p(y|M) using harmonic mean
%--------------------------------------------------------------------------
%
% Returns a harmonic mean estimate of the log-marginal likelihood or
% log-evidence and a sample from the posterior density of the free parameters
% of a model.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_MH.m 1143 2008-02-07 19:33:33Z spm $

% initialise parameters
%--------------------------------------------------------------------------
P(:,1) = spm_vec(B);

% MCMC - RMH
%--------------------------------------------------------------------------
n     = 2^8;                                          % number of burn in
N     = 2^16;                                         % number of samples
for i = 1:N

    % sample from proposal
    %----------------------------------------------------------------------
    pi = P(:,i);
    pp = pi + randn(size(P,1),1)/32;

    % compute importnace ratio
    %----------------------------------------------------------------------
    Lp = feval(L,spm_unvec(pp,B),y,M);
    Li = feval(L,spm_unvec(pi,B),y,M);
    r  = Lp/Li;

    % accept and marginal likelihood
    %----------------------------------------------------------------------
    if rand < r
        P(:,i + 1) = pp;
        F(i + 1)   = Lp;
    else
        P(:,i + 1) = pi;
        F(i + 1)   = Li;
    end

end

% burn in
%--------------------------------------------------------------------------
P = P(:,n:end);
F = F(n:end);
F = -log(mean(mean(F)./F)) + log(mean(F));


text:
% This script runs program rundivcrl
% for nine coordinate systems and puts
% the output in rundivcrl.tst
if exist('rundivcrl.tst','file')==2
  delete rundivcrl.tst
end
syms x y z
vxyz=[x*y,y*z,z*x];
diary rundivcrl.tst
% rundivcrl(vxyz,cordname,tn,noprint) call list
rundivcrl(vxyz,@cylin,[10*rand,rand*2*pi,rand*10]);
rundivcrl(vxyz,@sphr,[10*rand,rand*pi,rand*2*pi]);
rundivcrl(vxyz,@cone,[2*rand,rand*pi/2,rand*2*pi]);
rundivcrl(vxyz,@elipcyl,[2*rand,rand*2*pi,2*rand]);
rundivcrl(vxyz,@elipsod,[rand,1+rand,2+rand]);
rundivcrl(vxyz,@notort,[2*rand,rand*pi,rand*2*pi]);
rundivcrl(vxyz,@oblate,[2*rand,rand*pi,rand*2*pi]);
rundivcrl(vxyz,@parab,[2*rand,2*rand,2*rand]);
rundivcrl(vxyz,@toroid,[rand*2,rand*2*pi,rand*2*pi]);
diary off

author: Sable

repo: mcbench-benchmarks


text:
%IMM_PREDICT  Interacting Multiple Model (IMM) Filter prediction step
%
% Syntax:
%   [X_p,P_p,c_j,X,P] = IMM_PREDICT(X_ip,P_ip,MU_ip,p_ij,ind,dims,A,Q)
%
% In:
%   X_ip  - Cell array containing N^j x 1 mean state estimate vector for
%           each model j after update step of previous time step
%   P_ip  - Cell array containing N^j x N^j state covariance matrix for
%           each model j after update step of previous time step
%   MU_ip - Vector containing the model probabilities at previous time step
%   p_ij  - Model transition probability matrix
%   ind   - Indexes of state components for each model as a cell array
%   dims  - Total number of different state components in the combined system
%   A     - State transition matrices for each model as a cell array.
%   Q     - Process noise matrices for each model as a cell array.
%
% Out:
%   X_p   - Predicted state mean for each model as a cell array
%   P_p   - Predicted state covariance for each model as a cell array
%   c_j   - Normalizing factors for mixing probabilities
%   X     - Combined predicted state mean estimate
%   P     - Combined predicted state covariance estimate
%
% Description:
%   IMM filter prediction step.
%
% See also:
%   IMM_UPDATE, IMM_SMOOTH, IMM_FILTER

% History:
%   01.11.2007 JH The first official version.
%
% Copyright (C) 2007 Jouni Hartikainen
%
% $Id: imm_update.m 111 2007-11-01 12:09:23Z jmjharti $
%
% This software is distributed under the GNU General Public
% Licence (version 2 or later); please refer to the file
% Licence.txt, included with the software, for details.

function [X_p,P_p,c_j,X,P] = imm_predict(X_ip,P_ip,MU_ip,p_ij,ind,dims,A,Q)
    % Number of models
    m = length(X_ip);

    % Default values for state mean and covariance
    MM_def = zeros(dims,1);
    PP_def = diag(20*ones(dims,1));

    % Normalizing factors for mixing probabilities
    c_j = zeros(1,m);
    for j = 1:m
        for i = 1:m
            c_j(j) = c_j(j) + p_ij(i,j).*MU_ip(i);
        end
    end

    % Mixing probabilities
    MU_ij = zeros(m,m);
    for i = 1:m
        for j = 1:m
            MU_ij(i,j) = p_ij(i,j) * MU_ip(i) / c_j(j);
        end
    end

    % Calculate the mixed state mean for each filter
    X_0j = cell(1,m);
    for j = 1:m
        X_0j{j} = zeros(dims,1);
        for i = 1:m
            X_0j{j}(ind{i}) = X_0j{j}(ind{i}) + X_ip{i}*MU_ij(i,j);
        end
    end

    % Calculate the mixed state covariance for each filter
    P_0j = cell(1,m);
    for j = 1:m
        P_0j{j} = zeros(dims,dims);
        for i = 1:m
            P_0j{j}(ind{i},ind{i}) = P_0j{j}(ind{i},ind{i}) + MU_ij(i,j)*(P_ip{i} + (X_ip{i}-X_0j{j}(ind{i}))*(X_ip{i}-X_0j{j}(ind{i}))');
        end
    end

    % Space for predictions
    X_p = cell(1,m);
    P_p = cell(1,m);

    % Make predictions for each model
    for i = 1:m
        [X_p{i}, P_p{i}] = kf_predict(X_0j{i}(ind{i}),P_0j{i}(ind{i},ind{i}),A{i},Q{i});
    end

    % Output the combined predicted state mean and covariance, if wanted.
    if nargout > 3
        % Space for estimates
        X = zeros(dims,1);
        P = zeros(dims,dims);

        % Predicted state mean
        for i = 1:m
            X(ind{i}) = X(ind{i}) + MU_ip(i)*X_p{i};
        end

        % Predicted state covariance
        for i = 1:m
            P(ind{i},ind{i}) = P(ind{i},ind{i}) + MU_ip(i)*(P_p{i} + (X_ip{i}-X(ind{i}))*(X_ip{i}-X(ind{i}))');
        end
    end

author: EEA-sensors

repo: ekfukf

save_path: github-repos/MATLAB/EEA-sensors-ekfukf


text:
function X = hlp_insertSingletonDim(X,dim)
% insert a singleton dimension into dimension 'dim' of matrix X
% author: Tim Mullen, 2011

    nd = ndims(X);
    X = permute(X,unique_bc([1:dim-1 nd+1 dim:nd]));

repo: REST


text:
function mu_all = calc_mu_all(mu_jk, K_all)
% mu_all is B by M by K1
[~,B] = size(mu_jk{1});
[K1,M] = size(K_all);

mu_all = zeros(B,M,K1);
for i = 1:K1
    for j = 1:M
        mu_all(:,j,i) = mu_jk{j}(K_all(i,j),:)';
    end
end

end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function robot_disp = LoadRobotDisplay(robot, varargin)

    ip = inputParser;
    ip.addParameter('UseExported',true,@(x) isequal(x,true) || isequal(x,false));
    ip.addParameter('ExportPath','',@(x)ischar(x));
    ip.addParameter('SkipExporting',false,@(x) isequal(x,true) || isequal(x,false));

    ip.parse(varargin{:});

    opts = ip.Results;
    if isempty(opts.ExportPath)
        root_path = utils.get_root_path();
        export_path = fullfile(root_path,'gen','animator');
        opts.ExportPath = export_path;
    else
        export_path = opts.ExportPath;
    end
    if ~exist(export_path,'dir')
        mkdir(export_path);
    end
    addpath(export_path);

    f = figure(1000);clf;
    robot_disp = frost.Animator.Display(f, robot, opts);frame  = robot.Joints(6);

    pelvis = robot.Links(getLinkIndices(robot,'mtorso'));
    offset = pelvis.Offset;
    torso_bar = frost.Animator.Cylinder(robot_disp.axs, robot, frame, offset, 'TorsoBar', opts);
    robot_disp.addItem(torso_bar);

    torso = frost.Animator.LinkSphere(robot_disp.axs, robot, pelvis, 'TorsoCoM', opts);
    torso.radius = 0.1;
    robot_disp.addItem(torso);
    % add feet
    left_foot = sys.frames.LeftSole(robot);
    param = sys.GetExtraParams;
    name = 'LeftFoot';
    top_offset = [0,0,-param.hf];
    wf = param.wf;
    lt = param.lt;
    lh = param.lh;
    bottom_offset = [lt,-wf/2,0;
        lt,wf/2,0;
        -lh,wf/2,0;
        -lh,-wf/2,0
        ];
    left_foot = frost.Animator.Pyramid(robot_disp.axs, robot, left_foot, top_offset, bottom_offset, name, opts);
    robot_disp.addItem(left_foot);

    right_foot = sys.frames.RightSole(robot);
    param = sys.GetExtraParams;
    name = 'RightFoot';
    top_offset = [0,0,-param.hf];
    wf = param.wf;
    lt = param.lt;
    lh = param.lh;
    bottom_offset = [lt,-wf/2,0;
        lt,wf/2,0;
        -lh,wf/2,0;
        -lh,-wf/2,0
        ];
    right_foot = frost.Animator.Pyramid(robot_disp.axs, robot, right_foot, top_offset, bottom_offset, name, opts);
    robot_disp.addItem(right_foot);

end


text:
function Xi = NonlinearSparseRegression(x,u,dt,options_method,lambda)

Ninputs = min(size(u));
%% Compute Derivative
% compute derivative using fourth order central difference
% use TVRegDiff if more error
dx = zeros(length(x)-5,3);
for i=3:length(x)-3
    for k=1:size(x,2)
        dx(i-2,k) = (1/(12*dt))*(-x(i+2,k)+8*x(i+1,k)-8*x(i-1,k)+x(i-2,k));
    end
end
% concatenate
xaug = [x(3:end-3,:) u(3:end-3,:)];
dx(:,size(x,2)+1:size(x,2)+Ninputs) = repmat(0*dx(:,size(x,2)),[1,Ninputs]);

n = size(dx,2);

%% Sparse regression
clear Theta Xi
Theta = poolData(xaug,n,options_method.order,options_method.usesine);
Xi = sparsifyDynamics(Theta,dx,lambda,n);

author: eurika-kaiser

repo: SINDY-MPC


text:
function  [yz, dy] = pzextr1(iest, xest, yest)
% pzextr.m  uses extrapolation to evaluate nv functions at x=0 by fitting a
% polynomial to a sequence of estimates with progressively smaller values x= xest(1..np), and
% corresponding function vectors yest(1..nv,1..np). This call is number iest in the sequence
% of calls, see also bsstep1.m and pzextr.m
%
% output : the extrapolated function values, yz(1..nv,1..np) and their estimated error, dy(1..nv,1..np).

% D Vangheluwe 8 mrt 2005

global x_bulirsch_stoer  d_bulirsch_stoer

[nv np] = size(yest);

if np > 1
  x_bulirsch_stoer(iest,:) = xest;
% save current independent variable
  yz = yest;
  dy = yz;

%  x_bulirsch_stoer
%  d_bulirsch_stoer

  if iest == 1
% store first estimate in first column
     d_bulirsch_stoer(:,:,1) = yest;
  else
     c = yest;
     for k1 = 1:iest-1
        delta = 1 ./ (x_bulirsch_stoer(iest-k1,:) - xest);
        f1 = xest .* delta;
        f2 = x_bulirsch_stoer(iest-k1,:) .* delta;
        q = d_bulirsch_stoer(:,:,k1);
        d_bulirsch_stoer(:,:,k1) = dy;
        deltay = c - q;
        dy = repmat(f1, nv, 1) .* deltay;
        c = repmat(f2, nv, 1) .* deltay;
        yz = yz + dy;
     end
     d_bulirsch_stoer(:,:,iest) = dy;
  end

else  %np == 1

  x_bulirsch_stoer(iest) = xest';
% save current independent variable
  yz = yest;
  dy = yz;

%  x_bulirsch_stoer
%  d_bulirsch_stoer
  if iest == 1
% store first estimate in first column
     d_bulirsch_stoer(:,1) = yest;
  else
     c = yest;
     for k1 = 1:iest-1
        delta = 1 ./ (x_bulirsch_stoer(iest-k1) - xest);
        f1 = xest .* delta;
        f2 = x_bulirsch_stoer(iest-k1) .* delta;
        q = d_bulirsch_stoer(:,k1);
        d_bulirsch_stoer(:,k1) = dy;
        deltay = c - q;
        dy = f1 .* deltay;
        c = f2 .* deltay;
        yz = yz + dy;
     end
     d_bulirsch_stoer(:,iest) = dy;
  end

end  %np > 1


text:
function [windows]=mvg_makeSuperpixelsWindows(superPix,maxNumConnectedSuperpix)

%% Default settings
if nargin<2
    maxNumConnectedSuperpix=3; % Defaults to triplets
end

%% Make sure superpxels are in labeled format
if size(superPix,3)>1
    superPix=mvg_numerizeLabels(superPix);
end
superPix=double(superPix);

%% Make single superpixel windows
superpixLabels=unique(superPix(:));
windows=zeros(length(superpixLabels),4);
for i=1:length(superpixLabels)
    [rw,cl]=find(superPix==superpixLabels(i));
    windows(i,:)=[min(cl),min(rw),max(cl),max(rw)];
end

%% Make tuplets if needed
if maxNumConnectedSuperpix>1

    %% Find neighboring superpixels (4-neighbors)
    neighborLabels=zeros(size(superPix,1),size(superPix,2),4);
    neighborLabels(:,:,1)=[superPix(:,1),superPix(:,1:end-1)]; % left
    neighborLabels(:,:,2)=[superPix(:,2:end),superPix(:,end)]; % right
    neighborLabels(:,:,3)=[superPix(1,:);superPix(1:end-1,:)]; % top
    neighborLabels(:,:,4)=[superPix(2:end,:);superPix(end,:)]; % bottom

    %% Find superpixel border pixels
    borderPairs=[neighborLabels(:),repmat(superPix(:),[4,1])];
    superPixNeighbors=unique(sort(borderPairs,2),'rows');

    %% Remove single superpixels
    superPixNeighbors(superPixNeighbors(:,1)==superPixNeighbors(:,2),:)=[];

    %% Make bounding boxes for all pairs
    numWindowsTuplets=size(superPixNeighbors,1);
    windowsAdd=zeros(numWindowsTuplets,4);
    for i=1:numWindowsTuplets
        [rw,cl]=find(superPix==superPixNeighbors(i,1) | superPix==superPixNeighbors(i,2));
        windowsAdd(i,:)=[min(cl),min(rw),max(cl),max(rw)];
    end

    %% Remove multiple windows
    windowsAdd=unique(windowsAdd,'rows');

    %% Add windows to existing ones
    windows=[windows;windowsAdd];

end

%% Add more connected components if needed
if maxNumConnectedSuperpix>2
    % Make superpixel pair matrix
    superpixelPairsMatrix=zeros(length(superpixLabels),length(superpixLabels));
    superpixelPairsMatrix(superPixNeighbors(:,1)+(superPixNeighbors(:,2)-1)*length(superpixLabels))=1;
    superpixelPairsMatrix=superpixelPairsMatrix+superpixelPairsMatrix';

    % Initialize
    superpixComb=superPixNeighbors;

    % Continue to add windows
    for i=3:maxNumConnectedSuperpix
        % Extend combination by one
        superpixComb=extendSupepixComb_(superpixelPairsMatrix,superpixComb);

        % Make new windows
        windowsAdd=zeros(size(superpixComb,1),4);
        for j=1:size(superpixComb,1)
            [rw,cl]=find(ismember(superPix,superpixComb(j,:)));
            windowsAdd(j,:)=[min(cl),min(rw),max(cl),max(rw)];
        end

        % Remove multiple windows
        windowsAdd=unique(windowsAdd,'rows');

        % Add windows to existing ones
        windows=[windows;windowsAdd];

    end

end

%% Remove any dublicate windows
windows=unique(windows,'rows');

%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Additional functions %%
%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Add new superpixel combination to the existing ones %%%
function [superpixExtComb]=extendSupepixComb_(superpixelPairsMatrix,superpixComb)

%% Initialize
superpixExtComb=[];

%% Loop over existing combination and extend
for i=1:size(superpixComb,1)
    % Find connected superpixels
    connSuperpix=superpixelPairsMatrix(:,superpixComb(i,:));
    connSuperpix=find(sum(connSuperpix,2)>eps);
    connSuperpix=setdiff(connSuperpix,superpixComb(i,:));

    % Add connected components to
    if ~isempty(connSuperpix)
        superpixExtComb=[superpixExtComb; [repmat(superpixComb(i,:),[length(connSuperpix),1]), connSuperpix(:)]];
    end

end

%% Remove possible duplicates
superpixExtComb=unique(sort(superpixExtComb,2),'rows');

%%% Numerize superpixel labels %%%
function N = numerizeSuperpixLabels_(S)

% converts the segmentation image S from true-color
% to ordered integer labels
%

N = zeros(size(S,1),size(S,2),'uint16');
col2n = zeros(256,256,256);
totcol = 0;
for x = 1:size(S,2)
    for y = 1:size(S,1)
        p = reshape(S(y,x,:),1,3)+1;
        cix = col2n(p(1),p(2),p(3));
        if cix > 0
          N(y,x) = cix;
        else
          totcol = totcol+1;
          col2n(p(1),p(2),p(3)) = totcol;
          N(y,x) = totcol;
        end
    end
end

% MVG version (faster)
% S=double(S);
% temp=S(:,:,1)+S(:,:,2)*256+S(:,:,3)*256^2;
% lbl=unique(temp(:));
% %lblCode=[lbl,(1:length(lbl))'];
% N=zeros(size(S,1),size(S,2),'uint16');
% for i=1:length(lbl)
%     N(temp==lbl(i))=i;
% end

author: Cloud-CV

repo: object-proposals


text:
function analyse_eval_lin_and_qual_fusion(key_dev,ndx_eval,dev_scores_obj_array,eval_scores_obj_array,dev_scr_lin,eval_scr_lin,dev_scr_qual,eval_scr_qual)
% Does some analyses on evaluation scores for all systems and all
% fusions.  For each system, it calculates the estimated percentage
% of targets in the database.  It also calculates a KL divergence
% matrix of all systems against all systems.
% Inputs:
%   key_dev: The Key for the dev trials.
%   ndx_eval: The Ndx describing the eval trials.
%   dev_scores_obj_array: An array of Scores objects containing
%     system scores for the dev trials.
%   eval_scores_obj_array: An array of Scores objects containing
%     system scores for the eval trials.
%   dev_scr_lin: The dev scores resulting from the linear fusion.
%   eval_scr_lin: The eval scores resulting from the linear fusion.
%   dev_scr_qual: The dev scores resulting from the quality fusion.
%   eval_scr_qual: The eval scores resulting from the quality fusion.

assert(nargin==8)
assert(isa(key_dev,'Key'))
assert(isa(ndx_eval,'Ndx'))
assert(isa(dev_scr_lin,'Scores'))
assert(isa(eval_scr_lin,'Scores'))
assert(isa(dev_scr_qual,'Scores'))
assert(isa(eval_scr_qual,'Scores'))
assert(key_dev.validate())
assert(ndx_eval.validate())
assert(dev_scr_lin.validate())
assert(eval_scr_lin.validate())
assert(dev_scr_qual.validate())
assert(eval_scr_qual.validate())

eval_stacked_scores = stackScores(ndx_eval,eval_scores_obj_array);
trialmask = ndx_eval.trialmask;
eval_lin_fusion = eval_scr_lin.scoremat(trialmask(:))';
eval_qfusion = eval_scr_qual.scoremat(trialmask(:))';
eval_scores = [eval_stacked_scores;eval_lin_fusion;eval_qfusion];

ptar = estimate_target_proportions(eval_scores,0.08);
logprint(Logger.Info,'\nSystem Ptar estimates:\n');
logprint(Logger.Info,'%s',sprintfmatrix(100*ptar'));

ptar_best = ptar(end-1);

dev_stacked_scores = stackScores(key_dev,dev_scores_obj_array);
trialmask = key_dev.to_ndx().trialmask;
dev_lin_fusion = dev_scr_lin.scoremat(trialmask(:))';
dev_qfusion = dev_scr_qual.scoremat(trialmask(:))';
dev_scores = [dev_stacked_scores;dev_lin_fusion;dev_qfusion];

devKL = 100*kldivergence_of_systems(dev_scores,ptar_best,key_dev);
evalKL = 100*kldivergence_of_systems(eval_scores,ptar_best);
logprint(Logger.Info,'\nDev score divergence:\n%s',sprintfmatrix(devKL));
logprint(Logger.Info,'\nEval score divergence:\n%s',sprintfmatrix(evalKL));
end

author: nesl


text:
function [wnew,freqnew]=updateGlyph(P,X,C,w,a,s,freq)
%Syntax: [wnew,freqnew]=updateGlyph(P,X,C,w,a,s,freq)
%____________________________________________________
%
% Weight matrix update of a Spherical Self-Organized Feature Map.
%
% wnew is the updated weights' matrix.
% freqnew is the updated count-dependent parameter.
% P is the A-by-B matrix with the patterns to be classified.
% X is the N-by-3 matrix with the cartesian coordinates of the points on
%   the sphere.
% C is a cell array with the neighbors. The {i,j} cell represents the
%   neighbors of radius i of the j-th sphere point.
% w is the weights matrix to be updated.
% a is the training parameter.
% s is the neighborhood size parameter.
% freq is the count-dependent parameter.
%
%
% References:
%
% Sangole A., Knopf G. K. (2002): Representing high-dimensional data sets
% as close surfaces. Journal of Information Visualization 1: 111-119
%
% Sangole A., Knopf G. K. (2003): Geometric representations for
% high-dimensional data using a spherical SOFM. International Journal of
% Smart Engineering System Design 5: 11-20
%
% Sangole A., Knopf G. K. (2003): Visualization of random ordered numeric
% data sets using self-organized feature maps. Computers and Graphics 27:
% 963-976
%
% Sangole A. P. (2003): Data-driven Modeling using Spherical
% Self-organizing Feature Maps. Doctor of Philosophy (Ph.D.) Thesis.
% Department of Mechanical and Materials Engineering. Faculty of
% Engineering. The University of Western Ontario, London, Ontario, Canada.
%
%
% Archana P. Sangole, PhD., P.E. (TX chapter)
% School of Physical & Occupational Therapy
% McGill University
% 3654 Promenade Sir-William-Osler
% Montreal, PQ, H3G 1Y5
% e-mail: archana.sangole@mail.mcgill.ca
%
% CRIR, Rehabilitation Institute of Montreal
% 6300 Ave Darlington
% Montreal, PQ, H3S 2J5
% Tel: 514.340.2111 x2188
% Fax: 514.340.2154
%
%
% Alexandros Leontitsis, PhD
% Department of Education
% University of Ioannina
% 45110- Dourouti
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% 23-Mar-2006

if nargin<2 | isempty(P)==1 | isempty(X)==1
    error('More input arguments needed. Both P and X are rquired.');
else
    % P must be an A-by-B matrix
    if ndims(P)~=2
        error('P must be an A-by-B matrix.');
    end
    % X must be an M-by-3 matrix
    if size(X,2)~=3 | ndims(X)~=2
        error('X must be an N-by-3 matrix.');
    end
end

if isempty(C)==0
    % C and X must be of the same length
    if length(C)~=length(X)
        error('C and X must be of the same length.');
    end
    % Set the neighborhood radious
    r=size(C,2);
    % Set the h function
    h=exp(-(1:r).^2/(s*r));
end

if nargin<4 | isempty(w)==1
    w=randn(size(X,1),size(P,2));
    for i=1:size(X,2)
        w(:,i)=w(:,i)*std(P(:,i));
    end
else
    % w must be an N-by-B matrix
    if size(w,1)~=size(X,1) | size(w,2)~=size(P,2) | ndims(X)~=2
        error('w must be an N-by-B matrix.');
    end
end

if nargin<5 | isempty(a)==1
    a=0.1;
else
    % a must be a scalar
    if sum(size(a))>2
        error('a must be a scalar.');
    end
end

if nargin<6 | isempty(s)==1
    s=2;
else
    % s must be a scalar
    if sum(size(s))>2
        error('s must be a scalar.');
    end
    % s must be positive
    if s<=0
        error('s must be positive.');
    end
end

if nargin<7 | isempty(freq)==1
    freq=zeros(length(X),1);
else
    % freq must be a vector
    if min(size(freq))>1
        error('freq must be a vector.');
    end
    % freq must have the same length with X
    if length(freq)~=length(X)
        error('freq must have the same length with X.');
    end
end

% Initialize the weights
wnew=w;

% Initialize freqnew
freqnew=freq;

% Randomize the order of the patterns
q=randperm(size(P,1));

% For each uniformly random pattern
for j=1:size(P,1)
    % Measure the difference between the current P and all the weights
    A=q(j)*ones(size(wnew,1),1);
    d=P(A,:)-wnew;
    % Calculate the norm and apply the count-dependent parameter
    d=sqrt(sum((d.^2)'))'.*(freqnew+1);
    % Obtain the minimum distance and the corresponding index
    [D,I]=min(d);
    % Calculate the new weights
    wnew(I,:)=wnew(I,:)+a*(P(q(j),:)-wnew(I,:));
    % If the neighborhood matrix is given
    if isempty(C)==0
        % For each neghborhood radius
        for k=1:r
            A=q(j)*ones(size(C{I,k}));
            % Calculate the new weights of the neighbors
            wnew(C{I,k},:)=wnew(C{I,k},:)+a*h(k)*(P(A,:)-wnew(C{I,k},:));
            % Update the count-dependent parameter for the neighbors
            freqnew(C{I,k})=freqnew(C{I,k})+h(k);
        end
    end
    % Update the count-dependent parameter
    freqnew(I)=freqnew(I)+1;
end


text:
function G=constructMatrixOfMonomials(g,order)
%-fanDTasia ToolBox------------------------------------------------------------------
% This Matlab script is part of the fanDTasia ToolBox: a Matlab library for Diffusion
% Weighted MRI (DW-MRI) Processing, Diffusion Tensor (DTI) Estimation, High-order
% Diffusion Tensor Analysis, Tensor ODF estimation, Visualization and more.
%
% A Matlab Tutorial on DW-MRI can be found in:
% http://www.cise.ufl.edu/~abarmpou/lab/fanDTasia/tutorial.php
%
%-CITATION---------------------------------------------------------------------------
% If you use this software please cite the following work:
% A. Barmpoutis and B.C. Vemuri, "A Unified Framework for Estimating Diffusion Tensors
% of any order with Symmetric Positive-Definite Constraints",
% In the Proceedings of ISBI, 2010
%
%-DESCRIPTION------------------------------------------------------------------------
% This function computes all possible monomials in 3 variables of a certain order.
% The 3 variables are given in the form of 3-dimensional vectors.
%
%-USE--------------------------------------------------------------------------------
% G=constructMatrixOfMonomials(g,order);
%
% g: is a list of N 3-dimensional vectors stacked in a matrix of size Nx3
% order: is the order of the computed monomials
% G: contains the computed monomials. It is a matrix of size N x (2+order)!/2(order)!
%
%-DISCLAIMER-------------------------------------------------------------------------
% You can use this source code for non commercial research and educational purposes
% only without licensing fees and is provided without guarantee or warrantee expressed
% or implied. You cannot repost this file without prior written permission from the
% authors. If you use this software please cite the following work:
% A. Barmpoutis and B.C. Vemuri, "A Unified Framework for Estimating Diffusion Tensors
% of any order with Symmetric Positive-Definite Constraints", In Proc. of ISBI, 2010.
%
%-AUTHOR-----------------------------------------------------------------------------
% Angelos Barmpoutis, PhD
% Computer and Information Science and Engineering Department
% University of Florida, Gainesville, FL 32611, USA
% abarmpou at cise dot ufl dot edu
%------------------------------------------------------------------------------------

fprintf(1,'Constructing matrix of monomials G...');
for k=1:length(g)
    c=1;
    for i=0:order
		for j=0:order-i
			G(k,c)=(g(k,1)^i)*(g(k,2)^j)*(g(k,3)^(order-i-j));
			c=c+1;
        end
    end
end
fprintf(1,'Done\n');


text:
function geometry_test0366 ( )

%*****************************************************************************80
%
%% TEST0366 tests SEGMENT_POINT_DIST_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
  test_num = 3;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0366\n' );
  fprintf ( 1, '  SEGMENT_POINT_DIST_3D computes the distance\n' );
  fprintf ( 1, '    between a line segment and point in 3D.\n' );

  for test = 1 : test_num

    [ p1(1:3,1), seed ] = r8vec_uniform_01 ( 3, seed );
    [ p2(1:3,1), seed ] = r8vec_uniform_01 ( 3, seed );
    [ p(1:3,1), seed ] = r8vec_uniform_01 ( 3, seed );

    dist = segment_point_dist_3d ( p1, p2, p );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  TEST = %d', test );
    fprintf ( 1, '  P1 =   %12f  %12f  %12f\n', p1(1:3,1) );
    fprintf ( 1, '  P2 =   %12f  %12f  %12f\n', p2(1:3,1) );
    fprintf ( 1, '  P =    %12f  %12f  %12f\n', p(1:3,1) );
    fprintf ( 1, '  DIST = %12f\n', dist );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Script demonstrating usage of the bpdndl function.
%
% Author: Brendt Wohlberg <brendt@lanl.gov>  Modified: 2015-07-30
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'Copyright' and 'License' files
% distributed with the library.

% Training images
S0 = zeros(512, 512, 5);
S0(:,:,1) = single(stdimage('lena.grey')) / 255;
S0(:,:,2) = single(stdimage('barbara.grey')) / 255;
S0(:,:,3) = single(stdimage('kiel.grey')) / 255;
S0(:,:,4) = single(rgb2gray(stdimage('mandrill'))) / 255;
tmp = single(stdimage('man.grey')) / 255;
S0(:,:,5) = tmp(101:612, 101:612);

%Reduce images size to speed up demo script
tmp = zeros(256, 256, 5);
for k = 1:size(S0,3),
  tmp(:,:,k) = imresize(S0(:,:,k), 0.5);
end
S0 = tmp;

% Extract all 8x8 image blocks, reshape, and subtract block means
SB = imageblocks(S0, [8 8]);
SB = reshape(SB, size(SB,1)*size(SB,2), size(SB,3));
S = bsxfun(@minus, SB, mean(SB, 1));

% Construct initial dictionary
D0 = randn(size(S,1), 64);

% Set up bpdndl parameters
lambda = 0.2;
opt = [];
opt.Verbose = 1;
opt.MaxMainIter = 1000;
opt.rho = 50*lambda + 0.5;
opt.sigma = size(S,2)/200;
opt.AutoRho = 1;
opt.AutoRhoPeriod = 10;
opt.RhoRsdlRatio = 2;
opt.RhoScaling = 5;
opt.AutoRhoScaling = 1;
opt.AutoSigma = 1;
opt.AutoSigmaPeriod = 10;
opt.XRelaxParam = 1.8;
opt.DRelaxParam = 1.8;

% Do dictionary learning
[D, X, optinf] = bpdndl(D0, S, lambda, opt);

% Display learned dictionary
figure;
imdisp(tiledict(D, [8 8]));

% Plot functional value evolution
figure;
plot(optinf.itstat(:,2));
xlabel('Iterations');
ylabel('Functional value');

author: thfylsty


text:
function popDist = totaldistance(pop,dis)
% TOTALDISTANCE
% popDist = TOTALDISTANCE(pop, dis) calculate total distance of pop(routes)
% with the distance matrix dis. Evaluate Each Population Member (Calculate
% Total Distance)

[popSize, numberofcities] = size(pop);
for i = 1:popSize
    d = dis(pop(i,end),pop(i,1)); % Closed Path
    for k = 2:numberofcities
        d = d + dis(pop(i,k-1),pop(i,k));
    end
    popDist(i) = d;
end

author: HuangCongQing


text:
function [ a, b] = p02_interval ( )

%*****************************************************************************80
%
%% P02_INTERVAL returns a starting interval for optimization for problem 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2002
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, two points defining an interval in which
%    the local minimizer should be sought.
%
  a = 0.0;
  b = 1.0;

  return
end


text:
function x = rosser1_null_right ( )

%*****************************************************************************80
%
%% ROSSER1_NULL_RIGHT returns a right null vector of the ROSSER1 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real X(8,1), the null vector.
%
  x = [ ...
    1.0; ...
    2.0; ...
   -2.0; ...
   -1.0; ...
   14.0; ...
   14.0; ...
    7.0; ...
    7.0 ];

  return
end


text:
function cvx_expert_check( fname, varargin )

if ~isempty( varargin ) && ~any(cellfun('isclass',varargin,'cvx')), return; end
global cvx___
if cvx___.expert, return; end

url = [ 'file:///', cvx___.where, cvx___.fs, 'doc', cvx___.fs, 'advanced.html#the-successive-approximation-method' ];
fprintf( 1, [ 'CVX Warning:\n', ...
'   Models involving "%s" or other functions in the log, exp, and entropy\n', ...
'   family are solved using an experimental successive approximation method.\n', ...
'   This method is slower and less reliable than the method CVX employs for\n', ...
'   other models. Please see the section of the user''s guide entitled\n', ...
'       The successive approximation method\n', ...
'   for more details about the approach, and for instructions on how to\n', ...
'   suppress this warning message in the future.\n' ], fname, url );

cvx___.expert = true;

% Note that we do *not* call cvx_save_prefs here. We only save the
% preferences if an explicit setting of "cvx_expert true" is made.

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function s = stdmv(X, dim);
% handles missing values (NaN).

error(nargchk(1,2,nargin));

if nargin < 2
  dim = min(find(size(X)~=1));
  if isempty(dim), dim = 1; end
end

sz = ones(length(size(X)), 1);
sz(dim) = size(X,dim);

% NaNs to zero
Imv = isnan(X);
X0 = X - repmat(meanmv(X,dim), sz);
X0(Imv) = 0;

% Sum and divide by the number of non-missing values
s = sum(X0.^2, dim) ./ sum(~Imv, dim);
%mu(isinf(xmean)) = NaN;

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function projmat_art = comb_projmat(proj_mat1,proj_mat2)%,angles1,angles2)
projmat_art = [proj_mat1,fliplr(proj_mat2)];
end


text:
function toolboxSfmCompile(doSba,gccVer)
% Compiles all the private routines
%
% This function can recompile all that you need for the toolbox to run.
% That includes:
%  - the SBA static library that is needed to build the ...
%  - ... SBA mex file
%  - 3 normal mex files
% By default, the toolbox comes with precompiled binaries for windows 32
% It is easy to compile on Linux and I am developing on Linux so it works
% for sure there.
%
% USAGE
%  toolboxCompile()
%
% INPUTS
%  doSba      - [false] if true, it will compile the dll projection files
%               for SBA and that requires matlab to be launched from the
%               visual c++ terminal on Windows
%  gccVer     - on Linux, specify your gcc version here (e.g. 4.2 if
%               gcc-4.2 is a supported compiler on Matlab), which can be
%               useful if matlab thinks your gcc is too old/new
%
% OUTPUTS
%
% EXAMPLE
%
% See also
%
% Vincent's Structure From Motion Toolbox      Version 3.1
% Copyright (C) 2008-2011 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the GPL [see external/gpl.txt]

if nargin==0; doSba=false; end
if nargin<=1; gccVer=num2str(4.2); else gccVer=num2str(gccVer); end

if strcmp(computer,'GLNX86') || strcmp(computer,'GLNXA64')
  gcc_extra = [ 'CXX=g++-' gccVer ' CC=g++-' gccVer ' LD=g++-' gccVer ];
end

disp('Compiling.......................................');
disp('sba');
savepwd=pwd; cd(fileparts(mfilename('fullpath'))); cd('../');

dir_external = [ pwd '/external/' ];
dir_sba = [ dir_external 'sba/'];

% build the SBA mex file
cd(dir_sba);
% delete previous object files
fileList = {'sba_levmar.o', 'sba_levmar_wrap.o', 'sba_lapack.o', ...
  'sba_crsm.o', 'sba_chkjac.o', 'libsba.a', 'sba.lib'};
for i=1:length(fileList)
  if exist([dir_sba, '/', fileList{i}], 'file'), delete(fileList{i}); end
end

sba_files = [ 'sba.c ../sba_chkjac.c ../sba_levmar_wrap.c ../sba_levmar.c ' ...
  '../sba_lapack.c ../sba_crsm.c' ];

cd matlab
switch computer
  case {'PCWIN','PCWIN64'},
    ext = mexext;
    ext = ext(5:6);
    library_flags = [ '-D_WIN32 -I' dir_sba ' -L' dir_sba ...
      ' -L' dir_external 'blasLapack/win' ext ' -output sba.' mexext ' '];
    % Matlab on Windows 32/64
    if strcmp(computer, 'PCWIN')
      eval([ 'mex LINKFLAGS="$LINKFLAGS /NODEFAULTLIB:libcmt.lib" ' ...
        library_flags '-lclapack -lblas -llibF77 -llibI77 ' sba_files ]);
    else
      eval([ 'mex LINKFLAGS="$LINKFLAGS /NODEFAULTLIB:libcmt.lib" ' ...
        library_flags '-lclapack_nowrap -lf2c -lBLAS_nowrap ' sba_files ]);
    end
  case {'GLNX86'},
    % Matlab on Linux
    eval([ 'mex ' gcc_extra ' -I../ -O -ldl ' sba_files ...
      ' /usr/lib/liblapack.a /usr/lib/atlas/libblas.a ' ...
      '/usr/lib/libgfortran.so.3' ]);
  case {'GLNXA64'},
    % Matlab on Linux 64
    cd ../
    system('make cleanall');
    system('make libsba.a');
    cd matlab
    eval([ 'mex -lgfortran ' gcc_extra ' -I../ -O -ldl sba.c ' ...
      '../libsba.a /usr/lib/liblapack_pic.a ' ...
      '../../blasLapack/linux64/libf77blas.a ' ...
      '../../blasLapack/linux64/libatlas.a '
      ]);
  case {'i686-pc-linux-gnu', 'i686-unknown-linux-gnu', ...
      'x86_64-pc-linux-gnu', 'x86_64-unknown-linux-gnu'},
    % Octave on Linux
    mkoctfile --mex ./sba.c ../sba_chkjac.c ../sba_levmar_wrap.c ../sba_levmar.c ...
    ../sba_lapack.c ../sba_crsm.c -I../
end
cd ../../..

cd sfm/private/sba
disp('sba projection routines');
switch computer
  case {'PCWIN'},
    if doSba
      % Matlab on Windows 32
      system('cl /nologo /O2 sbaProjection.c /link /dll /out:sbaProjection32.dll');
    end
  case {'PCWIN64'},
    if doSba
      % Matlab on Windows 64
      cc = mex.getCompilerConfigurations('C');
      system([ '"' cc.Location '\VC\bin\amd64\cl.exe" /nologo /O2 ' ...
        '/GS- sbaProjection.c /link /dll /NOENTRY ' ...
        '/out:sbaProjection64.dll']);
    end
  case {'GLNX86','GLNXA64','i686-pc-linux-gnu', 'i686-unknown-linux-gnu', ...
      'x86_64-pc-linux-gnu', 'x86_64-unknown-linux-gnu'},
    % Matlab and Octave on Linux
    system('gcc -Wall -fPIC -O3 -shared -o sbaProjection.so sbaProjection.c');
end
cd ../../..

% compile the remaining mex files
rd=fileparts(mfilename('fullpath')); rd=rd(1:end-9);

% general compile options (can make architecture specific)
optsAfter={};
switch computer
  case {'PCWIN','PCWIN64'},
    if strcmp(computer, 'PCWIN')
      win_folder = 'win32';
    else
      win_folder = 'win64';
    end
    % Matlab on Windows
    lapacklib = fullfile(matlabroot, 'extern', 'lib', win_folder, ...
      'microsoft', 'libmwlapack.lib');
    blaslib = fullfile(matlabroot, 'extern', 'lib', win_folder, ...
      'microsoft', 'libmwblas.lib');
    opts={'-output'};
    optsAfter = {lapacklib, blaslib};
  case {'GLNX86', 'GLNXA64'},
    % Matlab on Linux
    opts = { ['CXX=g++-' gccVer] ['CC=g++-' gccVer] ['LD=g++-' ...
      gccVer] '-largeArrayDims' '-l' 'mwlapack' '-l' 'mwblas' '-output'};
  case {'i686-pc-linux-gnu', 'i686-unknown-linux-gnu', ...
      'x86_64-pc-linux-gnu', 'x86_64-unknown-linux-gnu'},
    % Octave on Linux
    opts = {'-o'};
end

% general compile options
fs={'computeCsfmInfimumMex','refineExteriorOrientationMex', ...
  'computeH'};%, ...
%  'msfmRotationDerivativeMex'};
ds={'nrsfm/private/', 'sfm/private/', 'external/em-sfm/', ...
  'msfm/private/'};
for i=1:length(fs)
  disp(fs{i});
  mex(opts{:},[rd '/' ds{i} '/' fs{i} '.' mexext], ...
    [rd '/' ds{i} '/' fs{i} '.c'], optsAfter{:});
end

cd(savepwd); disp('..................................Done Compiling');

end

author: vrabaud


text:
function [ EMGwhisk ] = GetWhiskFromPiezo( basePath,varargin )
%[ EMGwhisk ] = GetWhiskFromEMG(basePath )
%This is a detector that extracts whisking/nonwhisking epochs from
%implanted EMG in the whisker pad. Extracts also the EMG and EMG envelope.
%
%INPUT
%   Assumes presence of the following files:
%       basePath/baseName.abf    (whisking/camera pulses from clampex)
%       basePath/analogin.dat    (camera pulses to intan)
%   where basePath is a folder of the form:
%       whateverPath/baseName/
%If basePath not specified, tries the current path.
%
%   (options)
%       'PulseChannel'
%       'EMGChannel'
%
%OUTPUT
%Creates file:
%   basePath/baseName.EMGwhisk.states.mat
%
%
%
%To Add: option for no camera pulses, which just align time to intan/pupil
%
%
%DLevenstein 2017
%% DEV
% basePath = '/mnt/proraidDL/Database/WMProbeData/';
% baseName = 'Layers_LFP_Test02_170323_151411';
%%

if ~exist('basePath','var')
    basePath = pwd;
end
baseName = bz_BasenameFromBasepath(basePath);
%%
abfname = fullfile(basePath,[baseName,'.abf']);
analogName = fullfile(basePath,['analogin.dat']);
savefile = fullfile(basePath,[baseName,'.EMGwhisk.states.mat']);
figfolder = fullfile(basePath,'DetectionFigures');

if ~exist(abfname,'file')
    display('No .abf file named basePath/baseName/baseName.abf')
end
if ~exist(analogName,'file')
    display('No analogin.dat in basePath/baseName/')
end

%% Clampex File

[abffile,si,file_info] = abfload(abfname);
sf_abf = 1./(si.*10^-6);

%Prompt user for channel
numabfchans = length(abffile(1,:));
if ~exist('chanNums','var')
    piezochan = listdlg('ListString',file_info.recChNames,...
        'PromptString',['Which channels is Whisker? ']);
    %Replace this with prompt file_info.recChNames
end
%%
EMG = abffile(:,piezochan);
t_abf = [1:length(EMG)]'./sf_abf;

%%
% figure
% plot(t_abf,EMG,'k')
%% Get the whisking envelope
%EMGrange = [400 3000];
%EMG = FiltNPhase(EMG,EMGrange,sf_abf);

downsamplefactor = 16; %Downsample to same as the LFP;
EMG = downsample(EMG,downsamplefactor);
t_EMG = downsample(t_abf,downsamplefactor);
sf_down = sf_abf./downsamplefactor;
%%

EMGparms.gausswidth = 0.05;  %Gaussian width for smoothing (s)
EMGparms.Whthreshold = 3;    %EMG Threshold for Whisking (modSTDs)
EMGparms.NWhthreshold = 0.5;    %EMG Threshold for Whisking (modSTDs)
%EMGparms.threshold = 1;    %EMG Threshold for Whisking (modSTDs)
EMGparms.minwhisk = 0.1;     %Minimum whisking duration (s)
EMGparms.minNWh = 0.1;       %Minimum nonwhisking duration (s)
EMGparms.whiskmerge = 0.1;     %Minimum interwhisking duration (s)
EMGparms.NWhmerge = 0.02;       %Minimum internonwhisking duration (s)

%% Z-Score the EMGZ and get EMG envelope with RMS
EMGz = NormToInt(EMG,'modZ'); %Modified Z score - robust to outliers
EMGsm = RMSEnvelope(EMGz,EMGparms.gausswidth,1/sf_down);
EMGsm = EMGsm-min(EMGsm);
EMGwhisk.EMGenvelope = EMGsm;

%% Set the thresholds by whisking troughs -
%find by "gradient descent"(ish) from initial guess
EMGbins = linspace(-1.5,2,100);
EMGhist = hist(log10(EMGsm),EMGbins);
EMGgrad = smooth(gradient(EMGhist),4);

%Find troughs (gradient crossing from - to +)
troughidx = find(diff(EMGgrad>0)==1);
troughs = 10.^EMGbins(troughidx);

%Get sign of gradient at each of the thresholds and use that to pick trough
Whsign = sign(interp1(EMGbins,EMGgrad,log10(EMGparms.Whthreshold),'nearest'));
if Whsign==-1
    EMGparms.Whthreshold = troughs(find(troughs>EMGparms.Whthreshold,1,'first'));
    if isempty(EMGparms.Whthreshold)
        EMGparms.Whthreshold = troughs(end);
    end
elseif Whsign==1
    EMGparms.Whthreshold = troughs(find(troughs<EMGparms.Whthreshold,1,'last'));
end

NWhsign = sign(interp1(EMGbins,EMGgrad,log10(EMGparms.NWhthreshold),'nearest'));
if NWhsign==-1
    EMGparms.NWhthreshold = troughs(find(troughs>EMGparms.NWhthreshold,1,'first'));
elseif NWhsign==1
    EMGparms.NWhthreshold = troughs(find(troughs<EMGparms.NWhthreshold,1,'last'));
end
%%

% figure
% bar(EMGbins,EMGhist)
% hold on
% plot(EMGbins,EMGgrad)
% plot([1 1].*log10(EMGparms.Whthreshold),get(gca,'ylim'),'g','LineWidth',2)
% plot([1 1].*log10(EMGparms.NWhthreshold),get(gca,'ylim'),'r','LineWidth',2)
%
% axis tight
% xlabel('EMG Envelope (modZ)');
% LogScale('x',10)
% xlim([-1.5 max(log10(EMGsm))])

%% Identify Whisking on/offsets: EMG envelope crosses threshold
wh_thresh = EMGsm > EMGparms.Whthreshold;
wh_on = find(wh_thresh(2:end)>wh_thresh(1:end-1))+1; %whisking onsets (si)
wh_off = find(wh_thresh(2:end)< wh_thresh(1:end-1))+1;%whisking offsets (si)

% If data starts/ends in the middle of an epoch, drop first/last trigger
if wh_off(1)<wh_on(1)
    wh_off = wh_off(2:end);
end
if wh_off(end) < wh_on(end)
    wh_on = wh_on(1:end-1);
end

%% Identify NonWhisking on/offsets: EMG envelope crosses below threshold
nwh_thresh = EMGsm < EMGparms.NWhthreshold;
nwh_on = find(nwh_thresh(2:end)>nwh_thresh(1:end-1))+1; %whisking onsets (si)
nwh_off = find(nwh_thresh(2:end)< nwh_thresh(1:end-1))+1;%whisking offsets (si)

% If data starts/ends in the middle of an epoch, drop first/last trigger
if nwh_off(1)<nwh_on(1)
    nwh_off = nwh_off(2:end);
end
if nwh_off(end) < nwh_on(end)
    nwh_on = nwh_on(1:end-1);
end

%% Merge and Duration Thresholds
%Convert to seconds
wh_on = wh_on./sf_down;
wh_off = wh_off./sf_down;
nwh_on = nwh_on./sf_down;
nwh_off = nwh_off./sf_down;

%Merge brief interruptions
[ NWhints ] = MergeSeparatedInts( [nwh_on,nwh_off],EMGparms.NWhmerge );
[ Whints ] = MergeSeparatedInts( [wh_on,wh_off],EMGparms.whiskmerge );

%Drop nonwhisk epochs smaller than a minimum
[nwh_on,nwh_off] = MinEpochLength(NWhints(:,1),NWhints(:,2),EMGparms.minNWh,1);
%Drop whisking epochs smaller than a minimum
[wh_on,wh_off] = MinEpochLength(Whints(:,1),Whints(:,2),EMGparms.minwhisk,1);

%% Durations
Whdur = wh_off-wh_on;
NWhdur = nwh_off-nwh_on;

numbinst = 25;
durhist.bins = linspace(-1,2,numbinst);
durhist.Wh = hist(log10(Whdur),durhist.bins);
durhist.NWh = hist(log10(NWhdur),durhist.bins);

%% Reset time to align to LFP

%t_EMG = t_EMG-firstpulstime_abf+firstpulstime_lfp;
Whints = [wh_on wh_off];
NWhints = [nwh_on nwh_off];

%% Figure

figure
subplot(4,1,1)
    plot(t_EMG,EMGz,'k')

    hold on
    plot(t_EMG,EMGsm,'b','linewidth',2)
    plot(Whints',EMGparms.Whthreshold.*ones(size(Whints))','g','linewidth',2)
    plot(NWhints',EMGparms.NWhthreshold.*ones(size(NWhints))','r','linewidth',2)
    axis tight
    ylim([-100 100])
    ylabel('EMG (modZ)');

subplot(4,1,2)
    plot(t_EMG,EMGz,'k')

    hold on
    plot(t_EMG,EMGsm,'b','linewidth',2)
    plot(Whints',EMGparms.Whthreshold.*ones(size(Whints))','g','linewidth',2)
    plot(NWhints',EMGparms.NWhthreshold.*ones(size(NWhints))','r','linewidth',2)
    xlim([100 160])
    ylim([-20 40])
    ylabel('EMG (modZ)');

subplot(4,2,6)
hist(log10(EMGsm),100)
hold on
plot([1 1].*log10(EMGparms.Whthreshold),get(gca,'ylim'),'g')
plot([1 1].*log10(EMGparms.NWhthreshold),get(gca,'ylim'),'r')

axis tight
xlabel('EMG Envelope (modZ)');
LogScale('x',10)
xlim([-1.5 max(log10(EMGsm))])

subplot(4,2,8)
plot(durhist.bins,durhist.NWh,'r','linewidth',2)
hold on
plot(durhist.bins,durhist.Wh,'g','linewidth',2)
LogScale('x',10)
xlabel('Duration (s)')
ylabel('# Epochs')
legend('NWh','Wh')

NiceSave('WhiskingDetection',figfolder,baseName)

%%

EMGwhisk.ints.Wh = Whints;
EMGwhisk.ints.NWh = NWhints;
EMGwhisk.detectorparms = EMGparms;
EMGwhisk.detectorname = 'GetWhiskFromEMG';
EMGwhisk.detectiondate = today('datetime');
EMGwhisk.EMG = EMGz;
EMGwhisk.EMGsm = EMGsm;
EMGwhisk.t = t_EMG;

save(savefile,'EMGwhisk')
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
classdef Max < Algorithm

    methods (Access = public)

        function obj = Max()
            obj.name = 'Max';
            obj.inputPort = DataType.kSignal;
            obj.outputPort = DataType.kFeature;
        end

        function dataOut = compute(~,signal)
            dataOut = max(signal);
        end

        function metrics = computeMetrics(~,input)
            n = size(input,1);
            flops = 1 * n;
            memory = 1;
            outputSize = Constants.kFeatureBytes;
            metrics = Metric(flops,memory,outputSize);
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
% OtTableEdit Put the table into edit mode.
%
% Usage:
%    OtTableEdit(htable)

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% 2006-04-11 Created (Johan Meijdam, VORtech)

function iret = OtTableEdit(htable)
   iret = Ot('table.edit', htable);
end


text:
function soln = cannon_multipleShooting(guess,target,param)
% soln = cannon_multipleShooting(guess,target,param)
%
% This function uses single shooting to solve the cannon problem.
%
% INPUTS:
%   guess.initSpeed
%   guess.initAngle
%   target.x
%   target.y
%   param.c = quadratic drag coefficient
%   param.nGrid = number of grid points for the integration method
%   param.nSegment = number of trajectory segments
%
% OUTPUTS:
%   soln.t
%   soln.x
%   soln.y
%   soln.dx
%   soln.dy
%

global ITER_LOG_MULTIPLESHOOTING;  %Used for diagnostics and visualization only
ITER_LOG_MULTIPLESHOOTING = [];

P.nSegment = param.multipleShooting.nSegment;
P.nSubStep = param.multipleShooting.nSubStep;
P.nGrid = P.nSegment*P.nSubStep;
P.c = param.dynamics.c;

%%% Run a simulation to get the initial guess:
init.speed = guess.initSpeed;
init.angle = guess.initAngle;
traj = simulateCannon(init,P);
guess.T = traj.t(end);  %Trajectory duration

%%% Break the guess trajectory at segment bounds:
guess.t = linspace(0,guess.T,P.nSegment+1); guess.t(end) = [];
guess.z = interp1(traj.t', [traj.x; traj.y; traj.dx; traj.dy]', guess.t', 'spline')';

%%% Store the initial guess for the problem:
nState = 4;  %Number of states in the problem (x,y,dx,dy)
problem.x0 = [guess.T, reshape(guess.z,1,nState*P.nSegment)];
problem.lb = [];    % Lower bound on decision variables
problem.ub = [];   % Upper bound on decision variables

%%% Set up the linear constraints (there are none);
problem.Aineq = [];  problem.Aeq = [];
problem.bineq = [];  problem.beq = [];

%%% Set up the user-defined functions:
problem.objective = @(decVar)objective(decVar(4),decVar(5));  %Objective (cost) function
problem.nonlcon = @(decVar)nonLinCst(decVar,target,P);   %NonLinear constraints

%%% Set up the options for the solver:
problem.solver = 'fmincon';
if param.diagnostics.enable  %Then record full diagnostics:
    problem.options = optimset(...
        'Display','iter',...
        'MaxFunEvals',1e4,...
        'MaxIter',100,...
        'OutputFcn',@(decVar,optimVal,state)outFun(decVar,optimVal,state));
else %Run things quickly!
    problem.options = optimset(...
        'MaxFunEvals',1e4,...
        'MaxIter',100,...
        'Display','off');
end

%Use FMINCON to solve the constrained optimization problem:
[zSoln, fVal, exitFlag] = fmincon(problem);

%Call the constraint function one final time to get the trajectory:
[~, ~, t, zTraj] = nonLinCst(zSoln,target,P);

%%% Store the trajectory in a nice format:
soln.t = t;
soln.x = zTraj(1,:);
soln.y = zTraj(2,:);
soln.dx = zTraj(3,:);
soln.dy = zTraj(4,:);
soln.success = exitFlag == 1;
soln.cost = fVal;
soln.method = 'Multiple Shooting';

%%% Run diagnostics on the solution if desired:
if param.diagnostics.enable
    diagnostics_multipleShooting(target,param,soln)
else
    figure(param.diagnostics.figNum.multipleShooting); clf;
    plotSoln(soln, target, param);
end

end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%           Non-linear constraint function for multiple shooting          %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

function [C, Ceq,tTraj,zTraj] = nonLinCst(decVar,target,P)
% This is the key function!
% There are two key constraints here: Boundary Conditions & Dynamics
% The dynamics are implemented using a multiple simulations
% The decision variables are time, followed by the states at the start of
% each trajectory segment.
% There are constraints placed on both the initial and final states
% The simulation is implemented using 4th-order runge-kutta*.

nState = 4;  %Number of states in the problem (x,y,dx,dy)
nSegment = P.nSegment;
tEnd = decVar(1);
z0 = reshape(decVar(2:end),nState,P.nSegment);

% Run a simulation from the start of each segment, in parallel
nSub = P.nSubStep;  %Number of sub-steps for the integration method
tSim = linspace(0,tEnd/nSegment, nSub+1);
z = rk4_cannon(tSim,z0,P.c); %Simulate the trajectory
% Index Info:  z(nState,nSegment,nSubStep)

%%% Boundary Value Constraints:
BoundaryInit = [z(1,1,1); z(2,1,1)]; %Initial Position
BoundaryFinal = [z(1,end,end); z(2,end,end)]...
    - [target.x; target.y]; %Final Position

%%% Defect Constraints:
zEnd = z(:, 1:(end-1), end);  %States at the end of a segment
zStart = z(:, 2:end, 1);   %States at the beginning of a segment
Defects = reshape(zStart-zEnd,nState*(nSegment-1),1);

C = [];  %No inequality constraints
Ceq = [BoundaryInit; Defects; BoundaryFinal];  %Boundary Condition

if nargout==4  %Only used for post-processing  --  return trajectory

    tTraj = zeros(1,nSegment*nSub+1);
    zTraj = zeros(nState,nSegment*nSub+1);

    %Stitch together the trajectory
    idx = 0;

    tSegment = linspace(0,tEnd,nSegment+1); tSegment(end) = [];
    for i=1:nSegment  %Slow looping, but only run once, so not too bad
        for j=1:nSub
            idx = idx+1;
            tTraj(idx) = tSim(j) + tSegment(i);
            zTraj(:,idx) = z(:,i,j);
        end
    end
    tTraj(end) = decVar(1);
    zTraj(:,end) = z(:,end,end);
end

end

% * I use a fixed-order method (rather than ode45) because of the improved
% consistency in the evaluation (the fixed-order method performs the exact
% same arithmetic operations on every call, where a variable order might
% not - for example, by adjusting the grid spacing, which can cause noise
% in the gradient estimates in the optimization method).

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                     Data-logging and Diagnostics:                       %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

function stop = outFun(decVar,optimVal,state)
% This function is used for logging the progress of fmincon throughout the
% optimization run. It does not affect the optimization process.

global ITER_LOG_MULTIPLESHOOTING;  %Keeping track of iteration details

stop = false;

switch state
    case 'init'
    case 'iter'
        iter = optimVal.iteration+1;
        ITER_LOG_MULTIPLESHOOTING(iter).optimVal = optimVal;
        ITER_LOG_MULTIPLESHOOTING(iter).decVar = decVar;
    case 'done'
    otherwise
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
%% circlefit
% Below is a demonstration of the features of the |circlefit| function

%%
clear; close all; clc;

%% Syntax
% |[Vc,R]=circlefit(V);|

%% Description
% This function returns the centre Vc and radius R for a circle fitted to
% the input point set V. The input can be 2D or 3D but the output centre is
% always a point in 3D space.

%% Examples
%

%%
% Plot settings
markerSize=30;
lineWidth=3;

%%
% Create input circle

n=50; % Number of points on the circle
r=2; % True radius
t=linspace(0,2*pi,n+1)'; t=t(1:end-1); % Angles
V_true=r.*[cos(t) sin(t) zeros(size(t))]; %Circle true coordinates
V=V_true+r/20*randn(size(V_true));

%%
% Fit a circle

[Vc,R]=circlefit(V);

nf=250;
r=2;
t=linspace(0,2*pi,nf+1)'; t=t(1:end-1);
Vf=Vc+R.*[cos(t) sin(t) zeros(size(t))];

%%
% Visualize result

cFigure; hold on;
hp1=plotV(V,'k.','MarkerSize',markerSize);
hp2=plotV(V_true,'g.-','MarkerSize',markerSize,'LineWidth',lineWidth);
hp3=plotV(Vf,'r-','LineWidth',lineWidth);
legend([hp1 hp2 hp3],{'Noisy input data','True circle','Fitted circle'},'Location','NorthEastOutside')
axisGeom; view(2);
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function h=entropy(p,dim,cond,arg,step)
%ENTROPY calculates the entropy of discrete and sampled continuous distributions H=(P,DIM,COND,ARG,STEP)
%
%  Inputs:  P        is a vector or matrix of probabilities - one dimension per variable
%           DIM      lists dimensions along which to evaluate the entropy [default: 1st non singleton dimension]
%           COND     lists dimensions to use as conditional variables [default - none]
%           ARG      lists dimensions to use as parameters in the ouput [default - none]
%           STEP     for continuous distributions STEP gives the sample increment for each dimension of P
%                    if STEP is a scalar, the increment is assumed to be the same for each dimension
%
% Outputs:  H        is the entropy. It will have the same number of dimensions as the length of the ARG input.
%                    If the STEP argument is specified then this will be the differential entropy.
%
% Example: Suppose P(W,X,Y,Z) represents the joint probability of four correlated random variables
%
%               (a) H(W,X,Y,Z) = entropy(P,[1 2 3 4]).
%               (b) H(W) = entropy(P), or equivalently entropy(P,1)
%               (c) H(W,Z | X,Y) = entropy(P,[1 4],[2 3])
%               (d) H(W | X, Z=z) = entropy(P,1,2,4); this is a function of z and will be a column vector
%
% As a special case, if the dimensions included in DIM are all singletons, the entries in P are treated
% as Bernoulli variable probabilities.

%	   Copyright (C) Mike Brookes 2006
%      Version: $Id: entropy.m,v 1.3 2007/05/04 07:01:38 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<5
    stp=zeros(ndims(p),1);
else
    stp=repmat(step(1),ndims(p),1);
    stp(1:length(step))=step(:);
    stp=log2(stp);
end
if nargin<4
    arg=[];
else
    arg=arg(arg>0);
end
if nargin<3
    cond=[];
else
    cond=cond(cond>0);
end
if ~length(cond)
    s=size(p);
    if nargin<2
        dim=find(s>=min(2,max(s)));
        dim=dim(1);
    else
        dim=dim(dim>0);
    end
    st=prod(s);
    sd=prod(s(dim));
    sa=prod(s(arg));
    marg=1:length(s);
    marg(arg)=0;
    marg(dim)=0;
    marg=marg(marg>0);
    sm=st/sd/sa;
    if sm>1
        ip=[arg dim(:)' marg(:)'];
        sp=[s([arg dim(:)']) prod(s(marg))];
        q=sum(reshape(permute(p,[arg(:)' dim(:)' marg]),sa,sd,sm),3);
    else
        q=reshape(permute(p,[arg(:)' dim(:)' marg]),sa,sd);
    end
    if sd==1
        h=-log2(q+(q==0)).*q-log2(1-q+(q==1)).*(1-q);   % special treatment for bernoulli variables
    else
        sq=sum(q,2);
        h=sum(-log2(q+(q==0)).*q,2)./sq+log2(sq);
    end
    if length(arg)>1
        h=reshape(h,s(arg));
    end
    h=h+sum(stp(dim));
else
    % we could probably make this more efficient by avoiding the recursive call
    h=entropy(p,[dim(:); cond(:)],0,arg)-entropy(p,cond,0,arg);
end

author: decouples


text:
function A11 = synA11min(A10, A01, cmax)
%-----------------------------------------------------------------------------
%
% For each point of colour 11 this function assigns the minimum value at the
% neighbouring gridpoints of colours 10 and 01.
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: December 7, 2001.
% (c) 1998-2002 Stichting CWI, Amsterdam
%-----------------------------------------------------------------------------
[n10, m10]=size(A10);
[n01, m01]=size(A01);
n11=n10;
m11=m01;
%[n11, m11]=size(A11);
if     m11 == m10
  S=min(stripL(extR(A10, cmax)), A10);
elseif m11 == m10-1
  S=min(stripL(A10), stripR(A10));
else
  disp([' size A10 = ' int2str(size(A10)) ' size A01 = ' int2str(size(A01))]);
  error(' synA11min - A10 and A01 do not match ');
end
if     n11 == n01
  T=min(A01, stripU(extD(A01, cmax)));
elseif n11 == n01-1
  T=min(stripD(A01), stripU(A01));
else
  disp([' size A10 = ' int2str(size(A10)) ' size A01 = ' int2str(size(A01))]);
  error(' synA11min - A10 and A01 do not match ');
end
%Note: all(size(S) == size(T)) & all(size(S) == [n11 m11]) always holds.
A11=min(S, T);
%-----------------------------------------------------------------------------

author: Sable


text:
classdef io

methods(Static)

function result = yes_no(flag)
    if flag == 0
        result = 'No';
    else
        result = 'Yes';
    end
end

function result = true_false(flag)
    if flag == 0
        result = 'false';
    else
        result = 'true';
    end
end

function result = true_false_short(flag)
    if flag == 0
        result = 'F';
    else
        result = 'T';
    end
end

function result = tf_pair(a, b)
    a = spx.io.true_false_short(a);
    b = spx.io.true_false_short(b);
    result = [a b];
end

function result = success_failure(flag)
    if flag == 0
        result = 'Failure';
    else
        result = 'Success';
    end
end

function result = success_failure_short(flag)
    if flag == 0
        result = 'F';
    else
        result = 'S';
    end
end

function result = bytes2str(bytes)
    % Identify the scale of number of bytes
    scale = floor(log(bytes)/log(1024));
    switch scale
        case 0
            result = sprintf('%d Bytes', bytes);
        case 1
            result = sprintf('%.2f KB', bytes / (1024^1));
        case 2
            result = sprintf('%.2f MB', bytes / (1024^2));
        case 3
            result = sprintf('%.2f GB', bytes / (1024^3));
        case 4
            result = sprintf('%.2f TB', bytes / (1024^4));
        case -inf
            result = 'Not available';
        otherwise
            result = 'Too big';
    end

end

end

end

author: indigits

repo: sparse-plex


text:
function exact = p44_exact ( )

%*****************************************************************************80
%
%% P44_EXACT returns the exact integral for problem 44.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the value of the integral.
%
  alpha = p44_param_get ( );

  exact = ( 20.0 * sin ( 2.0^alpha ) ...
    - 2.0^alpha * cos ( 2.0^alpha ) ...
    + 2.0^alpha * exp ( -20.0 ) ) / ( 400.0 + 4.0^alpha );

  return
end


text:
function [ret]=realTime_startDirectServoCartesian( t_Kuka )
%% Syntax:
% realTime_startDirectServoCartesian( t_Kuka )

%% About:
% This function is used to turn on the direct servo function on the robot
% This function is for direct-servoing in Cartesian space, i.e. controlling the
% robot at EEF level.
% After starting the direct servo, using this function, you have to send
% the EEF target postitions using the function,
% ((sendEEfPosition)), refer also to the function
% ((realTime_stopDirectServoCartesian)).

%% Arreguments
% t_Kuka: is the TCP/IP connection object

% Copy right, Mohammad SAFEEA, 1st of April 2018

theCommand='stDcEEf_';
fprintf(t_Kuka, theCommand);
message=fgets(t_Kuka);

[ret]=checkAcknowledgment(message);
if ret==true
    disp('Acknowledge, realtime control turned on');
end

delay(0.5); % introduce some time delay, so the robot turns on the direct servo before starting the motion
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function [r, g, b] = imSplitChannels(rgb)
% Split the 3 channels of a 2D or 3D image.
%
%   [R, G, B] = imSplitChannels(RGB);
%
%   Example
%   img = imread('peppers.png');
%   [r, g, b] = imSplitChannels(img);
%   figure; imshow(r); title('red');
%   figure; imshow(g); title('green');
%
%   See also
%     imMergeChannels

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2010-02-02,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

if ndims(rgb) == 3
    % process planar image
    r = rgb(:,:,1);
    g = rgb(:,:,2);
    b = rgb(:,:,3);

elseif ndims(rgb) == 4
    % process 3D image
    r = squeeze(rgb(:,:,1,:));
    g = squeeze(rgb(:,:,2,:));
    b = squeeze(rgb(:,:,3,:));
else
    error('unprocessed image dimension');
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function adjacency_structure = gr_adjacency_structure ( node_num, ...
  node_coordinates, edge_num, edge_nodes )

%*****************************************************************************80
%
%% GR_ADJACENCY_STRUCTURE returns the adjacency structure of a graph.
%
%  Discussion:
%
%    Since we are using MATLAB, we return an actual structure.
%    That is, a cell array.  Entry {I} of the cell array is a
%    vector containing all the nodes which share an edge with
%    node I.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_COORDINATES(2,NODE_NUM), the coordinates of the nodes.
%
%    Input, integer EDGE_NUM, the number of edges.
%
%    Input, integer EDGE_NODES(2,EDGE_NUM), the indices of the two nodes
%    that form each edge.
%
%    Output, cell array ADJACENCY_STRUCTURE{NODE_NUM}, the
%    adjacency structure.
%
  adjacency_list = gr_adjacency_list ( node_num, ...
    node_coordinates, edge_num, edge_nodes );

  adjacency_pointer = gr_adjacency_pointer ( node_num, ...
    node_coordinates, edge_num, edge_nodes );

  adjacency_structure = cell ( node_num );

  for i = 1 : node_num
    i1 = adjacency_pointer(i);
    i2 = adjacency_pointer(i+1) - 1;
    adjacency_structure{i} = adjacency_list(i1:i2);
  end

  return
end


text:
function calpak_test502 ( )

%*****************************************************************************80
%
%% CALPAK_TEST502 tests WEEKDAY_TO_NAME_HEBREW.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 April 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST502\n' );
  fprintf ( 1, '  For the Hebrew calendar,\n' );
  fprintf ( 1, '  WEEKDAY_TO_NAME_HEBREW names the days of the week.\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 7
    s = weekday_to_name_hebrew ( i );
    fprintf ( 1, '  %2d  %s\n', i, s );
  end

  return
end


text:
% 2D Lattice Boltzmann (BGK) model of a fluid.
%  c4  c3   c2    D2Q9 model. At each timestep, particle densities propagate
%    \  |  /      outwards in the directions indicated in the figure. An
%  c5 -c9 - c1    equivalent 'equilibrium' density is found, and the densities
%    /  |  \      relax towards that state, in a proportion governed by omega.
%  c6  c7   c8
%                  Iain Haslam, April 2006.

omega=1.0; density=1.0; t1=4/9; t2=1/9; t3=1/36; c_squ=1/3; nx=1; ny=101;
F=repmat(density/9,[nx ny 9]); F_EQ=F; msize=nx*ny; qi=[0:msize:msize*7];
bound=zeros(nx,ny); bound(:,[1 ny])=1;%bound([14:19],[14:16])=1;
obs=find(bound);
bounds_to_bounce=[obs+qi(1) obs+qi(2) obs+qi(3) obs+qi(4) ...
	obs+qi(5) obs+qi(6) obs+qi(7) obs+qi(8)];
bounds_bounced=[obs+qi(5) obs+qi(6) obs+qi(7) obs+qi(8) ...
	obs+qi(1) obs+qi(2) obs+qi(3) obs+qi(4)];
avu=1; prevavu=1; ts=0; deltaU=1e-9; numactivenodes=sum(sum(1-bound));

while (ts<100000 & 1e-10<abs((prevavu-avu)/avu)) | ts<100
    ts=ts+1;

    % Propagate
	F(:,:,4)=F([2:nx 1],[ny 1:ny-1],4);    F(:,:,3)=F(:,[ny 1:ny-1],3);
	F(:,:,2)=F([nx 1:nx-1],[ny 1:ny-1],2); F(:,:,5)=F([2:nx 1],:,5);
	F(:,:,1)=F([nx 1:nx-1],:,1);           F(:,:,6)=F([2:nx 1],[2:ny 1],6);
	F(:,:,7)=F(:,[2:ny 1],7);              F(:,:,8)=F([nx 1:nx-1],[2:ny 1],8);
	BOUNCEDBACK=F(bounds_to_bounce); %Densities bouncing back at next timestep

    % Relax; calculate equilibrium state with equivalent speed and density to F
	DENSITY = sum(F,3);
	UX=(sum(F(:,:,[1 2 8]),3)-sum(F(:,:,[4 5 6]),3))./DENSITY;
	UY=(sum(F(:,:,[2 3 4]),3)-sum(F(:,:,[6 7 8]),3))./DENSITY;

	%UX(1,2:ny-1)=UX(1,2:ny-1)+deltaU;
	UX=UX+deltaU;
	UX(obs)=0; UY(obs)=0; DENSITY(obs)=0;
	U_SQU=UX.^2+UY.^2; U_C2=UX+UY; U_C4=-UX+UY; U_C6=-U_C2; U_C8=-U_C4;

	% stationary
	F_EQ(:,:,9)=t1*DENSITY.*(1-U_SQU/(2*c_squ));

	% nearest-neighbours
	F_EQ(:,:,1)=t2*DENSITY.*(1+UX/c_squ+0.5*(UX/c_squ).^2-U_SQU/(2*c_squ));
	F_EQ(:,:,3)=t2*DENSITY.*(1+UY/c_squ+0.5*(UY/c_squ).^2-U_SQU/(2*c_squ));
	F_EQ(:,:,5)=t2*DENSITY.*(1-UX/c_squ+0.5*(UX/c_squ).^2-U_SQU/(2*c_squ));
	F_EQ(:,:,7)=t2*DENSITY.*(1-UY/c_squ+0.5*(UY/c_squ).^2-U_SQU/(2*c_squ));

	% next-nearest neighbours
	F_EQ(:,:,2)=t3*DENSITY.*(1+U_C2/c_squ+0.5*(U_C2/c_squ).^2-U_SQU/(2*c_squ));
	F_EQ(:,:,4)=t3*DENSITY.*(1+U_C4/c_squ+0.5*(U_C4/c_squ).^2-U_SQU/(2*c_squ));
	F_EQ(:,:,6)=t3*DENSITY.*(1+U_C6/c_squ+0.5*(U_C6/c_squ).^2-U_SQU/(2*c_squ));
	F_EQ(:,:,8)=t3*DENSITY.*(1+U_C8/c_squ+0.5*(U_C8/c_squ).^2-U_SQU/(2*c_squ));
	F=omega*F_EQ+(1-omega)*F;
	F(bounds_bounced)=BOUNCEDBACK;
	prevavu=avu;avu=sum(sum(UX))/numactivenodes; ts=ts+1;
end
%The bounceback boundary condition means that the actual boundary lies
%mid-way between the open and closed nodes. The width of the channel is
%therefore ny-2.
L=(ny-2)/2; arrindices=[-(L-0.5):(L-0.5)];     EXACT=-(arrindices.^2-L^2);
tau = 1/omega;

%Calculate the analytical solution's peak velocity
temp_F         = density * deltaU / tau;
temp_viscosity = (tau-0.5)/3;
calcfactor     = temp_F/(2*temp_viscosity);
CALC           = [UX(1,2:ny-1)];
EXACT          = EXACT*calcfactor;
figure;
plot(arrindices,abs(1-(EXACT./CALC)));
maxerror=max(abs(1-(EXACT./CALC)));
title('Relative error of LBM solution across channel width');
fprintf('RESULTS\n')
fprintf('ts=%d, maxerror=%g\n',ts, maxerror);

%Calculate relative proportional error in centre of channel
[maxcalc, maxindex] = max(CALC);
maxerrorcentre = (EXACT(maxindex)-CALC(maxindex))/EXACT(maxindex)
figure;plot(arrindices,EXACT);
hold on;
plot(arrindices,CALC,'o');
title('Comparison of analytical with LBM results');
xlabel('location in channel'); ylabel('speed')

repo: Aero-matlab