text:
function info = Demo()

varargin=cell(1,2);

varargin(1,1)={'train'};
varargin(1,2)={struct('gpus', 1)};

run ../matconvnet/matlab/vl_setupnn ;
addpath ../matconvnet/examples ;

opts.expDir = 'exp/' ;
opts.dataDir = 'exp/data/' ;
opts.modelType = 'tracking' ;
opts.sourceModelPath = 'exp/models/' ;
[opts, varargin] = vl_argparse(opts, varargin) ;

% experiment setup
opts.imdbPath = fullfile(opts.expDir, 'imdb.mat') ;
opts.imdbStatsPath = fullfile(opts.expDir, 'imdbStats.mat') ;
opts.vocEdition = '11' ;
opts.vocAdditionalSegmentations = false ;

global resize;
display=1;

g=gpuDevice(1);
clear g;

test_seq='Skiing';
[config]=config_list(test_seq);

result=CREST_tracking(opts,varargin,config,display);

repo: CREST-Release


text:
function d=drectangle(p,x1,x2,y1,y2)

%   Copyright (C) 2004-2006 Per-Olof Persson. See COPYRIGHT.TXT for details.

d=-min(min(min(-y1+p(:,2),y2-p(:,2)),-x1+p(:,1)),x2-p(:,1));

repo: ifem


text:
function dt = DistanceTransform(b)

%DistanceTransform - Compute distance transform of a binary matrix (distance to nearest 1 value).
%
%  This function assigns to every point (x,y) in a binary matrix B the distance
%  to the nearest point in B with value 1. It uses the Euclidean metric for
%  computing distances. The code is based on the algorithm by Meijster et al.
%  (2000).
%
%  USAGE
%
%    dt = DistanceTransform(b)
%
%    b              binary matrix
%

% Copyright (C) 2014 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% The code is not self explanatory, but is a direct transcription of the algorithms provided
% in Meijster et al. (2000)

[m,n] = size(b);
g = nan(m,n);

for x = 1:n,
	% Scan 1
	if b(1,x),
		g(1,x) = 0;
	else
		g(1,x) = m + n;
	end
	for y = 2:m,
		if b(y,x),
			g(y,x) = 0;
		else
			g(y,x) = 1 + g(y-1,x);
		end
	end
	% Scan 2
	for y = m-1:-1:1,
		if g(y+1,x) < g(y,x),
			g(y,x) = 1 + g(y+1,x);
		end
	end
end

for y = 1:n,
	q = 1;
	s(1) = 1;
	t(1) = 1;
	for u = 2:n,
		% Scan 3
		while true,
			x = t(q);
			i = s(q);
			f1 = (x-i)^2 + g(y,i)^2;
			i = u;
			f2 = (x-i)^2 + g(y,i)^2;
			if f1 <= f2, break; end
			q = q - 1;
			if q < 1, break; end
		end
		if q < 1,
			q = 1;
			s(1) = u;
		else
			i = s(q);
			a = (u^2-i^2+g(y,u)^2-g(y,i)^2);
			b = 2*(u-i);
			w = 1 + floor(a/b);
			if w <= n,
				q = q + 1;
				s(q) = u;
				t(q) = w;
			end
		end
	end
	for u = n:-1:1,
		% Scan 4
		x = u;
		i = s(q);
		dt(y,u) = sqrt((x-i)^2 + g(y,i)^2);
		if u == t(q),
			q = q - 1;
		end
	end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function print_eps(outfile)
%PRINT_EPS stores the figure as eps file
%
%   Usage: print_eps(outfile)
%
%   Input parameters:
%       outfile     - file name to store the figure
%
%   See also: plot_sound_field, print_png

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input parameter =====================================
nargmin = 1;
nargmax = 1;
narginchk(nargmin,nargmax);

%% ===== Main ============================================================
% Font type
set_font_type('Arial');
print(outfile,'-depsc2','-r150');
close;

author: sfstoolbox

repo: sfs-matlab


text:
 function R = Robject(kappa, varargin)
%function R = Robject(kappa, [options])
%|
%| Build roughness penalty regularization "object" based on C = Cdiff() object,
%| for regularized solutions to inverse problems.
%|
%| General form of nonquadratic penalty function:
%|	R(x) = \sumk w_k \pot([Cx]_k), where [Cx]_k = \sum_j c_{kj} x_j.
%|
%| For quadratic case, \potk(t) = t^2/2 in which case
%|	R(x) = x' C' W C x / 2, where W depends on beta and edge_type.
%|
%| Penalty gradient is C' D C x, where D = diag{\wpot_k([Cx]_k)}.
%|
%| in
%|	kappa	[nx,ny[,nz]]	kappa array, or logical support mask
%|
%| options
%|	edge_type '?'		how to handle mask edge conditions (see Cdiff)
%|		'none'		no roughness penalty (NOT DONE)
%|		'tight'		only penalize within-mask differences (default)
%|		'leak'		penalize between mask and neighbors
%|				(its primary use for consistency with ASPIRE)
%|	'order', ?		1st-order or 2nd-order differences (see Cdiff)
%|	'offsets', [?]		offsets to neighboring pixels
%|					(see Cdiff for the defaults)
%|				use '3d:26' to penalize all 13 pairs of nbrs
%|				use '0' for C = I (identity matrix)
%|	'beta', ?		global regularization parameter
%|					default: 2^0
%|	'delta', ?		potential parameter, see potential_func()
%|					or {delta, param}.  default: inf
%|	'potential', '?'	e.g., 'huber', see potential_func()
%|				default: 'quad' for quadratic regularization.
%|	'type_denom', '?'	type of "denominator"
%|					(for quadratic surrogates like SPS)
%|		'matlab'	denominator for SPS
%|		'aspire'	denominator for SPS that matches aspire
%|		'none'		no denominator (default)
%|				(because R.E and R.denom needed only for SPS)
%|	'distance_power', ?	1 classical (default), 2 possibly improved
%|					see penalty_mex()
%|	'user_wt', [?]		User-provided array of penalty weight values
%|					size: [nx,ny[,nz],#offsets]
%|				of dimension [size(mask) length(offsets)].
%|				These are .* the usual wt values for edge_type.
%|	'mask'			Override default: mask = (kappa ~= 0)
%|
%| out
%|	R structure has the following "methods"
%|	R.penal(R, x)	evaluates R(x)
%|	R.cgrad(R, x)	evaluates \cgrad R(x) (column gradient)
%|	R.denom(R, x)	evaluates denominator for separable surrogate
%|	[pderiv pcurv] = feval(R.dercurv, R, C1*x) derivatives and curvatures
%|			for non-separable parabola surrogates
%|	R.diag(R)	diagonal of Hessian of R (at x=0), for preconditioners.
%|	R.C1		differencing matrix, with entries 1 and -1,
%|			almost always should be used in conjunction with R.wt
%|
%| Typical use:	mask = true([128 128]); % or something more conformal
%|		R = Robject(mask, 'beta', 2^7);
%|
%| Copyright 2004-11-14, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if streq(kappa, 'test'), run_mfile_local 'Robject_test', return, end

% option defaults
R.potential = 'quad';
R.beta = 2^0;
R.delta = inf;
R.edge_type = 'tight';
R.type_denom = 'none';
R.distance_power = 1;
R.order = 1; % 1st-order differences, used in call to Cdiff
R.user_wt = []; % user-provided wt values
R.mask = [];
R.offsets = [];

% parse name/value option pairs
R = vararg_pair(R, varargin);

% dimensions, and default offsets
if isempty(R.offsets)
	if ndims(kappa) == 2
		[nx ny] = size(kappa);
		R.offsets = [1 nx nx+1 nx-1];
	elseif ndims(kappa) == 3
		[nx ny nz] = size(kappa);
		R.offsets = [1 nx nx*ny];
	else
		error 'only 2D and 3D done'
	end
elseif streq(R.offsets, '3d:26') % all 26 neighbors (13 pairs)
	[nx ny nz] = size(kappa);
	R.offsets = [1 nx+[0 1 -1] nx*ny+col(outer_sum([-1:1],[-1:1]*nx))'];
end

if ~iscell(R.delta)
	R.delta = {R.delta};
end

R.offsets = int32(R.offsets);

if streq(R.edge_type, 'none')
	error 'todo: unpenalized'
end

if R.beta < 0, warn('Negative beta? This is probably wrong!'), end

R.isquad = streq(R.potential, 'quad');

if isempty(R.mask)
	R.mask = kappa ~= 0; % default is to infer from kappas
end
if ~islogical(R.mask), error 'mask must be logical', end

%
% build plain sparse differencing object "C1", containing 1 and -1 values
% (or identity matrix)
%
R.C_is_I = any(R.offsets == 0);
if R.C_is_I
	if any(R.offsets), error 'identity is offsets=0', end
	R.C1 = diag_sp(ones(sum(R.mask(:)),1)); % identity matrix
	R.wt = R.beta;

else
	R.C1 = Cdiff(R.mask, 'offsets', R.offsets, ...
			'edge_type', 'none', 'order', R.order);

	% wk factors
	R.wt_string = sprintf('wk,%s,%d', R.edge_type, R.order);
	R.wt = penalty_mex(R.wt_string, single(kappa), R.offsets, R.distance_power);
	R.wt = R.beta * single(R.wt(:)); % absorb beta into wk factors
end

% incorporate user provided wk values
if ~isempty(R.user_wt)
	R.wt = single(R.wt .* R.user_wt(:));
	R.user_wt = []; % save memory
end

% trick: for quadratic case, provide a R.C object where C = sqrt(W) * C1
if R.isquad
	R.C = R.C1;
	R.C.cascade_after = diag_sp(sqrt(R.wt));
end

%
% desired potential function
%
R.pot = potential_func(R.potential, R.delta{:});

%
% functions
%
R.dercurv = @Robject_dercurv;
R.handle_denom = @Robject_denom;
R.handle_diag = @Robject_diag;
% trick: the following form of R.penal allows multiple realizations of x
R.penal = @(R, x) ... % 2014-04-27 added 'double' to better match new version
	sum(repmat(R.wt, ncol(x)) .* R.pot.potk(R.pot, R.C1 * x), 'double');
R.cgrad = @(R, x) R.C1' * (diag_sp(R.wt) * R.pot.dpot(R.pot, R.C1 * x));
R.diag = @(R, x) R.handle_diag(R);
R.denom = @(R, x) R.handle_denom(R, x);

%
% precompute denominator for separable quadratic case if needed
%
if ~streq(R.type_denom, 'none') && R.isquad
	if streq(R.type_denom, 'matlab') || streq(R.type_denom, 'aspire')
		t = R.wt .* R.pot.wpot(R.pot, 0);
		if R.C_is_I
			R.denom_max0 = t;
		else
			t = reshape(t, [size(R.mask) length(R.offsets)]);
			if R.order == 1 % fix: order=2 denom?
				t = 2 * t; % "2" because (1,-1) differences
				R.denom_max0 = penalty_mex('diff1,back2', ...
						single(t), R.offsets);
			elseif R.order == 2 % (-1,2,1) so ?
				t = 4 * t; % 4 = |-1| + |2| + |-1|
				R.denom_max0 = penalty_mex('diff2,backA', ...
						single(t), R.offsets);
			else
				error 'order not done for denom'
			end
			R.denom_max0 = single(R.denom_max0(R.mask));
		end
	else
		error(['Unknown type_denom: ' R.type_denom])
	end
end

%
% Robject_denom()
% jth penalty separable surrogate curvature is d_j = \sumk |\ckj| \ck \wpotk
% where \ck = \sumj |\ckj|.
% Here, ck = 2 since there is a +1 and a -1 per row of C1 (for 1st-order)
% Also, |\ckj| = |\ckj|^2 since \ckj = +/- 1, so we can use 'diff1,back2'
%
function denom = Robject_denom(R, x)
if streq(R.type_denom, 'none')
	error 'denom not initialized'
end

if R.isquad
	denom = R.denom_max0;
else
	t = single(R.wt .* R.pot.wpot(R.pot, R.C1 * x));
	if R.C_is_I
		denom = single(t);
	else
		t = reshape(t, [size(R.mask) length(R.offsets)]);
		if R.order == 1
			t = 2 * t; % because (1,-1)
			denom = penalty_mex('diff1,back2', t, R.offsets);
		elseif R.order == 2 % (-1,2,1) so ck = 4
			t = 4 * t; % 4 = |-1| + |2| + |-1|
			denom = penalty_mex('diff2,backA', ...
					single(t), R.offsets);
		else
			error 'order not done'
		end
		denom = single(denom(R.mask));
	end
end

%
% evaluate diagonal of Hessian of (quadratic surrogate for) R at x=0
% \sum_k w_k |c_kj|^2 \wpot(0)
%
function rjj = Robject_diag(R)
t = single(R.wt .* R.pot.wpot(R.pot, 0));
t = reshape(t, [size(R.mask) length(R.offsets)]);
if R.C_is_I
	error 'not done, ask jeff'
else
	rjj = penalty_mex('diff1,back2', t, R.offsets);
	rjj = double(rjj(R.mask));
end

%
% evaluate \dpoti and \wpoti
%
function [deriv, curv] = Robject_dercurv(R, C1x)
deriv = R.wt .* R.pot.wpot(R.pot, C1x) .* C1x;
curv = R.wt .* R.pot.wpot(R.pot, C1x);

%
% Compute both cgrad and denom (of separable surrogate) efficiently.
% Unused for now, but could be used if anonymous functions are too inefficient
% since both R.cgrad and R.denom use C*x so there is redundancy.
% What we really need is an anonymous function that has two output arguments.
% No, the feval with a function_handle will suffice!
%
%function [cgrad, denom] = Robject_cgrad_denom(R, x)
%if R.isquad
%	cgrad = R.C' * (R.C * x);
%	denom = R.denom;
%else
%	Cx = R.C * x;
%	wx = R.wt .* R.pot.wpot(R.pot, Cx);
%	cgrad = R.C' * (wx .* Cx);
%	denom = R.E * wx;
%end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% Reeds Shepp path planner sample code
%
% based on python code from Python Robotics by Atsushi Sakai(@Atsushi_twi)
%
% Peter 3/18
%
% Finds the shortest path between 2 configurations:
% - robot can move forward or backward
% - the robot turns at zero or maximum curvature
% - there are discontinuities in velocity and steering commands (cusps)
% to see what it does run
%
% >> ReedsShepp.test
%
% References::
% - Reeds, J. A.; Shepp, L. A.
%   Optimal paths for a car that goes both forwards and backwards.
%   Pacific J. Math. 145 (1990), no. 2, 367--393.
%   https://projecteuclid.org/euclid.pjm/1102645450

% each path is described by a 3-letter word.
% the algorithm finds a bunch of possible paths, then chooses the shortest
% one.  Each word is represented by a structure with fields:
% - word      a 3-letter sequence drawn from the letters LRLS
% - L         total path length
% - lengths   a 3-vector of lengths, signed to indicate the direction of
%             curvature
% - traj      a cell array of 3xN matrices giving the path for each segment
% - dir       the direction of travel: +1 or -1
%
% TODO: display all the solutions in one figure, as subplots

classdef ReedsShepp < handle
    properties
        best  % the best path
        words
        maxc
    end

    methods
        function obj = ReedsShepp(q0, qf, maxcurv, dl)

            obj.maxc = maxcurv;

            % return the word describing the shortest path
            obj.words = generate_path(q0, qf, maxcurv);

            if isempty(obj.words)
                error('no path');
            end

            % find shortest path
            [~,k] = min( [obj.words.L] );

            obj.best = obj.words(k);

            % add the trajectory
            obj.best = generate_trajectories(obj.best, maxcurv, dl, q0);
        end

        function p = path(obj)
            p = [obj.best.traj{:}]';
        end

        function show(obj)
            for w=obj.words
                fprintf('%s (%g): [%g %g %g]\n', w.word, w.L, w.lengths);
            end
        end

        function plot(obj, varargin)

            opt.circles = [];
            opt.join = [];

            opt = tb_optparse(opt, varargin);

            if ~ishold
                clf
            end
            hold on
            word = obj.best;

            for i=1:3

                if word.dir(i) > 0
                    color = 'b';
                else
                    color = 'r';
                end
                if i == 1
                    x = word.traj{i}(1,:);
                    y = word.traj{i}(2,:);
                else
                    % ensure we join up the lines in the plot
                    x = [x(end) word.traj{i}(1,:)];
                    y = [y(end) word.traj{i}(2,:)];
                end

                if ~isempty(opt.join) && i<3
                    plot(x(end), y(end), opt.join{:});
                end
                if ~isempty(opt.circles)
                    T = SE2(word.traj{i}(:,1));
                    R = 1/obj.maxc;
                    c = T*[0; word.dir(i)*R];

                    plot_circle(c, R, opt.circles)
                    plot_point(c, 'k+')
                end

                plot(x, y, color, 'LineWidth', 2);
            end
            grid on; xlabel('X'); ylabel('Y')
            hold off
            axis equal
            title('Reeds-Shepp path');
        end

        function s = char(obj)
            s = '';
            s = strvcat(s, sprintf('Reeds-Shepp path:  %s, length %f', obj.best.word, obj.best.L));
            s = strvcat(s, sprintf(' segment lengths:  %f %f %f', obj.best.lengths));
        end

        function display(obj)
            disp( char(obj) );
        end
    end

    methods(Static)
        function test()
            maxcurv = 1;
            dl = 0.05;
            q0 = [0 0 pi/4]'; qf = [0 0 pi]';
            p = ReedsShepp(q0, qf, maxcurv, dl)

            p.plot('circles', 'k--', 'join', {'Marker', 'o', 'MarkerFaceColor', 'k'});
        end
    end
end % class ReedsShepp

function out = generate_trajectories(word, maxc, d, q0)

    % initialize the configuration
    p0 = q0;

    % output struct is same as input struct, but we will add:
    %  - a cell array of trajectories
    %  - a vector of directions -1 or +1
    out = word;

    for i=1:3
        m = word.word(i);
        l = word.lengths(i);

        x = [0:d:abs(l) abs(l)];

        p = pathseg(x, sign(l), m, maxc, p0);

        % add new fields to the struct

        if i == 1
            out.traj{i} = p;
        else
            % for subsequent segments skip the first point, same as last
            % point of previous segment
            out.traj{i} = p(:,2:end);
        end
        out.dir(i) = sign(l);

        % initial state for next segment is last state of this segment
        p0 = p(:,end);
    end
end

function q = pathseg(l, dir, m, maxc, p0)
    q0 = p0(:);
    switch m
        case 'S'
            f = @(t,q) dir*[cos(q(3)), sin(q(3)), 0]';
        case {'L', 'R'}
            f = @(t,q) dir*[cos(q(3)), sin(q(3)), dir*maxc]';
    end
    [t,q] = ode45(f, l, q0);
    q = q';  % points are column vectors
end

function words = generate_path(q0, q1, maxc)
    % return a list of all possible words
    q0 = q0(:); q1 = q1(:);
    dq = q1 - q0;
    dth = dq(3);

    xy = rot2(q0(3))' * dq(1:2) * maxc;
    x = xy(1); y = xy(2);

    words = [];
    words = SCS(x, y, dth, words);
    words = CSC(x, y, dth, words);
    words = CCC(x, y, dth, words);

    % account for non-unit curvature
    for i=1:numel(words)
        words(i).lengths = words(i).lengths / maxc;
        words(i).L = words(i).L / maxc;
    end

end

%%
function owords = SCS(x, y, phi, words)

    words = SLS([ x  y  phi], 1, 'SLS', words);
    words = SLS([ x -y -phi], 1, 'SRS', words);

    owords = words;
end

function owords = CCC(x, y, phi, words)

    words = LRL([ x  y  phi],  1, 'LRL', words);
    words = LRL([-x  y -phi], -1, 'LRL', words);
    words = LRL([ x -y -phi],  1, 'RLR', words);
    words = LRL([-x -y  phi], -1, 'RLR', words);

    % backwards
    xb = x * cos(phi) + y * sin(phi);
    yb = x * sin(phi) - y * cos(phi);

    flip = [0 1 0; 1 0 0; 0 0 1];  % flip u and v

    words = LRL([ xb  yb  phi],  flip, 'LRL', words);
    words = LRL([-xb  yb -phi], -flip, 'LRL', words);
    words = LRL([ xb -yb -phi],  flip, 'RLR', words);
    words = LRL([-xb -yb  phi], -flip, 'RLR', words);

    owords = words;
end

function owords = CSC(x, y, phi, words)

    words = LSL([ x  y  phi],  1, 'LSL', words);
    words = LSL([-x  y -phi], -1, 'LSL', words);
    words = LSL([ x -y -phi],  1, 'RSR', words);
    words = LSL([-x -y  phi], -1, 'RSR', words);
    words = LSR([ x  y  phi],  1, 'LSR', words);
    words = LSR([-x  y -phi], -1, 'LSR', words);
    words = LSR([ x -y -phi],  1, 'RSL', words);
    words = LSR([-x -y  phi], -1, 'RSL', words);

    owords = words;
end

% requires LSL, LSR, SLS, LRL

%%

function owords = SLS(q, sign, word, words)
    x = q(1); y = q(2); phi = mod(q(3), 2*pi);

    if y > 0.0 && phi > 0.0 && phi < pi * 0.99
        xd = - y / tan(phi) + x;
        t = xd - tan(phi / 2.0);
        u = phi;
        v = norm( [(x - xd) y]) - tan(phi / 2.0);
        owords = addpath(words, sign*[t, u, v], word);
    elseif y < 0.0 && phi > 0.0 && phi < pi * 0.99
        xd = - y / tan(phi) + x;
        t = xd - tan(phi / 2.0);
        u = phi;
        v = -norm([(x - xd) y]) - tan(phi / 2.0);
        owords = addpath(words, sign*[t, u, v], word);
    else
        owords = words;
    end
end

function owords = LSL(q, sign, word, words)
    x = q(1); y = q(2); phi = mod(q(3), 2*pi);

    [t,u] = cart2pol(x - sin(phi), y - 1.0 + cos(phi));
    if t >= 0.0
        v = angdiff(phi - t);
        if v >= 0.0
            owords = addpath(words, sign*[t, u, v], word);
            return
        end
    end

    owords = words;
end

function owords = LRL(q, sign, word, words)
    x = q(1); y = q(2); phi = mod(q(3), 2*pi);

    [t1,u1] = cart2pol(x - sin(phi), y - 1.0 + cos(phi));

    if u1 <= 4.0
        u = -2.0 * asin(0.25 * u1);
        t = angdiff(t1 + 0.5 * u + pi);
        v = angdiff(phi - t + u);

        if t >= 0.0 && u <= 0.0
            owords = addpath(words, [t, u, v]*sign, word);
            return
        end
    end

    owords = words;
end

function owords = LSR(q, sign, word, words)
    x = q(1); y = q(2); phi = mod(q(3), 2*pi);

    [t1,u1] = cart2pol(x + sin(phi), y - 1.0 - cos(phi));
    u1 = u1^2;
    if u1 >= 4.0
        u = sqrt(u1 - 4.0);
        theta = atan2(2.0, u);
        t = angdiff(t1 + theta);
        v = angdiff(t - phi);

        if t >= 0.0 && v >= 0.0
            owords = addpath(words, sign*[t, u, v], word);
            return
        end
    end

    owords = words;
end

%%
function owords = addpath(words, lengths, ctypes)

    % create a struct to represent this segment
    word.word = ctypes;
    word.lengths = lengths;

    % check same path exist
    for p = words
        if strcmp(p.word, word.word)
            if sum(p.lengths) - sum(word.lengths) <= 0.01
                owords = words;
                return %not insert path
            end
        end
    end

    word.L = sum(abs(lengths));

    % long enough to add?
    if word.L >= 0.01
        owords = [words word];
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function n = numel(varargin)
%NUMEL        Overloaded functions, necessary for Matlab V6.1f
%

% written  10/03/02     S.M. Rump
%

  n = 1;

repo: OpenMAS


text:
function [Forecasteval]=bvarfeval(data_endo_c,data_endo_c_lags,data_exo_c,stringdates3,Fstartdate,Fcenddate,Fcperiods,Fcomp,const,n,p,k,It,Bu,beta_gibbs,sigma_gibbs,forecast_record,forecast_estimates,names,endo,pref)

% function []=bvarfeval(data_endo_c,data_endo_c_lags,data_exo_c,stringdates3,Fstartdate,Fcenddate,Fcperiods,Fcomp,const,n,p,k,It,Bu,beta_gibbs,sigma_gibbs,forecast_estimates,names,endo,datapath)
% calculates, display and saves forecast evaluation results for a BVAR model
% inputs:  - matrix 'data_endo_c': matrix of endogenous data for the forecast evaluation period (i.e. period for which forecast is estimated and actual data exists)
%          - matrix 'data_endo_c_lags': lagged endogenous values to serve as initial conditions to the forecast evaluation period
%          - matrix 'data_exo_c': matrix of predicted data for the forecast period
%          - cell 'stringdates3': date strings for the forecast evaluation period (i.e. period for which forecast is estimated and actual data exists)
%          - string 'Fstartdate': start date of the forecasts
%          - string 'Fcenddate': end date of the forecat evaluation (i.e. period for which forecast is estimated and actual data exists)
%          - integer 'Fcperiods': number of periods for which forecast evaluation can be conducted (i.e.for which forecast is estimated and actual data exists)
%          - integer 'Fcomp': 0-1 value indicating if forecast evaluation is possible
%          - integer 'const': 0-1 value to determine if a constant is included in the model
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'p': number of lags included in the model (defined p 7 of technical guide)
%          - integer 'k': number of coefficients to estimate for each equation in the BVAR model (defined p 7 of technical guide)
%          - integer 'It': total number of iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - integer 'Bu': number of burn-in iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - matrix 'beta_gibbs': record of the gibbs sampler draws for the beta vector
%          - matrix'sigma_gibbs': record of the gibbs sampler draws for the sigma matrix (vectorised)
%          - cell 'forecast_estimates': lower bound, point estimates, and upper bound for the unconditional forecasts
%          - cell 'names': cell containing the excel spreadsheet labels (names and dates)
%          - cell 'endo': list of endogenous variables of the model
%          - string 'datapath': user-supplied path to excel data spreadsheet
% outputs: none

% changed ad 2017-01-06 to store results
Forecasteval.RMSE=[];
Forecasteval.MAE=[];
Forecasteval.MAPE=[];
Forecasteval.Ustat=[];
Forecasteval.CRPS_estimates=[];
Forecasteval.S1_estimates=[];
Forecasteval.S2_estimates=[];

% first, note that forecast evaluation can only be conducted if there is some observable data after the beginning of the forecast
if Fcomp==1

   % preliminary task: obtain a matrix of forecasts over the common periods
   for ii=1:n
   forecast_c(:,ii)=forecast_estimates{ii,1}(2,1:Fcperiods)';
   end
   % then compute the matrix of forecast errors
   ferrors=data_endo_c-forecast_c;

   % compute first the sequential RMSE, from (1.8.11)

   % square the forecast error matrix entrywise
   sferrors=ferrors.^2;
   % sum entries sequentially
   sumsferrors=sferrors(1,:);
   for ii=2:Fcperiods
   sumsferrors(ii,:)=sumsferrors(ii-1,:)+sferrors(ii,:);
   end
   % divide by the number of forecast periods and take square roots to obtain RMSE
   for ii=1:Fcperiods
   Forecasteval.RMSE(ii,:)=((1/ii)*sumsferrors(ii,:)).^0.5;
   end

   % compute then the sequential MAE, from (1.8.12)

   % take the absolute value of the forecast error matrix
   absferrors=abs(ferrors);
   % sum entries sequentially
   sumabsferrors=absferrors(1,:);
   for ii=2:Fcperiods
   sumabsferrors(ii,:)=sumabsferrors(ii-1,:)+absferrors(ii,:);
   end
   % divide by the number of forecast periods to obtain MAE
   for ii=1:Fcperiods
   Forecasteval.MAE(ii,:)=(1/ii)*sumabsferrors(ii,:);
   end

   % compute the sequential MAPE, from (1.8.13)

   % divide entrywise by actual values and take absolute values
   absratioferrors=abs(ferrors./data_endo_c);
   % sum entries sequentially
   sumabsratioferrors=absratioferrors(1,:);
   for ii=2:Fcperiods
   sumabsratioferrors(ii,:)=sumabsratioferrors(ii-1,:)+absratioferrors(ii,:);
   end
   % divide by 100*(number of forecast periods) to obtain MAPE
   for ii=1:Fcperiods
   Forecasteval.MAPE(ii,:)=(100/ii)*sumabsratioferrors(ii,:);
   end

   % compute the Theil's inequality coefficient, from (1.8.14)

   % first compute the left term of the denominator
   % square entrywise the matrix of actual data
   sendo=data_endo_c.^2;
   % sum entries sequentially
   sumsendo=sendo(1,:);
   for ii=2:Fcperiods
   sumsendo(ii,:)=sumsendo(ii-1,:)+sendo(ii,:);
   end
   % divide by the number of forecast periods and take square roots
   for ii=1:Fcperiods
   leftterm(ii,:)=((1/ii)*sumsendo(ii,:)).^0.5;
   end
   % then compute the right term of the denominator
   % square entrywise the matrix of forecast values
   sforecasts=forecast_c.^2;
   % sum entries sequentially
   sumsforecasts=sforecasts(1,:);
   for ii=2:Fcperiods
   sumsforecasts(ii,:)=sumsforecasts(ii-1,:)+sforecasts(ii,:);
   end
   % divide by the number of forecast periods and take square roots
   for ii=1:Fcperiods
   rightterm(ii,:)=((1/ii)*sumsforecasts(ii,:)).^0.5;
   end
   % finally, compute the U stats
   Forecasteval.Ustat=Forecasteval.RMSE./(leftterm+rightterm);

   % now compute the continuous ranked probability score

   % create the cell storing the results
   Forecasteval.CRPS=cell(n,1);
   % loop over endogenous variables
   for ii=1:n
      % loop over forecast periods on which actual data is known
      for jj=1:Fcperiods
      % compute the continuous ranked probability score
      score=bear.crps(forecast_record{ii,1}(:,jj),forecast_estimates{ii,1}(2,jj));
      Forecasteval.CRPS{ii,1}(1,jj)=score;
      end
   end
   Forecasteval.CRPS_estimates=(cell2mat(Forecasteval.CRPS))';

   % finally, consider the big piece: the computation of the log predictive score
   % this part implements algorithm a.8.1, described p 115 of technical guide
   % details of the procedure can be found p 111-115 of the technical guide

   % preliminary task: if there is a constant in the model, concatenate a column of ones to data_exo_c
   if const==1
   data_exo_c=[ones(Fcperiods,1) data_exo_c];
   end

   % create the cells storing the results
   S1=cell(n,1);
   S2=cell(n,1);

   for ll=1:It-Bu
      % step 2: draw beta and sigma
      beta=beta_gibbs(:,ll);
      sigma=reshape(sigma_gibbs(:,ll),n,n);

      % step 3: obtain mu, the mean vector, from (a.8.19)
      temp=bear.forecastsim(data_endo_c_lags,data_exo_c,beta,n,p,k,Fcperiods);
      mu=reshape(temp',n*Fcperiods,1);

      % step 4: obtain upsilon, the covariance matrix, from (a.8.22)
      % recover the A1,A2,... matrices by reshaping beta
      temp=reshape(beta,k,n);
      temp=(temp(1:p*n,:))';
      Amatrices=cell(1,p);
      for ii=1:p
      Amatrices{1,ii}=temp(:,n*(ii-1)+1:n*ii);
      end
      % initiate upsilon with upsilon_1,1
      upsilon=cell(Fcperiods,Fcperiods);
      upsilon{1,1}=sigma;
      % complete the first column using (a.8.23)
      for ii=2:Fcperiods
      upsilon_ij=Amatrices{1,1}*upsilon{ii-1,1};
      summ=min(ii-1,p);
         for jj=2:summ
         upsilon_ij=upsilon_ij+Amatrices{1,jj}*upsilon{ii-jj,1};
         end
      upsilon{ii,1}=upsilon_ij;
      end
      % complete the first row, using symmetry
      for ii=2:Fcperiods
      upsilon{1,ii}=upsilon{ii,1}';
      end
      % now take care of the other columns
      for jj=2:Fcperiods
      % first complete the variance matrix (the diagonal one), using (a.8.24)
      upsilon_jj=sigma;
      summ=min(jj-1,p);
         for kk=1:summ
         upsilon_jj=upsilon_jj+Amatrices{1,kk}*upsilon{jj-kk,jj};
         end
      upsilon{jj,jj}=upsilon_jj;
      % then complete the rest of the column, using (a.8.23)
         for ii=jj+1:Fcperiods
         % initiate upsilon_ij
         upsilon_ij=Amatrices{1,1}*upsilon{ii-1,jj};
         % continue the summation
         summ=min(ii-1,p);
            for kk=2:summ
            upsilon_ij=upsilon_ij+Amatrices{1,kk}*upsilon{ii-kk,jj};
            end
         upsilon{ii,jj}=upsilon_ij;
         end
         % complete the corresponding row, using symmetry
         for kk=jj+1:Fcperiods
         upsilon{jj,kk}=upsilon{kk,jj}';
         end
      end
      upsilon=cell2mat(upsilon);

      % step 5: obtain the predictive density
      % first consider scenario 1 (forecast at period T+i)
      % loop over variables
      for ii=1:n
         % loop over periods
         for jj=1:Fcperiods
         % define R
         R=zeros(1,n*Fcperiods);
         R(1,n*(jj-1)+ii)=1;
         % define the mean
         mean=R*mu;
         % define the variance
         covar=R*upsilon*(R');
         % define the actual value
         actual=data_endo_c(jj,ii);
         % determine the density
         [~,density]=bear.mndensity(actual,mean,covar,1);
         % record the result
         S1{ii,1}(ll,jj)=density;
         end
      end
      % then consider scenario 2 (forecast from beginning until period T+i)
      % loop over variables
      for ii=1:n
      R=[];
         % loop over periods
         for jj=1:Fcperiods
         % define R
         R=[R;zeros(1,n*Fcperiods)];
         R(jj,n*(jj-1)+ii)=1;
         % define the mean vector
         mean=R*mu;
         % define the covariance matrix
         covar=R*upsilon*(R');
         % define the vector of actual values
         actual=data_endo_c(1:jj,ii);
         % determine the density
         [~,density]=bear.mndensity(actual,mean,covar,jj);
         % record the result
         S2{ii,1}(ll,jj)=density;
         end
      end
   end

   % step 6: compute the log predictive score from (a.8.26)
   Forecasteval.S1_estimates=cell(n,1);
   Forecasteval.S2_estimates=cell(n,1);
   for ii=1:n
      for jj=1:Fcperiods
      Forecasteval.S1_estimates{ii,1}(1,jj)=log((1/(It-Bu))*sum(S1{ii,1}(:,jj)));
      Forecasteval.S2_estimates{ii,1}(1,jj)=log((1/(It-Bu))*sum(S2{ii,1}(:,jj)));
      end
   end
   Forecasteval.S1_estimates=(cell2mat(Forecasteval.S1_estimates))';
   Forecasteval.S2_estimates=(cell2mat(Forecasteval.S2_estimates))';

% if forecast evaluation is not possible, do not do anything
elseif Fcomp==0
end

% now, print the results and display them

filelocation=fullfile(pref.results_path, [pref.results_sub '.txt']);
fid=fopen(filelocation,'at');

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

Fevalinfo='Forecast evaluation:';
fprintf('%s\n',Fevalinfo);
fprintf(fid,'%s\n',Fevalinfo);

fprintf('%s\n','');
fprintf(fid,'%s\n','');

% if forecast evaluation is not possible, return a message to signal it

if Fcomp==0

finfo1=['Forecast evaluation cannot be conducted.'];
fprintf('%s\n',finfo1);
fprintf(fid,'%s\n',finfo1);
finfo2=['Forecasts start in ' Fstartdate ', while observable data is available only until ' names{end,1} '.'];
fprintf('%s\n',finfo2);
fprintf(fid,'%s\n',finfo2);
finfo3=['To obtain forecast evaluation, the forecast start date must be anterior to the end of the data set.'];
fprintf('%s\n',finfo3);
fprintf(fid,'%s\n',finfo3);

% if forecast evaluation is possible, display the results
elseif Fcomp==1

finfo1=['Evaluation conducted over ' num2str(Fcperiods) ' periods (from ' Fstartdate ' to ' Fcenddate ').'];
fprintf('%s\n',finfo1);
fprintf(fid,'%s\n',finfo1);

   % loop over endogenous variables
   for ii=1:n

   fprintf('%s\n','');
   fprintf(fid,'%s\n','');

   endoinfo=['Endogenous: ' endo{ii,1}];
   fprintf('%s\n',endoinfo);
   fprintf(fid,'%s\n',endoinfo);

   temp='fprintf(''%12s';
      for jj=1:Fcperiods-1
      temp=[temp ' %10s'];
      end
   temp=[temp ' %10s\n'','''''];
      for jj=1:Fcperiods
      temp=[temp ',''' stringdates3{jj,1} ''''];
      end
   temp=[temp ');'];
   eval(temp);
   temp='fprintf(fid,''%12s';
      for jj=1:Fcperiods-1
      temp=[temp ' %10s'];
      end
   temp=[temp ' %10s\n'','''''];
      for jj=1:Fcperiods
      temp=[temp ',''' stringdates3{jj,1} ''''];
      end
   temp=[temp ');'];
   eval(temp);

   label='RMSE:       ';
   values=Forecasteval.RMSE(1:Fcperiods,ii)';
   temp='fprintf(''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);
   temp='fprintf(fid,''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);

   label='MAE:        ';
   values=Forecasteval.MAE(1:Fcperiods,ii)';
   temp='fprintf(''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);
   temp='fprintf(fid,''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);

   label='MAPE:       ';
   values=Forecasteval.MAPE(1:Fcperiods,ii)';
   temp='fprintf(''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);
   temp='fprintf(fid,''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);

   label='Theil''s U:  ';
   values=Forecasteval.Ustat(1:Fcperiods,ii)';
   temp='fprintf(''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);
   temp='fprintf(fid,''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);

   label='CRPS:       ';
   values=Forecasteval.CRPS_estimates(1:Fcperiods,ii)';
   temp='fprintf(''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);
   temp='fprintf(fid,''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);

   label='Log score 1:';
   values=Forecasteval.S1_estimates(1:Fcperiods,ii)';
   temp='fprintf(''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);
   temp='fprintf(fid,''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);

   label='Log score 2:';
   values=Forecasteval.S2_estimates(1:Fcperiods,ii)';
   temp='fprintf(''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);
   temp='fprintf(fid,''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);

    end

end

fclose(fid);

author: european-central-bank

repo: BEAR-toolbox


text:
function c = correlation_power ( n, rho, rho0, e )

%*****************************************************************************80
%
%% CORRELATION_POWER evaluates the power correlation function.
%
%  Discussion:
%
%    The power correlation is
%
%      C(rho) = ( 1 - |rho| )^e  if 0 <= |rho| <= 1
%             = 0                otherwise
%
%      The constraint on the exponent is 2 <= e.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 October 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of arguments.
%
%    Input, real RHO(N,1), the arguments.
%    0.0 <= RHO.
%
%    Input, real RHO0, the correlation length.
%    0.0 < RHO0.
%
%    Input, real E, the exponent.
%    E has a default value of 2;
%    2 <= E.
%
%    Output, real C(N,1), the correlations.
%
  if ( nargin < 4 )
    e = 2.0;
  end

  c = zeros ( n, 1 );

  rho1 = abs ( rho ( : ) ) / rho0;

  i = find ( rho1 <= 1 );

  c(i) = ( 1.0 - rho1(i) ) .^ e;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function updateLineseries(obj,plotIndex)

%----SCATTER FIELDS----%

% x - [DONE]
% y - [DONE]
% r - [HANDLED BY SCATTER]
% t - [HANDLED BY SCATTER]
% mode - [DONE]
% name - [NOT SUPPORTED IN MATLAB]
% text - [DONE]
% error_y - [HANDLED BY ERRORBAR]
% error_x - [NOT SUPPORTED IN MATLAB]
% connectgaps - [NOT SUPPORTED IN MATLAB]
% fill - [HANDLED BY AREA]
% fillcolor - [HANDLED BY AREA]
% opacity --- [TODO]
% textfont - [NOT SUPPORTED IN MATLAB]
% textposition - [NOT SUPPORTED IN MATLAB]
% xaxis [DONE]
% yaxis [DONE]
% showlegend [DONE]
% stream - [HANDLED BY PLOTLYSTREAM]
% visible [DONE]
% type [DONE]

% MARKER
% marler.color - [DONE]
% marker.size - [DONE]
% marker.line.color - [DONE]
% marker.line.width - [DONE]
% marker.line.dash - [NOT SUPPORTED IN MATLAB]
% marker.line.opacity - [NOT SUPPORTED IN MATLAB]
% marker.line.smoothing - [NOT SUPPORTED IN MATLAB]
% marker.line.shape - [NOT SUPPORTED IN MATLAB]
% marker.opacity --- [TODO]
% marker.colorscale - [NOT SUPPORTED IN MATLAB]
% marker.sizemode - [NOT SUPPORTED IN MATLAB]
% marker.sizeref - [NOT SUPPORTED IN MATLAB]
% marker.maxdisplayed - [NOT SUPPORTED IN MATLAB]

% LINE

% line.color - [DONE]
% line.width - [DONE]
% line.dash - [DONE]
% line.opacity --- [TODO]
% line.smoothing - [NOT SUPPORTED IN MATLAB]
% line.shape - [NOT SUPPORTED IN MATLAB]

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

%-AXIS INDEX-%
axIndex = obj.getAxisIndex(obj.State.Plot(plotIndex).AssociatedAxis);

%-PLOT DATA STRUCTURE- %
plot_data = get(obj.State.Plot(plotIndex).Handle);

%-CHECK FOR MULTIPLE AXES-%
[xsource, ysource] = findSourceAxis(obj,axIndex);

%-AXIS DATA-%
eval(['xaxis = obj.layout.xaxis' num2str(xsource) ';']);
eval(['yaxis = obj.layout.yaxis' num2str(ysource) ';']);

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

%-scatter xaxis-%
obj.data{plotIndex}.xaxis = ['x' num2str(xsource)];

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

%-scatter yaxis-%
obj.data{plotIndex}.yaxis = ['y' num2str(ysource)];

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

%-scatter type-%
obj.data{plotIndex}.type = 'scatter';

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

%-scatter visible-%
obj.data{plotIndex}.visible = strcmp(plot_data.Visible,'on');

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

%-scatter x-%
obj.data{plotIndex}.x = plot_data.XData;

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

%-scatter y-%
obj.data{plotIndex}.y = plot_data.YData;

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

if isfield(plot_data,'ZData')
    if any(plot_data.ZData)
        %-scatter z-%
        obj.data{plotIndex}.z = plot_data.ZData;

        %overwrite type
        obj.data{plotIndex}.type = 'scatter3d';
    end
end

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

%-scatter name-%
obj.data{plotIndex}.name = plot_data.DisplayName;

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

%-scatter mode-%
if ~strcmpi('none', plot_data.Marker) ...
        && ~strcmpi('none', plot_data.LineStyle)
    mode = 'lines+markers';
elseif ~strcmpi('none', plot_data.Marker)
    mode = 'markers';
elseif ~strcmpi('none', plot_data.LineStyle)
    mode = 'lines';
else
    mode = 'none';
end

obj.data{plotIndex}.mode = mode;

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

%-scatter line-%
obj.data{plotIndex}.line = extractLineLine(plot_data);

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

%-scatter marker-%
obj.data{plotIndex}.marker = extractLineMarker(plot_data);

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

%-scatter showlegend-%
leg = get(plot_data.Annotation);
legInfo = get(leg.LegendInformation);

switch legInfo.IconDisplayStyle
    case 'on'
        showleg = true;
    case 'off'
        showleg = false;
end

obj.data{plotIndex}.showlegend = showleg;

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

end

author: goGPS-Project


text:
function handle = xyzhumanevaVisualise3d(pos,fid)

% XYZHUMANEVAVISUALISE3D
%
% COPYRIGHT : Carl Henrik Ek and Neil Lawrence, 2008

% MOCAP

if(nargin<2)
  fid = 2;
end

figure(fid);

% determine parametrization
if(max(max(abs(pos)))<100)
  type = 'jon';
else
  type = 'raquel';
end

joint = xyzhumaneva2joint(pos);
figure(fid);
handle = xyzhumanevaDraw(joint);axis equal;view([1 1 1]);
switch type
 case 'raquel'
  set(gca,'XLim',[-300 300],'YLim',[-300 300],'ZLim',[-800 800]);
 case 'jon'
  set(gca,'XLim',[-0.35 0.35],'YLim',[-0.35 0.35],'ZLim',[-1 0.3]);
end

return

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%% THE AGENT BASE CLASS (agent.m) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This class is designed to define a generic agent and import this variables
% into the simulation space for the purpose of multi-vehicle control simulation.
% The agent object is a child of the objectDefintion; the prinicple
% distinctions being:
% sensorRange      - The agent is assumed capable of observing its
%                    surroundings.
% controlFrequency - The frequency at which the control cycle is computed.

% Author: James A. Douthwaite

classdef agent < objectDefinition & agent_tools
    %% AGENT BASE CLASS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % This class contains the basic properties of a generic agent, neither
    % aerial or ground based.
    properties
        % DEFAULT BEHAVIOUR
        v_nominal = 2;                  % Default nominal linear speed (m/s)
        v_max  = 4;                     % Default maximal linear speed (m/s)
        w_max  = 0.5;                   % Default maximal angular speed (rad/s)
        detectionRadius = inf;
        % WAYPOINTS
        targetWaypoint;                 % The current waypoint target
        achievedWaypoints;              % The agents list of locally achieved waypoints
        % DYNAMIC PARAMETERS
        DYNAMICS;
        % SENSOR MEASUREMENT PARAMETERS
        SENSORS = struct('range',inf);  % Only generic parameter is a visual horizon
        % AGENT-SIDE OUTPUT DATA
        DATA;                           % The output container for agent-side data.
    end

    %% ///////////////////////// MAIN METHODS /////////////////////////////
    methods
        % Constructor
        function [this] = agent(varargin)
            % This function is to construct the agent object using the
            % object defintions held in the 'objectDefinition' base class.
            % INPUTS:
            % namestr - Assigned name string
            % OUTPUTS:
            % this     - The constructed object

            % Call the super class
            this@objectDefinition(varargin);

            % //////////////////// SENSOR PARAMETERS //////////////////////
            [this.SENSORS] = this.GetDefaultSensorParameters();               % Default sensing
            % /////////////////////////////////////////////////////////////
            this.SetGLOBAL(this.CreateGLOBAL());    % Assign (agent) GLOBAL structure
            this.DYNAMICS = this.CreateDYNAMICS();  % Assign (agent) DYNAMICS structure
            this.SetBufferSize(1);                  % Assign default buffer size

            % Assign defaults
            this.radius = 0.5;                      % Default radius (m)
            this.localState = zeros(12,1);          % Assign state

            % //////////////// Check for user overrides ///////////////////
            % - It is assumed that overrides to the properties are provided
            %   via the varargin structure.
            this = this.ApplyUserOverrides(varargin);                         % Recursive overrides
            % /////////////////////////////////////////////////////////////
        end
        % Setup (agent default)
        function [this] = setup(this,localXYZvelocity,localXYZrotations)   % [x y z phi theta psi]
            % This function is called in order to build the initial state
            % vector for the generic agent class 'objectDefinition'.
            % INPUTS:

            % ASSUMPTIONS:
            % - The designed state is described in the ENU axes.
            % - The state is of the form:
            % - [x y z phi theta psi dx dy dz dphi dtheta dpsi]

            switch this.Is3D()
                case false
                    % Initialises the state vector [x y psi ; dx dy dpsi]
                    this.setup_2DVelocities(localXYZvelocity,localXYZrotations);
                case true
                    % Initialises the state vector [x y z phi theta psi ; dx dy dz dphi dtheta dpsi]
                    this.setup_3DVelocities(localXYZvelocity,localXYZrotations)
                otherwise
                    error('Object dimensionality unknown');
            end
        end
        % Main
        function [this] = main(this,ENV,varargin)
            % This function is designed to house a generic agent process
            % cycle that results in an acceleration vector in the global axis.
            % INPUTS:
            % ENV     - The TIME simulations structure
            % varargin - Cell array of inputs
            % OUTPUTS:
            % this      - The updated project

            % PLOT AGENT FIGURE
            visualiseProblem = 0;
            visualiseAgent = 1;
            if this.objectID == visualiseAgent && visualiseProblem == 1
                %overHandle = figure('name','testFigure');
                overHandle = gcf;
                ax = gca;
                hold on; grid on;
                axis equal;
                xlabel('x_{m}'); ylabel('y_{m}'); zlabel('z_{m}');
            end

            % //////////// CHECK FOR NEW INFORMATION UPDATE ///////////////
            % UPDATE THE AGENT WITH THE NEW ENVIRONMENTAL INFORMATION
            [this,obstacleSet,agentSet] = this.GetAgentUpdate(ENV,varargin{1});

            % /////////////////// WAYPOINT TRACKING ///////////////////////
            % Design the current desired trajectory from the waypoint.
            [headingVector] = this.GetTargetHeading();
            desiredVelocity = headingVector*this.v_nominal;

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %                                                             %
            %                         DO NOTHING                          %
            %                                                             %
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            % PASS THE DESIRED VELOCITY TO THE DEFAULT CONTROLLER
            [this] = this.Controller(ENV.dt,desiredVelocity);
        end
    end

    %% ////////////////////// BASIC UPDATE FUNCTIONS //////////////////////
    methods
        % AGENT UPDATE (IDEAL) - Update from perfect environmental knowledge
        function [this,obstacleSet,agentSet,waypointSet] = GetAgentUpdate(this,ENV,observedobjects)
            % This function computes the agents process for updating its
            % knowledge of the environment.
            % INPUTS:
            % ENV.dt              - The unit timstep
            % observedobjects - The full observed object structure
            % OUTPUTS:
            % this             - The updated object
            % obstacleSet     - The list of obstacle structures
            % waypointSet     - The list of waypoint structures

            agentSet = []; obstacleSet = []; waypointSet = [];             % No observed objects

            % Check #1 - Nothing is observed..
            if isempty(observedobjects)
                return
            end

            % Input sanity check
            assert(isstruct(ENV),'Expecting environment update structure.');
            assert(isnumeric(ENV.dt),'The time-step must be a numeric timestep.');
            assert(isstruct(observedobjects),'Second parameter is a vector of observation structures.');

            % Update agent memory structure
            for entry = 1:numel(observedobjects)
                % Apply sensor model if there is one
                sensedobject = this.SensorModel(ENV.dt,observedobjects(entry));
                % Update memory structure from measurements
                this = this.UpdateMemoryFromObject(ENV.currentTime,sensedobject);
            end

            % SORT object SET BY PRIORITY
            [this] = this.UpdateMemoryOrderByField('priority');

            % DECERN object TYPES
            % this.MEMORY now contains an updated understanding of all the
            % objects in its visual horizon. From this, the waypoints and
            % obstacles must be parsed.

            % PULL MEMORY ITEMS FOR LATER MANIPULATION
            agentSet    = this.MEMORY([this.MEMORY.type] == OMAS_objectType.agent);
            waypointSet = this.MEMORY([this.MEMORY.type] == OMAS_objectType.waypoint);
            obstacleSet = this.MEMORY([this.MEMORY.type] == OMAS_objectType.obstacle); % Differenciate between the different object types

            % UPDATE THE TARGET WAYPOINT, PRIORTIES
            [this] = this.UpdateTargetWaypoint(waypointSet);                    % Update the target waypoint and heading vector
        end
    end

    %% //////////////////////// SENSING FUNCTIONS /////////////////////////
    methods
        % SENSOR MODEL - TOP LEVEL (Default)
        function [observedobject] = SensorModel(this,dt,observedobject)
            % This function provides an overridable method resembling the
            % sensor model through which all objects are processed.

            % Input sanity check
            assert(isnumeric(dt) && numel(dt) == 1,'The time step value is invalid.');
            assert(numel(observedobject.position) <= 3,'Expecting a position value to be [2x1].');
            assert(numel(observedobject.velocity) <= 3,'Expecting a velocity value to be [2x1].');

            % Get measurements from a camera
            [psi_j,theta_j,alpha_j] = this.GetCameraMeasurements(observedobject);
            % Override ideal parameters with camera model
            observedobject.heading = psi_j;
            observedobject.elevation = theta_j;
            observedobject.width = alpha_j;
            % Get the range estimate
            [observedobject.range] = this.GetRangeFinderMeasurements(observedobject);
        end
        % SENSOR MODEL - CAMERA & RANGE FINDER
        function [psi_j,theta_j,alpha_j] = GetCameraMeasurements(this,observedobject)
            % This function takes the simulation data and calculates the
            % spherical position and radius that would otherwise be sensed
            % by the system.
            % INPUTS:
            % obstacleData - The object observation structure
            % OUTPUTS:
            % range        - The obstacles apparent range
            % azimuth      - The obstacles angular position in the XY plane
            % elevation    - The obstacles angular position in the XZ plane
            % angularWidth - The obstacles angular width in the azimuth

            % Observed pixel coordinates
            psi_j   = observedobject.heading   + this.SENSORS.sigma_camera*randn(1);
            theta_j = observedobject.elevation + this.SENSORS.sigma_camera*randn(1);
            % Observed width
            alpha_j = observedobject.width + this.SENSORS.sigma_camera*randn(1);      % Uncertainty in the angular measurement
        end
        % SENSOR MODEL - RANGE FINDER
        function [d_j] = GetRangeFinderMeasurements(this,observedobject)
            % This function takes a simulation data resembling an object
            % and calculates the apparent range to the agent.

            % EMULATE MEASURED POSITIONAL VARIABLES (measured = range(true) + distortion)
            d_j = observedobject.range + this.SENSORS.sigma_rangeFinder*randn(1);
        end
        % SENSOR MODEL - LOCAL GPS & PITOT TUBE
        function [p_i,v_i,r_i] = GetAgentMeasurements(this)
            % This function makes estimates of the agent's current position
            % and velocity states (defined in this.localState).

            if this.Is3D
                positionIndices = 1:3;
                velocityIndices = 7:9;
            else
                positionIndices  = 1:2;
                velocityIndices = 4:5;
            end
            % GET THE LOCAL ABSOLUTE VARIABLES
            p_i =  this.localState(positionIndices,1) + this.SENSORS.sigma_position*rand(numel(positionIndices),1); % Get the absolute position measurement
            v_i =  this.localState(velocityIndices,1) + this.SENSORS.sigma_velocity*rand(numel(velocityIndices),1); % Get the absolute velocity measurement
            % RADIUS IS ASSUMED KNOWN
            r_i = this.radius();
        end
        % CALCULATE THE NEW STATE ESTIMATE
        function [position,velocity] = linearStateEstimation(this,dt,p0,v0,p1)
            % This function takes the previous known state of the obstacle
            % and estimates its new state.

            % GET THE POSITION
            dX = (p1 - p0);
            velocity = dX/dt;    % Defines the average velocity
            position = p1;
            % NO PREVIOUS VELOCITY RECORDED
            if any(isnan(v0))
                return
            end
        end
    end
    methods (Static)
        % SENSOR MODEL - REPRESENTATIVE SENSING
        function [SENSORS] = GetCustomSensorParameters()
            % This function is designed to populate the SENSOR structure
            % with representative sensing parameters for the assembly of
            % the sensing intervals.
            % BUILD THE SENSOR FACTOR
            SENSORS = struct(...
                'range',inf,...             % Assume the agent has perfect environmental knowledge (m)
                'sigma_position',0.5,...    % Accurate to within 0.5m
                'sigma_velocity',0.1,...    % Accurate to within 0.1m/s
                'sigma_rangeFinder',0.1,... % Accurate to within 0.1m
                'sigma_camera',5.208E-5,... % One pixel in a 1080p image
                'sampleFrequency',inf);     % object has perfect precision
        end
        % SENSOR MODEL - PERFECT SENSING
        function [SENSORS] = GetDefaultSensorParameters()
            % This function is designed to populate the SENSOR field with
            % perfect sensor capabilities.
            % BUILD THE SENSOR FACTOR
            SENSORS = struct(...
                'range',inf,...             % Assume the agent has perfect environmental knowledge (m)
                'sigma_position',0.0,...    % Perfect position measurement
                'sigma_velocity',0.0,...    % Perfect velocity measurement
                'sigma_rangeFinder',0.0,... % Perfect range acquistion
                'sigma_camera',0.0,...      % Infinte resolution
                'sampleFrequency',inf);     % object has perfect precision
        end
    end

    %% //////////////////////// DYNAMICS & CONTROL ////////////////////////
    methods
        % Simple controller
        function [this] = Controller(this,dt,desiredVelocity)
            % This function computes the agents change in state as a result
            % of the desired velocity vector

            % Input sanity check #1 - Is feasible
            assert(isnumeric(dt) && numel(dt) == 1,'The time step must be a numeric scalar.');
            assert(isnumeric(desiredVelocity) && numel(desiredVelocity) == 3,'Requested velocity vector is not a 2D numeric vector');
            assert(~any(isnan(desiredVelocity)),'The desired velocity vector contains NaNs.');

            % ///////////// UPDATE object GLOBAL PROPERTIES ///////////////
            % Input sanity check #2 - Zero vector
            [heading,speed] = this.nullVelocityCheck(desiredVelocity);

            % Get the relative heading
            [dPsi,dTheta]   = this.GetVectorHeadingAngles([1;0;0],heading); % Relative heading angles
            % The desired heading rate
            omega = ([0;dTheta;-dPsi] - zeros(3,1))/dt;

            % Apply kinematic constraints to state changess
            [omega_actual,speed_actual] = this.ApplyKinematicContraints(dt,omega,speed);

            % OMIT TRAJECTORY CHANGES IF IDLE
            if this.IsIdle()
                omega_actual = zeros(3,1);
                speed_actual = 0;
                this.v_nominal = 0;
            end

            dX = [speed_actual;0;0;omega_actual];

            % /////////// SIMPLE DYNAMICS + PURE TRANSLATION //////////////
            [dX] = this.SimpleDynamics(this.localState(1:6),[speed_actual;0;0],omega_actual);
            this.localState(1:6)  = this.localState(1:6) + dt*dX;
            this.localState(7:12) = dX;

            % ////// GLOBAL UPDATE FOR STATE WITH RETAINED VELOCITES //////
            [this] = this.GlobalUpdate_3DVelocities(dt,this.localState);
        end
        % SPEED AND HEADING (4D) PID CONTROLLER
        function [this] = Controller_PID(this,dt,desiredVelocity)
            % This function is desiged to generate feedack from a desired
            % local velocity vector.
            % INPUTS:
            % targetHeading - The unit heading vector
            % targetSpeed   - The target speed in that direction
            % OUTPUTS:
            % heading_fb    - The feedback on the agent heading
            % speed_fb      - The feedback on the agent speed
            % this           - The error-updated agent object

            % Input sanity check
            assert(isnumeric(dt) && numel(dt) == 1,'The time step must be a numeric scalar.');
            assert(isnumeric(desiredVelocity) && numel(desiredVelocity) == 3,'Requested velocity vector must be a 3D local vector.');
            assert(~any(isnan(desiredVelocity)),'The requested velocity contains NaNs.');

            % NUMERICAL SANITY CHECK TWO
            [unitDirection,desiredSpeed] = this.nullVelocityCheck(desiredVelocity);

            % APPLY SPEED CONSTRAINT
            if abs(desiredSpeed) > this.v_max
                desiredSpeed = sign(desiredSpeed)*this.v_max;
            end

            % GET THE EQUIVALENT HEADING ANGLE
            [dPsi,dTheta] = this.GetVectorHeadingAngles([1;0;0],unitDirection); % Relative heading angles
            dHeading = [0;dTheta;-dPsi];

            % RELATIVE SPEED
            e_speed = desiredSpeed - norm(this.localState(7:9));            % The speed error
            controlError = [e_speed;dHeading];                             % -ve in the NED control frame of reference

            % TUNING PARAMETERS
            Kp_linear = 0.8;
            Kd_linear = 0;
            Kp_angular = 1;
            Kd_angular = 0;

            % CALCULATE THE CONTROL FEEDBACK
            control_fb = diag([Kp_linear Kp_angular Kp_angular Kp_angular])*controlError + ...
                diag([Kd_linear Kd_angular Kd_angular Kd_angular])*(controlError - this.priorError);
            % REMEMBER PREVIOUS ERROR
            this.priorError = controlError;

            % PARSE THE CONTROL FEEDBACK FOR EXTERNAL USE
            speedFeedback = control_fb(1);                                 % Absolute speed input
            headingFeedback = control_fb(2:4);
            omega = headingFeedback/dt;

            % /////////// SIMPLE DYNAMICS + PURE TRANSLATION //////////////
            [dX] = this.SimpleDynamics(this.localState(1:6),[speedFeedback;0;0],omega);
            this.localState(1:6)  = this.localState(1:6) + dt*dX;
            this.localState(7:12) = dX;

            % ///////////// UPDATE object GLOBAL PROPERTIES ///////////////
            if this.GLOBAL.idleStatus
                this.localState(7:12) = zeros(6,1);
            end

            % ////// GLOBAL UPDATE FOR STATE WITH RETAINED VELOCITES //////
            [this] = this.GlobalUpdate_3DVelocities(dt,this.localState);
        end
        % Apply the DYNAMICS constraints
        function [omega_actual,speed_actual] = ApplyKinematicContraints(this,dt,omega,speed)
            % This function is designed to take a desired heading and
            % velocity and compare them to the kinematic contraints of the
            % agent.

            % Determine agent dimension
            if this.Is3D
                currentVelocity = this.localState(7:9);
                currentOmega    = this.localState(10:12);
            else
                currentVelocity = this.localState(4:5);
                currentOmega    = this.localState(6);
            end

            % For clarity
            v_lim      = this.DYNAMICS.maxLinearVelocity;
            dv_lim     = this.DYNAMICS.maxLinearAcceleration;
            omega_lim  = this.DYNAMICS.maxAngularVelocity;
            dOmega_lim = this.DYNAMICS.maxAngularAcceleration;

            % /////////////// Constrained heading change //////////////////
            % The implied angular acceleration
            omega_dot = (omega - currentOmega)/dt;
            % Constrain against heading acceleration limits
            bounded_omega_dot = boundValue(omega_dot,-dOmega_lim,dOmega_lim);
            % Calculate actual heading velocity limits
            omega = currentOmega + bounded_omega_dot*dt;
            % Constrain against heading rate limits
            omega_actual = boundValue(omega,-omega_lim,omega_lim);

            % ///////////////// Constrained speed change //////////////////
            % The current speed
            currentSpeed = norm(currentVelocity);
            % The proposed acceleration
            s_dot = (speed - currentSpeed)/dt;
            % Constrain against linear acceleration limits
            bounded_s_dot = boundValue(s_dot,-dv_lim(1),dv_lim(1));
            % Calcullate the actual speed
            speed = currentSpeed + bounded_s_dot*dt;
            % Interesect the velocity with permissible velocities
            speed_actual = boundValue(speed,-v_lim(1),v_lim(1));
        end
        % Create the DYNAMICS structure
        function [DYNAMICS] = CreateDYNAMICS(this)
            % Define basic DYNAMIC container
            [DYNAMICS] = this.CreateDYNAMICS_default();
        end
        % Create the default DYNAMICS constraint structure
        function [DYNAMICS] = CreateDYNAMICS_default(this)
            % Define an infinite "throw" range by default
            DYNAMICS = struct();
            % Define kinematic constraints
            if this.Is3D
                % States defined as [x,y,z,phi,theta,psi,dx,dy,dz,dphi,dtheta,dpsi]^T
                DYNAMICS.maxLinearVelocity      = ones(3,1)*this.v_max;    % Limits on the agents linear velocity
                DYNAMICS.maxLinearAcceleration  = inf(3,1);                % Limits on the agents linear acceleration
                DYNAMICS.maxAngularVelocity     = ones(3,1)*this.w_max;    % Limits on the agents angular velocity
                DYNAMICS.maxAngularAcceleration = inf(3,1);                % Limits on the agents angular acceleration
            else
                % States defined as [x,y,psi,dx,dy,dpsi]^T
                DYNAMICS.maxLinearVelocity      = ones(2,1)*this.v_max;    % Limits on the agents linear velocity
                DYNAMICS.maxLinearAcceleration  = inf(2,1);                % Limits on the agents linear acceleration
                DYNAMICS.maxAngularVelocity     = this.w_max;              % Limits on the agents angular velocity
                DYNAMICS.maxAngularAcceleration = inf(1,1);                % Limits on the agents angular acceleration
            end
        end
    end

    %% ////////////////////// GENERAL STATIC METHODS //////////////////////
    methods (Static)
        % DEFINE THE PITCH AND YAW TO MEET TARGET HEADING (2D & 3D)
        function [lambda,theta] = GetVectorHeadingAngles(V,U)
            % This function calculates the Line Of Sight (LOS) [horezontal]
            % angle which aids in defining the bank angle as a control input to the
            % dynamics.
            % INPUTS:
            % V - The current unity velocity vector
            % U - The unit correction vector
            % OUTPUTS:
            % lambda - The azimuth angle (2D)
            % theta  - The pitch angle   (3D)

            % NORMALISE THE ELEMENTS
            V = V/norm(V);
            U = U/norm(U);
            % GET THE HORIZONTAL ELEMENTS
            Vh = [V(1);V(2);0];
            Uh = [U(1);U(2);0];                                        % Reject the vertical elements
            % GET THE LINE OF SIGHT ANGLE
            rotationAxis = cross(Vh,Uh);

            % HANDLE 3D CASE (ELEVATION & LOS)
            if numel(V) == 3 && numel(U) == 3
                % GET THE LINE OF SIGHT ANGLE
                lambda = sign(rotationAxis(3))*acos(dot(Vh,Uh)/norm(Vh));  % Get the angle, signed by the direction of its cross product
                % GET THE ELEVATION ANGLE
                theta = atan2(U(3),norm(Uh));
            else
                % HANDLE 2D CASE (LOS)
                % GET THE LINE OF SIGHT ANGLE
                lambda = sign(rotationAxis(3))*acos(dot(Vh,Uh)/norm(Vh));
                % GET THE SUDO ELEVATION ANGLE
                theta = 0;
            end
        end
        % CALCULATE THE RADIUS
        function [r]    = GetRadiusFromAngularWidth(d,alpha)
            % Calculate the radius of the object
            r = (sin(alpha/2)/(1-sin(alpha/2)))*d;
        end
        % VALIDATE THE OBSTACLE
        function [tau]  = validateCollision(p,v)
            % CONFIRM WE KNOW THE HEADING OF THE OBSTACLE
            if any(isnan(v))
                tau = -inf;
                return
            end
            % Define the time to closest approach (+ve converging)
            tau = -(dot(p,v)/dot(v,v));
        end
        % BASIC VELOCITY VECTOR CHECK (2D & 3D)
        function [v_unit,v_mag] = nullVelocityCheck(v)
            v_mag = norm(v);
            if v_mag == 0
                v_unit = zeros(numel(v),1);
                v_unit(1) = 1;      % [1 0 0] % Default to current local forward direction
            else
                v_unit = v/v_mag;
            end
        end
    end

    %% ///////////////// VISUALISATION AND GENERAL TOOLS //////////////////
    methods
        % PLOT ALL THE OBSERVABLE OBSTACLES IN THE LOCAL FRAME
        function [figureHandle] = GetobjectScene(this,figureHandle)
            % FIGURE PREPERATION
            if ~exist('figureHandle','var')
                figureHandle = figure();
            else
                cla reset
            end

            ax = get(figureHandle,'CurrentAxes');
            set(ax,'NextPlot','replacechildren')
            hold on; grid on;
            axis equal;
            xlabel('X (m)'); ylabel('Y (m)'); zlabel('Z (m)');
            view([-120 35]);
            title(ax,'The agents perspective of its neighbourhood');

            % PLOT THE AGENT FOR PERSPECTIVE
            if size(this.GEOMETRY.vertices,1) < 1
                geometry = OMAS_geometry.defineSphere(zeros(3,1),this.radius);
                % REPRESENT GEOMETRY AS A PATCH
                entityHandle = patch(ax,...
                    'Vertices',geometry.vertices,...
                    'Faces',geometry.faces,...
                    'FaceColor',[0.4940, 0.1840, 0.5560]);
            else
                % PLOT THE AGENTS VERTEX DATA
                entityHandle = patch(ax,...
                    'Vertices',this.GEOMETRY.vertices,...
                    'Faces',this.GEOMETRY.faces,...
                    'FaceColor',this.GLOBAL.colour);
            end
            % ADD ANNOTATION
            annotationText = sprintf(' \t%s [ID-%d]',this.name,this.objectID);
            text(0,0,0,char(annotationText));

            % SET THE ENTITY DATA
            set(entityHandle,...
                'EdgeColor','k',...
                'EdgeAlpha',0.3,...
                'FaceLighting','gouraud',...
                'FaceAlpha',0.8,...
                'LineWidth',0.1,...
                'EdgeAlpha',0.7);

            % PLOT AGENT VELOCITY
            q = quiver3(0,0,0,1,0,0,'b');
            q.AutoScaleFactor = 1;
            % MOVE THROUGH THE OBSTACLE SET AND PLOT THE OBSTACLE POSITIONS
            for item = 1:numel(this.MEMORY)
                % Get the second objects data
                [p_j] = this.GetLastMeasurementByobjectID(this.MEMORY(item).objectID,'position');
                [v_j] = this.GetLastMeasurementByobjectID(this.MEMORY(item).objectID,'velocity');
                [r_j] = this.GetLastMeasurementByobjectID(this.MEMORY(item).objectID,'radius');

                if numel(this.MEMORY(item).geometry.vertices) < 1
                    % DEFINE GEOMETRY AS A SPHERE
                    if this.Is3D
                        [geometry] = OMAS_graphics.defineSphere(p_j,r_j);
                    else
                        [geometry] = OMAS_graphics.defineSphere([p_j;0],r_j);
                    end
                else
                    % DEFINE FROM OWN GEOMETRY
                    geometry = this.MEMORY(item).geometry;                 % Should be in the current frame
                end
                % REPRESENT GEOMETRY AS A PATCH
                entityHandle = patch(ax,...
                    'Vertices',geometry.vertices,...
                    'Faces',geometry.faces,...
                    'EdgeColor','k',...
                    'EdgeAlpha',0.2,...
                    'FaceLighting','gouraud',...
                    'FaceAlpha',0.2,...
                    'LineWidth',1);
                % PLOT REPRESENTATION
                switch this.MEMORY(item).type
                    case OMAS_objectType.agent
                        set(entityHandle,'FaceColor','b');
                    case OMAS_objectType.obstacle
                        set(entityHandle,'FaceColor','r');
                    case OMAS_objectType.waypoint
                        set(entityHandle,'FaceColor','g');
                    otherwise
                        set(entityHandle,'FaceColor','m');
                end
                % ADD ANNOTATION
                annotationText = sprintf(' \t%s [ID-%d]',this.MEMORY(item).name,this.MEMORY(item).objectID);
                if this.Is3D
                    text(p_j(1),p_j(2),p_j(3),char(annotationText));       % PLOT AGENT VELOCITY
                    q = quiver3(ax,p_j(1),p_j(2),p_j(3),...
                        v_j(1),v_j(2),v_j(3),'k');              % The velocity vector
                else
                    text(p_j(1),p_j(2),0,char(annotationText));
                    q = quiver3(ax,p_j(1),p_j(2),0,...
                        v_j(1),v_j(2),0,'k');                   % The velocity vector
                end
                q.AutoScaleFactor = 1;
            end
            drawnow;
            hold off;
        end
        % APPEND AN FIGURE FRAME TO ANIMATION
        function [figureHandle] = GetAnimationFrame(this,ENV,figureHandle,fileName)
            % INPUT HANDLING
            if nargin < 4
                fileName = sprintf('bodyAxes - %s [ID-%.0f].gif',this.name,this.objectID);
            end

            % FIGURE PREPERATION
            filePath = strcat(ENV.outputPath,fileName);
            drawnow;                                                       % Ensure all items are plotted before taking the image
            % CAPTURE FIGURE AS FRAME
            frame = getframe(figureHandle);
            im = frame2im(frame);                                          % Capture the figure as an image
            [imind,cm] = rgb2ind(im,256);
            % Write to the GIF File
            if ENV.currentStep == 1
                imwrite(imind,cm,filePath,'gif',...
                    'Loopcount',inf,...
                    'DelayTime',ENV.dt,...
                    'Location',[0 0],...
                    'Comment',sprintf('%s[ID-%.0f]',this.name,this.objectID));
            else
                imwrite(imind,cm,filePath,'gif',...
                    'WriteMode','append',...
                    'DelayTime',ENV.dt,...
                    'Location',[0 0]);
            end
            pause(0.1);                                                    % Pause for stability
        end
        % AGENT DATA STORAGE
        function [this] = writeAgentData(this,TIME,loopIndicator,loopDuration)
            % This function writes a regular DATA structure to the agent
            % .DATA structure to allow it to be easily interpretted.

            % THE AGENT COMPUTATION TIMES
            this.DATA.indicator(TIME.currentStep) = loopIndicator;
            this.DATA.dt(TIME.currentStep) = loopDuration;
            this.DATA.steps(TIME.currentStep) = TIME.currentStep;
            this.DATA.time(TIME.currentStep) = TIME.currentTime;
        end
    end
    %% ///////////////////////// GET/SET METHODS //////////////////////////
    methods
        % Set the maximum linear speed
        function set.v_max(this,v)
            assert(numel(v) == 1,'Expecting a scalar linear rate.');
            this.v_max = v;
            this.DYNAMICS.maxLinearVelocity(:) = v;
        end
        % Set the maximum angular speed
        function set.w_max(this,w)
            assert(numel(w) == 1,'Expecting a scalar angular rate.');
            this.w_max = w;
            this.DYNAMICS.maxAngularVelocity(:) = w;
        end
        % Set the detection radius
        function set.detectionRadius(this,r)
            assert(isnumeric(r) && numel(r) == 1,'Detection radius must be a scalar value.');
            % Override the SENSORS parameter to reflect the change
            if isstruct(this.SENSORS) && isfield(this.SENSORS,'range')
                this.SENSORS.range = r;                                    % Set the sensor range value
            end
            % Assign the value to the GLOBAL structure
            this.SetGLOBAL('detectionRadius',r);                           % Set the detection radius
        end
        % User update (Override)
        function [this] = ApplyUserOverrides(this,pairArray)
            % This function is designed to parse a generic set of user
            % inputs and allow them to be compared to a default input
            % structure. This should be called in the class constructor

            % Input sanity check #1
            if nargin < 2 || numel(pairArray) == 0
                return
            end

            % Call the all parameter overrider
            [this] = GetParameterOverrides_recursive(this,pairArray);

            % Preform dependant value checks
            this.MEMORY = this.CreateMEMORY(this.maxSamples);        % Initialise memory structure (with 2D varient)

            if this.detectionRadius ~= this.SENSORS.range
                % The detection radius is determined by SENSORS.range
                this.SENSORS.range = this.detectionRadius;
            end
        end
    end
end


text:
function [ n_data, a, x, fx ] = poisson_cdf_values ( n_data )

%*****************************************************************************80
%
%% POISSON_CDF_VALUES returns some values of the Poisson CDF.
%
%  Discussion:
%
%    CDF(X)(A) is the probability of at most X successes in unit time,
%    given that the expected mean number of successes is A.
%
%    In Mathematica, the function can be evaluated by:
%
%      Needs["Statistics`DiscreteDistributions`]
%      dist = PoissonDistribution [ a ]
%      CDF [ dist, x ]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%    Daniel Zwillinger,
%    CRC Standard Mathematical Tables and Formulae,
%    30th Edition, CRC Press, 1996, pages 653-658.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real A, the parameter of the function.
%
%    Output, integer X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 21;

  a_vec = [ ...
     0.02E+00, ...
     0.10E+00, ...
     0.10E+00, ...
     0.50E+00, ...
     0.50E+00, ...
     0.50E+00, ...
     1.00E+00, ...
     1.00E+00, ...
     1.00E+00, ...
     1.00E+00, ...
     2.00E+00, ...
     2.00E+00, ...
     2.00E+00, ...
     2.00E+00, ...
     5.00E+00, ...
     5.00E+00, ...
     5.00E+00, ...
     5.00E+00, ...
     5.00E+00, ...
     5.00E+00, ...
     5.00E+00 ];

  fx_vec = [ ...
     0.9801986733067553E+00, ...
     0.9048374180359596E+00, ...
     0.9953211598395555E+00, ...
     0.6065306597126334E+00, ...
     0.9097959895689501E+00, ...
     0.9856123220330293E+00, ...
     0.3678794411714423E+00, ...
     0.7357588823428846E+00, ...
     0.9196986029286058E+00, ...
     0.9810118431238462E+00, ...
     0.1353352832366127E+00, ...
     0.4060058497098381E+00, ...
     0.6766764161830635E+00, ...
     0.8571234604985470E+00, ...
     0.6737946999085467E-02, ...
     0.4042768199451280E-01, ...
     0.1246520194830811E+00, ...
     0.2650259152973617E+00, ...
     0.4404932850652124E+00, ...
     0.6159606548330631E+00, ...
     0.7621834629729387E+00 ];

  x_vec = [ ...
     0, 0, 1, 0, ...
     1, 2, 0, 1, ...
     2, 3, 0, 1, ...
     2, 3, 0, 1, ...
     2, 3, 4, 5, ...
     6 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    a = 0.0;
    x = 0;
    fx = 0.0;
  else
    a = a_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function [mssim, ssim_map] = ssim(img1, img2, K, window, L)

% ========================================================================
% SSIM Index with automatic downsampling, Version 1.0
% Copyright(c) 2009 Zhou Wang
% All Rights Reserved.
%
% ----------------------------------------------------------------------
% Permission to use, copy, or modify this software and its documentation
% for educational and research purposes only and without fee is hereby
% granted, provided that this copyright notice and the original authors'
% names appear on all copies and supporting documentation. This program
% shall not be used, rewritten, or adapted as the basis of a commercial
% software or hardware product without first obtaining permission of the
% authors. The authors make no representations about the suitability of
% this software for any purpose. It is provided "as is" without express
% or implied warranty.
%----------------------------------------------------------------------
%
% This is an implementation of the algorithm for calculating the
% Structural SIMilarity (SSIM) index between two images
%
% Please refer to the following paper and the website with suggested usage
%
% Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, "Image
% quality assessment: From error visibility to structural similarity,"
% IEEE Transactios on Image Processing, vol. 13, no. 4, pp. 600-612,
% Apr. 2004.
%
% http://www.ece.uwaterloo.ca/~z70wang/research/ssim/
%
% Note: This program is different from ssim_index.m, where no automatic
% downsampling is performed. (downsampling was done in the above paper
% and was described as suggested usage in the above website.)
%
% Kindly report any suggestions or corrections to zhouwang@ieee.org
%
%----------------------------------------------------------------------
%
%Input : (1) img1: the first image being compared
%        (2) img2: the second image being compared
%        (3) K: constants in the SSIM index formula (see the above
%            reference). defualt value: K = [0.01 0.03]
%        (4) window: local window for statistics (see the above
%            reference). default widnow is Gaussian given by
%            window = fspecial('gaussian', 11, 1.5);
%        (5) L: dynamic range of the images. default: L = 255
%
%Output: (1) mssim: the mean SSIM index value between 2 images.
%            If one of the images being compared is regarded as
%            perfect quality, then mssim can be considered as the
%            quality measure of the other image.
%            If img1 = img2, then mssim = 1.
%        (2) ssim_map: the SSIM index map of the test image. The map
%            has a smaller size than the input images. The actual size
%            depends on the window size and the downsampling factor.
%
%Basic Usage:
%   Given 2 test images img1 and img2, whose dynamic range is 0-255
%
%   [mssim, ssim_map] = ssim(img1, img2);
%
%Advanced Usage:
%   User defined parameters. For example
%
%   K = [0.05 0.05];
%   window = ones(8);
%   L = 100;
%   [mssim, ssim_map] = ssim(img1, img2, K, window, L);
%
%Visualize the results:
%
%   mssim                        %Gives the mssim value
%   imshow(max(0, ssim_map).^4)  %Shows the SSIM index map
%========================================================================

if (nargin < 2 || nargin > 5)
   mssim = -Inf;
   ssim_map = -Inf;
   return;
end

if (size(img1) ~= size(img2))
   mssim = -Inf;
   ssim_map = -Inf;
   return;
end

[M N] = size(img1);

if (nargin == 2)
   if ((M < 11) || (N < 11))
	   mssim = -Inf;
	   ssim_map = -Inf;
      return
   end
   %window = fspecial('gaussian', 11, 1.5);	%
   window = Gaussian_window(11,1.5);
   K(1) = 0.01;					% default settings
   K(2) = 0.03;					%
   L = 255;                                     %
end

if (nargin == 3)
   if ((M < 11) || (N < 11))
	   mssim = -Inf;
	   ssim_map = -Inf;
      return
   end
   %window = fspecial('gaussian', 11, 1.5);
   window = Gaussian_window(11,1.5);
   L = 255;
   if (length(K) == 2)
      if (K(1) < 0 || K(2) < 0)
		   mssim = -Inf;
   		ssim_map = -Inf;
	   	return;
      end
   else
	   mssim = -Inf;
   	ssim_map = -Inf;
	   return;
   end
end

if (nargin == 4)
   [H W] = size(window);
   if ((H*W) < 4 || (H > M) || (W > N))
	   mssim = -Inf;
	   ssim_map = -Inf;
      return
   end
   L = 255;
   if (length(K) == 2)
      if (K(1) < 0 || K(2) < 0)
		   mssim = -Inf;
   		ssim_map = -Inf;
	   	return;
      end
   else
	   mssim = -Inf;
   	ssim_map = -Inf;
	   return;
   end
end

if (nargin == 5)
   [H W] = size(window);
   if ((H*W) < 4 || (H > M) || (W > N))
	   mssim = -Inf;
	   ssim_map = -Inf;
      return
   end
   if (length(K) == 2)
      if (K(1) < 0 || K(2) < 0)
		   mssim = -Inf;
   		ssim_map = -Inf;
	   	return;
      end
   else
	   mssim = -Inf;
   	ssim_map = -Inf;
	   return;
   end
end

img1 = double(img1);
img2 = double(img2);

% automatic downsampling
f = max(1,round(min(M,N)/256));
%downsampling by f
%use a simple low-pass filter
if(f>1)
    lpf = ones(f,f);
    lpf = lpf/sum(lpf(:));
    img1 = imfilter(img1,lpf,'symmetric','same');
    img2 = imfilter(img2,lpf,'symmetric','same');

    img1 = img1(1:f:end,1:f:end);
    img2 = img2(1:f:end,1:f:end);
end

C1 = (K(1)*L)^2;
C2 = (K(2)*L)^2;
window = window/sum(sum(window));

mu1   = filter2(window, img1, 'valid');
mu2   = filter2(window, img2, 'valid');
mu1_sq = mu1.*mu1;
mu2_sq = mu2.*mu2;
mu1_mu2 = mu1.*mu2;
sigma1_sq = filter2(window, img1.*img1, 'valid') - mu1_sq;
sigma2_sq = filter2(window, img2.*img2, 'valid') - mu2_sq;
sigma12 = filter2(window, img1.*img2, 'valid') - mu1_mu2;

if (C1 > 0 && C2 > 0)
   ssim_map = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))./((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2));
else
   numerator1 = 2*mu1_mu2 + C1;
   numerator2 = 2*sigma12 + C2;
	denominator1 = mu1_sq + mu2_sq + C1;
   denominator2 = sigma1_sq + sigma2_sq + C2;
   ssim_map = ones(size(mu1));
   index = (denominator1.*denominator2 > 0);
   ssim_map(index) = (numerator1(index).*numerator2(index))./(denominator1(index).*denominator2(index));
   index = (denominator1 ~= 0) & (denominator2 == 0);
   ssim_map(index) = numerator1(index)./denominator1(index);
end

mssim = mean(ssim_map(:));

return

author: andrewssobral


text:
%% Economic Parameters

% Preferences
coefrra = 1;
frisch = 0.5;
adjfricshgridfrac = 0.85;
meanlabeff = 3; 	%so that at h=1/3 outout will be approximately = 1;
labdisutil = meanlabeff ./ ( (0.75 .^(-coefrra)) * ((1.0/3.0).^(1.0/frisch))); %guess labdisutil so that at average wages and average consumption hours =1/3 (sets C/Y = 0.75);
maxhours = 1;

% Production
ceselast = 10;							 % elasticity of substitution / demand
priceadjust = 100;

% Policy parameters
taylor_inflation = 1.25;	%taylor rule coefficient on inflation
taylor_outputgap = 0;		%taylor rule coefficient on output
labtax = 0.2;			% marginal tax rate
govbondtarget = 6;		%multiple of quarterly GDP
lumptransferpc = 0.06;	 %6% of quarterly GDP in steady state
govbcrule_fixnomB = 0.0;

% Some aggregates
Y_SS = 1;
B_SS = govbondtarget .* Y_SS;
m_SS = (ceselast-1)/ceselast;
w_SS = m_SS;
N_SS = 1/3;						    % steady state hours: so that quartelry GDP = 1 in s.s
profit_SS = (1 - m_SS) * Y_SS;
lumptransfer = lumptransferpc .* Y_SS;

% Aggregate shocks
ssigma_MP = sqrt(0.05);
ttheta_MP = 0.25;

%% Approximation Parameters
% HJB
maxit_HJB = 500;
tol_HJB  = 1.0e-8;
Delta_HJB = 1.0e6;

maxit_KFE = 1000;
tol_KFE  = 1.0e-12;
Delta_KFE = 1.0e6;

crit  = 10^(-6);

% steady state r
r0 = 0.005;
rmax = 0.08;
rmin = 0.001;
Ir = 100;

% steady state rho iteration
rrho0 = 0.02;
rhomax = 0.05;
rhomin = 0.005;
crit_S = 10^(-5);

% hours loop
niter_hours = 10;

%% Set Grids

%% Asset
I          = 100;
agridparam = 1; % 1 for linear
amin       = 0;
amax       = 40;

a          = linspace(0,1,I)';
a          = a.^(1/agridparam);
a          = amin + (amax-amin).*a;

%% Income
J = 2;

ygrid_combined = [0.2; 1];
ymarkov_combined = [-0.5 0.5; 0.0376 -0.0376];

% Compute stationary income distribution
AT = ymarkov_combined';
g_z = ones(J,1)/J;
for n = 1 : 50
	g_z_new = (speye(J) - AT * 1000) \ g_z;
	diff = max(abs(g_z_new - g_z));
	if diff < 10e-6
		break
	end
	g_z = g_z_new;
end

z = exp(ygrid_combined);

% scale so that mean = meanlabeff
z  	  = meanlabeff.*z./sum(z .* g_z);
z_bar = sum(z .* g_z);
z     = z';
zz = ones(I,1)*z;

%% Kronecker to matrix
aa = a*ones(1,J);

daf = [a(2:I)-a(1:I-1); a(I)-a(I-1)];
dazf = repmat(daf,1,J);
dab = [a(2)-a(1); a(2:I)-a(1:I-1)];
dazb = repmat(dab,1,J);

adelta = zeros(size(a));
adelta(1) = 0.5.*daf(1);
adelta(2:I-1) = 0.5.*(daf(1:I-2)+daf(2:I-1));
adelta(I) = 0.5*daf(I-1);

azdelta = repmat(adelta,1,J);
azdelta_mat = spdiags(azdelta(:),0,I*J,I*J);

Aswitch = kron(ymarkov_combined,speye(I));

%% Others
IterateR = 0;
IterateRho = 1;

author: gregkaplan

repo: phact

save_path: github-repos/MATLAB/gregkaplan-phact


text:
function LMthumbnailsbar(folderlist, filelist, pagename, HOMEIMAGES);
%
% LMthumbnailsbar(folderlist, filelist, pagename, HOMEIMAGES);
%
% HOMEIMAGES should be the path for the images in your local copy of the
% database.
%
% This function allows building an interface that communicates with the
% LabelMe web annotation tool.
%
% It can be used to label specific images. It adds a bar with thumbnails at
% the left side of the annotation tool.
%
% You can use this to label images that have some characteristic that you
% want. You can use this function in combination with the LMquery function.
%
% For instance, if you want to create a web page with images only of
% kitchens so that the thumbnails are connected to the LabelMe web
% annotation tool online, you can do the next thing:
%
% [D,j] = LMquery(database, 'folder', 'kitchen');
% for i = 1:length(D);
%    folderlist{i} = D(i).annotation.folder;
%    filelist{i} = D(i).annotation.filename;
% end
% LMthumbnailsbar(folderlist, filelist, 'myphotoalbum.html', HOMEIMAGES);
%
%
% See also LMphotoalbum
%

Nimages = length(folderlist);
webpage = 'http://labelme.csail.mit.edu/tool.html?collection=LabelMe&mode=i'
thumbpage = 'http://people.csail.mit.edu/brussell/research/LabelMe/Thumbnails/'

% Header
page = {};

% Create links for each image
for i = 1:Nimages
    imageline = sprintf('', thumbpage, folderlist{i}, filelist{i});
    page = addline(page, sprintf('%d
',i));
    page = addline(page, ...
        [sprintf('%s', ...
        webpage, folderlist{i}, filelist{i}, imageline)]);
end

% Close web page
page = addline(page, '');

% write web page
fid = fopen(['bar' pagename],'w');
for i = 1:length(page);
    fprintf(fid, [page{i} '\n']);
end
fclose(fid);

% Header
page = {};
page = addline(page, '<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Frameset//EN"');
page = addline(page, '"http://www.w3.org/TR/html4/frameset.dtd">');
page = addline(page, '<title>LabelMe thumbnail bar</title>');
page = addline(page, '<FRAMESET cols="154, *">');
page = addline(page, sprintf('<FRAME src="%s">', ['bar' pagename]));
page = addline(page, sprintf('<FRAME src = "%s&folder=%s&image=%s" name="labelme">', webpage, folderlist{1}, filelist{1}));
page = addline(page, '</FRAMESET>');

% write web page
fid = fopen(pagename,'w');
for i = 1:length(page);
    fprintf(fid, [page{i} '\n']);
end
fclose(fid);

function page = addline(page, line)

page = [page {line}];

author: CSAILVision

repo: LabelMeToolbox


text:
function T = logm(q,varargin)
% the logarithmic map that translates a rotation into a spin tensor
%
% Syntax
%   T = logm(q) % spin tensor with reference to the identical rotation
%   T = logm(q,q_ref) % spin tensor with reference q_ref
%
% Input
%  q - @quaternion
%  q_ref - @quaternion
%
% Output
%  T - @spinTensor
%
% See also
% spinTensor/exp

tq = log(q,varargin{:});

M = zeros([3,3,size(q)]);

M(2,1,:) =  tq.z;
M(3,1,:) = -tq.y;
M(3,2,:) =  tq.x;

M(1,2,:) = -tq.z;
M(1,3,:) =  tq.y;
M(2,3,:) = -tq.x;

% make it a spinTensor
T = spinTensor(M);

author: mtex-toolbox

repo: mtex

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


text:
function M=ImpliedExpRets(S,w)

M_=S*w;

s=sqrt(mean(diag(S)));
M=M_/mean(M_)*.5*s;

author: Sable

repo: mcbench-benchmarks


text:
function value = f2 ( x )

%*****************************************************************************80
%
%% F2 evaluates a function that can be used for Chebyshev interpolation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 September 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(), one or more points where the function is to be evaluated.
%
%    Output, real VALUE(), the function values.
%
  value = exp ( x );

  return
end


text:
function [C_df,Df,S_df] = extract_DF_F(Y,A,C,S,i,options)

% extract DF/F signals after performing NMF
% inputs:  Y raw data (d X T matrix, d # number of pixels, T # of timesteps)
%          A matrix of spatial components (d x K matrix, K # of components)
%          C matrix of temporal components (K x T matrix)
%          S matrix of deconvolved activity ((K-1) x T matrix) (optional)
%          i index of component that represent the background (optional, if not
%          given it's estimated)
%          options structure used for specifying method for determining DF
%           default method is the median of the trace. By changing
%           options.df_prctile an arbitray percentile can be used (between 0 and 100).
%           a moving window can also be established by specifying options.df_window

% outputs:  C_df temporal components in the DF/F domain
%           Df   background for each component to normalize the filtered raw data
%           S_df deconvolved activity/spikes in the DF/F domain

% Written by:
% Eftychios A. Pnevmatikakis, Simons Foundation, 2015

defoptions = CNMFSetParms;
if nargin < 6 || isempty(options)
    options = defoptions;
end
if ~isfield(options,'df_prctile') || isempty(options.df_prctile)
    options.df_prctile = defoptions.df_prctile;
end
if ~isfield(options,'df_window') || isempty(options.df_window)
    options.df_window = defoptions.df_window;
end

nA = sqrt(sum(A.^2))';
[K,~] = size(C);
A = A/spdiags(nA,0,K,K);    % normalize spatial components to unit energy
C = spdiags(nA,0,K,K)*C;

if nargin < 5 || isempty(i)
    [~,i] = min(sum(A.^6)); % identify background component
end

non_bg = true(1,K);
non_bg(i) = false;      % non-background components
Yf = A'*Y - (A'*A(:,non_bg))*C(non_bg,:);

if isempty(options.df_window) || (options.df_window > size(C,2))
    if options.df_prctile == 50
        Df = median(Yf,2);
    else
        Df = prctile(Yf,options.df_prctile,2);
    end
    C_df = spdiags(Df,0,K,K)\C;
else
    if options.df_prctile == 50
        Df = medfilt1(Yf,options.df_window,[],2,'truncate');
    else
        Df = zeros(size(Yf));
        for i = 1:size(Df,1);
            df_temp = running_percentile(Yf(i,:), options.df_window, options.df_prctile);
            Df(i,:) = df_temp(:)';
        end
    end
    C_df = C./Df;
end

C_df(i,:) = 0;

if nargin < 4 || isempty(S) || nargout < 3
    S_df = [];
    if nargout == 3
        warning('Merged spikes matrix is returned as empty because the original matrix was not provided.');
    end
else
    if isempty(options.df_window) || (options.df_window > size(C,2))
        S_df = spdiags(Df(non_bg(:)),0,sum(non_bg),sum(non_bg))\S;
    else
        S_df = S./Df(non_bg,:);
    end
end

author: zhoupc


text:
function [ f, d, s, t ] = cubic_value ( x )

%*****************************************************************************80
%
%% CUBIC_VALUE evaluates a cubic function.
%
%  Discussion:
%
%    f(x) =   x^3 -  7 x^2 + 10 x
%    d(x) = 3 x^2 - 14 x   + 10
%    s(x) = 6 x   - 14
%    t(x) = 6
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real F, D, S, T, the value and first three derivatives of
%    the cubic function.
%

%
%  If X is a vector, force it to be a column vector.
%
  x = x ( : );

  f = 0.0 + x .* ( 10.0 + x .* (  - 7.0 + x * 1.0 ) );
  d =              10.0 + x .* ( - 14.0 + x * 3.0 );
  s =                            - 14.0 + x * 6.0;
  t =                                         6.0;

  return
end