text:
function out = tn(FUN,x0,varargin)
%TN   Truncated Newton minimization.
%
%   OUT = TN(FUN,X0) minimizes FUN starting at the point X0 using a
%   Hessian-free truncated Newton method. The outer loop is Newton's
%   method, and the inner loop uses a conjugate gradient method to compute
%   an approximation to the Newton direction. Furthermore, the Hessian
%   vector product is approximated using forward finite differences using
%   HESSVEC_FD. FUN is a handle for a function that takes a single vector
%   input and returns two arguments --- the scalar function value and the
%   vector-valued gradient. See POBLANO_OUT for details of the output
%   parameters.
%
%   OUT = TN(FUN,X0,'param',value,...) specifies a
%   parameters and its value. See POBLANO_PARAMS for further details on
%   standard parameters. Additionally, TN requires
%
%   'CGSolver' - Matlab CG method to use {'symmlq'}
%     'symmlq' : symmlq (designed for symmetric indefinite systems)
%     'pcg' :    pcg (designed for symmetric positive definite systems)
%
%   'CGIters' - maximum number of conjugate gradient iterations allowed {5}
%
%   'CGTolType' - CG stopping tolerance type used ('quadratic')
%     'quadratic' :    || R || / || G || <  min(0.5,|| G ||)
%     'superlinear' :  || R || / || G || <  min(0.5,sqrt(|| G ||))
%     'fixed' :        || R || < CGTol
%   where R is the residual and G is the gradient of FUN at X.
%
%   'CGTol' - CG stopping tolerance when CGTolType is 'fixed' {1e-6}
%
%   'HessVecFDStep' - Hessian vector product finite difference step {1e-10}
%     0 - use the default step given in HESSVEC_FD
%     >0 - fixed value to use at the difference step
%
%   PARAMS = TN('defaults') returns a structure containing the
%   default parameters for the particular Poblano method.
%
%
%   Examples
%
%   Suppose the function and gradient of the objective function are
%   specified in an mfile named mysin.m:
%
%     function [f,g]=example1(x,a)
%     if nargin < 2, a = 1; end
%     f = sin(a*x);
%     g = a*cos(a*x);
%
%   We can call the optimization method (using its default
%   parameters) using the command:
%
%     out = tn(@(x) example1(x,3), pi/4);
%
%   To change a parameter, we can specify a param/value input pair
%   as follows:
%
%     out = tn(@(x) example1(x,3), pi/4, 'Display', 'final');
%
%   Alternatively, we can use a structure to define the parameters:
%
%     params.MaxIters = 2;
%     out = tn(@(x) example1(x,3), pi/4, params);
%
%   See also POBLANO_OUT, POBLANO_PARAMS, HESSVEC_FD, FUNCTION_HANDLE.
%
%MATLAB Poblano Toolbox.
%Copyright 2009-2012, Sandia Corporation.

%% Parse parameters

% Create parser
params = inputParser;

% Set Poblano parameters
params = poblano_params(params);

% Set parameters for this method
params.addParamValue('CGIters',5,@(x) x > 0);
params.addParamValue('CGTolType','quadratic',@(x) ismember(x,{'quadratic','superlinear','fixed'}));
params.addParamValue('CGTol',1e-6, @(x) x > 0);
params.addParamValue('HessVecFDStep',1e-10, @(x) x >= 0);
params.addParamValue('CGSolver','symmlq',@(x) ismember(x,{'symmlq','pcg'}));

% Parse input
params.parse(varargin{:});

%% Check input arguments
if (nargin == 1) && isequal(FUN,'defaults') && (nargout == 1)
    out = params.Results;
    return;
elseif (nargin < 2)
    error('Error: invalid input arguments');
end

%% Initialize
xk = x0;
[fk,gk] = feval(FUN,xk);
out = poblano_out(xk,fk,gk,1,params);
ak = 1.0;

cgIters = params.Results.CGIters;
tolType = params.Results.CGTolType;
cgTol = params.Results.CGTol;
hessvecFDstep = params.Results.HessVecFDStep;
solver = params.Results.CGSolver;

%% Main loop
while out.ExitFlag == -1

    % Compute step direction, pk
    ngk = norm(gk);
    switch tolType
        case 'quadratic'
            cgTol = min(1-eps,min(0.5,ngk)*ngk);       % quadratic convergence
        case 'superlinear'
            cgTol = min(1-eps,min(0.5,sqrt(ngk))*ngk); % superlinear convergence
        otherwise
    end

    % Setup Hessian vector product finite difference approximation
    if hessvecFDstep > 0
        hv = @(v) hessvec_fd(v,FUN,xk,gk,hessvecFDstep);
    else
        hv = @(v) hessvec_fd(v,FUN,xk,gk);
    end

    % Keep track of number of function calls are made using hessvec_fd
    global nfev_hessvec_fd;
    nfev_hessvec_fd = 0;

    % Compute step direction, pk
    [pk,cg_flag,cg_relres,ncgfev]= feval(solver, hv, -gk, cgTol, cgIters);

    % Compute step length
    [xk,fk,gk,ak,lsinfo,nfev] = poblano_linesearch(FUN,xk,fk,gk,ak,pk,params.Results);
    if (lsinfo ~= 1)
        if strcmp(params.Results.Display, 'iter')
            fprintf(1,[mfilename,': line search warning = %d\n'],lsinfo);
        end
        pk = -gk;
        [xk,fk,gk,ak,lsinfo,nfev] = poblano_linesearch(FUN,xk,fk,gk,ak,pk,params.Results);
    end

    % Update counts, check exit conditions, etc.
    out = poblano_out(xk,fk,gk,nfev+nfev_hessvec_fd,params,out);
end

clear nfev_hessvec_fd;

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
 function [xs, info] = qpwls_pcg1(x, A, W, yi, C, varargin)
%function [xs, info] = qpwls_pcg1(x, A, W, yi, C, [options])
%|
%| quadratic penalized weighted least squares (QPWLS) via
%| preconditioned conjugate gradients (PCG) algorithm
%| cost(x) = (y-Ax)'W(y-Ax) / 2 + x'C'Cx / 2
%|
%| in
%|	x	[np 1]		initial estimate
%|	A	[nd np]		system matrix
%|	W	[nd nd]		data weighting matrix, usually Gdiag(wi)
%|	yi	[nd 1]		noisy data
%|	C	[nc np]		penalty 'differencing matrix' (0 for unregularized)
%|
%| options
%|	niter			# total iterations (default: 1)
%|					(max # if tol used)
%|	isave	[]		list of iterations to archive (default: 'last')
%|	userfun	@		user defined function handle (see default below)
%|					taking arguments (x, userarg{:})
%|	userarg {}		user arguments to userfun (default {})
%|	precon	[np np]		preconditioner (matrix or object) (or 1)
%|	dircheck 0|1		check descent direction? (default: 1)
%|				set to 0 to save time, if you dare...
%|	stop_diff_tol		stop iterations if norm(xnew-xold)/norm(xnew)
%|				is less than this unitless value.  default: 0
%|	stop_diff_norm		use norm(.,type) for stop rule
%|				choices: 1 | 2 (default) | inf
%|	stop_grad_tol		stop if norm(grad) / y'W y < tol; default: 0
%|	stop_grad_norm		which norm(grad) to use.  default: 2
%|	chat	0|1		verbosity (default 0)
%|
%| out
%|	xs	[np niter]	estimates each iteration
%|	info	[niter 3]	gamma, step size, time each iteration
%|
%| Copyright Jan 1998, Jeff Fessler, University of Michigan

if nargin == 1 && streq(x, 'test'), qpwls_pcg1_test0, return, end
if nargin < 5, ir_usage, end

% defaults
arg.precon = 1;
arg.niter = 1;
arg.isave = [];
arg.userfun = @userfun_default;
arg.userarg = {};
arg.key = 1;
arg.dircheck = true; % default is to check descent direction
arg.stop_diff_tol = 0;
arg.stop_diff_norm = 2;
arg.stop_grad_tol = 0;
arg.stop_grad_norm = 2;
arg.chat = 0;

arg = vararg_pair(arg, varargin, 'subs', ...
{'stop_threshold', 'stop_diff_tol'; 'stop_norm_type', 'stop_diff_norm'});

arg.isave = iter_saver(arg.isave, arg.niter);
if arg.stop_diff_tol
	norm_diff = @(x) norm(x, arg.stop_diff_norm);
end
if arg.stop_grad_tol
	norm_grad = @(g) norm(g, arg.stop_grad_norm) / reale(yi' * (W * yi));
end

if ~isreal(yi) && ~isequal(arg.precon, 1)
	persistent warned
	if isempty(warned), warned = 0; end
	if ~warned
		warning 'not 100% sure about the complex preconditioned case'
		warned = 1;
	end
end

cpu etic
if isempty(x), x = zeros(ncol(A),1); end

x = x(:);
np = length(x);
xs = zeros(np, length(arg.isave), 'single');
if any(arg.isave == 0)
	xs(:, arg.isave == 0) = single(x);
end

%info = zeros(arg.niter, ?); % trick: do not initialize because size may change

% initialize projections
ticker(mfilename, 1, arg.niter)
Ax = A * x;
Cx = C * x;

% iterate
for iter = 1:arg.niter
	ticker(mfilename, iter, arg.niter)

	% (negative) gradient
	ngrad = A' * (W * (yi-Ax)) - C' * Cx;

	if arg.stop_grad_tol && norm_grad(ngrad) < arg.stop_grad_tol
		if arg.chat
			printm('stop at iteration %d with grad %g < %g', ...
				iter, norm_grad(ngrad), arg.stop_grad_tol)
		end
		if isequal(arg.isave, arg.niter) % saving last iterate only?
			xs = single(x); % save 'final' iterate
		else % saving many iterates?
			xs(:, arg.isave > iter) = []; % clear out unused
		end
	return
	end

	% preconditioned gradient
	pregrad = arg.precon * ngrad;

	% search direction
	newinprod = ngrad' * pregrad;
	% fix: should i take the real part?
	newinprod = reale(newinprod, 'warn', 'inprod');
	if iter == 1
		ddir = pregrad;
		gamma = 0;
	else
		if oldinprod == 0
			warn 'inprod=0. going nowhere!'
			gamma = 0;
		else
			gamma = newinprod / oldinprod;	% Fletcher-Reeves
%			gamma = (newinprod - oldgrad' * pregrad) / oldinprod;
		end
		ddir = pregrad + gamma * ddir;
	end
	oldgrad = ngrad;
	oldinprod = newinprod;

	% check if descent direction
	if arg.dircheck && real(ddir' * ngrad) < 0
		warn 'wrong direction; try using stop_grad_tol'
		ratio = norm(ngrad(:), arg.stop_grad_norm) / (yi'*W*yi);
		pr ratio % see how small it is
		if arg.key, keyboard, end
	end

	% step size in search direction
	Adir = A * ddir;
	Cdir = C * ddir;

	denom = Adir'*(W*Adir) + Cdir'*Cdir;
	denom = reale(denom, 'error', 'denom');
	if denom == 0
		warning 'found exact solution??? step=0 now!?'
		step = 0;
	else
		step = (ddir' * ngrad) / denom;
%		step = reale(step, 'warn', 'step');
		step = real(step); % real step sizes seems only logical
	end

	if step < 0
		warning 'downhill?'
		if arg.key, keyboard, end
	end

	% update
	Ax = Ax + step * Adir;
	Cx = Cx + step * Cdir;
	x = x + step * ddir;

	if any(arg.isave == iter)
		xs(:, arg.isave == iter) = single(x);
	end
	info(iter,:) = arg.userfun(x, arg.userarg{:});

	% check norm(xnew-xold) / norm(xnew) vs threshold
	if arg.stop_diff_tol && ...
		norm_diff(step * ddir) / norm_diff(x) < arg.stop_diff_tol
		if arg.chat
			ratio = norm_diff(step * ddir) / norm_diff(x);
			printm('stop at iteration %d with diff %g < %g', ...
				iter, ratio, arg.stop_diff_tol)
		end
		if isequal(arg.isave, arg.niter) % saving last iterate only?
			xs = single(x); % save the 'final' iterate
		else % saving many iterates?
			xs(:, arg.isave > iter) = []; % clear out unused
		end
	return
	end
end

% default user function.
% using this evalin('caller', ...) trick, one can compute anything of interest
function out = userfun_default(x, varargin)
gamma = evalin('caller', 'gamma');
step = evalin('caller', 'step');
out = [gamma step cpu('etoc')];

% qpwls_pcg1_test0()
function qpwls_pcg1_test0
mask = true([8 7]); mask(1) = false;
A = Gblur(mask, 'psf', ones(3)/9);
tmp = ones(size(mask));
A = Gdiag(tmp(mask), 'mask', mask);
xtrue = zeros(size(mask), 'single');
xtrue(end/2, round(end/2)) = 1;
y = A * xtrue(mask);
beta = 2^-7;
beta = 2^-2;
R = Reg1(mask, 'beta', beta, 'order', 1);
qpwls_psf(A, R.C, 1, mask, 1, 'loop', 0);
hess = full(A' * A + R.C' * R.C);
xhat = hess \ (A' * y);
xhat = embed(xhat, mask);
pr fwhm2(xhat)
im clf, im(xhat)

% user functions for tracking time and distance to a reference image
f.userfun = @(x, xref) [cpu('etoc') norm(x(:) - xref(:))];
f.userarg = {xhat(mask)}; % reference image just for testing

xinit = 0 * mask;
xpcg = qpwls_pcg1(xinit(mask(:)), A, 1, y, R.C, 'niter', 100, ...
	'userfun', f.userfun, 'userarg', f.userarg, ...
        'stop_grad_tol', 1e-8, 'stop_grad_norm', 2, ...
        'stop_diff_tol', 0e-6, 'stop_diff_norm', 2, 'chat', 1);
xpcg = embed(xpcg, mask);

im plc 2 2
im(1, xtrue)
im(2, xhat)
im(3, xpcg)
im(4, xpcg - xhat)
equivs(xpcg, xhat)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
x = load('ex3x.dat');
y = load('ex3y.dat');
m = length(y);
x = [ones(m, 1), x];

# figure;
# plot(x(:,2), y, 'x');
# ylabel('Price in $');
# xlabel('Living area in square feet');

# figure;
# plot(x(:,3), y, 'x');
# ylabel('Price in $');
# xlabel('Number of bedrooms');

sigma = std(x);
mu = mean(x);
x(:,2) = (x(:,2) - mu(2)) ./ sigma(2);
x(:,3) = (x(:,3) - mu(3)) ./ sigma(3);

num_iterations = 50;

alphas = [0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1, 1.3, 3, 10];
J = zeros(num_iterations, length(alphas));

for iter = 1:length(alphas)
  theta = zeros(size(x(1,:)))';
  alpha = alphas(iter);
  for num_iterations = 1:50
    J(num_iterations, iter) = (1/2*m) * (x * theta - y)' * (x * theta - y);
    theta = theta - alpha * (1/m) * x' * (x * theta - y);
  end
end

figure;
hold on;

for iter = 1:length(alphas)-2
  alpha = alphas(iter);

  % now plot J
  % technically, the first J starts at the zero-eth iteration
  % but Matlab/Octave doesn't have a zero index
  plot(0:49, J(1:50, iter), '-');
  xlabel(['Number of iterations, alpha = ' num2str(alpha)]);
  ylabel('Cost J');
end

alpha = 1;
theta = zeros(size(x(1,:)))';
for num_iterations = 1:100
  theta = theta - alpha * (1/m) * x' * (x * theta - y);
end

theta

x1 = ([2, 1650, 3] - mu) ./ (sigma + [1,0,0]);
y_predicted_1 = x1 * theta

% Normal Equations

% Reload data
x = load('ex3x.dat');
y = load('ex3y.dat');
m = length(y);
x = [ones(m, 1), x];

theta = (x' * x) \ (x' * y)

x2 = [1 1650 3];
y_predicted_2 = x2 * theta

author: zellyn


text:
function obj = configureDynamics(obj, varargin)
    % Set the symbolic expressions of the dynamics of the rigid body system
    %
    % Parameters:
    % varargin: extra options @type varargin

    opts = struct(varargin{:});
    if isfield(opts, 'DelayCoriolisSet')
        delay_set = opts.DelayCoriolisSet;
    else
        delay_set = false;
    end

    if isfield(opts, 'OmitCoriolisSet')
        omit_set = opts.OmitCoriolisSet;
    else
        omit_set = false;
    end

    n_link = length(obj.Links);
    links = getTwists(obj.Links);
    for i=1:n_link
        links{i}.Mass = obj.Links(i).Mass;
        links{i}.Inertia = obj.Links(i).Inertia;
    end

    % set the inertia matrix
    tic
    %     De = eval_math_fun('InertiaMatrix',[links,{obj.numState}]);
    n_link = length(links);
    De = cell(n_link,1);
    m = 100;
    fprintf('Evaluating the inertia matrix D(q): \t');
    bs = '\b';
    sp = ' ';
    msg = '%d percents completed.\n'; % carriage return included in case error interrupts loop
    msglen = length(msg)-3; % compensate for conversion character and carriage return
    fprintf(1,sp(ones(1,ceil(log10(m+1))+msglen)));

    for i=1:n_link
        k = floor(i*100/n_link);
        fprintf(1,[bs(mod(0:2*(ceil(log10(k+1))+msglen)-1,2)+1) msg],k);
        De{i} = eval_math_fun('InertiaMatrixByLink',[links(i),{obj.numState}]);
    end

    De_motor = zeros(obj.numState);
    for i=1:obj.numState
        if ~isempty(obj.Joints(i).Actuator)
            actuator = obj.Joints(i).Actuator;

            if ~isempty(actuator.Inertia) && ~isempty(actuator.Ratio)
                % reflected motor inertia: I*r^2
                De_motor(i,i) = actuator.Inertia * actuator.Ratio^2;
            end
        end
    end

    %     fprintf(1,sp(ones(1,40)));
    Mmat_full = [De; {De_motor}];
    obj.setMassMatrix(Mmat_full);

    toc

    Qe = obj.States.x;
    dQe = obj.States.dx;

    %|@note !!! The minus sign !!!
    % Fvec appears at the right hand side (typically -C(q,dq)dq - G(q)
    % appears at the left hand side of the Euler-Lagrangian Equation)
    if omit_set
        tic
        fprintf('Evaluating the gravity vector G(q): \t');
        Ge = SymFunction(['Ge_vec_',obj.Name],-eval_math_fun('GravityVector',[links,{Qe}]),{Qe,dQe});
        toc
        vf = {Ge};
    else
        tic
        n_mass_mat = numel(obj.Mmat);
        n_dof      = obj.numState;
        Ce1 = cell(n_mass_mat*n_dof,1);
        Ce2 = cell(n_mass_mat*n_dof,1);
        Ce3 = cell(n_mass_mat*n_dof,1);
        m = 100;
        fprintf('Evaluating the coriolis term C(q,dq)dq: \t');
        bs = '\b';
        sp = ' ';
        msg = '%d percents completed.\n'; % carriage return included in case error interrupts loop
        msglen = length(msg)-3; % compensate for conversion character and carriage return
        fprintf(1,sp(ones(1,ceil(log10(m+1))+msglen)));
        for i=1:n_mass_mat
            for j=1:n_dof
                k = floor(((i-1)*3+(j-1)+1)*100/(3*n_mass_mat*n_dof));
                fprintf(1,[bs(mod(0:2*(ceil(log10(k+1))+msglen)-1,2)+1) msg],k);
                Ce1{(i-1)*n_dof+j} = SymFunction(...
                    ['Ce1_vec_L',num2str(i),'_J',num2str(j),'_', obj.Name],...
                    eval_math_fun('InertiaToCoriolisPart1New',{Mmat_full{i},Qe,dQe,j},[],'DelayedSet',delay_set),...
                    {Qe,dQe});
                k = floor(((i-1)*3+(j-1)+2)*100/(3*n_mass_mat*n_dof));
                fprintf(1,[bs(mod(0:2*(ceil(log10(k+1))+msglen)-1,2)+1) msg],k);
                Ce2{(i-1)*n_dof+j} = SymFunction(...
                    ['Ce2_vec_L',num2str(i),'_J',num2str(j),'_', obj.Name],...
                    eval_math_fun('InertiaToCoriolisPart2New',{Mmat_full{i},Qe,dQe,j},[],'DelayedSet',delay_set),...
                    {Qe,dQe});
                k = floor(((i-1)*3+(j-1)+3)*100/(3*n_mass_mat*n_dof));
                fprintf(1,[bs(mod(0:2*(ceil(log10(k+1))+msglen)-1,2)+1) msg],k);
                Ce3{(i-1)*n_dof+j} = SymFunction(...
                    ['Ce3_vec_L',num2str(i),'_J',num2str(j),'_', obj.Name],...
                    eval_math_fun('InertiaToCoriolisPart3New',{Mmat_full{i},Qe,dQe,j},[],'DelayedSet',delay_set),...
                    {Qe,dQe});
            end
        end
        fprintf(1,sp(ones(1,40)));
        toc
        % We add dQe as dependent variables to the gravity vector in order to have a same structure
        % as the corilios forces
        tic
        fprintf('Evaluating the gravity vector G(q): \t');
        Ge = SymFunction(['Ge_vec_',obj.Name],-eval_math_fun('GravityVector',[links,{Qe}]),{Qe,dQe});
        toc
        vf = [Ce1;Ce2;Ce3;{Ge}];
    end

    obj.setDriftVector(vf);

end

author: ayonga

repo: frost-dev

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


text:
function [segmentation] = ft_datatype_segmentation(segmentation, varargin)

% FT_DATATYPE_SEGMENTATION describes the FieldTrip MATLAB structure for segmented
% voxel-based data and atlasses. A segmentation can either be indexed or probabilistic
% (see below).
%
% A segmentation is a volumetric description which is usually derived from an anatomical
% MRI, which describes for each voxel the tissue type. It for example distinguishes
% between white matter, grey matter, csf, skull and skin. It is mainly used for masking
% in visualization, construction of volume conduction models and for construction of
% cortical sheets. An volume-based atlas is basically a very detailed segmentation with
% an anatomical label for each voxel.
%
% For example, the AFNI TTatlas+tlrc segmented brain atlas (which can be created
% with FT_READ_ATLAS) looks like this
%
%              dim: [161 191 141]        the size of the 3D volume in voxels
%        transform: [4x4 double]         affine transformation matrix for mapping the voxel coordinates to head coordinate system
%         coordsys: 'tal'                the transformation matrix maps the voxels into this (head) coordinate system
%             unit: 'mm'                 the units in which the coordinate system is expressed
%           brick0: [161x191x141 uint8]  integer values from 1 to N, the value 0 means unknown
%           brick1: [161x191x141 uint8]  integer values from 1 to M, the value 0 means unknown
%      brick0label: {Nx1 cell}
%      brick1label: {Mx1 cell}
%
%
% An example segmentation with binary values that can be used for construction of a
% BEM volume conduction model of the head looks like this
%
%           dim: [256 256 256]         the dimensionality of the 3D volume
%     transform: [4x4 double]          affine transformation matrix for mapping the voxel coordinates to head coordinate system
%      coordsys: 'ctf'                 the transformation matrix maps the voxels into this (head) coordinate system
%          unit: 'mm'                  the units in which the coordinate system is expressed
%         brain: [256x256x256 logical] binary map representing the voxels which belong to the brain
%         scalp: [256x256x256 logical] binary map representing the voxels which belong to the scalp
%         skull: [256x256x256 logical] binary map representing the voxels which belong to the skull
%
% An example of a whole-brain anatomical MRI that was segmented using FT_VOLUMESEGMENT
% looks like this
%
%         dim: [256 256 256]         the size of the 3D volume in voxels
%   transform: [4x4 double]          affine transformation matrix for mapping the voxel coordinates to head coordinate system
%    coordsys: 'ctf'                 the transformation matrix maps the voxels into this (head) coordinate system
%        unit: 'mm'                  the units in which the coordinate system is expressed
%        gray: [256x256x256 double]  probabilistic map of the gray matter
%       white: [256x256x256 double]  probabilistic map of the white matter
%         csf: [256x256x256 double]  probabilistic map of the cerebrospinal fluid
%
% The examples above demonstrate that a segmentation can be indexed, i.e. consisting of
% subsequent integer numbers (1, 2, ...) or probabilistic, consisting of real numbers
% ranging from 0 to 1 that represent probabilities between 0% and 100%. An extreme case
% is one where the probability is either 0 or 1, in which case the probability can be
% represented as a binary or logical array.
%
% The only difference to the volume data representation is that the segmentation
% structure contains the additional fields xxx and xxxlabel. See FT_DATATYPE_VOLUME for
% further details.
%
% Required fields:
%   - dim, transform
%
% Optional fields:
%   - coordsys, unit
%
% Deprecated fields:
%   - none
%
% Obsoleted fields:
%   - none
%
% Revision history:
% (2012/latest) The explicit distunction between the indexed and probabilistic
% representation was made. For the indexed representation the additional
% xxxlabel cell-array was introduced.
%
% (2005) The initial version was defined.
%
% See also FT_DATATYPE, FT_DATATYPE_VOLUME, FT_DATATYPE_PARCELLATION

% Copyright (C) 2012, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% get the optional input arguments, which should be specified as key-value pairs
version           = ft_getopt(varargin, 'version', 'latest');
segmentationstyle = ft_getopt(varargin, 'segmentationstyle');  % can be indexed or probabilistic
hasbrain          = ft_getopt(varargin, 'hasbrain', 'no');     % no means that it is not required, if present it won't be removed

% convert from string into boolean
hasbrain = istrue(hasbrain);

if strcmp(version, 'latest')
  segversion = '2012';
  volversion = 'latest';
  clear version
else
  segversion = version;
  volversion = version;
  clear version
end

if isempty(segmentation)
  return;
end

switch segversion
  case '2012'
    % determine whether the style of the input fields is probabilistic or indexed
    fn = fieldnames(segmentation);
    %fn = setdiff(fn, 'inside'); % exclude the inside field from any conversions
    fn(strcmp(fn, 'inside')) = [];
    [indexed, probabilistic] = determine_segmentationstyle(segmentation, fn, segmentation.dim);

    % ignore the fields that do not contain a segmentation
    sel = indexed | probabilistic;
    fn            = fn(sel);
    indexed       = indexed(sel);
    probabilistic = probabilistic(sel);

    % convert from an exclusive to cumulative representation
    % this is only only for demonstration purposes
    % for i=1:length(sel)
    %   segmentation.(fn{sel(i)}) = volumefillholes(segmentation.(fn{sel(i)}));
    % end

    [dum, i] = intersect(fn, {'scalp', 'skull', 'brain'});
    if numel(i)==3
      % put them in the preferred order
      fn(i) = {'scalp', 'skull', 'brain'};
    end
    [dum, i] = intersect(fn, {'skin', 'skull', 'brain'}); % this is not likely
    if numel(i)==3
      % put them in the preferred order
      fn(i) = {'skin', 'skull', 'brain'};
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % ensure that the segmentation is internally consistent
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if any(probabilistic)
      segmentation = fixsegmentation(segmentation, fn(probabilistic), 'probabilistic');
    end
    if any(indexed)
      segmentation = fixsegmentation(segmentation, fn(indexed), 'indexed');
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % convert the segmentation to the desired style
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if isempty(segmentationstyle)
      % keep it as it is
    elseif strcmp(segmentationstyle, 'indexed') && any(probabilistic)
      segmentation  = convert_segmentationstyle(segmentation, fn(probabilistic), segmentation.dim, 'indexed');
      indexed(probabilistic)       = true;  % these are now indexed
      probabilistic(probabilistic) = false; % these are now indexed
    elseif strcmp(segmentationstyle, 'probabilistic') && any(indexed)
      segmentation  = convert_segmentationstyle(segmentation, fn(indexed), segmentation.dim, 'probabilistic');
      probabilistic(indexed) = true;  % these are now probabilistic
      indexed(indexed)       = false; % these are now probabilistic
    end % converting between probabilistic and indexed

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % add the brain if requested
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if hasbrain
      if all(indexed)
        fn = fieldnames(segmentation);
        sel = false(size(fn));
        for i=1:numel(fn)
          sel(i) = any(strcmp(fn, [fn{i} 'label']));
        end
        fn = fn(sel);

        if numel(fn)>1
          ft_error('cannot construct a brain mask on the fly; this requires a single indexed representation');
        else
          seg      = segmentation.(fn{1});
          seglabel = segmentation.([fn{1} 'label']);
          if ~any(strcmp(seglabel, 'brain'))
            threshold = 0.5;
            smooth    = 5;
            % ensure that the segmentation contains the brain mask, if not then construct it from gray+white+csf
            if length(intersect(seglabel, {'gray' 'white' 'csf'}))~=3
              ft_error('cannot construct a brain mask on the fly; this requires gray, white and csf');
            end
            gray  = seg==find(strcmp(seglabel, 'gray'));
            white = seg==find(strcmp(seglabel, 'white'));
            csf   = seg==find(strcmp(seglabel, 'csf'));
            brain = gray + white + csf;
            clear gray white csf seg
            brain = volumesmooth(brain,    smooth,    'brain');
            brain = volumethreshold(brain, threshold, 'brain');
            % store it in the output
            segmentation.brain = brain;
          end % try to construct the brain
        end

      elseif all(probabilistic)
        if ~isfield(segmentation, 'brain')
          if ~all(isfield(segmentation, {'gray' 'white' 'csf'}))
            ft_error('cannot construct a brain mask on the fly; this requires gray, white and csf');
          end
          threshold = 0.5;
          smooth    = 5;
          % ensure that the segmentation contains the brain mask, if not then construct it from gray+white+csf tissue probability maps
          gray  = segmentation.gray;
          white = segmentation.white;
          csf   = segmentation.csf;
          brain = gray + white + csf;
          clear gray white csf
          brain = volumesmooth(brain,    smooth,    'brain');
          brain = volumethreshold(brain, threshold, 'brain');
          % store it in the output
          segmentation.brain = brain;
        end
      else
        ft_error('cannot construct a brain mask on the fly; this requires a uniquely indexed or a uniquely probabilitic representation');
      end
    end % if hasbrain

  case '2005'
    % the only difference is that the indexed representation for xxx did not have the xxxlabel field prior to the 2012 version
    fn = fieldnames(segmentation);
    sel = ~cellfun(@isempty, regexp(fn, 'label$'));
    segmentation = rmfield(segmentation, fn(sel));
    % furthermore it corresponds to the oldest version of the volume representation
    volversion = '2003';

  otherwise
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ft_error('unsupported version "%s" for segmentation datatype', segversion);
end

% the segmentation is a speciat type of volume structure, so ensure that it also fulfills the requirements for that
segmentation = ft_datatype_volume(segmentation, 'version', volversion);

% For the pass through ft_datatype_volume it is perhaps necessary to remove
% the fields that are specific for the segmentation and add them later again.
% At this moment ft_datatype_volume nicely passes all fields, so there is no
% special handling of the segmentation fields needed.

author: spm


text:
% STD_SAVEDAT - save measure for computed data
%
% Usage: std_savedat( filename, structure);
%
% Authors:  Arnaud Delorme, SCCN, INC, UCSD, 2006-

% Copyright (C) Arnaud Delorme, SCCN, INC, UCSD, October 11, 2004, arno@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function std_savedat( tmpfile, structure)

    delims = find( tmpfile == '.');
    if ~isfield(structure, 'datafile') && ~isfield(structure, 'datafiles')
        structure.datafile = [ tmpfile(1:delims(end)-1) '.set' ];
    end

    % fix reading problem (bug 764)
    tmpfile2  = which(tmpfile);
    if isempty(tmpfile2), tmpfile2 = tmpfile; end;
    tmpfile = tmpfile2;

    eeglab_options;
    if option_saveversion6
        try
            save('-v6' , tmpfile, '-struct', 'structure');
        catch
            fields = fieldnames(structure);
            for i=1:length(fields)
                eval([ fields{i} '=structure.'  fields{i} ';']);
            end
            save('-mat', tmpfile, fields{:});
        end
    else
        save('-v7.3' , tmpfile, '-struct', 'structure');
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
%
%% PRIMES_SCRIPT.M
%
fprintf ( 1, '\n' );
fprintf ( 1, 'PRIMES_SCRIPT:\n' );
fprintf ( 1, '  MATLAB version.\n' );
tic;
s4 = sum ( isprime(3:2:10^4) ) + 1
t4 = toc
%
tic;
s5 = sum ( isprime(3:2:10^5) ) + 1
t5 = toc
%
%tic;
%s6 = sum ( isprime(3:2:10^6) ) + 1
%t6 = toc
%
tic;
p4 = length ( primes(10^4) )
t4 = toc;
%
tic;
p5 = length ( primes(10^5) )
t5 = toc;
%
tic;
p6 = length ( primes(10^6) )
t6 = toc;
output = fopen ( 'primes_report.txt', 'w' );
fprintf ( output, 'There are %d prime numbers up to 10^4.\n', s4 );
fprintf ( output, 'There are %d prime numbers up to 10^5.\n', s5 );
%fprintf ( output, 'There are %d prime numbers up to 10^6.\n', s6 );
fclose ( output );
%
%  Let's find out what files exist now.
%
ls
%
%  Terminate.
%
fprintf ( 1, '\n' );
fprintf ( 1, 'PRIMES_SCRIPT:\n' );
fprintf ( 1, '  Normal end of execution.\n' );

author: johannesgerer

repo: jburkardt-m


text:
N = 3;
x = rand(N,2); % each row is a feature vector
m = mean(x,1);
xc = x-repmat(m, N, 1);

C = eye(N) - (1/N)*ones(N,N);
xc2 = C*x;
assert(approxeq(xc, xc2))

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%% 2D GEOMETRIC COLLISION AVOIDANCE AGENT (agent_2D_VO.m) %%%%%%%%%%%%%%%%%

% Author: James A. Douthwaite

classdef agent_2D_RVO < agent_2D_VO & agent_RVO
%% INITIALISE THE AGENT SPECIFIC PARAMETERS
    properties
        % PROPERTIES UNIQUE TO THE 2D RECIPROCAL VO METHOD
    end
%%  CLASS METHODS
    methods
        % CONSTRUCTOR
        function obj = agent_2D_RVO(varargin)
            % Construct the agent object and initialise with the following
            % specific parameters.

            % CALL THE SUPERCLASS CONSTRUCTOR
            obj@agent_2D_VO(varargin);

            % //////////////// Check for user overrides ///////////////////
            % - It is assumed that overrides to the properties are provided
            %   via the varargin structure.
            [obj] = obj.ApplyUserOverrides(varargin);
            % /////////////////////////////////////////////////////////////
        end

        % MAIN CYCLE IS INHERITED FROM THE SUPERCLASS
        % The waypoint handling, dynamics and update procedures must also
        % be the same as the super-class.

        % CALCULATE THE NECESSARY 2D AVOIDANCE VELOCITY
        function [headingVector,speed] = GetAvoidanceCorrection(obj,dt,desiredVelocity,visualiseProblem)
            % This function calculates the 2D avoidance velocity command
            % and returns it to be achieved by the controller.

            % Input sanity check
            assert(numel(desiredVelocity) == 2,'Desired velocity must be 2D');

            % AGENT KNOWLEDGE (2D)
            [p_i,v_i,r_i] = obj.GetAgentMeasurements();

            % Define the obstacle list
            obstacleIDs = [obj.MEMORY([obj.MEMORY.type] ~= OMAS_objectType.waypoint).objectID];

            % MOVE THROUGH THE PRIORITISED OBSTACLE SET
            VO = [];
            for item = 1:numel(obstacleIDs)
                % Get object data from memory structure
                p_j = obj.GetLastMeasurementByID(obstacleIDs(item),'position');
                v_j = obj.GetLastMeasurementByID(obstacleIDs(item),'velocity');
                r_j = obj.GetLastMeasurementByID(obstacleIDs(item),'radius');

                % Neighbour conditions
                neighbourConditionA = item < obj.maxNeighbours;            % Maximum number of neighbours
                neighbourConditionB = norm(p_j) < obj.neighbourDist;       % [CONFIRMED]
                neighbourConditionC = ~any(isnan(v_j));                    % Wait for a valid velocity reading
                if ~neighbourConditionB || ~neighbourConditionC
                    continue
                end

                % OBSTACLE KNOWLEDGE
                p_j = p_j + p_i;
                v_j = v_j + v_i;                                           % Convert relative parameters to absolute
                tau_j = 0;

                % OBSTACLE TYPE BEHAVIOUR
                type_j = obj.GetLastMeasurementByID(obstacleIDs(item),'type');
                if OMAS_objectType.agent == type_j
                    % DEFINE RVO AS AGENT BEHAVIOUR
                    [VO_j] = obj.define2DReciprocalVelocityObstacle(p_i,v_i,r_i,p_j,v_j,r_j,tau_j,visualiseProblem);
                else
                    % OBSTACLE BEHAVIOUR
                    [VO_j] = obj.define2DVelocityObstacle(p_i,v_i,r_i,p_j,v_j,r_j,tau_j,visualiseProblem);
                end
                % CONCATINATE THE VO SET
                VO = [VO,VO_j];
            end

            % THE CLEAR PATH STRATEGY
            [avoidanceVelocity] = obj.strategy_clearPath(v_i,desiredVelocity,VO,visualiseProblem);

            % SPECIAL CASE- VELOCITY MAGNITUDE IS ZERO
            speed = norm(avoidanceVelocity);
            headingVector = avoidanceVelocity/speed;
            if isnan(headingVector)
                headingVector = [1;0];  % Retain previous heading
            end

            % PLOT THE VECTOR CONSTRUCT
            if visualiseProblem && obj.objectID == visualiseProblem
                % PLOT THE LEADING TANGENT VECTOR
                OMAS_axisTools.drawTriad([p_i;0],eye(3));
                % CURRENT VELOCITY
                q = quiver(gca,p_i(1),p_i(2),v_i(1),v_i(2),'m');
                q.AutoScaleFactor = 1;
                % DESIRED VELOCITY
                q = quiver(gca,p_i(1),p_i(2),desiredVelocity(1),desiredVelocity(2),'g');
                q.AutoScaleFactor = 1;
                % AVOIDANCE VELOCITY
                q = quiver(gca,p_i(1),p_i(2),avoidanceVelocity(1),avoidanceVelocity(2),'b');
                q.AutoScaleFactor = 1;
            end
        end
    end
    % /////////////// RECIPRICOL VELOCITY OBSTACLE METHODS ////////////////
    methods (Access = public)
        % DEFINE THE 2D RECIPROCAL VELOCITY OBSTACLE (RVO)
        function [RVO] = define2DReciprocalVelocityObstacle(obj,p_i,v_i,r_i,p_j,v_j,r_j,tau_j,plotOn)
            % This function assembles the reciprocal velocity obstacle in 2D.

            % MAP THE 2D INPUTS TO 3D
            p_i = [p_i;0]; v_i = [v_i;0];
            p_j = [p_j;0]; v_j = [v_j;0];
            % CALL THE 3D VO GENERATION FUNCTION
            [RVO] = obj.define3DReciprocalVelocityObstacle(p_i,v_i,r_i,p_j,v_j,r_j,tau_j,plotOn);
            % MAP THE 3D INPUTS BACK TO 2D INPUTS
            RVO.apex = RVO.apex(1:2,1);
            RVO.axisUnit = RVO.axisUnit(1:2,1);
            RVO.leadingEdgeUnit = RVO.leadingEdgeUnit(1:2,1);
            RVO.trailingEdgeUnit = RVO.trailingEdgeUnit(1:2,1);
        end
    end
end
% AGENT STATE VECTOR [x;y;phi;xdot;ydot;phidot]

repo: OpenMAS


text:
% This is an example of using this software to train a binary
% decision tree classifier, and use this classifier to classify
% testing data.

% Copyright (C) 2012 Quan Wang <wangq10@rpi.edu>,
% Signal Analysis and Machine Perception Laboratory,
% Department of Electrical, Computer, and Systems Engineering,
% Rensselaer Polytechnic Institute, Troy, NY 12180, USA
%
% You are free to use this software for academic purposes if you cite our paper:
% Q. Wang, Y. Ou, A.A. Julius, K.L. Boyer, M.J. Kim,
% Tracking tetrahymena pyriformis cells using decision trees,
% in: 2012 International Conference on Pattern Recognition, Tsukuba Science City, Japan.
%
% For commercial use, please contact the authors.

clear;clc;close all;

%% settings
Depth=5; % maximal depth of decision tree
Splits=100; % number of candidate thresholds at each node
MinNode=10; % minimal size of a non-leaf node

%% training
load TrainingData.mat;

tic;
T=create01Tree(X,Y,Depth,Splits,MinNode);
t1=toc;

clear X Y;

%% testing
load TestingData.mat;

tic;
y=[];
for i=1:size(X,1)
    x=X(i,:);
    y(i,1)=decide01Tree(x,T);
end
t2=toc;

%% evaluation
errorRate=sum(abs(y-Y))/max(size(Y));
fprintf('Error rate = %.4f\n',errorRate);
fprintf('Training time = %.4f seconds\n',t1);
fprintf('Testing time = %.4f seconds\n',t2);

author: Sable

repo: mcbench-benchmarks


text:
function definput = arg_fwt(definput)

%definput.flags.ext=  {'per','zpd','sym','symw','asym','asymw','ppd','sp0'};
definput.flags.ext = {'per','zero','even','odd'};

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [ SigmaInv] = CalcSigmaCCNFflat(alphas, betas, n, precalc_B_without_beta, precalc_eye, precalc_zeros)
%CALCSIGMACCNFflat Computing SigmaInv matrices (represented as a vector as
%it is a symmetric matrix)

    A = sum(alphas) .* precalc_eye;

    % calculating the B + C from the paper (here referred to as B)
    Btmp = precalc_B_without_beta * betas;

    B = precalc_zeros;
    on = tril(true(n,n));
    B(on) = Btmp;
    B = B';
    B(on) = Btmp;

    SigmaInv = 2 * (A + B);

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function acosh_test ( )

%*****************************************************************************80
%
%% ACOSH_TEST tests R4_ACOSH and R8_ACOSH.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../test_values' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'ACOSH_TEST:\n' );
  fprintf ( 1, '  Test ARCCOSH_VALUES, R4_ACOSH, R8_ACOSH\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '             X      ARCCOSH(X)\n' );
  fprintf ( 1, '                   R4_ACOSH(X)        Diff\n' );
  fprintf ( 1, '                   R8_ACOSH(X)        Diff\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx1 ] = arccosh_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = r4_acosh ( single ( x ) );
    fx3 = r8_acosh ( x );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %14.4f  %14.6g\n', x, fx1 );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx2, abs ( fx1 - fx2 ) );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx3, abs ( fx1 - fx3 ) );

  end

  rmpath ( '../test_values' );

  return
end


text:
function rects = get_rect_from_pts(allpoints)
%from 5 points to rect -- l, t, r, b
if false
    rects = fast_get_rect_from_pts(allpoints);
else
    std_points = [
        0.2 0.2
        0.8 0.2
        0.5 0.5
        0.3 0.75
        0.7 0.75];

    rects = nan(size(allpoints, 1), 4);

    tic;
    for i = 1:size(allpoints, 1)
        try
            points = allpoints(i, :);
            points = reshape(points, 2, [])';
            points = points(1:5, :);
            t = cp2tform(double(points), std_points, 'similarity');

            [xc, yc] = tforminv(t, 0.5, 0.5);
            [xtl, ytl] = tforminv(t, 0, 0);
            [xtr, ytr] = tforminv(t, 1, 0);

            w = sqrt((xtl-xtr).^2+(ytl-ytr).^2);
            rect = [xc-w/2, yc-w/2, xc+w/2, yc+w/2];
            rect = round(rect);
            rects(i, :) = rect;
            if (mod(i, 10000) == 0)
                fprintf('%d/%d... ', i, size(allpoints, 1));
                toc;
            end
        catch
            fprintf('Fail to process: %d!\n', i);
            continue;
        end
    end
end

end

author: liuyuisanai

repo: RSA-for-object-detection


text:
function dx = lorenz(t,x,Beta)
dx = [
Beta(1)*(x(2)-x(1));
x(1)*(Beta(2)-x(3))-x(2);
x(1)*x(2)-Beta(3)*x(3);
];

author: dynamicslab


text:
function [h,Coord,Color,height] = som_dendrogram(Z,varargin)

%SOM_DENDROGRAM Visualize a dendrogram.
%
% [h,Coord,Color,height] = som_dendrogram(Z, [[argID,] value, ...])
%
%  Z = som_linkage(sM);
%  som_dendrogram(Z);
%  som_dendrogram(Z,sM);
%  som_dendrogram(Z,'coord',co);
%
%  Input and output arguments ([]'s are optional):
%   h        (vector) handle to the arc lines
%   Z        (matrix) size n-1 x 1, the hierarchical cluster matrix
%                     returned by functions like LINKAGE and SOM_LINKAGE
%                     n is the number of original data samples.
%   [argID,  (string) See below. The values which are unambiguous can
%    value]  (varies) be given without the preceeding argID.
%   Coord    (matrix) size 2*n-1 x {1,2}, the coordinates of the
%                     original data samples and cluster nodes used
%                     in the visualization
%   Color    (matrix) size 2*n-1 x 3, the colors of ...
%   height   (vector) size 2*n-1 x 1, the heights of ...
%
% Here are the valid argument IDs and corresponding values. The values
% which are unambiguous (marked with '*') can be given without the
% preceeding argID.
%   'data'  *(struct) map or data struct: many other optional
%                     arguments require this
%            (matrix) data matrix
%   'coord'  (matrix) size n x 1 or n x 2, the coordinates of
%                     the original data samples either in 1D or 2D
%            (matrix) size 2*n-1 x {1,2}, the coordinates of both
%                     original data samples and each cluster
%           *(string) 'SOM', 'pca#', 'sammon#', or 'cca#': the coordinates
%                     are calculated using the given data and the
%                     required projection algorithm. The '#' at the
%                     end of projection algorithms refers to the
%                     desired output dimension and can be either 1 or 2
%                     (2 by default). In case of 'SOM', the unit
%                     coordinates (given by SOM_VIS_COORDS) are used.
%   'color'  (matrix) size n x 3, the color of the original data samples
%            (matrix) size 2*n-1 x 3, the colors of both original
%                     data samples and each cluster
%            (string) color specification, e.g. 'r.', used for each node
%   'height' (vector) size n-1 x 1, the heights used for each cluster
%            (vector) size 2*n-1 x 1, the heights used for both original
%                     data samples and each cluster
%           *(string) 'order', the order of combination determines height
%                     'depth', the depth at which the combination
%                     happens determines height
%   'linecolor' (string) color specification for the arc color, 'k' by default
%               (vector) size 1 x 3
%
% See also SOM_LINKAGE, DENDROGRAM.

% Copyright (c) 2000 by Juha Vesanto
% Contributed to SOM Toolbox on June 16th, 2000 by Juha Vesanto
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 2.0beta juuso 160600

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% read the arguments

% Z
nd = size(Z,1)+1;
nc = size(Z,1);

% varargin
Coordtype = 'natural'; Coord = []; codim = 1;
Colortype = 'none'; Color = [];
height = [zeros(nd,1); Z(:,3)];
M = [];
linecol = 'k';

i=1;
while i<=length(varargin),
  argok = 1;
  if ischar(varargin{i}),
    switch varargin{i},
     case 'data', i = i + 1; M = varargin{i};
     case 'coord',
      i=i+1;
      if isnumeric(varargin{i}), Coord = varargin{i}; Coordtype = 'given';
      else
	if strcmp(varargin{i},'SOM'), Coordtype = 'SOM';
	else Coordtype = 'projection'; Coord = varargin{i};
	end
      end
     case 'color',
      i=i+1;
      if isempty(varargin{i}), Colortype = 'none';
      elseif ischar(varargin{i}), Colortype = 'colorspec'; Color = varargin{i};
      else Colortype = 'given'; Color = varargin{i};
      end
     case 'height', i=i+1; height = varargin{i};
     case 'linecolor', i=i+1; linecol = varargin{i};
     case 'SOM',
      Coordtype = 'SOM';
     case {'pca','pca1','pca2','sammon','sammon1','sammon2','cca','cca1','cca2'},
      Coordtype = 'projection'; Coord = varargin{i};
     case {'order','depth'}, height = varargin{i};
    end
  elseif isstruct(varargin{i}), M = varargin{i};
  else
    argok = 0;
  end
  if ~argok,
    disp(['(som_dendrogram) Ignoring invalid argument #' num2str(i+1)]);
  end
  i = i+1;
end

switch Coordtype,
 case 'SOM',
  if isempty(M) | ~any(strcmp(M.type,{'som_map','som_topol'})) ,
    error('Cannot determine SOM coordinates without a SOM.');
  end
  if strcmp(M.type,'som_map'), M = M.topol; end
 case 'projection',
  if isempty(M), error('Cannot do projection without the data.'); end
  if isstruct(M),
    if strcmp(M.type,'som_data'), M = M.data;
    elseif strcmp(M.type,'som_map'), M = M.codebook;
    end
  end
  if size(M,1) ~= nd,
    error('Given data must be equal in length to the number of original data samples.')
  end
 case 'given',
  if size(Coord,1) ~= nd & size(Coord,1) ~= nd+nc,
    error('Size of given coordinate matrix does not match the cluster hierarchy.');
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% initialization

% Coordinates
switch Coordtype,
 case 'natural', o = leavesorder(Z)'; [dummy,Coord] = sort(o); codim = 1;
 case 'SOM', Coord = som_vis_coords(M.lattice,M.msize); codim = 2;
 case 'projection',
  switch Coord,
   case {'pca','pca2'},       Coord = pcaproj(M,2);   codim = 2;
   case 'pca1',               Coord = pcaproj(M,1);   codim = 1;
   case {'cca','cca2'},       Coord = cca(M,2,20);    codim = 2;
   case 'cca1',               Coord = cca(M,1,20);    codim = 1;
   case {'sammon','sammon2'}, Coord = sammon(M,2,50); codim = 2;
   case 'sammon1',            Coord = sammon(M,1,50); codim = 1;
  end
 case 'given', codim = min(size(Coord,2),2); % nill
end

if size(Coord,1) == nd,
  Coord = [Coord; zeros(nc,size(Coord,2))];
  for i=(nd+1):(nd+nc),
    leaves = leafnodes(Z,i,nd);
    if any(leaves), Coord(i,:) = mean(Coord(leaves,:),1); else Coord(i,:) = Inf; end
  end
end

% Colors
switch Colortype,
 case 'colorspec', % nill
 case 'none', Color = '';
 case 'given',
  if size(Color,1) == nd,
    Color = [Color; zeros(nc,3)];
    for i=(nd+1):(nd+nc),
      leaves = leafnodes(Z,i,nd);
      if any(leaves), Color(i,:) = mean(Color(leaves,:),1);
      else Color(i,:) = 0.8;
      end
    end
  end
end

% height
if ischar(height),
  switch height,
   case 'order', height = [zeros(nd,1); [1:nc]'];
   case 'depth', height = nodedepth(Z); height = max(height) - height;
  end
else
  if length(height)==nc, height = [zeros(nd,1); height]; end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% draw

% the arcs
lfrom = []; lto = [];
for i=1:nd+nc,
  if i<=nd, ch = [];
  elseif ~isfinite(Z(i-nd,3)), ch = [];
  else ch = Z(i-nd,1:2)';
  end
  if any(ch),
    lfrom = [lfrom; i*ones(length(ch),1)];
    lto = [lto; ch];
  end
end

% the coordinates of the arcs
if codim == 1,
  Lx = [Coord(lfrom), Coord(lto), Coord(lto)];
  Ly = [height(lfrom), height(lfrom), height(lto)];
  Lz = [];
else
  Lx = [Coord(lfrom,1), Coord(lto,1), Coord(lto,1)];
  Ly = [Coord(lfrom,2), Coord(lto,2), Coord(lto,2)];
  Lz = [height(lfrom), height(lfrom), height(lto)];
end

washold = ishold;
if ~washold, cla; end

% plot the lines
if isempty(Lz),
  h = line(Lx',Ly','color',linecol);
else
  h = line(Lx',Ly',Lz','color',linecol);
  if ~washold, view(3); end
  rotate3d on
end

% plot the nodes
hold on
switch Colortype,
 case 'none', % nill
 case 'colorspec',
  if codim == 1, plot(Coord,height,Color);
  else plot3(Coord(:,1), Coord(:,2), height, Color);
  end
 case 'given',
  som_grid('rect',[nd+nc 1],'line','none','Coord',[Coord, height],...
	   'Markersize',10,'Markercolor',Color);
end
if ~washold, hold off, end

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% subfunctions

function depth = nodedepth(Z)

  nd = size(Z,1)+1;
  nc = size(Z,1);
  depth = zeros(nd+nc,1);
  ch = nc+nd-1;
  while any(ch),
    c  = ch(1); ch = ch(2:end);
    if c>nd & isfinite(Z(c-nd,3)),
      chc = Z(c-nd,1:2);
      depth(chc) = depth(c) + 1;
      ch = [ch, chc];
    end
  end
  return;

function inds = leafnodes(Z,i,nd)

  inds = [];
  ch = i;
  while any(ch),
    c  = ch(1); ch = ch(2:end);
    if c>nd & isfinite(Z(c-nd,3)), ch = [ch, Z(c-nd,1:2)]; end
    if c<=nd, inds(end+1) = c; end
  end
  return;

function order = leavesorder(Z)

  nd = size(Z,1)+1;
  order = 2*nd-1;
  nonleaves = 1;
  while any(nonleaves),
    j = nonleaves(1);
    ch = Z(order(j)-nd,1:2);
    if j==1, oleft = []; else oleft = order(1:(j-1)); end
    if j==length(order), oright = []; else oright = order((j+1):length(order)); end
    order = [oleft, ch, oright];
    nonleaves = find(order>nd);
  end
  return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


function y = autopilot(uu,P)
%
% autopilot for mavsim
%
% Modification History:
%   2/11/2010 - RWB
%   5/14/2010 - RWB
%   9/30/2014 - RWB
%

    % process inputs
    NN = 0;
%    pn       = uu(1+NN);  % inertial North position
%    pe       = uu(2+NN);  % inertial East position
    h        = uu(3+NN);  % altitude
    Va       = uu(4+NN);  % airspeed
%    alpha    = uu(5+NN);  % angle of attack
%    beta     = uu(6+NN);  % side slip angle
    phi      = uu(7+NN);  % roll angle
    theta    = uu(8+NN);  % pitch angle
    chi      = uu(9+NN);  % course angle
    p        = uu(10+NN); % body frame roll rate
    q        = uu(11+NN); % body frame pitch rate
    r        = uu(12+NN); % body frame yaw rate
%    Vg       = uu(13+NN); % ground speed
%    wn       = uu(14+NN); % wind North
%    we       = uu(15+NN); % wind East
%    psi      = uu(16+NN); % heading
%    bx       = uu(17+NN); % x-gyro bias
%    by       = uu(18+NN); % y-gyro bias
%    bz       = uu(19+NN); % z-gyro bias
    NN = NN+19;
    Va_c     = uu(1+NN);  % commanded airspeed (m/s)
    h_c      = uu(2+NN);  % commanded altitude (m)
    chi_c    = uu(3+NN);  % commanded course (rad)
    NN = NN+3;
    t        = uu(1+NN);   % time

    autopilot_version = 2;
        % autopilot_version == 1 <- used for tuning
        % autopilot_version == 2 <- standard autopilot defined in book
        % autopilot_version == 3 <- Total Energy Control for longitudinal AP
    switch autopilot_version
        case 1,
           [delta, x_command] = autopilot_tuning(Va_c,h_c,chi_c,Va,h,chi,phi,theta,p,q,r,t,P);
        case 2,
           [delta, x_command] = autopilot_uavbook(Va_c,h_c,chi_c,Va,h,chi,phi,theta,p,q,r,t,P);
        case 3,
               [delta, x_command] = autopilot_TECS(Va_c,h_c,chi_c,Va,h,chi,phi,theta,p,q,r,t,P);
    end
    y = [delta; x_command];
end

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% autopilot_tuning
%   - used to tune each loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [delta, x_command] = autopilot_tuning(Va_c,h_c,chi_c,Va,h,chi,phi,theta,p,q,r,t,P)

    mode = 5;
    switch mode
        case 1, % tune the roll loop
            phi_c = chi_c; % interpret chi_c to autopilot as course command
            delta_a = roll_hold(phi_c, phi, p, P);
            delta_r = 0; % no rudder
            % use trim values for elevator and throttle while tuning the lateral autopilot
            delta_e = P.u_trim(1);
            delta_t = P.u_trim(4);
            theta_c = 0;
        case 2, % tune the course loop
            if t==0,
                phi_c   = course_hold(chi_c, chi, r, 1, P);
            else
                phi_c   = course_hold(chi_c, chi, r, 0, P);
            end
            delta_a = roll_hold(phi_c, phi, p, P);
            delta_r = 0; % no rudder
            % use trim values for elevator and throttle while tuning the lateral autopilot
            delta_e = P.u_trim(1);
            delta_t = P.u_trim(4);
            theta_c = 0;
        case 3, % tune the throttle to airspeed loop and pitch loop simultaneously
            theta_c = 20*pi/180 + h_c;
            chi_c = 0;
            if t==0,
                phi_c   = course_hold(chi_c, chi, r, 1, P);
                delta_t = airspeed_with_throttle_hold(Va_c, Va, 1, P);
           else
                phi_c   = course_hold(chi_c, chi, r, 0, P);
                delta_t = airspeed_with_throttle_hold(Va_c, Va, 0, P);
            end
            delta_e = pitch_hold(theta_c, theta, q, P);
            delta_a = roll_hold(phi_c, phi, p, P);
            delta_r = 0; % no rudder
            % use trim values for elevator and throttle while tuning the lateral autopilot
        case 4, % tune the pitch to airspeed loop
            chi_c = 0;
            delta_t = P.u_trim(4);
            if t==0,
                phi_c   = course_hold(chi_c, chi, r, 1, P);
                theta_c = airspeed_with_pitch_hold(Va_c, Va, 1, P);
           else
                phi_c   = course_hold(chi_c, chi, r, 0, P);
                theta_c = airspeed_with_pitch_hold(Va_c, Va, 0, P);
            end
            delta_a = roll_hold(phi_c, phi, p, P);
            delta_e = pitch_hold(theta_c, theta, q, P);
            delta_r = 0; % no rudder
            % use trim values for elevator and throttle while tuning the lateral autopilot
        case 5, % tune the pitch to altitude loop
            chi_c = 0;
            if t==0,
                phi_c   = course_hold(chi_c, chi, r, 1, P);
                theta_c = altitude_hold(h_c, h, 1, P);
                delta_t = airspeed_with_throttle_hold(Va_c, Va, 1, P);
           else
                phi_c   = course_hold(chi_c, chi, r, 0, P);
                theta_c = altitude_hold(h_c, h, 0, P);
                delta_t = airspeed_with_throttle_hold(Va_c, Va, 0, P);
            end
            delta_a = roll_hold(phi_c, phi, p, P);
            delta_e = pitch_hold(theta_c, theta, q, P);
            delta_r = 0; % no rudder
            % use trim values for elevator and throttle while tuning the lateral autopilot
         case 6, % tune the pitch loop
            phi_c = 0;
            theta_c = 20*pi/180 + h_c;
            delta_a = P.u_trim(2);
            delta_r = 0; % no rudder
            % use trim values for elevator and throttle while tuning the lateral autopilot
            delta_e = pitch_hold(theta_c, theta, q, P);
            delta_t = P.u_trim(4);
    end
    %----------------------------------------------------------
    % create outputs

    % control outputs
    delta = [delta_e; delta_a; delta_r; delta_t];
    % commanded (desired) states
    x_command = [...
        0;...                    % pn
        0;...                    % pe
        h_c;...                  % h
        Va_c;...                 % Va
        0;...                    % alpha
        0;...                    % beta
        phi_c;...                % phi
        %theta_c*P.K_theta_DC;... % theta
        theta_c;
        chi_c;...                % chi
        0;...                    % p
        0;...                    % q
        0;...                    % r
        ];

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% autopilot_uavbook
%   - autopilot defined in the uavbook
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [delta, x_command] = autopilot_uavbook(Va_c,h_c,chi_c,Va,h,chi,phi,theta,p,q,r,t,P)

    %----------------------------------------------------------
    % lateral autopilot
    if t==0,
        % assume no rudder, therefore set delta_r=0
        delta_r = 0;%coordinated_turn_hold(beta, 1, P);
        phi_c   = course_hold(chi_c, chi, r, 1, P);

    else
        phi_c   = course_hold(chi_c, chi, r, 0, P);
        delta_r = 0;%coordinated_turn_hold(beta, 0, P);
    end
    delta_a = roll_hold(phi_c, phi, p, P);

    %----------------------------------------------------------
    % longitudinal autopilot

    % define persistent variable for state of altitude state machine
    persistent altitude_state;
    persistent initialize_integrator;
    % initialize persistent variable
    if h<=P.altitude_take_off_zone,
        altitude_state = 1;
    elseif h<=h_c-P.altitude_hold_zone,
        altitude_state = 2;
    elseif h>=h_c+P.altitude_hold_zone,
        altitude_state = 3;
    else
        altitude_state = 4;
    end
    initialize_integrator = 1;
%     display(altitude_state);
    % implement state machine
    switch altitude_state,
        case 1,  % in take-off zone
            delta_t = 0.5;
            theta_c = P.theta_max;
        case 2,  % climb zone
            delta_t = 0.5;
            theta_c = airspeed_with_pitch_hold(Va_c, Va, 0, P);
        case 3, % descend zone
            delta_t = 0;
            theta_c = airspeed_with_pitch_hold(Va_c, Va, 0, P);
        case 4, % altitude hold zone
            delta_t = airspeed_with_throttle_hold(Va_c, Va, flag, P);
            theta_c = altitude_hold(h_c, h, flag, P);
    end

    delta_e = pitch_hold(theta_c, theta, q, P);
    % artificially saturation delta_t
    delta_t = sat(delta_t,1,0);

    %----------------------------------------------------------
    % create outputs

    % control outputs
    delta = [delta_e; delta_a; delta_r; delta_t];
    % commanded (desired) states
    x_command = [...
        0;...                    % pn
        0;...                    % pe
        h_c;...                  % h
        Va_c;...                 % Va
        0;...                    % alpha
        0;...                    % beta
        phi_c;...                % phi
        %theta_c*P.K_theta_DC;... % theta
        theta_c;
        chi_c;...                % chi
        0;...                    % p
        0;...                    % q
        0;...                    % r
        ];

    y = [delta; x_command];

    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% autopilot_TECS
%   - longitudinal autopilot based on total energy control systems
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [delta, x_command] = autopilot_TECS(Va_c,h_c,chi_c,Va,h,chi,phi,theta,p,q,r,t,P)

    %----------------------------------------------------------
    % lateral autopilot
    if t==0,
        % assume no rudder, therefore set delta_r=0
        delta_r = 0;%coordinated_turn_hold(beta, 1, P);
        phi_c   = course_hold(chi_c, chi, r, 1, P);

    else
        phi_c   = course_hold(chi_c, chi, r, 0, P);
        delta_r = 0;%coordinated_turn_hold(beta, 0, P);
    end
    delta_a = roll_hold(phi_c, phi, p, P);

    %----------------------------------------------------------
    % longitudinal autopilot based on total energy control

    delta_e = 0;
    delta_t = 0;

    %----------------------------------------------------------
    % create outputs

    % control outputs
    delta = [delta_e; delta_a; delta_r; delta_t];
    % commanded (desired) states
    x_command = [...
        0;...                    % pn
        0;...                    % pe
        h_c;...                  % h
        Va_c;...                 % Va
        0;...                    % alpha
        0;...                    % beta
        phi_c;...                % phi
        %theta_c*P.K_theta_DC;... % theta
        theta_c;
        chi_c;...                % chi
        0;...                    % p
        0;...                    % q
        0;...                    % r
        ];

    y = [delta; x_command];

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Autopilot functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function delta_a = roll_hold(phi_c, phi, p, P)
    persistent integrator;
    persistent error_d1;
    if isempty(integrator) % reset (initialize) persistent variables when flag==1
        integrator = 0;
        error_d1 = 0;   % _d1 means delayed by one time step
    end
    error = phi_c - phi;    % compute the current error
    integrator = integrator + (P.Ts / 2) * (error + error_d1); % update integrator
    error_d1 = error; % update the error for next time through the loop
    delta_a = sat(...                   % implement PID control
        P.kp_phi * error + ...         % proportional term
        P.ki_phi * integrator - ...    % integral term
        P.kd_phi * p,...  % derivative term
        P.delta_a_max, -P.delta_a_max ... % ensure abs(u)<=limit
        );
    if P.ki_phi~=0        % implement integrator anti-windup
        u_unsat = P.kp_phi * error + P.ki_phi * integrator - P.kd_phi * p;
        integrator = integrator + P.Ts / P.ki_phi * (delta_a - u_unsat);
    end
end

function phi_c = course_hold(chi_c, chi, r, flag, P)
    persistent integrator;
    persistent error_d1;
    if isempty(integrator)        % reset (initialize) persistent variables when flag==1
        integrator = 0;
        error_d1 = 0;   % _d1 means delayed by one time step
    end
    error = chi_c - chi;    % compute the current error
    integrator = integrator + (P.Ts / 2) * (error + error_d1); % update integrator
    error_d1 = error; % update the error for next time through the loop
    phi_c = sat(...                   % implement PID control
        P.kp_chi * error + ...         % proportional term
        P.ki_chi * integrator, ...    % integral term
        P.roll_max, -P.roll_max ... % ensure abs(u)<=limit
        );
    if P.ki_chi~=0        % implement integrator anti-windup
        u_unsat = P.kp_chi * error + P.ki_chi * integrator;
        integrator = integrator + P.Ts / P.ki_chi * (phi_c - u_unsat);
    end
end

function delta_e = pitch_hold(theta_c, theta, q, P)
    persistent error_d1;
    if isempty(error_d1)       % reset (initialize) persistent variables when flag==1
        error_d1 = 0;   % _d1 means delayed by one time step
    end
    error = theta_c - theta;    % compute the current error\
    error_d1 = error; % update the error for next time through the loop
    delta_e = sat(...                   % implement PID control
        P.kp_theta * error - ...         % proportional term
        P.kd_theta * q,...  % derivative term
        P.delta_e_max, -P.delta_e_max ... % ensure abs(u)<=limit
        );
end

function delta_t = airspeed_with_throttle_hold(Va_c, Va, flag, P)
    persistent integrator;
    persistent error_d1;
    if isempty(integrator)        % reset (initialize) persistent variables when flag==1
        integrator = 0;
        error_d1 = 0;   % _d1 means delayed by one time step
    end
    error = Va_c - Va;    % compute the current error
    integrator = integrator + (P.Ts / 2) * (error + error_d1); % update integrator
    error_d1 = error; % update the error for next time through the loop
    delta_t = sat(...                   % implement PID control
        P.kp_v * error + ...         % proportional term
        P.ki_v * integrator, ...    % integral term
        P.delta_t_max, P.delta_t_min ... % ensure abs(u)<=limit
        );
    if P.ki_v~=0        % implement integrator anti-windup
        u_unsat = P.kp_v * error + P.ki_v * integrator;
        integrator = integrator + P.Ts / P.ki_v * (delta_t - u_unsat);
    end
end

function theta_c = airspeed_with_pitch_hold(Va_c, Va, flag, P)
    persistent integrator;
    persistent error_d1;
    if isempty(integrator)        % reset (initialize) persistent variables when flag==1
        integrator = 0;
        error_d1 = 0;   % _d1 means delayed by one time step
    end
    error = Va_c - Va;    % compute the current error
    integrator = integrator + (P.Ts / 2) * (error + error_d1); % update integrator
    error_d1 = error; % update the error for next time through the loop
    theta_c = sat(...                   % implement PID control
        P.kp_v2 * error + ...         % proportional term
        P.ki_v2 * integrator, ...    % integral term
        P.pitch_max, -P.pitch_max ... % ensure abs(u)<=limit
        );
    if P.ki_v2~=0        % implement integrator anti-windup
        u_unsat = P.kp_v2 * error + P.ki_v2 * integrator;
        integrator = integrator + P.Ts / P.ki_v2 * (theta_c - u_unsat);
    end
end

function theta_c = altitude_hold(h_c, h, flag, P)
    persistent integrator;
    persistent error_d1;
    if isempty(integrator)        % reset (initialize) persistent variables when flag==1
        integrator = 0;
        error_d1 = 0;   % _d1 means delayed by one time step
    end
    error = h_c - h;    % compute the current error
    integrator = integrator + (P.Ts / 2) * (error + error_d1); % update integrator
    error_d1 = error; % update the error for next time through the loop
    theta_c = sat(...                   % implement PID control
        P.kp_h * error + ...         % proportional term
        P.ki_h * integrator, ...    % integral term
        P.pitch_max, -P.pitch_max ...    % ensure abs(u)<=limit
        );
    if P.ki_h~=0        % implement integrator anti-windup
        u_unsat = P.kp_h * error + P.ki_h * integrator;
        integrator = integrator + P.Ts / P.ki_h * (theta_c - u_unsat);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sat
%   - saturation function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = sat(in, up_limit, low_limit)
  if in > up_limit,
      out = up_limit;
  elseif in < low_limit;
      out = low_limit;
  else
      out = in;
  end
end


%  Script file: timings.m
%
%  Purpose:
%    This program calculates the time required to
%    calculate the squares of all integers from 1 to
%    10,000 in four different ways:
%    1.  Using a for loop with an uninitialized output
%        array.
%    2.  Using a for loop with a pre-allocated output
%        array and NO JIT compiler.
%    3.  Using a for loop with a pre-allocated output
%        array and the JIT compiler.
%    4.  Using vectors.
%
%  Record of revisions:
%      Date       Programmer          Description of change
%      ====       ==========          =====================
%    01/29/07    S. J. Chapman        Original code
%
% Define variables:
%   ii, jj       -- Loop index
%   average1     -- Average time for calculation 1
%   average2     -- Average time for calculation 2
%   average3     -- Average time for calculation 3
%   average4     -- Average time for calculation 4
%   maxcount     -- Number of times to loop calculation
%   square       -- Array of squares

% Perform calculation with an uninitialized array
% "square".  This calculation is done only once
% because it is so slow.
maxcount = 1;               % Number of repetitions
tic;                        % Start timer
for jj = 1:maxcount
   clear square             % Clear output array
   for ii = 1:10000
     square(ii) = ii^2;     % Calculate square
   end
end
average1 = (toc)/maxcount;  % Calculate average time

% Perform calculation with a pre-allocated array
% "square", calling an external function to square
% the number.  This calculation is averaged over 10
% loops.
maxcount = 10;              % Number of repetitions
tic;                        % Start timer
for jj = 1:maxcount
   clear square             % Clear output array
   square = zeros(1,10000); % Pre-initialize array
   for ii = 1:10000
     square(ii) = sqr(ii);  % Calculate square
   end
end
average2 = (toc)/maxcount;  % Calculate average time

% Perform calculation with a pre-allocated array
% "square".  This calculation is averaged over 100
% loops.
maxcount = 100;             % Number of repetitions
tic;                        % Start timer
for jj = 1:maxcount
   clear square             % Clear output array
   square = zeros(1,10000); % Pre-initialize array
   for ii = 1:10000
     square(ii) = ii^2;     % Calculate square
   end
end
average3 = (toc)/maxcount;  % Calculate average time

% Perform calculation with vectors.  This calculation
% averaged over 1000 executions.
maxcount = 1000;            % Number of repetitions
tic;                        % Start timer
for jj = 1:maxcount
   clear square             % Clear output array
   ii = 1:10000;            % Set up vector
   square = ii.^2;          % Calculate square
end
average4 = (toc)/maxcount;  % Calculate average time

% Display results
fprintf('Loop / uninitialized array        = %8.4f\n', average1);
fprintf('Loop / initialized array / no JIT = %8.4f\n', average2);
fprintf('Loop / initialized array / JIT    = %8.4f\n', average3);
fprintf('Vectorized                        = %8.4f\n', average4);


function [ber, rate] = MO_method()

global  H  Vn W_mopt Nt  Nrf Nr;
i = 0;

W_equal = W_mopt;
w = trace (W_equal' * W_equal);
%random initialization
V_RF = exp( 1i*unifrnd(0,2*pi,Nt,Nrf));
W_RF = exp( 1i*unifrnd(0,2*pi,Nr,Nrf));
%iteration trigger, the normal initialization just for pass into functions
trigger = 1;
m_MSE_new = 100;

%limit the iterations number by i<10
while (trigger > 1e-5)

    % precoding
    H1 = H' * W_equal;
    Vn1 = Vn * w;
    [V_RF, V_U] = mo_algorithm(V_RF, Vn1, H1);
    V_equal = V_RF *V_U;
    v = trace (V_equal * V_equal');   %beta^(-2)

    %combining
    H2 = H * V_equal;
    Vn2 = Vn * v;
    [W_RF, W_B] = mo_algorithm(W_RF, Vn2, H2);
    W_equal = W_RF * W_B;
    w = trace (W_equal' * W_equal);
    %modified MSE
    H_equal = W_equal'*H2;

    m_MSE_old = m_MSE_new;
    m_MSE_new = trace(H_equal * H_equal' - H_equal - H_equal') + Vn * v * w;
    trigger = m_MSE_old - m_MSE_new;

    i = i + 1;
end

V_B = V_U / sqrt(v);

V = V_RF * V_B;
W = W_RF * W_B;

ber = get_ber(V, W);
rate = get_rate(V, W);


text:
function pp=lpcra2pp(ra)
%LPCAR2PP LPC: Convert ar filter autocorrelation to power spectrum polynomial in cos(w) PP=(RA)
% pp is a polynomial such that |polyval(ra,e^jw)| = polyval(pp,cos(w))

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

persistent tp;
[nf,p1]=size(ra);
% we check here if p is the same as before and reuse the tp matrix
if size(tp,1)~=p1
   p=p1-1;
   % chebyshev polynomials up to order p
   tp=zeros(p1,p1);
   tp(1,p1)=2;
   tp(2,p)=2;
   for i=3:p1
      tp(i,p+2-i:p)=2*tp(i-1,p+3-i:p1)-tp(i-2,p+2-i:p);
      tp(i,p1)=-tp(i-2,p1);
   end
   tp(1,p1)=1;
end
pp=ra*tp;

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function [tru] = computeTrueAnomFromHypAnom(HypA, ecc)
%computeTrueAnomFromHypAnom Summary of this function goes here
%   Detailed explanation goes here

    upper = sqrt(ecc+1) .* tanh(HypA/2);
    lower = sqrt(ecc-1);
    tru = atan2(upper, lower) * 2;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [declusCat,mNumDeclus] = MonteReasenberg(numSim,Cat)
    % Monte Carlo Simulation for Reasenberg-declustering parameters
    %
    % Wrapper function to do Monte Carlo simulations of the
    % input parameters into the Reasenberg Declustering algorithm
    %%

    % if ~exist('FileOut','var')
    %     disp(['You must provide the output results filename!'])
    %     return
    % end

    resFileOut = 'DeclusRes';
    parmFileOut = 'DeclusParms';
    mNumDeclus=[];

    %%
    % set Default ranges for Reasenberg input variables and find their range
    raTaumin = [.5 2.5];
    raTaumax = [3 20];
    raP = [.9 .999];
    raXk = [0.4 0.6];
    raXmeff = [2.7 2.9];
    raRfact = [5 20];
    raErr = [2 4];
    raDerr = [4 6];

    tauminDiff = (raTaumin(2) - raTaumin(1));
    taumaxDiff = (raTaumax(2) - raTaumax(1));
    pDiff = (raP(2) - raP(1));
    xkDiff = (raXk(2) - raXk(1));
    xmeffDiff = (raXmeff(2) - raXmeff(1));
    rfactDiff = (raRfact(2) - raRfact(1));

    %% add column for independence probability
    % actually will just be number of times the event has appeared in a catalogue
    % (will need to divide by simNum to get P)

    Cat(:,10) = 0;

    % set the rand number generator state
    rng('shuffle');

    % simulate parameter values and run the delcustering code

    lErr = raErr(1);
    lDerr = raDerr(1);
    rdc= ReasenbergDeclusterClass(Cat, 'err', lErr, 'derr', lDerr, ....
        'AutoShowPlots',false,'DelayProcessing',true,'InteractiveMode',false);

    for simNum = 1:numSim

        randNum = rand(1,8);
        rdc.taumin = raTaumin(1) + tauminDiff*randNum(1);
        rdc.taumax = raTaumax(1) + taumaxDiff*randNum(2);
        rdc.P = raP(1) + pDiff*randNum(3);
        rdc.xk = raXk(1) + xkDiff*randNum(4);
        rdc.xmeff = raXmeff(1) + xmeffDiff*randNum(5);
        rdc.rfact = raRfact(1) + rfactDiff*randNum(6);

        [declusCat, is_mainshock] = rdc.ReasenbergDeclus();

        % [declusCat,is_mainshock] = ReasenbergDeclus(lTaumin,lTaumax,lXk,lXmeff,lP,lRfact,lErr,lDerr,Cat);

        Cat(is_mainshock,10) = Cat(is_mainshock,10) + 1;
        nIst=zeros(length(Cat),1);
        nIst(is_mainshock)=1;
        nIst=logical(nIst);
        mNumDeclus=[mNumDeclus,(nIst==1)];
        save(resFileOut,'Cat');

        monteParms(simNum) = {[rdc.taumin;rdc.taumax;rdc.P;rdc.xk;rdc.xmeff;rdc.rfact;rdc.err;rdc.derr]};
        save(parmFileOut,'monteParms');
        disp(num2str(simNum));

    end
end

author: CelsoReyes


text:
function slope = calc_slope_grigsby(structNum,Thresholds)
%function slope = calc_slope_grigsby(structNum,Thresholds)
%
%This function returns shape features.
%
%IEN & APA, 07/06/2006

global planC stateS
indexS = planC{end};
for i=1:length(Thresholds)
    contourSUV(structNum,Thresholds(i));
    structVol(i) = getStructureVol(length(planC{indexS.structures}));
    %runCERRCommand('del','structure',length(planC{indexS.structures}))
    n = length(planC{indexS.structures});
    len = n;
    planC = delUniformStr(n, planC); %Update the uniform data.
    planC{indexS.structures}(n:len-1) = planC{indexS.structures}(n+1:len);
    planC{indexS.structures} = planC{indexS.structures}(1:len-1);
end
coeff = polyfit(Thresholds,structVol,1);
slope=coeff(1);
return;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function J=calcPolarConvJacob(zPolar,systemType,useHalfRange,lTx,lRx,M)
%%CALCPOLARCONVJACOB Calculate the Jacobian for a monostatic or bistatic
%            range and polar angle measurement in 2D with respect to
%            Cartesian position. Atmospheric effects are ignored. This
%            type of Jacobian is useful when performing tracking using
%            Cartesian-converted measurements where the clutter density is
%            specified in the measurement coordinate system, not the
%            converted measurement coordinate system.
%
%INPUTS: zPolar A 2X1 point in polar coordinates in the format
%          [range;azimuth], where the angle is given in radians and the
%          range can be bistatic.
% systemType An optional parameter specifying the axis from which the
%          azimuth angle is measured. It is assumed that the azimuth
%          angle is given in radians. Possible values are
%          0 (The default if omitted) The azimuth angle is
%             counterclockwise from the x axis.
%          1 The azimuth angle is measured clockwise from the y axis.
% useHalfRange A boolean value specifying whether the bistatic range value
%          should be divided by two. This normally comes up
%          when operating in monostatic mode, so that the range reported is
%          a one-way range. The default if this parameter is not provided
%          is false.
%      lTx The 2X1 [x;y] location vector of the transmitter in Cartesian
%          coordinates. If this parameter is omitted or an empty matrix is
%          passed, then the transmitter is assumed to be at the origin.
%      lRx The 2X1 [x;y] location vector of the receiver in Cartesian
%          coordinates. If this parameter is omitted or an empty matrix is
%          passed, then the receiver is assumed to be at the origin.
%        M A 2X2 rotation matrices to go from the alignment of the global
%          coordinate system to that at the receiver. If omitted or an
%          empty matrix is passed, then it is assumed that the local
%          coordinate system is aligned with the global and M=eye(2,2)
%          --the identity matrix is used.
%
%OUTPUTS: J The 2X2 Jacobian matrix. Each row is a components of range, and
%           azimuth (in that order by row) with derivatives taken with
%           respect to [x,y] by column.
%
%This function converts the measurement into Cartesian coordinates and then
%calls rangeGradient and polAngGradient.
%
%February 2017 David F.Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<6||isempty(M))
   M=eye(2,2);
end

if(nargin<5||isempty(lRx))
   lRx=zeros(2,1);
end

if(nargin<4||isempty(lTx))
   lTx=zeros(2,1);
end

if(nargin<3||isempty(useHalfRange))
    useHalfRange=true;
end

if(nargin<2||isempty(systemType))
    systemType=0;
end

x=pol2Cart(zPolar,systemType,useHalfRange,lTx,lRx,M);

J=zeros(2,2);
J(1,:)=rangeGradient(x,useHalfRange,lTx,lRx);
J(2,:)=polAngGradient(x,systemType,lRx);

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [t,options] = loadTensor_json(fname,varargin)

[~,~,ext] = fileparts(fname);
if ~strcmpi(ext,'.json'), interfaceError(fname), end

s = loadjson(fname);

type = char(fieldnames(s));

s = s.(type);

cs = crystalSymmetry(s.CS.pointGroup,s.CS.abc,s.CS.abg,s.CS.alignment{:},'mineral',s.CS.mineral);

t = feval(type,s.M,cs,'rank',s.rank);
try t.opt = s.opt; end

options = {};

author: mtex-toolbox

repo: mtex

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


text:
function jacobi_eigenvalue_test ( )

%*****************************************************************************80
%
%% JACOBI_EIGENVALUE_TEST tests the JACOBI_EIGENVALUE library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 July 2013
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'JACOBI_EIGENVALUE_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the JACOBI_EIGENVALUE library.\n' );

  jacobi_eigenvalue_test01 ( );
  jacobi_eigenvalue_test02 ( );
  jacobi_eigenvalue_test03 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'JACOBI_EIGENVALUE_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( )

  return
end


text:
% Init Script
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     08/11/2016  Royi Avital
%   *   First release.
%

%% General Parameters
close('all');
clear();
clc();

% set(0, 'DefaultFigureWindowStyle', 'docked');
% defaultLoosInset = get(0, 'DefaultAxesLooseInset');
% set(0, 'DefaultAxesLooseInset', [0.05, 0.05, 0.05, 0.05]);

figPosSmall     = [100, 100, 400, 300];
figPosMedium    = [100, 100, 800, 600];
figPosLarge     = [100, 100, 960, 720];
figPosDefault   = [100, 100, 560, 420];

fontSizeTitle   = 14;
fontSizeAxis    = 12;
fontSizeString  = 12;

lineWidthThin   = 1;
lineWidthNormal = 3;
lineWidthThick  = 4;

markerSizeSmall     = 4;
markerSizeNormal    = 8;
markerSizeLarge     = 10;

% https://www.mathworks.com/help/matlab/graphics_transition/why-are-plot-lines-different-colors.html
% https://www.mathworks.com/matlabcentral/answers/160332
mColorOrder = get(groot, 'DefaultAxesColorOrder');

randomNumberStream  = RandStream('mlfg6331_64', 'NormalTransform', 'Ziggurat');
subStreamNumber     = round(sum(clock()));
% subStreamNumber     = 162;
subStreamNumber     = 20;
set(randomNumberStream, 'Substream', subStreamNumber);
RandStream.setGlobalStream(randomNumberStream);

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

author: RoyiAvital

repo: StackExchangeCodes


text:
function [t,t1] = geyer_icse(x,maxlag)
% GEYER_ICSE - Compute autocorrelation time tau using Geyer's
%              initial convex sequence estimator
%
%   C = GEYER_ICSE(X) returns autocorrelation time tau.
%   C = GEYER_ICSE(X,MAXLAG) returns autocorrelation time tau with
%       MAXLAG . Default MAXLAG = M-1.
%
%   References:
%      [1] C. J. Geyer, (1992). "Practical Markov Chain Monte Carlo",
%          Statistical Science, 7(4):473-511
%
% Copyright (C) 2002 Aki Vehtari

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

% compute autocorrelation
if nargin > 1
  cc=acorr(x,maxlag);
else
  cc=acorr(x);
end
[n,m]=size(cc);

% acorr returns values starting from lag 1, so add lag 0 here
cc=[ones(1,m);cc];
n=n+1;

% now make n even
if mod(n,2)
  n=n-1;
  cc(end,:)=[];
end

% loop through variables
t=zeros(1,m);
t1=zeros(1,m);
opt=optimset('LargeScale','off','display','off');
for i1=1:m
  c=cc(:,i1);
  c=sum(reshape(c,2,n/2),1);
  ci=find(c<0);
  if isempty(ci)
    warning(sprintf('Inital positive could not be found for variable %d, using maxlag value',i1));
    ci=n/2;
  else
    ci=ci(1)-1; % initial positive
  end
  c=[c(1:ci) 0];    % initial positive sequence
  t1(i1)=-1+2*sum(c); % initial positive sequence estimator
  if ci>2
    ca=fmincon(@se,c,[],[],[],[],0*c,c,@sc,opt,c); % monotone convex sequence
  else
    ca=c;
  end
  t(i1)=-1+2*sum(ca); % monotone convex sequence estimator
end

function e = se(x,xx)
% SE - Error in monotone convex sequene estimator
e=sum((xx-x).^2);

function [c,ceq] = sc(x,xx)
% SE - Constraint in monotone convex sequene estimator
ceq=0*x;
c=ceq;
d=diff(x);
dd=-diff(d);
d(d<0)=0;d=d.^2;
dd(dd<0)=0;dd=dd.^2;
c(1:end-1)=d;c(2:end)=c(2:end)+d;
c(1:end-2)=dd;c(2:end-1)=c(2:end-1)+dd;c(3:end)=c(3:end)+dd;

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
clear all
close all
path(path,'..\..\..\FUZZCLUST')
%the data
load motorcycle.txt
data.X = motorcycle(:,[1 2]);

%parameters
param.c=4;
param.m=2;
param.e=1e-6;
param.ro=ones(1,param.c);
param.val=1;
%normalization
data=clust_normalize(data,'range');
%clustering
result = FCMclust(data,param);
plot(data.X(:,1),data.X(:,2),'b.',result.cluster.v(:,1),result.cluster.v(:,2),'ro');
hold on
%draw contour-map
new.X=data.X;
eval=clusteval(new,result,param);
%validation
result = validity(result,data,param);
result.validity


text:
function [conn,connnum,count]=neighborelem(elem,nn);
%
% [conn,connnum,count]=neighborelem(elem,nn)
%
% create node neighbor list from a mesh
%
% author: Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
% date: 2007/11/21
%
% input:
%    elem:  element table of a mesh
%    nn  :  total node number of the mesh
%
% output:
%    conn:  output, a cell structure of length nn, conn{n}
%           contains a list of all neighboring elem ID for node n
%    connnum: vector of length nn, denotes the neighbor number of each node
%    count: total neighbor numbers
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

conn=cell(nn,1);
dim=size(elem);
for i=1:dim(1)
  for j=1:dim(2)
    conn{elem(i,j)}=[conn{elem(i,j)},i];
  end
end
count=0;
connnum=zeros(1,nn);
for i=1:nn
    conn{i}=sort(conn{i});
    connnum(i)=length(conn{i});
    count=count+connnum(i);
end

repo: REST


text:
% pdfbtrainimagethmt.m
% written by: Duncan Po
% Date: August 24, 2002
% Using the EM algorithm, train models for the specified image
% Usage: [model, stateprob] = pdfbtrainimagethmt(imname, imformat, initmodel, mD)
% Inputs:   imname      - name of the image file
%           imformat    - format of the image file (e.g. 'gif')
%           initmodel   - optional. Give an intial model to speed up the training
%                         process. Input '' if not providing initial model.
%           mD          - convergence value
% Output:   model       - the model generated
%           stateprob   - state probabilities

function [model, stateprob] = pdfbtrainimagethmt(imname, imformat, initmodel, mD)

pyrfilter = '9-7';
dirfilter = 'pkva';
levndir = [2 2 3 3];
ns = 2;
zeromean = 'yes';

coef = contourlet(pyrfilter, dirfilter, levndir, imname, imformat);

for dir = 1:2.^levndir(1)
    [tree, scaling] = contourlet2tree(coef, dir);
    if isempty(initmodel)
        [tempmodel, tempstateprob] = pdfbtrainthmt(tree, levndir, mD, ns,zeromean);
    else
        [tempmodel, tempstateprob] = pdfbtrainthmt(tree, levndir, mD,...
            initmodel{dir});
    end;
    model{dir} = tempmodel;
    stateprob{dir} = tempstateprob;
end;


text:
function [fx,dF_dX,dF_dTheta] = f_Henon(Xt,Theta,ut,inF)
% Henon chaotic map evolution function

x       = Xt;
a       = Theta(1);
b       = Theta(2);

fx      = zeros(2,1);
fx(1)   = x(2) + 1 -a.*x(1).^2;
fx(2)   = b.*x(1);

J       = zeros(2,2);
J(1,:)  = [-2.*a*x(1),1];
J(2,:)  = [b,0];
dF_dX   = J';

dF_dTheta = zeros(2,2);
dF_dTheta(1,:)      = [-x(1).^2,0];
dF_dTheta(2,:)      = [0,x(1)];
dF_dTheta = dF_dTheta';

author: MBB-team

repo: VBA-toolbox


text:
function writeEnergy(N,T,U)

%This function writes the energy file for the n-link pendulum.

nStateStr = num2str(2*N);

filename = ['energy_' num2str(N) '_link'];

comments{1} = ['[ENERGY, KINETIC, POTENTIAL] = ' upper(filename) '(Z,P)'];
comments{2} = ' ';
comments{3} = 'FUNCTION:  This function computes the energy of a double';
comments{4} = '    pendulum.';
comments{5}  = 'INPUTS: ';
comments{10} = ['    z = [' nStateStr ' X nTime]  matrix of states'];
comments{11} = '    P = struct of parameters';
comments{8}  = 'OUTPUTS: ';
comments{9}  = '    energy = [1 X nTime] vector of total energy';
comments{10} = '    kinetic = [1 X nTime] vector of kinetic energy';
comments{11} = '    potential = [1 X nTime] vector of potential energy';
comments{12} = ' ';
comments{13} = 'NOTES:';
comments{14} = ['    This file was automatically generated by ' mfilename '.m'];

params{1} = {'g ','gravity'};
params{2} = {'m','mass'};
params{3} = {'l','length'};
params{4} = {'I','moment of inertia about its center of mass'};
params{5} = {'d','distance between center of mass and parent joint'};

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                               write file                                %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
fid = fopen([filename '.m'],'w');

fprintf(fid, ['function [energy, kinetic, potential] = ' filename '(z,P) \n']);

for i=1:length(comments)
    fprintf(fid,['%%' comments{i} '\n']);
end
fprintf(fid,'\n');

for i=1:length(params)
    if i==1  %Gravity, scalar
    fprintf(fid,[params{i}{1} ' = P.' params{i}{1} '; %%' params{i}{2} '\n']);
    else  %unpack vectors:
        for j=1:N
           fprintf(fid,[params{i}{1} num2str(j) ...
               ' = P.' params{i}{1} '(' num2str(j) '); %% Link ' num2str(j) ' ' params{i}{2} '\n']);
        end
    end
end
fprintf(fid,'\n');

for j = 1:N
   fprintf(fid,['th' num2str(j) ' = z(' num2str(j) ',:); \n']);
end
fprintf(fid,'\n');

for j = 1:N
   fprintf(fid,['dth' num2str(j) ' = z(' num2str(j+N) ',:); \n']);
end
fprintf(fid,'\n');

fprintf(fid,['kinetic = ' vectorize(char(T)) ';\n']);
fprintf(fid,['potential = ' vectorize(char(U)) ';\n']);
fprintf(fid,'energy = potential + kinetic;\n\n');
fprintf(fid,'end \n');

fclose(fid);

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [IN,JN,XN,II,JI,XI] = globMatrixIFE3D(fun,mesh,femI,fem1,d1,fem2,d2)

%% USAGE: generate global matrix on a 3D mesh
%
% INPUTS:
% fun --- coefficient function
% mesh --- a struct data contains very rich mesh information.
% fem1 --- global DoF for test function space
% fem2 --- global DoF for trial function space
% d1 --- derivative info for test function
% d2 --- derivative info for trial function
%            d = [0,0,0]: function value
%            d = [1,0,0]: Dx value
%            d = [0,1,0]: Dy value
%            d = [0,0,1]: Dz value
%
% OUTPUTS:
% [IN JN XN] --- triplets of the sparse matrix from regular elements.
% [II JI XI] --- triplets of the sparse matrix from interface elements.

% Last Modified: 08/07/2020 by Xu Zhang

%% 0. Initializaiton
if nargin == 5
    fem2 = fem1; d2 = d1;
end

if strcmp(fem1.type,'P1')||strcmp(fem1.type,'DGP1')||strcmp(fem1.type,'CR')
    feEvalBas1 = @evalP1Bas3D;
elseif strcmp(fem1.type,'P2')||strcmp(fem1.type,'DGP2')
    feEvalBas1 = @evalP2Bas3D;
end

if strcmp(fem2.type,'P1')||strcmp(fem2.type,'DGP1')||strcmp(fem2.type,'CR')
    feEvalBas2 = @evalP1Bas3D;
elseif strcmp(fem2.type,'P2')||strcmp(fem2.type,'DGP2')
    feEvalBas2 = @evalP2Bas3D;
end

%% 1. Matrix on noninterface elements
dof1 = fem1.ldof; dof2 = fem2.ldof; nloc = dof1*dof2;
ntID = find(mesh.tLoc > 0); ntN = length(ntID);
AN = fem1.area(ntID);
gxN = fem1.gx(ntID,:); gyN = fem1.gy(ntID,:); gzN = fem1.gz(ntID,:); gw = fem1.gw;
XN = zeros(nloc*ntN, 1);

coefN = feval(fun,gxN,gyN,gzN);
Ibas = cell(dof1,1); Jbas = cell(dof2,1);
for i = 1:dof1
    Ibas{i} = feEvalBas1(fem1.bas(ntID,:,i), gxN, gyN, gzN, d1);
end
for j = 1:dof2
    Jbas{j} = feEvalBas2(fem2.bas(ntID,:,j), gxN, gyN, gzN, d2);
end

IN = reshape(repmat(fem1.t(ntID,:),4,1),nloc*ntN,1);
JN = repmat(reshape(fem2.t(ntID,:),dof2*ntN,1),4,1);
ind = 0;
for i = 1:dof1
    for j = 1:dof2
        XN(ind+1:ind+ntN) = AN.*sum(((Ibas{i}.*(coefN.*Jbas{j})).*gw'),2);
        ind = ind + ntN;
    end
end

%% 2. Matrix on interface elements
AI = femI.area; gw = femI.gw; gxI = femI.gx; gyI = femI.gy; gzI = femI.gz;
ntI = size(femI.t,1); % not number of interface element, but quadrature element
XI = zeros(nloc*ntI, 1);

coefI = feval(fun,gxI,gyI,gzI);
Ibas = cell(dof1,1); Jbas = cell(dof2,1);
for i = 1:dof1
    Ibas{i} = feEvalBas1(femI.bas(:,:,i), gxI, gyI, gzI, d1);
end
for j = 1:dof2
    Jbas{j} = feEvalBas2(femI.bas(:,:,j), gxI, gyI, gzI, d2);
end

II = reshape(repmat(femI.t,4,1),nloc*ntI,1);
JI = repmat(reshape(femI.t,dof2*ntI,1),4,1);
ind = 0;
for i = 1:dof1
    for j = 1:dof2
        XI(ind+1:ind+ntI) = AI.*sum(((Ibas{i}.*(coefI.*Jbas{j})).*gw'),2);
        ind = ind + ntI;
    end
end

repo: ifem


text:
function [stepsize, newx, newkey, ssstats] = ...
                    stepsize_sg(problem, x, d, iter, options, storedb, key) %#ok<INUSD>
% Standard step size selection algorithm for the stochastic gradient method
%
% Given a problem structure, a point x on the manifold problem.d and a
% tangent vector d at x, produces a stepsize (a positive real number) and a
% new point newx obtained by retraction -stepsize*d at x. Additional inputs
% include iter (the iteration number of x, where 0 marks the initial
% guess), an options structure, a storedb database and the key of point x
% in that database. Additional outputs include the key of newx in the
% database, newkey, as well as a structure ssstats collecting statistics
% about the work done during the call to this function.
%
% See in code for the role of available options:
%    options.stepsize_type
%    options.stepsize_init
%    options.stepsize_lambda
%    options.stepsize_decaysteps
%
% This function may create and maintain a structure called sssgmem inside
% storedb.internal. This gives the function the opportunity to remember
% what happened in previous calls.
%
% See also: stochasticgradient

% This file is part of Manopt: www.manopt.org.
% Original authors: Bamdev Mishra and Nicolas Boumal, March 30, 2017.
% Contributors: Hiroyuki Kasai and Hiroyuki Sato.
% Change log:

    % Allow omission of the key, and even of storedb.
    if ~exist('key', 'var')
        if ~exist('storedb', 'var')
            storedb = StoreDB();
        end
        key = storedb.getNewKey(); %#ok<NASGU>
    end

    % Initial stepsize guess.
    default_options.stepsize_init = 0.1;
    % Stepsize evolution type. Options are 'decay', 'fix' and 'hybrid'.
    default_options.stepsize_type = 'decay';
    % If stepsize_type = 'decay' or 'hybrid', lambda is a weighting factor.
    default_options.stepsize_lambda = 0.1;
    % If stepsize_type = 'hybrid', decaysteps states for how many
    % iterations the step size decays before becoming constant.
    default_options.stepsize_decaysteps = 100;

    if ~exist('options', 'var') || isempty(options)
        options = struct();
    end
    options = mergeOptions(default_options, options);

    type = options.stepsize_type;
    init = options.stepsize_init;
    lambda = options.stepsize_lambda;
    decaysteps = options.stepsize_decaysteps;

    switch lower(type)

        % Step size decays as O(1/iter).
        case 'decay'
            stepsize = init / (1 + init*lambda*iter);

        % Step size is fixed.
        case {'fix', 'fixed'}
            stepsize = init;

        % Step size decays only for the few initial iterations.
        case 'hybrid'
            if iter < decaysteps
                stepsize = init / (1 + init*lambda*iter);
            else
                stepsize = init / (1 + init*lambda*decaysteps);
            end

        otherwise
            error(['Unknown options.stepsize_type. ' ...
                   'Should be ''fix'', ''decay'' or ''hybrid''.']);

    end

    % Store some information.
    ssstats = struct();
    ssstats.stepsize = stepsize;

    % Compute the new point and give it a key.
    newx = problem.M.retr(x, d, -stepsize);
    newkey = storedb.getNewKey();

end

author: NicolasBoumal

repo: manopt