text:
%% Toolbox Sparsity - A toolbox for sparse coding and sparse regularization
%
% Copyright (c) 2008 Gabriel Peyre
%

%%
% The toolbox can be downloaded from Matlab Central
% http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=16204&objectType=FILE

%%
% We first includes in the path some additional useful scripts.
path(path, 'toolbox/');

%% Overview
%
% This toolbox implements several algorithms to compute sparse expansion in redundant dictionaries and to solve inverse problems with sparse regularization
% (and also TV regularization).
%

%%
% Sparse expansion of a signal |f| of size |n| in a redundant dictionary |D| (which is a matrix of size |n x p| with |p>n|, or an implicit operator)
% corresponds to the computation of coefficients |x| of size |p| such that |f=D*x| (or approximate equality |norm(y-D*x)<epsilon|). Since many |x| are possible, one
% can assume that |x| is sparse, for instance supposing that it as a small L1 norm |sum(abs(x))|. This can be achieved by Lagrangian optimization
%

%%
% |min_x 1/2*norm(f-D*x)^2 + T*sum(abs(x))|  (1)

%%
% where |T| is increased if |epsilon| is increased.

%%
% Sparse regularization for the resolution of |y=U*f + noise| (where |U| is a rank defficientis equivalent) corresponds to the computation of a solution
% |f| that is sparse in a dictionary |D|, thus being |f=D*x| with small L1 norm. This is performed by solving

%%
% |min_x 1/2*norm(y-U*D*x)^2 + T*sum(abs(x))|   (2)

%%
% And the solution is computed from the optimal coefficients |x| as |f=D*x|.
%
% Problems (1) and (2) are equivalent if one replaces the dictionary |D| of (1) by the transformed dictionary |U*D| in (2). One can thus use the same
% software to solve both (2) and (1).

%%
% TV regularization corresponds to replacing (2) by

%%
% |min_f 1/2*norm(y-U*f)^2 + T*TV(f)|   (3)
%
% where |TV(f)| is the discrete TV norm of the signal/image |f|. Intuitively, this corresponds to assuming that the gradient of |f| is sparse, so (3) is
% quite close to (2) if |D| is chosen a translation invariant Haar wavelet dictionary. But optimizing (3) and (2) necessitate different algorithms.

%%
%
% * To solve (1) or (2) with |epsilon=0| (no noise), which corresponds to the minimum L1 norm solution to |f=D*x| or |y=U*D*x|, this toolbox implements the Douglas-Rachford iterative algorithm.
% * To solve (1) or (2), this toolbox implements the iterative soft thresholding algorithm. If |T| is not known but |epsilon| is known, then the implementation automatically computes the correct |T| during the iterations.
% * To solve (3), this toolbox implements Chambolle's JMIV 2002 algorithm. If |T| is not known but |epsilon| is known, then the implementation automatically computes the correct |T| during the iterations.

%% Compressed sensing of exactely sparse signals with Fourier measurements.

%%
% First, we create an |s|-sparse signal with random spikes locations and
% random coefficients values.

% dimensionality
n = 1024;
% sparsity
s = 25;
% generate random spikes
x0 = compute_rand_sparse(n, s, 'uniform');

%%
% We select an operator |U|, which is a pre-defined implicit operator.
% This is done by filling in a structure |options| with the correct
% parameters. Many other measurement matrices are available for Compressive
% Sensing (Gaussian, Bernouilli, Fourier, Sinus, Hadamard, etc).

% number of measurements
p = 100;
% type of matrix
clear options;
options.cs_type = 'fourier';
options.n = n;
options.p = p;

%%
% We perform the measurements (without noise) by applying the sensing
% operator.

y = callback_sensing_rand(x0, +1, options);

%%
% We can display the signal vector and the measurements

clf;
subplot(2,1,1);
plot_sparse_diracs(x); title('Signal x_0');
subplot(2,1,2);
h = plot(real(y)); axis tight;
set(h, 'LineWidth', 2);
set(gca, 'FontSize', 20);
title('Measurements y=U*x_0');

%%
% Since there is no noise, the recovery can be performed using
% Douglas-Rachford algorithm. Special thanks to Jalal Fadili for his help
% on this algorithm.

% set to 1 these options if you want cool displays
options.verb = 0;
options.draw_iter = 0;
% parameter of the algorithm
options.niter = 3000;
options.x = zeros(n,1);
[xlun,lun] = perform_douglas_rachford(@callback_sensing_rand, y, options);
xlun = real(xlun);

%%
% We can display the result. This is a case of perfect recovery, |xlun=x0|!
% However, if you increase slightly the value of |s|, then the recovery
% fails ...

clf;
options.title = {'Signal' 'Recovery'};
options.val = .01; % remove small residual
plot_sparse_diracs({x0 xlun}, options);

%%
% We can also display the decrease of the L1 norm of the solution during
% the iterations, to check that the algorithm really does a good job.

clf;
h = plot(lun); axis tight;
set(h, 'LineWidth', 2);
set(gca, 'FontSize', 20);
xlabel('#iteration'); ylabel('L_1 norm |x|_1');

%% Compressed sensing of a noisy sparse signals with Gaussian measurements.

%%
% To be more robust to noise, we need to decrease the number of
% measurements. Also, we are going to use Gaussian real measurements, which
% in fact perform twice less measurements than complex Fourier
% measurements.

% sparsity
s = 20;
% generate random spikes
x0 = compute_rand_sparse(n, s, 'uniform');

%%
% We perform noisy measurements.
% We also select an other kind of matrix, given in explicit form (an array |U|
% of double). This is possible only for small size problems.

% Generate a matrix for sensing
U = randn(p,n);
% noiseless measurements
y = U*x0;
% make some noise
sigma = .06*std(y(:)); % noise level
y = y + sigma*randn(size(y));

%%
% Since there is some noise, we need to use iterative thresholding.
% The regularization parameter is adapted automatically to fit the noise
% level |sigma|.

% amplification is >1 to further sparsify the solution
amplification = 2;
% norm of the noise |epsilon|
options.etgt = amplification * sigma*sqrt(p);
% initial regularization parameter
options.T = 1;
% perform iterative thresholding
options.niter = 4000;
[xlun,err,lun,lambda] = perform_iterative_thresholding( U,y,options);

%%
% We can display the evolution of the threshold |lambda| during the iterations
% and also display the decay of the Lagrangian to see what the algorithm is
% doing.

clf;
subplot(2,1,1);
h = plot(lambda); axis tight;
set(h, 'LineWidth', 2);
set(gca, 'FontSize', 20);
title('Evolution of \lambda');
subplot(2,1,2);
lagr = .5*err.^2 + lambda(end)*lun;
h = plot(lagr); axis tight;
set(h, 'LineWidth', 2);
set(gca, 'FontSize', 20);
title('Evolution of |y-U x|^2+\lambda |x|_1');
axis([1 options.niter min(lagr) min(lagr)*1.5]);

%%
% The solution of the L1 minimization is biased, because L1 tends to
% under-estimate the value of the true coefficients (this is very much
% similar to the issue with soft thresholding). To remediate to this
% issue, it is possible to "debiased" the solution by extracting the
% support and then performing an L2 best fit (which is not biased) of the
% measurements. Be aware that this might give bad results if the
% support is baddly estimated.

% remove small amplitude coefficients because they are not very reliable
xdeb = xlun .* (abs(xlun)>.05);
% perform L2 regression
xdeb = perform_debiasing(U,xdeb,y);

%%
% We can display the result. There is no more perfect recovery because of
% the noise, but the solution is close to the original signal.

clf;
options.title = {'Signal' 'Recovery' 'Debiased'};
options.val = .01; % remove small residual
plot_sparse_diracs({x0 xlun xdeb}, options);

%% Tomography inversion of an image with TV regularization.
% We can recover from partial Fourier measurements using TV regularization.
% We add some noise to the measurement, and perform the recovery using
% Chambolle's algorithm for TV regularization.

%%
% First we load the image.

n = 256;
M = load_image('phantom', n);
M = rescale(M,.05,.95);

%%
% Then we load the tomography mask.

nrays = 18; % Candes/Romberg/Tao set up
mask = compute_tomography_mask(n,nrays);
% number of measurements
p = sum(mask(:)==1);
% display
imageplot({M mask}, {'Signal' 'Frequencies'});

%%
% We set up |options| for the tomography operator (sub-sampling of the
% frequencies) and perform measurements.

options.size = size(M);
options.mask = mask;
y = callback_tomography(M,+1,options);
% add some noise
sigma = .01*std(y(:));
y = y + sigma*randn(size(y));

%%
% We can perform L2 pseudo-inversion (zero padding in Fourier) by just applying
% the transpose operator.

M2 = callback_tomography(y,-1,options); M2 = real(M2);
imageplot({M clamp(M2)}, {'Original' 'L2 inversion'});

%%
% We can perform inversion with TV regularization.

% set up parameter for the regularization
options.TVoperator = @callback_tomography;
options.TVrhs = y;
options.etgt = 4*sigma*sqrt(p); % target noise level
options.niter = 300;
options.niter_inner = 30; % iteration for inner loop of Chambolle's algorithm.
options.tol = 1e-4;
options.display = 0; % set to 1 for sexy display
options.verb = 0;
options.tau = 1/8;
% perform the regularization
[Mtv,err,tv,lambda] = perform_tv_denoising(zeros(n),options);

%%
% We can display the results.

clf;
imageplot({M clamp(Mtv)}, {'Original' 'Recovery'});

%%
% We can display the evolution of the regularization
% parameter and the TV energy.

clf;
subplot(2,1,1);
h = plot(lambda); axis tight;
set(h, 'LineWidth', 2);
set(gca, 'FontSize', 20);
title('Evolution of \lambda');
subplot(2,1,2);
lagr = .5*err.^2 + lambda(end)*tv;
h = plot(lagr); axis tight;
set(h, 'LineWidth', 2);
set(gca, 'FontSize', 20);
title('Evolution of |y-U x|^2+\lambda |x|_1');
axis([1 options.niter min(lagr) min(lagr)*1.5]);

%% Sparse spikes deconvolution with L1 regularization and Matching Pursuit
% sparse spikes deconvolution corresponds to inverting a band pass filter.
% This is very much used in seismic imaging where the filter is called a
% wavelet. Since the signal (reflecticity of the ground) is composed of
% isolated spikes, one uses L1 minimisation in the diract basis to recover
% the signal. This can also be performed by other sparse optimization
% technics such as greedy matching pursuit. Note that the operator |U| is
% equal to the dictionary |D| in this case, and is composed of translated
% copy of the filter. Since this results in a highly ill posed operator,
% one needs to sub-sample the dictionary a little.

%%
% First we load a filter which is a second derivative of a Gaussian

clear options;
n = 1024;
options.sigma = .02;
h = compute_sparse_spike_filter('dergauss',n,options);
h = h/max(h);

%%
% We can display the filter and its fourier transform (to see how much
% frequencies are removed by the filtering).

hf = real(fft(h)); hf = hf / max(hf);
clf;
subplot(2,1,1);
hh = plot(fftshift(h)); axis tight;
set(hh, 'LineWidth', 2);
set(gca, 'FontSize', 20);
title('Filter');
subplot(2,1,2);
hh = plot(fftshift(hf)); axis tight;
set(hh, 'LineWidth', 2);
set(gca, 'FontSize', 20);
title('Fourier Transform');

%%
% Then we build a dictionary of sub-sampled translated copies of |h|.

q = 2; % sub-sampling (distance between wavelets)
p = n/q;
[Y,X] = meshgrid(1:q:n,1:n);
D = h(mod(X-Y,n)+1);

%%
% We can display a sub-set of elements of the dictionary.

clf;
sel = round( linspace(1,p,5) ); sel(end) = [];
hh = plot(D(:,sel)); axis tight;
set(hh, 'LineWidth', 2);
set(gca, 'FontSize', 20);
title('Translated wavelets');

%%
% Joel Tropp and Jean-Jacques Fuchs have derived in several IEEE Tr. Info.
% Theory papers some powerfull bounds of the recovery of signals with L1
% minimization. This includes the Fuchs (which depends on the sign of the coefficients),
% ERC and WERC (which depend only on the support of the coefficients) criterions.
% We display here these three criterions to see wether a Dirac spikes train
% of spacing |Delta| between two consecutive spikes can be recovered stably
% by both matching
% pursuit and L1 minimization (basis pursuit).
% The dashed vertical lines display the critical spacing as detected by
% these three criterions.

options.display = 1;
options.delta_max = 100;
options.verb = 0;
[minscale,crit,lgd] = compute_minimum_scale(D, options);
set(gca, 'FontSize', 20);
xlabel('\Delta');

%%
% We create a sparse spikes train, with a default in the middle (two close
% spikes). Note that we select the spacing between the spikes slightly below the minimum scale detected
% by ERC<1, which causes sparsity methods to become unstable, and matching pursuit to fail.

options.delta = 23; % regular spacing
options.delta0 = 3; % create default that is not recovered by BP
x = compute_rand_sparse(p, 0, 'seismic', options);

%%
% The signal is the convolution with the filters plus noise.

sigma = .03;
y = D*x + randn(n,1)*sigma;
% display
clf;
subplot(2,1,1);
plot_sparse_diracs(x);
title('Sparse spikes x');
subplot(2,1,2);
hh = plot(y); axis tight;
set(hh, 'LineWidth', 2);
set(gca, 'FontSize', 20);
title('Signal y=Dx+w');

%%
% We can recover the signal using L1 minimisation plus debiasing.

options.T = 1;
options.niter = 1000;
options.etgt = 5*sigma*sqrt(n); % noise level
options.x = zeros(p,1); % initial solution
[xlun,err,lun,Tlist]  = perform_iterative_thresholding(D,y,options);
% perform debiasins
xlun(abs(xlun)<.08)=0;
xdeb = perform_debiasing(D,xlun,y);
% display
options.title = {'Signal' 'L1' 'Debiased'};
plot_sparse_diracs({x xlun xdeb}, options);

%%
% Another way to solve the problem is using greedy matching pursuit
% procedure. Most of the time, it gives results comparable or a little worse
% than basis pursuit, but it is faster. Here, we are in a critical set up
% where basis pursuit performs better than matching pursuit. The toolbox
% implements basis matching pursuit and its orthogonalized version. Only
% the basic pursuit is implemented with implicit operator, and orthogonalized
% pursuit requires the computation of many pseudo inverse (so it is quite slow).

% parameters for the pursuits
options.tol = 3*norm(w)/norm(y);
options.nbr_max_atoms = 80;
xmp = perform_mp(D,y,options);
xomp = perform_omp(D,y,options);
% display
options.title = {'Signal' 'Matching pursuit' 'Orthogonal Matching Pursuit'};
plot_sparse_diracs({x xmp xomp}, options);

%% Image inpainting with sparse wavelets regularization
% Image inpainting consists in filling missing pixels in an image. Many inpainting
% algorithm rely on a diffusion PDE that try to propagate information from
% the boundary of the missing region. An alternative methods consist in
% treating this problem as an inverse problem. Inpainting indeed
% corresponds to a diagonal operator |U| with 0 and 1 on the diagonal (0
% indicating missing pixels). This can be regularized using iterative
% thresholding. Since there is no noise, we want to use a very low
% threshold. This can be achieved by decaying the threshold during the
% iterations, and corresponds to the Morphological Component Analysis
% strategy of Jean-Luc Starck, Miki Elad and David Donoho.
% See also the homepage of Jalal Fadili for many example of inpainting.

%%
% First we load an image and a mask.

clear options;
n = 128;
M = load_image('lena');
M = rescale( crop(M,n) );
options.rho = .7; % remove 70% of pixels
mask = compute_inpainting_mask('rand',n,options);
options.mask = mask;

%%
% Apply the operator to remove pixels.

y = callback_inpainting(M, +1, options);
imageplot({M y}, {'Image' 'Data to inpaint'});

%%
% First we set up the dictionary: a translation invariant wavelet
% dictionary.

options.Jmin = 4;
options.wavelet_type = 'biorthogonal';
options.wavelet_vm = 3;
options.D = @callback_atrou;

%%
% Then we set up decaying thresholding parameter for the MCA solver.

options.thresh = 'soft';
options.Tmax = .1;
options.Tmin = 0;
options.niter = 200;
options.tau = 1;
options.drawiter = 0;
options.verb = 0;

%%
% Perform the inpainting using sparse wavelet expansion.

% do the resolution
[MW,err,lun,Tlist] = perform_iterative_thresholding(@callback_inpainting, y, options);
% retrieve the image from its coefficients
Mlun = callback_atrou(MW, +1, options);

%%
% Display the result

imageplot({M clamp(Mlun)},{'Original' 'Recovered'});

%% Dictionary Learning
% Instead of using a fixed dictionary |D| to compress or denoise a signal
% with sparse coding (or to solve an inverse problem), it is possible to
% optimize |D| in order to sparsify a set of given exemplar. In practice,
% these exemplar are taken as small patches extracted from tons of natural
% images. As first shown by Olshausen and Field, the optimized atoms then
% ressemble those of an oriented (steerable) wavelet frame.

%%
% The parameter of the dictionary

% size of the patches
w = 12; n = w^2;
% redundancy of the dictionary
redun = 1.5;
% number of atoms
p = round( n*redun );
% overtraining factor
overtraining = 6;
% number of example
m = round( overtraining*p );

%%
% First we load a set of images
Mlist = load_image({'lena' 'barb' 'boat'});

%%
% Then we extract a large set of patches from these images

Y = load_random_patches(Mlist,w,m);

%%
% Parameters for the learning of the dictionary (optimization of the atoms
% using an iterative method).

% number of atoms
options.K = p;
% number of iterations
options.niter_learning = 60;
% solver used for the sparse coding stage
options.sparse_coder = 'omp';
options.sparse_coder = 'mp';
% sparsity targeter for the sparse code
options.nbr_max_atoms = 5;
% initialization method
options.options.init_dico = 'input';

%%
% Perform the learning using an optimization procedure called the MOD
% (Method of Direction) algorithm.

options.learning_method = 'mod';
options.verb = 0;
[D1,X1,E1] = perform_dictionary_learning(Y,options);

%%
% Display the vectors of the dictionary as small images.
% You can notice that the texture of barbara is occupating many atoms !

clf;
options.ndim = 2;
options.normalization = 'clamp';
display_dictionnary(D1, X1, [10 14], options );

author: gpeyre

repo: matlab-toolboxes


text:
function test_issue856

% MEM 3gb
% WALLTIME 00:10:00
% DEPENDENCY ft_determine_coordsys ft_convert_coordsys ft_determine_units ft_convert_units ft_transform_geometry

%%

mri0 = ft_read_mri(dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/Subject01.mri'));
grad0 = ft_read_sens(dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/Subject01.ds'), 'senstype', 'meg');
headmodel0 = ft_read_headmodel(dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/Subject01.hdm'));

ft_hastoolbox('spm12', 1);

%%
% these have mixed units

mri1 = ft_determine_units(mri0);
grad1 = ft_determine_units(grad0);
headmodel1 = ft_determine_units(headmodel0);

%%
% convert to meter

target = 'm';
mri2 = ft_convert_units(mri1, target);
grad2 = ft_convert_units(grad1, target);
headmodel2 = ft_convert_units(headmodel1, target);

assert(~isequal(mri1.transform, mri2.transform));
assert(~isequal(grad1.coilpos, grad2.coilpos));
assert(~isequal(headmodel1.r, headmodel2.r));

%%
% convert to acpc

method = 0; % only approximate
mri3 = ft_convert_coordsys(mri2, 'acpc', method);
grad3 = ft_convert_coordsys(grad2, 'acpc', method);
headmodel3 = ft_convert_coordsys(headmodel2, 'acpc', method);

%%
% convert to millimeter

target = 'mm';
mri2 = ft_convert_units(mri1, target);
grad2 = ft_convert_units(grad1, target);
headmodel2 = ft_convert_units(headmodel1, target);

%%
% spm_affreg
mri4 = ft_convert_coordsys(mri2, 'acpc', 1);

% this round-trip is not exactly identical
mri4b = ft_convert_coordsys(mri4, 'ctf', 0);
assert(~isequal(mri2.transform, mri4b.transform));

%%
% smp_normalise
mri5 = ft_convert_coordsys(mri2, 'acpc', 2);

% this round-trip is not exactly identical
mri5b = ft_convert_coordsys(mri5, 'ctf', 0);
assert(~isequal(mri2.transform, mri5b.transform));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function pass = test_stringParser(pref)
% This test ensures that the STRINGPARSER class is doing the correct thing in a
% few different situations.

%% BVPS

% Strings we want to test.
str = { 'u(0) = 1'
        'u(-1) = 1'
        'u(0)=1, u(2) = 3'
        'u = 1,v = 0,w = 3'
        'u(end),u(1,left)'
        'diff(sin(u))'
        'diff(sin(u),2)'
        'diff(u^2) + sin(x)*u = 1'
        'feval(u,sqrt(2))'
        'feval(u,0,left)'
        'sum(u) = 0'
        'sum(u,0,.5) = 0'
        'volt(sin(x-y),u) = 0'
        'fred(sin(x-y),u) = fred(cos(x-z),u)'
        'u + fred(sin(2*pi*(y-x)),u) = feval(u,0,left)'};

% Store the number of strings
n = numel(str);
outBVP = cell(n, 1);

for k = 1:n
    outBVP{k} = stringParser.str2anon(str{k}, 'bvp');
end

correctBVP = {  '@(u) feval(u,0)-1'
                '@(u) feval(u,-1)-1'
                '@(u) feval(u,0)-1;feval(u,2)-3'
                '@(u,v,w) [u-1;v;w-3]'
                '@(u) feval(u,''end'');feval(u,1,''left'')'
                '@(u) diff(sin(u))'
                '@(u) diff(sin(u),2)'
                '@(x,u) diff(u.^2)+sin(x).*u-1'
                '@(u) feval(u,sqrt(2))'
                '@(u) feval(u,0,''left'')'
                '@(u) sum(u)'
                '@(u) sum(u,0,.5)'
                '@(x,u) volt(@(x,y)sin(x-y),u)'
                '@(x,u) fred(@(x,y)sin(x-y),u)-fred(@(x,z)cos(x-z),u)'
                '@(x,u) u+fred(@(x,y)sin(2.*pi.*(y-x)),u)-feval(u,0,''left'')'};

passBVP = zeros(n, 1);
for k = 1:n
    passBVP(k) = strcmp(outBVP{k}, correctBVP{k});
end

%% EIGS

str = {'u''''-lambda*u = 0'
       'u'' + sin(x)*u = lambda*u'
       'u''''+u'' = lambda*(u + u'') + u'};
n2 = numel(str);
outEIG = cell(n2, 1);

for k = 1:n2
    outEIG{k} = stringParser.str2anon(str{k}, 'eig');
end

correct = {{'@(u) diff(u,2)'
            '@(u) u'}
            {'@(x,u) diff(u)+sin(x).*u'
            '@(x,u) u'}
            {'@(u) diff(u,2)+diff(u)-u'
            '@(u) u+diff(u)'}};

passEIG = zeros(n2, 1);
for k = 1:n2
    for j = 1:numel(outEIG{k})
        passEIG(k) = strcmp(outEIG{k}{j}, correct{k}{j});
    end
end

%% PDES

str = {'u_t+x*u'''' = u'''};
n3 = numel(str);

outPDE = cell(n3, 1);
for k = 1:n3
    outPDE{k} = stringParser.str2anon(str{k}, 'pde');
end

correct = {'@(x,t,u) -x.*diff(u,2)+diff(u)'};

passPDE = zeros(n3, 1);
for k = 1:n3
    passPDE(k) = strcmp(outPDE{k}, correct{k});
end

%% Concatenate all types
pass = [passBVP; passEIG; passPDE];

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% LSHAPEAFEMQUADCURL quad curl equations on the Lshaped domain
%
%   LSHAPEAFEMQUADCURL computes ND0-(CR-P0)-(linear ND) nonconforming approximations
%   of the quad curl equations in the unit cube. The mesh is refined
%   adaptively guided by a residual-based estimator using a separate marking
%   strategy with two marking parameters for eta_1 which estimates the error
%   for w and eta_2 for phi.
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

close all; clear;

%% Set up
maxN = 6e5;
maxIt = 24;
theta = 0.3;
theta1 = 0.5;
nDofMaxwell = zeros(maxIt,1);
nDofStokes = zeros(maxIt,1);
h = zeros(maxIt,1);
erru = zeros(maxIt,1);
erruL2 = zeros(maxIt,1);
errw = zeros(maxIt,1);
errphiL2 = zeros(maxIt,1);
errphi = zeros(maxIt,1);
etaTotal = zeros(maxIt,1);
etaw = zeros(maxIt,1);
etaphi = zeros(maxIt,1);
etau = zeros(maxIt,1);

%% Generate initial mesh
load meshLshape3.mat % a minimal Lshape mesh
bdFlag = setboundary3(node,elem,'Dirichlet');
for i = 1
    [node,elem,bdFlag] = uniformrefine3(node,elem,bdFlag);
end
[elem,bdFlag,HB] = label3(node,elem,'all',bdFlag);

%% PDE and options
gamma = 3+1/6; alpha = 2/3;
% gamma = 8/3; alpha = 2/3;
% gamma = 5/2; alpha = 1/2;
% gamma = 7/3; alpha = 1/3; % extreme case
pde = quadCurlDataLshape1(gamma,alpha);
% pde = quadCurlDataLshape2(gamma,alpha);
% pde = quadCurlDataLshape(gamma,alpha);
option.printlevel = 0;

%% AFEM cycles using only 1st Maxwell+ Stokes, 2nd Maxwell solved only for ref
for k = 1:maxIt
    fprintf('\n\n\nRefinement level %d\n', k);
    NT = size(elem,1);

    %% SOLVE the 1st Maxwell problem

    pdeCurl1.J = pde.quadcurlu;
    pdeCurl1.g_D = @(p) zeros(size(p,1),3);
    option.solver = 'direct';
    if NT > 2e4; option.solver = 'diag'; end
    fprintf('\n**************Solving the first Maxwell problem for w**************');
    soln.w = Maxwellsaddle(node,elem,bdFlag,pdeCurl1,option,HB);
%     soln.w = Maxwellsaddle(node,elem,bdFlag,pdeCurl1,option);
    fprintf('*********************************************************************\n');
    %% SOLVE the Stokes problem

    % RHS of Stokes
    [curlwh, volume, curlbasis] = curlu3(node,elem,soln.w);
    pdeStokes.f = curlwh;
    pdeStokes.g_D = pde.curlu;
    option.solver = 'diag';
    option.printlevel = 1;
    fprintf('\n**************Solving the Stokes problem for phi*******************\n');
    solnStokes = Stokes3CRP0(node,elem,bdFlag,pdeStokes,option,HB);
    fprintf('*********************************************************************\n');
    %% SOLVE the 2nd Maxwell problem (not used in the indicator for mesh refinement)

    % RHS os the 2nd Maxwell
    [elem2face,face] = dof3face(elem);
    NF = size(face,1);
    soln.phi = reshape(solnStokes.u,NF,3);

    phihCenter = zeros(NT,3); % at each elem center
    for j = 1:3 % each component
        phihj = soln.phi(:,j);
        phihj2elem = phihj(elem2face);
        phihCenter(:,j) = sum(phihj2elem,2)/4;
        % Crouzeix-Raviart shape function 1-3\lambda = 1/4 at center
    end

    [elem2edge,edge] = dof3edge(sort(elem,2)); % sort for edge element
    NE = size(edge,1);
    bt = zeros(NT,6);
    for j = 1:6 % 6 basis of the edge element
        bt(:,j) = dot(curlbasis(:,:,j),phihCenter,2).*volume;
    end
    fCurl = zeros(2*NE,1);
    fCurl(1:NE,:) = accumarray(elem2edge(:),bt(:),[NE 1]);

    % solve the last Maxwell problem
    option.solver = 'direct';
    if NE > 2e4; option.solver = 'diag'; end
    pdeCurl2.J = fCurl;
    pdeCurl2.g_D = pde.exactu;
    pdeCurl2.g = pde.g;
    fprintf('\n*****Solving the second Maxwell problem for u (for reference)*****');
    soln.u = Maxwell1saddle(node,elem,bdFlag,pdeCurl2,option,HB);
    fprintf('*********************************************************************\n');
    %% ESTIMATE
    soln.curlw = curlwh;
    nDofMaxwell(k) = 2*NE;
    nDofStokes(k) = 3*NF;

    option.includeU = true; % only use the indicator by w_h and phi_h
    option.scaling = false; % no scaling and separate marking
    [etaK, est] = estimatequadcurl(node,elem,soln,pde,option);

    hK = est.hK;
%     h = sqrt(sum(hK.^2)/NT);
    Dlambda = gradbasis3(node,elem);
    curlphi = curlu3CR(elem2face,soln.phi,Dlambda);
    [errphi(k), errphiK] = getL2error3(node,elem,pde.curlcurlu,curlphi);
    [errphiL2(k), errKphiL2] = getL2error3(node,elem,pde.curlu,soln.phi);
    etaTotal(k) = sqrt(sum(etaK.^2));
    curluh = curlu3(node,elem,soln.u);
    [erru(k), errK] = getL2error3(node,elem,pde.curlu,curluh);
    erruL2(k) = getL2error3ND1(node,elem,pde.exactu,soln.u);
%     est.eta1 = est.eta1.*hK;
%     est.eta1 = est.eta1*min(hK);
    etaw(k) = sqrt(sum((est.eta1).^2));
    etaphi(k) = sqrt(sum((est.eta2).^2));
    etau(k) = sqrt(sum((est.eta3).^2));

    fprintf('Total eta          = %g \n', etaTotal(k));
%     fprintf('Mesh size          = %g \n', max(hK));
    fprintf('# face             = %g \n', NF);
    fprintf('# edge             = %g \n', NE);
    fprintf('||phi_h - curl u|| = %g \n', errphiL2(k));
    %% Visualize
    figure(1);
    subplot(1,3,1)
    etawE2V = accumarray(elem(:),repmat(est.eta1,[4,1]),[size(node,1), 1]);
    showsolution3(node,elem,etawE2V,'z<0.25 ','EdgeColor','k','facecolor','flat');
    colorbar;
    title('eta for w')
    view(120,45)

    subplot(1,3,2)

    % errphiL2E2V = accumarray(elem(:),repmat(errphiL2K,[4,1]),[size(node,1), 1]);
    % showsolution3(node,elem,errphiL2E2V,'z<0.25 ','EdgeColor','k','facecolor','flat');

    etaphiE2V = accumarray(elem(:),repmat(est.eta2,[4,1]),[size(node,1), 1]);
    showsolution3(node,elem,etaphiE2V,'z<0.25 ','EdgeColor','k','facecolor','flat');
    colorbar;
    title('eta for phi')
    view(120,45)

    subplot(1,3,3)
%     errU2v = accumarray(elem(:),repmat(erruK,[4,1]), [size(node,1), 1]);
%     showsolution3(node,elem,errU2v,'z<0.25','EdgeColor','k','facecolor','flat');
    % errL2U2v = accumarray(elem(:),repmat(abs(erruL2K),[4,1]),[size(node,1), 1]);
    % showsolution3(node,elem,errL2U2v,'z<0.25 ','EdgeColor','k','facecolor','flat');
    errphiE2V = accumarray(elem(:),repmat(errphiK,[4,1]),[size(node,1), 1]);
    showsolution3(node,elem,errphiE2V,'z<0.25 ','EdgeColor','k','facecolor','flat');
    title('H^1 error for phi')
    colorbar;
    view(120,45)

    set(gcf,'color','w','Position', [100 600 1400 300])
    drawnow;

    %% MARK and REFINE
    if nDofStokes(k) > maxN; break; end
%
%     markedElem = mark(elem,etaK,theta);
%     markedElemC = mark(elem,est.elemRescurlwh,2/3*theta);
%     markedElem = unique([markedElem; markedElemC]);

    % markedElem = mark(elem,errKphiL2,theta); % cheating
    markedElem = mark2(elem,est.eta1,theta1,est.eta2,theta);

    fprintf('eta1(marked)^2/eta1^2 = %g \n', sum(est.eta1(markedElem).^2)/sum(est.eta1.^2));
    fprintf('eta2(marked)^2/eta2^2 = %g \n', sum(est.eta2(markedElem).^2)/sum(est.eta2.^2));
    if k < maxIt
        [node,elem,bdFlag,HB] = bisect3(node,elem,markedElem,bdFlag,HB);
    end

end

%%
close all;
figure(1);
nDofMaxwell = nDofMaxwell(1:k);
nDofStokes = nDofStokes(1:k);
erru = erru(1:k);
errphiL2 = errphiL2(1:k);
etaTotal = etaTotal(1:k);
erruL2 = erruL2(1:k);
etaw = etaw(1:k);
etaphi = etaphi(1:k);
etau = etau(1:k);
errphi = errphi(1:k);
ITER_THRESH = ceil(k/2);
r1 = showrate(nDofMaxwell,erru,ITER_THRESH,'k-o');
r2 = showrate(nDofStokes,errphiL2,ITER_THRESH,'r-o');
% r3 = showrate(nDofMaxwell+nDofStokes,etaTotal,ITER_THRESH,'b-*');
r3 = showrate(nDofMaxwell,erruL2,ITER_THRESH,'b-*');
r4 = showrate(nDofMaxwell/2,etaw,ITER_THRESH,'color',[0.7,0.2,0.7],'marker','*');
r5 = showrate(nDofStokes,etaphi,ITER_THRESH,'color',[0.2,0.8,0.7],'marker','*');
r6 = showrate(nDofMaxwell,etau,ITER_THRESH,'color',[0.7,0.8,0.2],'marker','*');
r7 = showrate(nDofStokes,errphi,ITER_THRESH,'g-o');
set(gca,'xscale','log','yscale','log');
T = title('Convergence');
XL = xlabel('$\#$ DoFs');
YL = ylabel('Error Magnitudes');
L = legend('$\Vert  \nabla\times(u-u_h)\Vert$', ...
    ['$(\# \mbox{edge})^{' num2str(r1) '}$'],...
    '$\Vert  \nabla\times u - \phi_h\Vert$', ...
    ['$(\# \mbox{face})^{' num2str(r2) '}$'],...
        '$\Vert u - u_h\Vert$', ...
    ['$(\# \mbox{edge})^{' num2str(r3) '}$'],...
     '$\eta_1(w_h)$', ['$(\# \mbox{edge})^{' num2str(r4) '}$'],...
    '$\eta_2(w_h, \phi_h)$', ['$(\# \mbox{face})^{' num2str(r5) '}$'],...
    '$\eta_3(\phi_h, u_h)$', ['$(\# \mbox{edge})^{' num2str(r6) '}$'],...
    '$|\nabla\times u-\phi_h|_{1,h}$', ['$(\# \mbox{face})^{' num2str(r7) '}$'],...
'LOCATION','southwest');

%     '$\eta(w_h,\phi_h)$', ...
%     ['$(\# \mbox{DoF})^{' num2str(r3) '}$'],...

set([XL,YL,L],'Interpreter','latex','FontSize', 16);
set(T,'Interpreter','latex','FontSize',20);
set(gca,'TickLabelInterpreter', 'latex');
grid on;
% set(gcf,'color','w','Position', [100 200 300 600])
set(gcf,'color','w','Position',[613 598 820 820])
drawnow;

%% sanity check
if nDofMaxwell<1e4
    figure(2);
    subplot(1,2,1);
    option.scale = 2;
    option.plot = 'quiver3';
    center = (node(elem(:,1),:) + node(elem(:,2),:) + ...
        node(elem(:,3),:) + node(elem(:,4),:))/4;
    showvector3(center,pde.curlu(center),option);
    subplot(1,2,2);
    curluh = curlu3(node,elem,soln.u);
    showvector3(center,curluh,option);

    set(gcf,'color','w','Position', [500 600 750 300])
    drawnow;
end

repo: ifem


text:
function [stan,rtfr,hat] = tfrrstan(x,t,N,G,h,trace);
%TFRRSTAN Reassigned Stankovic distribution.
%	[TFR,RTFR,HAT] = TFRRSTAN(X,T,N,H,TRACE)
%	computes the Stankovic distribution and its reassigned version.
%
%	X     : analysed signal.
%	T     : the time instant(s)      (default : 1:length(X)).
%	N     : number of frequency bins (default : length(X)).
%       G     : frequency averaging window
%	h     : stft window              (default : Hamming(N/4)).
%	TRACE : if nonzero, the progression of the algorithm is shown
%	                                 (default : 0).
%	TFR,  : time-frequency representation and its reassigned
%	RTFR    version. When called without output arguments,
%	        TFRRSTAN runs TFRQVIEW.
%	HAT   : Complex matrix of the reassignment vectors.
%
%	Example :
%	 sig=fmlin(128,0.1,0.4); t=1:2:128;
%	 G=tftb_window(9,'hanning'); h=tftb_window(61,'hanning'); tfrrstan(sig,t,128,G,h,1);
%
%	See also  all the time-frequency representations listed in
%	 the file CONTENTS (TFR*)

%	F. Auger, August, September 1997.
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if (nargin == 0),
 error('At least 1 parameter required');
end;
[xrow,xcol] = size(x);
if (nargin <= 2),
 N=xrow;
end;

hlength=floor(N/4);
hlength=hlength+1-rem(hlength,2);

if (nargin == 1),
 t=1:xrow; G=[0.25; 0.5; 0.25]; h = tftb_window(hlength); trace=0;
elseif (nargin == 2)|(nargin == 3),
 G=[0.25; 0.5; 0.25]; h = tftb_window(hlength); trace=0;
elseif (nargin == 4),
 h = tftb_window(hlength); trace = 0;
elseif (nargin == 5),
 trace = 0;
end;

[Grow,Gcol]=size(G); LG=(Grow-1)/2;
if (Gcol~=1)|(rem(Grow,2)==0),
 error('G must be a smoothing window with odd length');
end;

if (N<0),
 error('N must be greater than zero');
end;
[trow,tcol] = size(t);
if (xcol~=1),
 error('X must have only one column');
elseif (trow~=1),
 error('T must only have one row');
elseif (2^nextpow2(N)~=N),
 fprintf('For a faster computation, N should be a power of two\n');
end;

[hrow,hcol]=size(h); Lh=(hrow-1)/2;
if (hcol~=1)|(rem(hrow,2)==0),
 error('H must be a smoothing window with odd length');
end;

if (tcol==1),
 Dt=1;
else
 Deltat=t(2:tcol)-t(1:tcol-1);
 Mini=min(Deltat); Maxi=max(Deltat);
 if (Mini~=Maxi),
  error('The time instants must be regularly sampled.');
 else
  Dt=Mini;
 end;
 clear Deltat Mini Maxi;
end;

stan = zeros(N,tcol);
rtfr = zeros(N,tcol);
hat  = zeros(N,tcol);

if trace, disp('Stankovic distribution (with reassignement)'); end;
Ex=mean(abs(x(min(t):max(t))).^2); Threshold=1.0e-3*Ex;
Dh=dwindow(h); Th=h.*[-Lh:Lh]';
for icol=1:tcol,
 if trace, disprog(icol,tcol,10); end;
 ti= t(icol);
 tau=-min([round(N/2)-1,Lh,ti-1]):min([round(N/2)-1,Lh,xrow-ti]);
 indices= rem(N+tau,N)+1;
 tfr= zeros(N,3);
 tfr(indices,1)=x(ti+tau).*conj( h(Lh+1-tau));
 tfr(indices,2)=x(ti+tau).*conj(Th(Lh+1-tau));
 tfr(indices,3)=x(ti+tau).*conj(Dh(Lh+1-tau));
 tfr=fft(tfr);
 stan(:,icol)=G(LG+1) * abs(tfr(:,1)).^2;

 for jcol=1:N,
  stanTh=G(LG+1)*tfr(jcol,1)*conj(tfr(jcol,2));
  stanDh=G(LG+1)*tfr(jcol,1)*conj(tfr(jcol,3));
  for kstan=1:min(N/2-1,LG),
   stanbefore=rem(rem(jcol-kstan-1,N)+N,N)+1;
   stanafter =rem(rem(jcol+kstan-1,N)+N,N)+1;
   stan(jcol,icol)= stan(jcol,icol) ...
                  + G(LG+1-kstan)*tfr(stanbefore,1)*conj(tfr(stanafter ,1)) ...
                  + G(LG+1+kstan)*tfr(stanafter ,1)*conj(tfr(stanbefore,1));
   stanTh= stanTh + G(LG+1-kstan)*tfr(stanbefore,1)*conj(tfr(stanafter ,2)) ...
                  + G(LG+1+kstan)*tfr(stanafter ,1)*conj(tfr(stanbefore,2));
   stanDh= stanDh + G(LG+1-kstan)*tfr(stanbefore,1)*conj(tfr(stanafter ,3)) ...
                  + G(LG+1+kstan)*tfr(stanafter ,1)*conj(tfr(stanbefore,3));
  end;
  stan(jcol,icol)=real(stan(jcol,icol));
  if abs(stan(jcol,icol))>Threshold,
   icolhat = round(icol - real(stanTh/stan(jcol,icol))/Dt);
   jcolhat = round(jcol - N*imag(stanDh/stan(jcol,icol)/(2.0*pi)));

   jcolhat= rem(rem(jcolhat-1,N)+N,N)+1;
   icolhat= min(max(icolhat,1),tcol);

   rtfr(jcolhat,icolhat)=rtfr(jcolhat,icolhat) + stan(jcol,icol) ;
   hat(jcol,icol)= jcolhat + j * icolhat;
   %fprintf('%12.3f %12.3f , %12.3f %12.3f \n',jcol,icol,jcolhat,icolhat);
  else
   rtfr(jcol,icol)=rtfr(jcol,icol) + stan(jcol,icol);
   hat(jcol,icol)= jcol + j * icol;
  end;
 end;

end ;

if trace, fprintf('\n'); end;

if (nargout==0),
 TFTBcontinue=1;
 while (TFTBcontinue==1),
  choice=menu ('Choose the representation:',...
               'stop',...
               'Stankovic distribution',...
               'reassigned Stankovic distribution');
  if (choice==1), TFTBcontinue=0;
  elseif (choice==2),
   tfrqview(stan,x,t,'type1');
  elseif (choice==3),
   tfrqview(rtfr,x,t,'type1');
  end;
 end;
end;

author: HeLiangHIT


text:
function runTracker_VOT(hp_name, hp_value)
%RUN_TRACKER_VOT executes the tracker within the evaluation of the VOT toolkit

    %% Read params.txt
	if nargin>0
		% hyperparams mode
	    params = readParams('params.txt', hp_name, hp_value);
	else
		% normal mode
		params = readParams('params.txt');
	end

    %% Read files
    fid = fopen('images.txt','r');
    inter = textscan(fid,'%[^\n]');
    params.img_files = inter{1,1};
    im = imread(params.img_files{1});
    fclose(fid);

    % check if sequence is grayscale
    if(size(im,3)==1)
        params.grayscale_sequence = true;
    end

    params.bb_VOT = csvread('region.txt');
    region = params.bb_VOT(1,:);
    params.img_path = '';
    % save sequence name
    C = strsplit(params.img_files{1},'/');
    if ~isempty(C)
        params.sequence = C{end-1};
    else
        params.sequence = 'error_sequence_not_detected';
    end
    %% Convert to axis-aligned bbox
    if(numel(region)==8)
        % polygon format
        [cx, cy, w, h] = getAxisAlignedBB(region);
    else
        x = region(1);
        y = region(2);
        w = region(3);
        h = region(4);
        cx = x+w/2;
        cy = y+h/2;
    end

    % init_pos is the centre of the initial bounding box
    params.init_pos = [cy cx];
    params.target_sz = round([h w]);

    [params, bg_area, fg_area, area_resize_factor] = initializeAllAreas(im, params);

    params.fout = fopen('output.txt','w');

    % start actual tracking
    trackerMain(params,im, bg_area, fg_area, area_resize_factor);
    fclose('all');
    % within VOT benchmark the tracker has to terminate at the end of each sequence
    exit
end

author: bertinetto

repo: staple

save_path: github-repos/MATLAB/bertinetto-staple


text:
   function y = autopilot_quad(drone, t)
    % AUTOPILOT_QUAD - autopilot used for the quadcopter

    p_drone = drone.p_drone;
    p_sim = drone.p_sim;
    autopilot_version = drone.autopilot_version;

    % State variables
    pn = drone.z_hat(1); % position north
    pe = drone.z_hat(2); % position east
    h = drone.z_hat(3); % altitude
    vx = drone.z_hat(4); % x velocity
    vy = drone.z_hat(5); % y velocity
    vz = drone.z_hat(6); % z velocity
    Va = drone.z_hat(7); % airspeed
    phi = drone.z_hat(10); % roll angle
    theta = drone.z_hat(11); % pitch angle
    psi = drone.z_hat(12); % yaw angle
    p = drone.z_hat(14); % body frame roll rate
    q = drone.z_hat(15); % body frame pitch rate
    r = drone.z_hat(16); % body frame yaw rate

    %%%%% Autopilot version %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % autopilot_version == 1 <- attitude controller
    % autopilot_version == 2 <- speed controller
    % autopilot_version == 3 <- acceleration controller
    % TO IMPLEMENT: autopilot_version == 4 <- position controller

    % Commanded values
    if autopilot_version == 1% attitude autopilot
        h_c = drone.command(1); % commanded altitude        [m]
        phi_c = drone.command(2); % commanded roll          [rad]
        theta_c = drone.command(3); % commanded pitch       [rad]
        psi_c = drone.command(4); % commanded yaw           [rad]

    elseif autopilot_version == 2% velocity autopilot
        psi_c = drone.command(1); % commanded yaw            [rad]
        vn_c = drone.command(2); % commanded north velocity  [m/s]
        ve_c = drone.command(3); % commanded east velocity   [m/s]
        vd_c = drone.command(4); % commanded down velocity   [m/s]

    elseif autopilot_version == 3% acceleration autopilot
        psi_c = drone.command(1); % commanded yaw            [rad]
        an_c = drone.command(2); % commanded north accell    [m/s^2]
        ae_c = drone.command(3); % commanded east accell     [m/s^2]
        ad_c = drone.command(4); % commanded down accell     [m/s^2]

    else %if autopilot_version == 4% position autopilot
        psi_c = drone.command(1); % commanded yaw            [rad]
        pn_c = drone.command(2); % commanded north position  [m]
        pe_c = drone.command(3); % commanded east position   [m]
        pd_c = drone.command(4); % commanded down position   [m]
    end

    %Forces
    F_xyz = drone.forces;

    a_ned = rotate_b2i(F_xyz, phi, theta, psi) / p_drone.mass;
    an = a_ned(1);
    ae = a_ned(2);
    ad = a_ned(3);

    switch autopilot_version
        case 1
            [delta, x_command] = autopilot_attitude(h_c, phi_c, theta_c, psi_c, ...
                h, phi, theta, psi, p, q, r, t, p_drone, p_sim, drone);

        case 2
            [delta, x_command] = autopilot_velocity(psi_c, vn_c, ve_c, vd_c, ...
                vx, vy, vz, phi, theta, psi, ...
                p, q, r, h, t, p_drone, p_sim, drone);
        case 3
            [delta, x_command] = autopilot_acceleration(psi_c, an_c, ae_c, ad_c, ...
                an, ae, ad, phi, theta, psi, ...
                p, q, r, h, t, p_drone, p_sim, drone);
        case 4
            [delta, x_command] = autopilot_position(psi_c, pn_c, pe_c, pd_c, ...
                vx, vy, vz, phi, theta, psi, ...
                p, q, r, h, pn, pe, t, p_drone, p_sim, drone);
    end

    y = [delta; x_command];

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           Autopilot versions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [delta, x_command] = autopilot_attitude(h_c, phi_c, theta_c, ...
        psi_c, h, phi, theta, psi, p, q, r, t, p_drone, p_sim, drone)
% AUTOPILOT_ATTITUDE

    if h <= p_drone.altitude_take_off_zone
        drone.altitude_state = 2;
    else
        drone.altitude_state = 2;
    end

    if t == 0 || (drone.altitude_state ~= drone.altitude_state_prev)
        drone.initialize_integrator = 1;
    else
        drone.initialize_integrator = 0;
    end

    % Implement state machine
    switch drone.altitude_state
        case 1% in take-off zone
            d1 = 1;
            d2 = 1;
            d3 = 1;
            d4 = 1;

        case 2% outside take-off zone
            [thrust, drone.P_h_thrust] = h_thrust_hold(h_c, h, ...
                drone.initialize_integrator, p_drone, p_sim, drone.P_h_thrust);
            [torque_phi, drone.P_roll_torque] = roll_torque_hold(phi_c, ...
                phi, p / 180 * pi, drone.initialize_integrator, p_drone, p_sim, ...
                drone.P_roll_torque);
            [torque_theta, drone.P_pitch_torque] = pitch_torque_hold(...
                theta_c, theta, q / 180 * pi, ...
                drone.initialize_integrator, p_drone, p_sim, drone.P_pitch_torque);
            [torque_psi, drone.P_psi_torque] = yaw_torque_hold(psi_c, ...
                psi, r / 180 * pi, drone.initialize_integrator, p_drone, p_sim, ...
                drone.P_psi_torque);

            [d1, d2, d3, d4] = mixer(thrust, torque_phi, torque_theta, ...
                torque_psi, p_drone, p_sim);

            % artificially saturation delta_t
            d1 = saturate(d1, p_drone.dt_max, 0);
            d2 = saturate(d2, p_drone.dt_max, 0);
            d3 = saturate(d3, p_drone.dt_max, 0);
            d4 = saturate(d4, p_drone.dt_max, 0);
    end

    % Control outputs
    delta = [d1; d2; d3; d4];

    % Commanded (desired) states
    x_command = [...
                0; ...% pn
            0; ...% pe
            h_c; ...% h
            0; ...% vn
            0; ...% ve
            0; ...% vd
            0; % an
            0; % ae
            0; % ad
            0; ...% Va
            0; ...% alpha
            0; ...% beta
            phi_c; ...% phi
            theta_c; % theta
            psi_c; ...% psi
            0; % chi
            0; ...% p
            0; ...% q
            0; ...% r
            ];

    y = [delta; x_command];
    drone.altitude_state_prev = drone.altitude_state;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [delta, x_command] = autopilot_velocity(psi_c, vn_c, ve_c, ...
        vd_c, vx, vy, vz, phi, theta, psi, p, q, r, h, t, p_drone, p_sim, drone)
% AUTOPILOT_VELOCITY

    Rbi = Rb2i(phi, theta, psi);
    v_ned = Rbi * [vx vy vz]';
    vn = v_ned(1);
    ve = v_ned(2);
    vd = v_ned(3);

    if h <= p_drone.altitude_take_off_zone
        drone.altitude_state = 2;
    else
        drone.altitude_state = 2;
    end

    if t == 0 || (drone.altitude_state ~= drone.altitude_state_prev)
        drone.initialize_integrator = 1;
    else
        drone.initialize_integrator = 0;
    end

    % Implement state machine
    switch drone.altitude_state
        case 1% in take-off zone
            d1 = 1;
            d2 = 1;
            d3 = 1;
            d4 = 1;
        case 2% outside take-off zone
            [thrust, drone.P_vd_thrust] = vd_thrust_hold(vd_c, vd, ...
                drone.initialize_integrator, p_drone, p_sim, drone.P_vd_thrust);
            [phi_c, drone.P_ve_roll] = ve_roll_hold(ve_c, ve, ...
                drone.initialize_integrator, p_drone, p_sim, drone.P_ve_roll);
            [torque_phi, drone.P_roll_torque] = roll_torque_hold(phi_c, ...
                phi, p / 180 * pi, drone.initialize_integrator, p_drone, p_sim, ...
                drone.P_roll_torque);
            [theta_c, drone.P_vn_pitch] = vn_pitch_hold(vn_c, vn, ...
                drone.initialize_integrator, p_drone, p_sim, drone.P_vn_pitch);
            [torque_theta, drone.P_pitch_torque] = pitch_torque_hold(...
                theta_c, theta, q / 180 * pi, ...
                drone.initialize_integrator, p_drone, p_sim, drone.P_pitch_torque);
            [torque_psi, drone.P_psi_torque] = yaw_torque_hold(psi_c, ...
                psi, r / 180 * pi, drone.initialize_integrator, p_drone, p_sim,...
                drone.P_psi_torque);

            [d1, d2, d3, d4] = mixer(thrust, torque_phi, torque_theta, ...
                torque_psi, p_drone, p_sim);

            % artificially saturation delta_t
            d1 = saturate(d1, p_drone.dt_max, 0);
            d2 = saturate(d2, p_drone.dt_max, 0);
            d3 = saturate(d3, p_drone.dt_max, 0);
            d4 = saturate(d4, p_drone.dt_max, 0);
    end

    % Control outputs
    delta = [d1; d2; d3; d4];

    % Commanded (desired) states
    x_command = [...
                0; ...% pn
            0; ...% pe
            0; ...% h
            vn_c; ...% vn
            ve_c; ...% ve
            vd_c; ...% vd
            0; % an
            0; % ae
            0; % ad
            0; ...% Va
            0; ...% alpha
            0; ...% beta
            phi_c; ...% phi
            theta_c; % theta
            psi_c; ...% psi
            0; % chi
            0; ...% p
            0; ...% q
            0; ...% r
            ];

    y = [delta; x_command];
    drone.altitude_state_prev = drone.altitude_state;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% // TODO: Check it
function [delta, x_command] = autopilot_acceleration(psi_c, an_c, ...
        ae_c, ad_c, an, ae, ad, phi, theta, psi, p, q, r, h, t, p_drone, p_sim, drone)
% AUTOPILOT_ACCELERATION

    if h <= p_drone.altitude_take_off_zone
        drone.altitude_state = 2;
    else
        drone.altitude_state = 2;
    end

    if t == 0 || (drone.altitude_state ~= drone.altitude_state_prev)
        drone.initialize_integrator = 1;
    else
        drone.initialize_integrator = 0;
    end

    % Implement state machine
    switch drone.altitude_state
        case 1% in take-off zone
            d1 = 1;
            d2 = 1;
            d3 = 1;
            d4 = 1;

        case 2% outside take-off zone
            [thrust, drone.P_ad_thrust] = ad_thrust_hold(ad_c, ad, ...
                drone.initialize_integrator, p_drone, drone.P_ad_thrust);
            %[thrust, drone.P_h_thrust] = h_thrust_hold(15, h, ...
                % drone.initialize_integrator, p_drone, drone.P_h_thrust);
            [phi_c, drone.P_ae_roll] = ae_roll_hold(ae_c, ae, ...
                drone.initialize_integrator, p_drone, drone.P_ae_roll);
            [torque_phi, drone.P_roll_torque] = roll_torque_hold(phi_c, ...
                phi, deg2rad(p), drone.initialize_integrator, p_drone, p_sim, ...
                drone.P_roll_torque);
            [theta_c, drone.P_an_pitch] = an_pitch_hold(an_c, an, ...
                drone.initialize_integrator, p_drone, drone.P_an_pitch);
            [torque_theta, drone.P_pitch_torque] = pitch_torque_hold(...
                theta_c, theta, deg2rad(q), ...
                drone.initialize_integrator, p_drone, drone.P_pitch_torque);
            [torque_psi, drone.P_psi_torque] = yaw_torque_hold(psi_c, ...
                psi, deg2rad(r), drone.initialize_integrator, p_drone, p_sim,  ...
                drone.P_psi_torque);

            [d1, d2, d3, d4] = mixer(thrust, torque_phi, torque_theta, ...
                torque_psi, p_drone, p_sim);

            % artificially saturation delta_t
            d1 = saturate(d1, p_drone.dt_max, 0);
            d2 = saturate(d2, p_drone.dt_max, 0);
            d3 = saturate(d3, p_drone.dt_max, 0);
            d4 = saturate(d4, p_drone.dt_max, 0);
    end

    % Control outputs
    delta = [d1; d2; d3; d4];

    % Commanded (desired) states
    x_command = [...
                0; ...% pn
            0; ...% pe
            0; ...% h
            0; ...% vn
            0; ...% ve
            0; ...% vd
            an_c;
            ae_c;
            ad_c;
            0; ...% Va
            0; ...% alpha
            0; ...% beta
            phi_c; ...% phi
            theta_c; % theta
            psi_c; ...% psi
            0; % chi
            0; ...% p
            0; ...% q
            0; ...% r
            ];

    y = [delta; x_command];
    drone.altitude_state_prev = drone.altitude_state;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [delta, x_command] = autopilot_position(psi_c, pn_c, pe_c, ...
        pd_c, vx, vy, vz, phi, theta, psi, p, q, r, h, pn, pe, t, p_drone, p_sim, drone)
% AUTOPILOT_POSITION
% // TODO: Check it
    Rbi = Rb2i(phi, theta, psi);
    v_ned = Rbi * [vx vy vz]';
    vn = v_ned(1);
    ve = v_ned(2);
    vd = v_ned(3);

    if h <= p_drone.altitude_take_off_zone
        drone.altitude_state = 2;
    else
        drone.altitude_state = 2;
    end

    if t == 0 || (drone.altitude_state ~= drone.altitude_state_prev)
        drone.initialize_integrator = 1;
    else
        drone.initialize_integrator = 0;
    end

    % Implement state machine
    switch drone.altitude_state
        case 1% in take-off zone
            d1 = 1;
            d2 = 1;
            d3 = 1;
            d4 = 1;
        case 2% outside take-off zone
            [vd_c, drone.P_pd_thrust] = pd_vd_hold(pd_c, h, ...
                drone.initialize_integrator, p_drone, drone.P_pd_thrust);
            [thrust, drone.P_vd_thrust] = vd_thrust_hold(vd_c, vd, ...
                drone.initialize_integrator, p_drone, p_sim, drone.P_vd_thrust);
            [ve_c, drone.P_pe_roll] = pe_ve_hold(pe_c, pe, ...
                drone.initialize_integrator, p_drone, drone.P_pe_thrust);
            [phi_c, drone.P_ve_roll] = ve_roll_hold(ve_c, ve, ...
                drone.initialize_integrator, p_drone, p_sim, drone.P_ve_roll);
            [torque_phi, drone.P_roll_torque] = roll_torque_hold(phi_c, ...
                phi, deg2rad(p), drone.initialize_integrator, p_drone, p_sim, ...
                drone.P_roll_torque);
            [vn_c, drone.P_pn_pitch] = pn_vn_hold(pn_c, pn, ...
                drone.initialize_integrator, p_drone, drone.P_pn_pitch);
            [theta_c, drone.P_vn_pitch] = vn_pitch_hold(vn_c, vn, ...
                drone.initialize_integrator, p_drone, p_sim, drone.P_vn_pitch);
            [torque_theta, drone.P_pitch_torque] = pitch_torque_hold(...
                theta_c, theta, deg2rad(q), ...
                drone.initialize_integrator, p_drone, drone.P_pitch_torque);
            [torque_psi, drone.P_psi_torque] = yaw_torque_hold(psi_c, ...
                psi, deg2rad(r), drone.initialize_integrator, p_drone, p_sim, ...
                drone.P_psi_torque);

            [d1, d2, d3, d4] = mixer(thrust, torque_phi, torque_theta, ...
                torque_psi, p_drone, p_sim);

            % artificially saturation delta_t
            d1 = saturate(d1, p_drone.dt_max, 0);
            d2 = saturate(d2, p_drone.dt_max, 0);
            d3 = saturate(d3, p_drone.dt_max, 0);
            d4 = saturate(d4, p_drone.dt_max, 0);
    end

    % Control outputs
    delta = [d1; d2; d3; d4];

    % Commanded (desired) states
    x_command = [...
                pn_c; ...% pn
            pe_c; ...% pe
            pd_c; ...% h
            vn_c; ...% vn
            ve_c; ...% ve
            vd_c; ...% vd
            0; % an
            0; % ae
            0; % ad
            0; ...% Va
            0; ...% alpha
            0; ...% beta
            phi_c; ...% phi
            theta_c; % theta
            psi_c; ...% psi
            0; % chi
            0; ...% p
            0; ...% q
            0; ...% r
            ];

    y = [delta; x_command];
    drone.altitude_state_prev = drone.altitude_state;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Autopilot functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [vd_c, P_pd_thrust] = pd_vd_hold(pd_c, pd, reinit, p_drone, p_sim, P_pd_thrust)

    if reinit || isempty(P_pd_thrust)
        P_pd_thrust.int = 0;
        P_pd_thrust.diff = 0;
        P_pd_thrust.error_prev = 0;
        P_pd_thrust.y_c_prev = 0;
    end

    [vd_c, P_pd_thrust] = pid_loop(pd_c, pd, Inf, p_drone.kp_vd, p_drone.ki_vd, p_drone.kd_vd, ...
        p_drone.vd_max, p_sim.dt, 0.1, P_pd_thrust);
end

function [ve_c, P_pe_roll] = pe_ve_hold(pe_c, pe, reinit, p_drone, P_pe_roll)

    if reinit || isempty(P_pe_roll)
        P_pe_roll.int = 0;
        P_pe_roll.diff = 0;
        P_pe_roll.error_prev = 0;
        P_pe_roll.y_c_prev = 0;
    end

    [ve_c, P_pe_roll] = pid_loop(pe_c, pe, Inf, p_drone.kp_ve, p_drone.ki_ve, p_drone.kd_ve, ...
        p_drone.v_max, p_sim.dt, 0.1, P_pe_roll);
end

function [vn_c, P_pn_pitch] = pn_vn_hold(pn_c, pn, reinit, p_drone, p_sim, P_pn_pitch)

    if reinit || isempty(P_pn_pitch)
        P_pn_pitch.int = 0;
        P_pn_pitch.diff = 0;
        P_pn_pitch.error_prev = 0;
        P_pn_pitch.y_c_prev = 0;
    end

    [vn_c, P_pn_pitch] = pid_loop(pn_c, pn, Inf, p_drone.kp_vn, p_drone.ki_vn, p_drone.kd_vn, ...
        p_drone.v_max, p_sim.dt, 0.1, P_pn_pitch);
end

function [thrust, P_h_thrust] = h_thrust_hold(h_c, h, reinit, p_drone, p_sim, P_h_thrust)

    if reinit || isempty(P_h_thrust)
        P_h_thrust.int = 0;
        P_h_thrust.diff = 0;
        P_h_thrust.error_prev = 0;
        P_h_thrust.y_c_prev = 0;
    end

    [thrust, P_h_thrust] = pid_loop(h_c, h, Inf, p_drone.kp_h, p_drone.ki_h, p_drone.kd_h, ...
        p_drone.thrust_max, p_sim.dt, 0.1, P_h_thrust);
end

function [torque_phi, P_roll_torque] = roll_torque_hold(phi_c, phi, p, ...
        reinit, p_drone, p_sim, P_roll_torque)

    if reinit || isempty(P_roll_torque)
        P_roll_torque.int = 0;
        P_roll_torque.diff = 0;
        P_roll_torque.error_prev = 0;
        P_roll_torque.y_c_prev = 0;
    end

    [torque_phi, P_roll_torque] = pid_loop(phi_c, phi, p, p_drone.kp_phi, 0, ...
        p_drone.kd_phi, p_drone.torque_max, p_sim.dt, 0.1, P_roll_torque);
end

function [torque_theta, P_pitch_torque] = pitch_torque_hold(theta_c, ...
        theta, q, reinit, p_drone, p_sim, P_pitch_torque)

    if reinit || isempty(P_pitch_torque)
        P_pitch_torque.int = 0;
        P_pitch_torque.diff = 0;
        P_pitch_torque.error_prev = 0;
        P_pitch_torque.y_c_prev = 0;
    end

    [torque_theta, P_pitch_torque] = pid_loop(theta_c, theta, q, ...
        p_drone.kp_theta, 0, p_drone.kd_theta, p_drone.torque_max, p_sim.dt, 0.1, ...
        P_pitch_torque);
end

function [torque_psi, P_psi_torque] = yaw_torque_hold(psi_c, psi, r, ...
        reinit, p_drone, p_sim, P_psi_torque)

    if reinit || isempty(P_psi_torque)
        P_psi_torque.int = 0;
        P_psi_torque.diff = 0;
        P_psi_torque.error_prev = 0;
        P_psi_torque.y_c_prev = 0;
    end

    [torque_psi, P_psi_torque] = pid_loop(psi_c, psi, r, p_drone.kp_psi, 0, ...
        p_drone.kd_psi, p_drone.torque_max, p_sim.dt, 0.1, P_psi_torque);
end

function [theta_c, P_vn_pitch] = vn_pitch_hold(vn_c, vn, reinit, p_drone, p_sim, ...
        P_vn_pitch)

    if reinit || isempty(P_vn_pitch)
        P_vn_pitch.int = 0;
        P_vn_pitch.diff = 0;
        P_vn_pitch.error_prev = 0;
        P_vn_pitch.y_c_prev = 0;
    end

    [theta_c, P_vn_pitch] = pid_loop(vn_c, vn, Inf, p_drone.kp_vn, p_drone.ki_vn, 0, ...
        p_drone.theta_max, p_sim.dt, 0.1, P_vn_pitch);
end

function [phi_c, P_ve_roll] = ve_roll_hold(ve_c, ve, reinit, p_drone, p_sim, P_ve_roll)

    if reinit || isempty(P_ve_roll)
        P_ve_roll.int = 0;
        P_ve_roll.diff = 0;
        P_ve_roll.error_prev = 0;
        P_ve_roll.y_c_prev = 0;
    end

    [phi_c, P_ve_roll] = pid_loop(ve_c, ve, Inf, p_drone.kp_ve, p_drone.ki_ve, 0, ...
        p_drone.phi_max, p_sim.dt, 0.1, P_ve_roll);
end

function [thrust, P_vd_thrust] = vd_thrust_hold(vd_c, vd, reinit, p_drone, p_sim, ...
        P_vd_thrust)

    if reinit || isempty(P_vd_thrust)
        % P_vd_thrust.int = -0.5;
        P_vd_thrust.int = 0;
        P_vd_thrust.diff = 0;
        P_vd_thrust.error_prev = 0;
        P_vd_thrust.y_c_prev = 0;
    end

    [thrust, P_vd_thrust] = pid_loop(vd_c, vd, Inf, p_drone.kp_vd, p_drone.ki_vd, 0, ...
        p_drone.thrust_max, p_sim.dt, 0.1, P_vd_thrust);
end

function [thrust, P_ad_thrust] = ad_thrust_hold(ad_c, ad, reinit, p_drone, p_sim, ...
        P_ad_thrust)

    if reinit || isempty(P_ad_thrust)
        P_ad_thrust.int = 0;
        P_ad_thrust.diff = 0;
        P_ad_thrust.error_prev = 0;
        P_ad_thrust.y_c_prev = 0;
    end

    [thrust, P_ad_thrust] = pid_loop(ad_c, ad, Inf, p_drone.kp_ad, p_drone.ki_ad, 0, ...
        p_drone.thrust_max, p_sim.dt, 0.1, P_ad_thrust);
end

function [theta_c, P_an_pitch] = an_pitch_hold(an_c, an, reinit, p_drone, ...
        P_an_pitch)

    if reinit || isempty(P_an_pitch)
        P_an_pitch.int = 0;
        P_an_pitch.diff = 0;
        P_an_pitch.error_prev = 0;
        P_an_pitch.y_c_prev = 0;
    end

    [theta_c, P_an_pitch] = pid_loop(an_c, an, Inf, p_drone.kp_an, p_drone.ki_an, 0, ...
        p_drone.theta_max, p_sim.dt, 0.1, P_an_pitch);
end

function [phi_c, P_ae_roll] = ae_roll_hold(ae_c, ae, reinit, p_drone, P_ae_roll)

    if reinit || isempty(P_ae_roll)
        P_ae_roll.int = 0;
        P_ae_roll.diff = 0;
        P_ae_roll.error_prev = 0;
        P_ae_roll.y_c_prev = 0;
    end

    [phi_c, P_ae_roll] = pid_loop(ae_c, ae, Inf, p_drone.kp_ae, p_drone.ki_ae, 0, ...
        p_drone.phi_max, p_sim.dt, 0.1, P_ae_roll);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mixer
%     - mixes the effects of thrust and torques
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [d1, d2, d3, d4] = mixer(thrust, torque_phi, torque_theta, ...
        torque_psi, p_drone, p_sim)

    omega1_sqd = (thrust / (4 * p_drone.C_prop) + ...
        torque_theta / (2 * p_drone.C_prop * p_drone.l_arm) - ...
        torque_psi / (4 * p_drone.CD));

    if omega1_sqd > 0
        d1 = sqrt(omega1_sqd) / p_drone.k_omega;
    else
        d1 = 0;
    end

    omega2_sqd = (thrust / (4 * p_drone.C_prop) - ...
        torque_phi / (2 * p_drone.C_prop * p_drone.l_arm) + ...
        torque_psi / (4 * p_drone.CD));

    if omega2_sqd > 0
        d2 = sqrt(omega2_sqd) / p_drone.k_omega;
    else
        d2 = 0;
    end

    omega3_sqd = (thrust / (4 * p_drone.C_prop) - ...
        torque_theta / (2 * p_drone.C_prop * p_drone.l_arm) - ...
        torque_psi / (4 * p_drone.CD));

    if omega3_sqd > 0
        d3 = sqrt(omega3_sqd) / p_drone.k_omega;
    else
        d3 = 0;
    end

    omega4_sqd = (thrust / (4 * p_drone.C_prop) + ...
        torque_phi / (2 * p_drone.C_prop * p_drone.l_arm) + ...
        torque_psi / (4 * p_drone.CD));

    if omega4_sqd > 0
        d4 = sqrt(omega4_sqd) / p_drone.k_omega;
    else
        d4 = 0;
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Saturate
%       - saturation function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = saturate(in, up_limit, low_limit)

    if in > up_limit
        out = up_limit;
    elseif in < low_limit
        out = low_limit;
    else
        out = in;
    end

end

author: lis-epfl

repo: swarmlab

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


text:
function [X_t,Mu_t,Sig_t]=OUstep(X_0,t,Mu,Th,Sig)

[NumSimul,N]=size(X_0);

% location
ExpM=expm(-Th*t);

Mu_t = repmat((Mu-ExpM*Mu)',NumSimul,1) +  X_0*ExpM';

% scatter
TsT=kron(Th,eye(N))+kron(eye(N),Th);

VecSig=reshape(Sig,N^2,1);
VecSig_t=inv(TsT)*(eye(N^2)-expm(-TsT*t))*VecSig;
Sig_t=reshape(VecSig_t,N,N);
Sig_t=(Sig_t+Sig_t')/2;

Eps=mvnrnd(zeros(N,1),Sig_t,NumSimul);

X_t=Mu_t+Eps;
Mu_t=mean(Mu_t,1)';

author: Sable

repo: mcbench-benchmarks


text:
function value = rhs_test02 ( x_num, x, t )

%*****************************************************************************80
%
%% RHS_TEST02 evaluates the right hand side for problem 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X_NUM, the number of nodes.
%
%    Input, real X(X_NUM,1), the node coordinates.
%
%    Input, real T, the current time.
%
%    Output, real VALUE(X_NUM,1), the source term.
%
  value(1:x_num,1) = 0.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ value_min, value_mean, value_max, value_var ] = ...
 tet_mesh_quality1 ( node_num, node_xyz, tetra_order, tetra_num, tetra_node )

%*****************************************************************************80
%
%% TET_MESH_QUALITY1 returns a tet mesh quality factor.
%
%  Discussion:
%
%    The tet mesh quality measure is the minimum of the
%    corresponding tetrahedron quality measure, over all tetrahedrons in the
%    tet mesh.
%
%    This routine is designed for an order-4 tet mesh.  Order 10 tet meshes
%    may be input, but the extra nodes are ignored.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, double NODE_XYZ(3,NODE_NUM), the nodes.
%
%    Input, integer TETRA_ORDER, the order of the tetrahedrons.
%
%    Input, integer TETRA_NUM, the number of tetrahedrons.
%
%    Input, integer TETRA_NODE(TETRA_ORDER,TETRA_NUM), the indices of the nodes
%    that make up the tetrahedrons.
%
%    Output, double VALUE_MIN, VALUE_MEAN, VALUE_MAX, VALUE_VAR,
%    the minimum, mean, maximum and variance of the quality measure.
%
  dim_num = 3;

  for tetra = 1 : tetra_num

    tetrahedron(1:dim_num,1:4) = node_xyz(1:dim_num,tetra_node(1:4,tetra));

    tetrahedron_quality(tetra) = tetrahedron_quality1_3d ( tetrahedron );

  end

  value_max  = max  ( tetrahedron_quality(1:tetra_num) );
  value_min  = min  ( tetrahedron_quality(1:tetra_num) );
  value_mean = mean ( tetrahedron_quality(1:tetra_num) );
  value_var  = var  ( tetrahedron_quality(1:tetra_num) );

  return
end


text:
function value = function_2d_name ( ifunc )

%*****************************************************************************80
%
%% FUNCTION_2D_NAME returns the name of the current 2D function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer IFUNC, the index of the function.
%
%    Output, character FUNCTION_2D_NAME(7), the name of the function.
%
  if ( ifunc == 1 )
    value = '      1';
  elseif ( ifunc == 2 )
    value = '      X';
  elseif ( ifunc == 3 )
    value = '    X^2';
  elseif ( ifunc == 4 )
    value = '    X^3';
  elseif ( ifunc == 5 )
    value = '    X^4';
  elseif ( ifunc == 6 )
    value = '    X^5';
  elseif ( ifunc == 7 )
    value = '    X^6';
  elseif ( ifunc == 8 )
    value = '      R';
  elseif ( ifunc == 9 )
    value = ' SIN(X)';
  elseif ( ifunc == 10 )
    value = ' EXP(X)';
  elseif ( ifunc == 11 )
    value = '1/(1+R)';
  elseif ( ifunc == 12 )
    value = 'SQRT(R)';
  else
    value = '???????';
  end

  return
end


text:
%
% [f, Covf, pseudoX, logtheta, Kpp, Kxp, Lp] =
%     gplearn(y, V, X, pseudoX, kX, logtheta, ...)
%
% V is inverse covariance (i.e. precision) matrix for observation noise -
% sparse format is recommended.
%
% In practice, inv(V) should be diagonal. However, inv(V) can be more
% complex sparse matrix as long as it has zeros where kX is non-zero. That
% is, observation noise can be correlated if the corresponding latent
% function values are independent a priori. This more general property is
% used in PCA model, where inv(V) and kX are block-diagonal (although rows
% and columns may be in mixed order!) but the blocks "overlap" only on the
% diagonal, that is, they both have non-zero elements at the same time only
% on the diagonal.
%
% y = inv(V)*c, that is, c = V*y
% N( inv(V)*c | 0, inv(V) + Kxp inv(Kpp) Kpx )

% $$$ function [f, Covf, pseudoX, logtheta, Kpp, Kxp, Lp] = gplearn(y, V, ...
% $$$                                                   X, pseudoX, kX, logtheta, ...
% $$$                                                   varargin)
function [f, Covf, Xp, logtheta, Kpp, Kxp, Lp] = gplearn(logtheta, covfunc, ...
                                                  X, y, V, Xp, varargin)

% function [] = gp_learn(covfunc, theta, y, V)

warning('Deprecated. Use gp_learn or gp_learn_pseudo');

if ~iscell(covfunc)
  covfunc = {covfunc};
end

opts = struct( ...
    'maxsearch', 10, ...
    'vy', [], ...
    'cholv', [], ...
    'checkgrad', false, ...
    'updatepseudo', true, ...
    'priorlogtheta', []);

[ opts, errmsg, wrnmsg ] = argschk( opts, varargin{:} );
if ~isempty(errmsg), error( errmsg ), end
if ~isempty(wrnmsg), warning( wrnmsg ), end

m = cols(Xp); % number of pseudoinputs
d = rows(X); % dimensionality of the input space

if isempty(opts.cholv)
  LV = chol(V, 'lower');
%  [LV,V] = safechol(V, 1e-16, 'lower');
else
  LV = opts.cholv;
end

if isempty(opts.vy)
  c = V*y;
else
  c = opts.vy;
end

if opts.maxsearch ~= 0
  if opts.updatepseudo && ~isempty(Xp)
    % Learn pseudo-inputs and hyperparameters
    x0 = [Xp(:); logtheta(:)];
    func = @(x) bound2cost(x, c, V, LV, X, covfunc, m, d, true);
    fiter = inf;
    while ~any(~isinf(fiter))
      [x, fiter] = minimize(x0, func, opts.maxsearch);
      x0 = x0 - 1; % adhoc trying to help ill conditioned matrices..
    end
    Xp = reshape(x(1:(m*d)), [d,m]);
    logtheta = x((m*d+1):end);
  else
    % Learn only hyperparameters
    x0 = logtheta(:);
    func = @(x) bound2cost(x, c, V, LV, X, covfunc, [], d, false, Xp);
    fiter = inf;
    while ~any(~isinf(fiter))
      [logtheta, fiter] = minimize(x0, func, opts.maxsearch);
      x0 = x0 - 1; % adhoc trying to help ill conditioned matrices..
    end
  end
  if opts.checkgrad
    mycheckgrad(func, x0, 1e-8)
  end
end

% THE FOLLOWING WOULD COME AS A SIDE PRODUCT FROM THE OPTIMIZATION BUT HOW
% TO GET THEM??

if isempty(Xp)
  % Full case (no pseudo-inputs)

  n = cols(X);
  In = diag(sparse(ones(n,1)));

  K = feval(covfunc{:}, logtheta, X, X);
% $$$   norm_K = normest(K);
% $$$   K = K + 1e-14*norm_K*In; % SET PROPERLY
%  K = K + 1e-6*In;

  % DEBUG: REGULARIZE
  K = regularize(K);
  Lp = chol(K, 'lower');
%  [Lp,K] = safechol(K, 1e-10, 'lower');
  %K = Lp*Lp';
  Kpp = K;
  Kxp = [];

  % We assume that V is DIAGONAL
  V = sparse(V);
  L = chol(inv(V) + K, 'lower');
%  L = safechol(inv(V) + K, 1e-16, 'lower');

  % Posterior mean and covariance
  if false
    % TODO: NOTE, THAT DIAGONAL COVARIANCE APPROXIMATION IS A BIT
    % ADHOC!!!!! CAUSES LOGLIKELIHOOD TO BE INCORRECT!!!
    f = K * c - K * solve_triu(L', solve_tril(L, K*c));
    Covf = spalloc(n,n,n);
    for i=1:n
      r = solve_tril(L,K(:,i));
      Covf(i,i) = K(i,i) - r'*r;
    end
  else
    invLK = full(solve_tril(L,K));
    Covf = K - invLK' * invLK;
    f = Covf * c;
  end

else
  % Pseudo input case

  Im = diag(sparse(ones(m,1)));

  % GP-prior covariances
  Kpp = feval(covfunc{:}, logtheta, Xp, Xp);
% $$$   norm_Kpp = normest(Kpp);
% $$$   Kpp = Kpp + 1e-14*norm_Kpp*Im; % SET PROPERLY
%  Kpp = Kpp + 1e-3*Im; % SET PROPERLY
  Kxp = feval(covfunc{:}, logtheta, X, Xp);

  % DEBUG: REGULARIZE
  %Kpp = regularize(Kpp);

  % Helpful variables
  reg = 1;
  while 1
    try
      Kpp = regularize(Kpp, reg);
      Lp = chol(Kpp, 'lower');
      break
    catch
      warning('Ill conditioned covariance matrix in learning phase!!');
      reg = reg * 10;
    end
  end
  %  [Lp, Kpp] = safechol(Kpp, 1e-10, 'lower');
  %Kpp = Lp*Lp';

  %szLp = size(Lp)
  %szKxp = size(Kxp)
  %szLV = size(LV)
  R = solve_tril(Lp, Kxp' * LV);
  S = Im + R*R';
  LS = chol(S, 'lower');
%  LS = safechol(S, 1e-16, 'lower');
  z = LV \ c;
  beta = solve_tril(LS, R*z);
  invLsLp = solve_tril(LS, Lp');

  % Posterior mean and covariance
  f = invLsLp' * beta;
  Covf = invLsLp' * invLsLp;

end

% $$$ if issparse(Lp)
% $$$   density_in_gplearn = nnz(Lp) / prod(size(Lp));
% $$$ end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [f, df] = bound2cost(x, c, V, LV, inputs, kX, m, d, updatepseudo, ...
                              pseudoinputs)

if updatepseudo
  pseudoinputs = reshape(x(1:(m*d)), [d,m]);
  logtheta = x((m*d+1):end);
  [bound, dlogtheta, dpseudoinputs] = pseudobound(c, V, LV, inputs, ...
                                                  pseudoinputs, logtheta, kX);
else
  logtheta = x(:);
  [bound, dlogtheta] = pseudobound(c, V, LV, inputs, pseudoinputs, logtheta, ...
                                   kX);
end

f = -bound;

if updatepseudo
  df = -[dpseudoinputs(:); dlogtheta(:)];
else
  df = -dlogtheta(:);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%
function [bound, dlogtheta, dpseudoinputs] = ...
    pseudobound(c, V, LV, inputs, pseudoinputs, logtheta, covfunc)

if ~isempty(pseudoinputs)
  %% PSEUDO INPUT CASE
  %disp('Pseudo case in gplearn > pseudobound');

  d = rows(pseudoinputs); % dimensionality of pseudo-inputs
  m = cols(pseudoinputs); % number of pseudo-inputs

  % GP-prior (co)variances and their gradients
  switch nargout
   case 1
    Kpp = feval(covfunc{:}, logtheta, pseudoinputs, pseudoinputs);
    Kxp = feval(covfunc{:}, logtheta, inputs, pseudoinputs);
    Kx = feval(covfunc{:}, logtheta, inputs, []);
   case 2
    [Kpp, dKpp_logtheta] = feval(covfunc{:}, logtheta, pseudoinputs, ...
                                   pseudoinputs);
    [Kxp, dKxp_logtheta] = feval(covfunc{:}, logtheta, inputs, pseudoinputs);
    [Kx, dKx_logtheta] = feval(covfunc{:}, logtheta, inputs, []);
   case 3
    [Kpp, dKpp_logtheta, dKpp_pseudo] = feval(covfunc{:}, logtheta, ...
                                                pseudoinputs, pseudoinputs);
    [Kxp, dKxp_logtheta, dKxp_pseudo] = feval(covfunc{:}, logtheta, ...
                                                inputs, pseudoinputs);
    [Kx, dKx_logtheta] = feval(covfunc{:}, logtheta, inputs, []);
   otherwise
    error('No defined action here')
  end

% $$$ if ~iscell(covfunc), covfunc = {covfunc}; end
% $$$ [K,tmp] = feval(covfunc{:}, logtheta, pseudoinputs, pseudoinputs);
% $$$ figure
% $$$ imagesc(K-Kpp)
% $$$ return

  Im = diag(sparse(ones(m,1)));
% $$$   norm_Kpp = normest(Kpp)
% $$$   Kpp = Kpp + 1e-14*norm_Kpp*Im; % SET PROPERLY
%  Kpp = Kpp + 1e-3*Im; % SET PROPERLY
% $$$   figure
% $$$   imagesc(Kpp);
% $$$   error('jou')

  % DEBUG: REGULARIZE
  %Kpp = regularize(Kpp);

  % Helpful variables
  Kpp = regularize(Kpp);
  try
    Lp = chol(Kpp, 'lower');
  catch
    warning('Covariance matrix ill conditioned');
    bound = -inf;
    dlogtheta = nan * zeros(length(logtheta),1);
    dpseudoinputs = nan * zeros(d,m);
    return
  end
%  [Lp,Kpp] = safechol(Kpp, 1e-10, 'lower');
  %Kpp = Lp*Lp';

  R = solve_tril(Lp, Kxp' * LV);
  S = Im + R*R';
  LS = chol(S, 'lower');
%  LS = safechol(S, 1e-16, 'lower');
  LA = Lp * LS;
  z = LV \ c;
  H = solve_tril(LS, R);
  T = solve_triu(LA', H);
  beta = H * z;
  b = T * z;
  invLpR = solve_triu(Lp', R);
  diagV = full(diag(V));
  dif = c - V*(Kxp*b);

  % Lower bound for loglikelihood
  bound = -logdettri(LS) - 0.5*z'*z + 0.5*beta'*beta ...
          - 0.5*Kx'*diagV ...
          + 0.5*traceprod(R,R,true);

  % General gradient coefficients
  Tpp = 0.5 * solve_triu(LA', solve_triu(LA', LS'*LS-eye(m))') ...
        - 0.5 * invLpR*invLpR';
  Txp = -LV * T' ...
        + LV*invLpR';

  % Gradient for hyperparameters
  if nargout >= 2
    n = length(logtheta);
    dlogtheta = zeros(n,1);
    for i=1:n
      dlogtheta(i) = traceprod(dKpp_logtheta(:,:,i), Tpp, true) ...
          - 0.5 * b' * dKpp_logtheta(:,:,i) * b ...
          + traceprod(dKxp_logtheta(:,:,i), Txp, true) ...
          + dif' * dKxp_logtheta(:,:,i) * b ...
          - 0.5 * dKx_logtheta(:,:,i)' * diagV;
    end
  end

  if nargout >= 3
    % Gradient for pseudo-inputs
    dpseudoinputs = zeros(d,m);
    Dpp = spalloc(m,m,2*m);
    for i=1:m
      for j=1:d
        Dpp(:) = 0;
        Dpp(i,:) = dKpp_pseudo(j,:,i);
        Dpp(:,i) = Dpp(:,i) + dKpp_pseudo(j,:,i)';

        dpseudoinputs(j,i) = traceprod(Dpp, Tpp, true) ...
            - 0.5 * b' * Dpp * b ...
            + dKxp_pseudo(j,:,i) * Txp(:,i) ...
            + dif' * dKxp_pseudo(j,:,i)' * b(i);

      end
    end
  end

else

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % FULL CASE (NO PSEUDO INPUTS)
% $$$   disp('Full case in gplearn > pseudobound');
  d = rows(inputs); % dimensionality of pseudo-inputs
  n = cols(inputs);
  %m = cols(pseudoinputs); % number of pseudo-inputs

  % GP-prior (co)variances and their gradients
  switch nargout
   case 1
    K = feval(covfunc{:}, logtheta, inputs, inputs);
   case 2
    [K, dK_logtheta] = feval(covfunc{:}, logtheta, inputs, inputs);
   otherwise
    error('No pseudo inputs given');
  end

  In = diag(sparse(ones(n,1)));
% $$$   norm_K = normest(K);
% $$$   K = K + 1e-14*norm_K*In; % SET PROPERLY
%  K = K + 1e-6*In;

  % DEBUG: REGULARIZE
  %K = regularize(K);

  % ASSUME THAT V is DIAGONAL
  V = sparse(V);
  Ky = inv(V) + K;
  try
    L = chol(Ky, 'lower');
  catch
    warning('Covariance matrix ill conditioned');
    bound = -inf;
    dlogtheta = nan * zeros(length(logtheta),1);
    return
  end
%  [L,Ky] = safechol(Ky, 1e-10, 'lower');
  % DEBUG:
  %Ky = L*L';
  y = sparse(V) \ c;

  z = solve_tril(L, y);

  % Lower bound for loglikelihood
  bound = -logdettri(L) - 0.5*z'*z;

% $$$   % General gradient coefficients
% $$$   Tpp = 0.5 * solve_triu(LA', solve_triu(LA', LS'*LS-eye(m))') ...
% $$$         - 0.5 * invLpR*invLpR';
% $$$   Txp = -LV * T' ...
% $$$         + LV*invLpR';

  % Gradient for hyperparameters
  if nargout >= 2
    b = solve_triu(L',z);
    n = length(logtheta);
    dlogtheta = zeros(n,1);
    if issparse(Ky)
      invKy = sinv(Ky);
    else
      invKy = inv(Ky);
    end
    for i=1:n
      dlogtheta(i) = 0.5 * b' * dK_logtheta(:,:,i) * b ...
          - 0.5 * traceprod( invKy, dK_logtheta(:,:,i), true );
% $$$       dlogtheta(i) = traceprod(dKpp_logtheta(:,:,i), Tpp, true) ...
% $$$           - 0.5 * b' * dKpp_logtheta(:,:,i) * b ...
% $$$           + traceprod(dKxp_logtheta(:,:,i), Txp, true) ...
% $$$           + dif' * dKxp_logtheta(:,:,i) * b ...
% $$$           - 0.5 * dKx_logtheta(:,:,i)' * diagV;
    end
  end

end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
%--- help for recreate_parameters ---
%
%  INTERNAL FUNCTION
%
%

author: jmaih


text:
classdef TestDrawMarker
    %TestDrawMarker

    methods (Static)
        function test_1
            img = 255*ones(64,64,3,'uint8');
            out = cv.drawMarker(img, [32,32]);
            validateattributes(out, {class(img)}, {'size',size(img)});
        end

        function test_options
            img = zeros([100 400 3],'uint8');
            markers = '+x*ds^v';
            for i=1:numel(markers)
                img = cv.drawMarker(img, [50*i 50], ...
                    'Color',[255 0 0], 'MarkerType',markers(i), ...
                    'MarkerSize',30, 'Thickness',2, 'LineType','AA');
            end
            validateattributes(img, {'uint8'}, {'size',[100 400 3]});
        end

        function test_error_argnum
            try
                cv.drawMarker();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function success = testPrepareJAABAData(datatype,rootdir)

success = false;
if nargin < 2,
  rootdir = '/groups/branson/bransonlab/projects/JAABA/sampledata';
end

%% create the JAABA GUI
hfig = PrepareJAABAData;

% get handles
handles = guidata(hfig);

%% parameters

params = struct;
params.datatype = datatype;

switch datatype,

  case 'Ctrax'

    params.inputdir = fullfile(rootdir,'flies','GMR_14C07_AE_01_TrpA_Rig1Plate15BowlA_20120404T141155');
    params.inputmoviefilestr = 'movie.ufmf';
    params.inputmoviefile = fullfile(params.inputdir,params.inputmoviefilestr);
    params.inputfiledict = {'intrxfile','ctrax_results.mat',fullfile(params.inputdir,'ctrax_results.mat')
      'annfile','movie.ufmf.ann',fullfile(params.inputdir,'movie.ufmf.ann')};
    [~,params.experiment_name] = fileparts(params.inputdir);
    params.experiment_name = [params.experiment_name,'__Ctrax'];
    params.outputdir = '../tests/PrepareJAABADataTests';
    params.outputmoviefilestr = params.inputmoviefilestr;
    params.outputtrxfilestr = 'registered_trx.mat';
    params.outputperframedir = 'perframe';

    % sample jab file
    params.jabfile = fullfile(rootdir,'flies','ChaseAR_new_v10p0.jab');

    % options
    params.options = struct;
    params.options.SoftLinkFiles = true;
    params.options.fliplr = false;
    params.options.flipud = false;
    params.options.dotransposeimage = false;
    params.options.fps = 30.344385106161738;
    params.options.pxpermm = 7.888175824175824;
    params.options.OverRideFPS = true;
    params.options.OverRideArena = false;
    params.options.CropFirstFrame = 1;
    params.options.CropEndFrame = inf;
    params.options.ArenaType = 'circle';

    % read arena or fps?
    params.readarena = true;
    params.readfps = true;

  case 'CtraxPlusWings',

    params.inputdir = fullfile(rootdir,'flies','GMR_14C07_AE_01_TrpA_Rig1Plate15BowlA_20120404T141155');
    params.inputmoviefilestr = 'movie.ufmf';
    params.inputmoviefile = fullfile(params.inputdir,params.inputmoviefilestr);
    params.inputfiledict = {'intrxfile','wingtracking_results.mat',fullfile(params.inputdir,'wingtracking_results.mat')
      'annfile','movie.ufmf.ann',fullfile(params.inputdir,'movie.ufmf.ann')
      'inperframedir','perframe',fullfile(params.inputdir,'perframe')};
    [~,params.experiment_name] = fileparts(params.inputdir);
    params.experiment_name = [params.experiment_name,'__CtraxPlusWings'];
    params.outputdir = '../tests/PrepareJAABADataTests';
    params.outputmoviefilestr = params.inputmoviefilestr;
    params.outputtrxfilestr = 'wingtracking_results.mat';
    params.outputperframedir = 'perframe';

    % sample jab file
    params.jabfile = fullfile(rootdir,'flies','WingExtension_AR_v1p0.jab');

    % options
    params.options = struct;
    params.options.SoftLinkFiles = true;
    params.options.fliplr = false;
    params.options.flipud = false;
    params.options.dotransposeimage = false;
    params.options.fps = 30.344385106161738;
    params.options.pxpermm = 7.888175824175824;
    params.options.OverRideFPS = true;
    params.options.OverRideArena = false;
    params.options.CropFirstFrame = 1;
    params.options.CropEndFrame = inf;
    params.options.ArenaType = 'circle';

    % read arena or fps?
    params.readarena = true;
    params.readfps = true;

  case 'SimpleTwoFlies',

    params.inputdir = fullfile(rootdir,'courtship_bowls','FCF_pBDPGAL4U_1500437_TrpA_Rig2Plate17BowlD_20121121T152832');
    params.inputmoviefilestr = 'movie.ufmf';
    params.inputmoviefile = fullfile(params.inputdir,params.inputmoviefilestr);
    params.inputfiledict = {'intrxfile','registered_trx.mat',fullfile(params.inputdir,'registered_trx.mat')
      'inperframedir','perframe',fullfile(params.inputdir,'perframe')};
    [~,params.experiment_name] = fileparts(params.inputdir);
    params.experiment_name = [params.experiment_name,'__SimpleTwoFlies'];
    params.outputdir = '../tests/PrepareJAABADataTests';
    params.outputmoviefilestr = params.inputmoviefilestr;
    params.outputtrxfilestr = 'registered_trx.mat';
    params.outputperframedir = 'perframe';

    % sample jab file
    params.jabfile = fullfile(rootdir,'courtship_bowls','copulation_cb_ARv1p8.jab');

    % options
    params.options = struct;
    params.options.SoftLinkFiles = true;
    params.options.fliplr = false;
    params.options.flipud = false;
    params.options.dotransposeimage = false;
    params.options.OverRideFPS = true;
    params.options.OverRideArena = false;
    params.options.CropFirstFrame = 1;
    params.options.CropEndFrame = inf;
    params.options.ArenaType = 'None';

    % read arena or fps?
    params.readarena = false;
    params.readfps = true;

  case 'LarvaeRiveraAlba'

    params.inputdir = fullfile(rootdir,'larvae','OregonRA');
    params.inputmoviefilestr = 'movie.mov';
    params.inputmoviefile = fullfile(params.inputdir,params.inputmoviefilestr);
    params.inputfiledict = {'intrxfile','trx.mat',fullfile(params.inputdir,'trx.mat')
       'inperframedir','perframe',fullfile(params.inputdir,'perframe')};
    [~,params.experiment_name] = fileparts(params.inputdir);
    params.experiment_name = [params.experiment_name,'__LarvaeRiveraAlba'];
    params.outputdir = '../tests/PrepareJAABADataTests';
    params.outputmoviefilestr = params.inputmoviefilestr;
    params.outputtrxfilestr = 'trx.mat';
    params.outputperframedir = 'perframe';

    % sample jab file
    params.jabfile = fullfile(rootdir,'larvae','TODO.jab');

    % options
    params.options = struct;
    params.options.SoftLinkFiles = true;
    params.options.fliplr = false;
    params.options.flipud = false;
    params.options.dotransposeimage = false;
    params.options.OverRideFPS = false;
    params.options.OverRideArena = false;
    params.options.CropFirstFrame = 1;
    params.options.CropEndFrame = inf;
    params.options.ArenaType = 'circle';

    % read arena or fps?
    params.readarena = false;
    params.readfps = false;

  case 'Qtrax',

    params.inputdir = fullfile(rootdir,'qtrax','EH091202_15A01_p1');
    params.inputmoviefilestr = 'EH091202_15A01_p1.mp4';
    params.inputmoviefile = fullfile(params.inputdir,params.inputmoviefilestr);
    params.inputfiledict = {'featfile','EH091202_15A01_p1_1_feat.mat',fullfile(params.inputdir,'EH091202_15A01_p1_1_feat.mat')
      'roifile','EH091202_15A01_p1_1_roi.mat',fullfile(params.inputdir,'EH091202_15A01_p1_1_roi.mat')};
    [~,params.experiment_name] = fileparts(params.inputdir);
    params.experiment_name = [params.experiment_name,'__Qtrax'];
    params.outputdir = '../tests/PrepareJAABADataTests';
    [~,~,ext] = fileparts(params.inputmoviefilestr);
    params.outputmoviefilestr = ['movie',ext];
    params.outputtrxfilestr = 'trx.mat';
    params.outputperframedir = 'perframe';

    % sample jab file
    params.jabfile = fullfile(rootdir,'qtrax','qtrax_dummy20130819.jab');

    % options
    params.options = struct;
    params.options.SoftLinkFiles = true;
    params.options.fliplr = false;
    params.options.flipud = false;
    params.options.dotransposeimage = false;
    params.options.OverRideFPS = false;
    params.options.OverRideArena = false;
    params.options.CropFirstFrame = 1;
    params.options.CropEndFrame = inf;
    params.options.ArenaType = 'Circle';

    % read arena or fps?
    params.readarena = false;
    params.readfps = false;

  case 'MAGATAnalyzer',

    params.inputdir = fullfile(rootdir,'larvae_samuel','20120206_100300');
    params.inputmoviefilestr = '20120206_100300@FCF_attp2_1500062@UAS_TNT_2_0003@t7@n#n#n#n@30@.mmf';
    params.inputmoviefile = fullfile(params.inputdir,params.inputmoviefilestr);
    params.inputfiledict = {'expfile','expfile.mat',fullfile(params.inputdir,'FCF_attp2_1500062@UAS_TNT_2_0003_experiment_20120206_100300.mat')};
    % experiment name is encoded in movie name
    [~,params.experiment_name,ext] = fileparts(params.inputmoviefile);
    params.experiment_name = [params.experiment_name,'__MAGATAnalyzer'];
    params.outputdir = '../tests/PrepareJAABADataTests';
    params.outputmoviefilestr = ['movie',ext];
    params.outputtrxfilestr = 'trx.mat';
    params.outputperframedir = 'perframe';

    % sample jab file
    params.jabfile = fullfile(rootdir,'larvae_samuel','MAGATAnalyzer_Test.jab');

    % options
    params.options = struct;
    params.options.SoftLinkFiles = true;
    params.options.fliplr = false;
    params.options.flipud = false;
    params.options.dotransposeimage = false;
    params.options.OverRideFPS = false;
    params.options.OverRideArena = false;
    params.options.CropFirstFrame = 1;
    params.options.CropEndFrame = inf;
    params.options.ArenaType = 'None';

    % read arena or fps?
    params.readarena = false;
    params.readfps = false;

  case 'MWT',

    params.inputdir = fullfile(rootdir,'larvae_mwt','MZZ_ppk1d9GAL4_t2_TrpA_20120612_110553');
    params.inputmoviefilestr = '';
    params.inputmoviefile = '';
    tmp = dir(fullfile(params.inputdir,'*.blob*'));
    blobspaths = cellfun(@(x) fullfile(params.inputdir,x),{tmp.name},'UniformOutput',false);
    tmp = dir(fullfile(params.inputdir,'*.dat'));
    datpaths = cellfun(@(x) fullfile(params.inputdir,x),{tmp.name},'UniformOutput',false);
    params.inputfiledict = {'blobsfile','trx.blobs',blobspaths
      'datfiles','chor.dat',datpaths};
    [~,params.experiment_name] = fileparts(params.inputdir);
    params.experiment_name = [params.experiment_name,'__MWT'];
    params.outputdir = '../tests/PrepareJAABADataTests';
    params.outputmoviefilestr = params.inputmoviefilestr;
    params.outputtrxfilestr = 'trx.mat';
    params.outputperframedir = 'perframe';

    % sample jab file
    params.jabfile = fullfile(rootdir,'larvae_mwt','LarvaeMWT_Roll_v18TO.jab');

    % options
    params.options = struct;
    params.options.SoftLinkFiles = true;
    params.options.fliplr = false;
    params.options.flipud = false;
    params.options.dotransposeimage = false;
    params.options.fps = 30.319466231015564;
    params.options.pxpermm = 17.857142387323400;
    params.options.OverRideFPS = false;
    params.options.OverRideArena = false;
    params.options.CropFirstFrame = 1;
    params.options.CropEndFrame = inf;
    params.options.ArenaType = 'None';

    % read arena or fps?
    params.readarena = true;
    params.readfps = false;

  case 'LarvaeReid',

    name1 = 'FCF_attP2_1500062_50';
    name2 = 'o_ethylbutyrate015m8ul_0s1x360s0s#n#n#n@20';
    name3 = '20120427_134403';
    inputdir1 = fullfile(rootdir,'larvae_reid',name1,name2);
    params.inputdir = fullfile(inputdir1,name3);
    params.inputmoviefilestr = '';
    params.inputmoviefile = '';
    tmp = dir(fullfile(params.inputdir,'*.blob*'));
    blobspaths = cellfun(@(x) fullfile(params.inputdir,x),{tmp.name},'UniformOutput',false);
    params.inputfiledict = {'blobsfile','trx.blobs',blobspaths
      'kinmatfile','kinVariables.mat',fullfile(inputdir1,'kinVariables.mat')
      'eventmatfile','eventVariables.mat',fullfile(inputdir1,'eventVariables.mat')};
    params.experiment_name = sprintf('%s__%s__%s__LarvaeReid',name1,name2,name3);
    params.outputdir = '../tests/PrepareJAABADataTests';
    params.outputmoviefilestr = params.inputmoviefilestr;
    params.outputtrxfilestr = 'trx.mat';
    params.outputperframedir = 'perframe';

    % sample jab file
    params.jabfile = fullfile(rootdir,'larvae_reid','LarvaeReid_dummy.jab');

    % options
    params.options = struct;
    params.options.SoftLinkFiles = true;
    params.options.fliplr = false;
    params.options.flipud = false;
    params.options.dotransposeimage = false;
    params.options.fps = 30.319466231015564;
    params.options.pxpermm = 17.857142387323400;
    params.options.OverRideFPS = false;
    params.options.OverRideArena = false;
    params.options.CropFirstFrame = 1;
    params.options.CropEndFrame = inf;
    params.options.ArenaType = 'None';

    % read arena or fps?
    params.readarena = true;
    params.readfps = false;

  case 'LarvaeLouis',

    name1 = 'rawdata';
    name2 = '20111212-120514';
    inputdir1 = fullfile(rootdir,'larvae_louis',name1);
    params.inputdir = fullfile(inputdir1,name2);
    params.inputmoviefilestr = '20111212-120514_video.avi';
    params.inputmoviefile = fullfile(params.inputdir,params.inputmoviefilestr);
    params.inputfiledict = {'indatafile','20111212-120514_data.txt',fullfile(params.inputdir,'20111212-120514_data.txt')
      'inkinmatfile','kinVariables.mat',fullfile(inputdir1,'kinVariables.mat')};
    params.experiment_name = sprintf('%s__LarvaeLouis',name2);
    params.outputdir = '../tests/PrepareJAABADataTests';
    [~,~,ext] = fileparts(params.inputmoviefilestr);
    params.outputmoviefilestr = ['movie',ext];
    params.outputtrxfilestr = 'trx.mat';
    params.outputperframedir = 'perframe';

    % sample jab file
    params.jabfile = fullfile(rootdir,'larvae_louis','LarvaeLouis_dummy.jab');

    % options
    params.options = struct;
    params.options.SoftLinkFiles = true;
    params.options.fliplr = false;
    params.options.flipud = false;
    params.options.dotransposeimage = false;
    params.options.OverRideFPS = false;
    params.options.OverRideArena = false;
    params.options.CropFirstFrame = 1;
    params.options.CropEndFrame = inf;
    params.options.ArenaType = 'None';

    % read arena or fps?
    params.readarena = false;
    params.readfps = false;

end

%% set PrepareJAABAData fields to match

%% set data type
handles = guidata(hfig);
handles.InputDataTypeIndex = find(strcmp(params.datatype,handles.InputDataTypeNames));
handles.InputDataType = params.datatype;
handles = PrepareJAABAData('UpdateGUI',handles);
guidata(hfig,handles);

%% set input table entries
handles = guidata(hfig);

% movie
handles.inputmoviefilestr = params.inputmoviefilestr;
handles.InputVideoFile = params.inputmoviefile;

% rest of files
InputDataType = handles.InputDataTypes.(handles.InputDataType);
for i = 1:size(params.inputfiledict,1),
  j = find(strcmp({InputDataType.files.code},params.inputfiledict{i,1}));
  if numel(j) ~= 1,
    error('Could not match input file type %s',params.inputfiledict{i,1});
  end
  handles.inputfilestrs.(handles.InputDataType){j} = params.inputfiledict{i,2};
  handles.InputFiles.(handles.InputDataType){j} = params.inputfiledict{i,3};
end

handles = PrepareJAABAData('UpdateGUI',handles);
guidata(hfig,handles);

%% set output table entries

handles = guidata(hfig);

% experiment directory
handles.ExperimentDirectory = fullfile(params.outputdir,params.experiment_name);

% output file strings
handles.moviefilestr = params.outputmoviefilestr;
handles.trxfilestr = params.outputtrxfilestr;
handles.perframedirstr = params.outputperframedir;

handles = PrepareJAABAData('UpdateGUI',handles);
guidata(hfig,handles);

%% options

handles = guidata(hfig);

fns = fieldnames(params.options);
for i = 1:numel(fns),
  handles.(fns{i}) = params.options.(fns{i});
end

handles = PrepareJAABAData('UpdateGUI',handles);
guidata(hfig,handles);

%% read arena parameters

if params.readarena,

  PrepareJAABAData('pushbutton_ReadArenaParameters_Callback',hfig,[],guidata(hfig));
  handles = guidata(hfig);

end

%% read fps parameters

if params.readfps,

  PrepareJAABAData('pushbutton_ReadFPS_Callback',hfig,[],guidata(hfig));
  handles = guidata(hfig);

end

%% convert

outexpdir = fullfile(params.outputdir,params.experiment_name);
if exist(outexpdir,'dir'),
  rmdir(outexpdir,'s');
end

PrepareJAABAData('pushbutton_Convert_Callback',hfig,[],guidata(hfig));
handles = guidata(hfig);

%% check that this worked right

filesneeded = {
  fullfile(outexpdir,params.outputtrxfilestr)
  fullfile(outexpdir,params.outputperframedir)
  };
if ~isempty(params.outputmoviefilestr),
  filesneeded{end+1} = fullfile(outexpdir,params.outputmoviefilestr);
end

ismissing = false(size(filesneeded));
for i = 1:numel(filesneeded),
  ismissing(i) = ~exist(filesneeded{i},'file');
  if ismissing(i) && ispc,
    [~,didfind] = GetPCShortcutFileActualPath(filesneeded{i});
    ismissing(i) = didfind == 0;
  end
  fprintf('File %s exists? %d\n',filesneeded{i},~ismissing(i));
end
if any(ismissing),
  error(['Files missing:',sprintf(' %s',filesneeded{ismissing})]);
end

%% close the gui

PrepareJAABAData('figure1_CloseRequestFcn',hfig,[],guidata(hfig));

%% see if we can run JAABADetect on this experiment

JAABADetect(outexpdir,'jabfiles',{params.jabfile},'forcecompute',true);

success = true;

author: kristinbranson

repo: JAABA


text:
function [c] = kron(a,b)
%Kronecker product of two QTT_Tuckers
%   [C]=KRON(A,B) computes Kronecker product of A and B, where A and B are
%   in the QTT-Tucker format
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

if ( isempty(a) )
  c=b;
  return
elseif ( isempty(b) )
  c=a;
  return
end

c = qtt_tucker;
c.dphys = a.dphys+b.dphys;
c.tuck = cell(c.dphys, 1);
c.core = kron(a.core, b.core);
c.tuck(1:a.dphys) = a.tuck;
c.tuck(a.dphys+1:c.dphys) = b.tuck;

end

author: oseledets

repo: TT-Toolbox