text:
%CLEVALB Classifier evaluation (learning curve), bootstrap version
%
% 	E = CLEVALB(A,CLASSF,TRAINSIZES,NREPS)
%
% INPUT
%   A           Training dataset
%   CLASSF      Classifier to evaluate
%   TRAINSIZES  Vector of class sizes, used to generate subsets of A
%               (default [2,3,5,7,10,15,20,30,50,70,100])
%   NREPS       Number of repetitions (default 1)
%
% OUTPUT
%   E           Error structure (see PLOTE)
%
% DESCRIPTION
% Generates at random, for all class sizes defined in TRAINSIZES, training
% sets out of the dataset A and uses these for training the untrained
% classifier CLASSF. CLASSF may also be a cell array of untrained
% classifiers; in this case the routine will be run for all of them. The
% resulting trained classifiers are tested on all objects in A. This
% procedure is then repeated N times.
%
% Training set generation is done "with replacement" and such that for each
% run the larger training sets include the smaller ones and that for all
% classifiers the same training sets are used.
%
% If CLASSF is fully deterministic, this function uses the RAND random
% generator and thereby reproduces if its seed is reset (see RAND).
% If CLASSF uses RANDN, its seed may have to be set as well.
%
% Use FID = 1 to report progress to the command window.
%
% EXAMPLES
% See PREX_CLEVAL.
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, CLEVALB, TESTC, PLOTE

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: clevalb.m,v 1.4 2008/07/03 09:05:50 duin Exp $

function e = clevalb(a,classf,learnsizes,nreps,fid)

  % use of fid is outdated
  if (nargin < 4) || isempty(nreps)
		prwarning(2,'number of repetitions not specified, assuming NREPS = 1');
		nreps = 1;
	end;
  if (nargin < 3) || isempty(learnsizes)
		prwarning(2,'vector of training set class sizes not specified, assuming [2,3,5,7,10,15,20,30,50,70,100]');
 		learnsizes = [2,3,5,7,10,15,20,30,50,70,100];
	end;

	% If a single mapping is given, convert it to a 1 x 1 cell array.

  if (ismapping(classf)), classf = {classf}; end

	% Correct for old argument order.

  if (isdataset(classf)) & (ismapping(a))
  	tmp = a; a = classf; classf = {tmp};
  end
  if (isdataset(classf)) & (iscell(a)) & (ismapping(a{1}))
  	tmp = a; a = classf; classf = tmp;
  end
  if ~iscell(classf), classf = {classf}; end

	% Assert that all is right.
  isdataset(a); ismapping(classf{1});

	% Remove requested class sizes that are larger than the size of the
	% smallest class.

  mc = classsizes(a); [m,k,c] = getsize(a);
	toolarge = find(learnsizes >= min(mc));
	if (~isempty(toolarge))
		prwarning(2,['training set class sizes ' num2str(learnsizes(toolarge)) ...
								 ' larger than the minimal class size in A; removed them']);
	  learnsizes(toolarge) = [];
	end
  learnsizes = learnsizes(:)';

	% Fill the error structure.

  nw = length(classf(:));
  datname = getname(a);

  e.n       = nreps;
  e.error   = zeros(nw,length(learnsizes));
  e.std     = zeros(nw,length(learnsizes));
  e.xvalues = learnsizes(:)';
  e.names   = [];
  e.xlabel  = 'Training set size';
  if (nreps > 1)
  	e.ylabel= ['Averaged error (' num2str(nreps) ' experiments)'];
  elseif (nreps == 1)
  	e.ylabel = 'Error';
  else
		error('Number of repetitions NREPS should be >= 1.');
	end;
  if (~isempty(datname))
  	e.title = ['Bootstrapped learning curve on ' datname];
  end
  if (learnsizes(end)/learnsizes(1) > 20)
  	e.plot = 'semilogx'; 				% If range too large, use a log-plot for X.
  end

	% Report progress.

	s1 = sprintf('clevalb: %i classifiers: ',nw);
	prwaitbar(nw,s1);

	% Store the seed, to reset the random generator later for different
	% classifiers.

  seed = randreset;

	% Loop over all classifiers (with index WI).

  for wi = 1:nw

  	isuntrained(classf{wi});
    name = getname(classf{wi});
    e.names = char(e.names,name);
    prwaitbar(nw,wi,[s1 name]);

		% E1 will contain the error estimates.

  	e1 = zeros(nreps,length(learnsizes));

		% Take care that classifiers use same training set.

  	randreset(seed); seed2 = seed;

		% For NREPS repetitions...
		s2 = sprintf('cleval: %i repetitions: ',nreps);
		prwaitbar(nreps,s2);

  	for i = 1:nreps

			prwaitbar(nreps,i,[s2 int2str(i)]);
			% Store the randomly permuted indices of samples of class CI to use in
			% this training set in JR(CI,:).

  		JR = zeros(c,max(learnsizes));
  		for ci = 1:c
				JC = findnlab(a,ci);

				% Necessary for reproducable training sets: set the seed and store
				% it after generation, so that next time we will use the previous one.
  			randreset(seed2);

  			R = ceil(rand(1,max(learnsizes))*length(JC));
  			JR(ci,:) = JC(R)';

  			seed2 = randreset;
  		end

  		li = 0;										% Index of training set.
			nlearns = length(learnsizes);
			s3 = sprintf('cleval: %i sizes: ',nlearns);
			prwaitbar(nreps,s3);
  		for j = 1:nlearns
				nj = learnsizes(j);
        prwaitbar(nlearns,j,[s3 int2str(j) ' (' int2str(nj) ')']);
        li = li + 1;

				% J will contain the indices for this training set.

				J = [];
  			for ci = 1:c
  				J = [J;JR(ci,1:nj)'];
  			end;

				% Train classifier CLASSF{WI} on this training set and calculate
				% error.

  			W = a(J,:)*classf{wi};
    		e1(i,li) = testc(a,W);

  		end
      prwaitbar(0);

  	end
    prwaitbar(0);

		% Calculate average error and standard deviation for this classifier
		% (or set the latter to zero if there's been just 1 repetition).

  	e.error(wi,:) = mean(e1,1);
  	if (nreps == 1)
  		e.std(wi,:) = zeros(1,size(e.std,2));
  	else
  		e.std(wi,:) = std(e1)/sqrt(nreps);
  	end

  end
  prwaitbar(0);

  % The first element is the empty string [], remove it.
  e.names(1,:) = [];

return

author: marianux

repo: ecg-kit

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


text:
function [U S V info] = truncated_svd(A, p)
% Returns an SVD decomposition of A truncated to rank p.
%
% function [U S V] = truncated_svd(A, p)
%
% Input: A real matrix A of size mxn and an integer p <= min(m, n).
% Output: An orthonormal matrix U of size mxp, an orthonormal matrix Y of
%         size nxp and a diagonal matrix S of size pxp with nonnegative and
%         decreasing diagonal entries such that USV.' is the best rank p
%         approximation of A according to the Frobenius norm. All real.
%         This function produces an output akin to svds.
%
% The decomposition is obtained by maximizing
%   f(U, V) = .5*norm(U'*A*V, 'fro')^2
% where U, V are orthonormal. Notice that f(U*Q, V*R) = f(U, V) for all
% Q, R orthogonal pxp matrices. Hence, only the column spaces of U and V
% matter and we may perform the optimization over a product of two
% Grassmannian manifolds.
%
% It is easy to show that maximizing f is equivalent to minimizing g with
%   g(U, V) = min_S norm(U*S*V' - A, 'fro')^2,
% which confirms that we are going for a best low-rank approximation of A.
%
% The inner workings of the Grassmann manifold use the built-in svd
% function of Matlab but only for matrices of size mxp and nxp to
% re-orthonormalize them.
%
% Notice that we are actually chasing a best fixed-rank approximation of a
% matrix, which is best obtained by working directly over a manifold of
% fixed-rank matrices. This is simply an example script to demonstrate some
% functionalities of the toolbox.
%
% The code can be modified to accept a function handle for A(x) = A*x
% instead of a matrix A, which is often useful. This would further require
% a function handle At for the transpose of A, such that At(x) = A.'*x.

% This file is part of Manopt and is copyrighted. See the license file.
%
% Main author: Nicolas Boumal, July 5, 2013
% Contributors:
%
% Change log:
%

    % Generate some random data to test the function if none is given.
    if ~exist('A', 'var') || isempty(A)
        A = randn(42, 60);
    end
    if ~exist('p', 'var') || isempty(p)
        p = 5;
    end

    % Retrieve the size of the problem and make sure the requested
    % approximation rank is at most the maximum possible rank.
    [m n] = size(A);
    assert(p <= min(m, n), 'p must be smaller than the smallest dimension of A.');

    % Define the cost and its derivatives on the Grassmann manifold
    tuple.U = grassmannfactory(m, p);
    tuple.V = grassmannfactory(n, p);
    % All of the code will work just as well if we ignore the invariance
    % property of the cost function indicated above and thus place U and V
    % on the Stiefel manifold (orthonormal matrices) instead of the
    % Grassmann manifold. Working on Stiefel is expected to be slower
    % though, partly because de search space is higher dimensional and
    % partly because the optimizers are not isolated.
    % tuple.U = stiefelfactory(m, p);
    % tuple.V = stiefelfactory(n, p);
    M = productmanifold(tuple);

    % Define a problem structure, specifying the manifold M, the cost
    % function and its derivatives. Here, to demonstrate the rapid
    % prototyping capabilities of Manopt, we directly define the Euclidean
    % gradient and the Euclidean Hessian egrad and ehess instead of the
    % Riemannian gradient and Hessian grad and hess. Manopt will take care
    % of the conversion. This automatic conversion is usually not
    % computationally optimal though, because much of the computations
    % involved in obtaining the gradient could be reused to obtain the
    % Hessian. After the prototyping stage, when efficiency becomes
    % important, it makes sense to define grad and hess rather than egrad
    % an ehess, and to use the caching system (the store structure).
    problem.M = M;
    problem.cost  = @cost;
    problem.egrad = @egrad;
    problem.ehess = @ehess;

    % The functions below make many redundant computations. This
    % performance hit can be alleviated by using the caching system.

    % Cost function
    function f = cost(X)
        U = X.U;
        V = X.V;
        f = -.5*norm(U'*A*V, 'fro')^2;
    end
    % Euclidean gradient of the cost function
    function g = egrad(X)
        U = X.U;
        V = X.V;
        AV = A*V;
        AtU = A'*U;
        g.U = -AV*(AV'*U);
        g.V = -AtU*(AtU'*V);
    end
    % Euclidean Hessian of the cost function
    function h = ehess(X, H)
        U = X.U;
        V = X.V;
        Udot = H.U;
        Vdot = H.V;
        AV = A*V;
        AtU = A'*U;
        AVdot = A*Vdot;
        AtUdot = A'*Udot;
        h.U = -(AVdot*AV'*U + AV*AVdot'*U + AV*AV'*Udot);
        h.V = -(AtUdot*AtU'*V + AtU*AtUdot'*V + AtU*AtU'*Vdot);
    end

    % Execute some checks on the derivatives for early debugging.
    % These things can be commented out of course.
    checkgradient(problem);
    pause;
    checkhessian(problem);
    pause;

    % Issue a call to a solver. A random initial guess will be chosen and
    % default options are selected. Here, we specify a maximum trust
    % region radius (which in turn induces an initial trust region radius).
    % Note that this is not required: default values are used if we omit
    % this. The diameter of the manifold scales like sqrt(2*p), hence the
    % form of our (empirical) choice.
    options.Delta_bar = 4*sqrt(2*p);
    [X Xcost info] = trustregions(problem, [], options); %#ok<ASGLU>
    U = X.U;
    V = X.V;

    % Finish the job by rotating U and V such that the middle matrix S can
    % be diagonal with nonnegative, decreasing entries. This requires a
    % small svd of size pxp.
    Spp = U'*A*V;
    [Upp Spp Vpp] = svd(Spp);
    U = U*Upp;
    S = Spp;
    V = V*Vpp;

    % For our information, Manopt can also compute the spectrum of the
    % Riemannian Hessian on the tangent space at (any) X. Computing the
    % spectrum at the solution gives us some idea of the conditioning of
    % the problem. If we were to implement a preconditioner for the
    % Hessian, this would also inform us on its performance.
    %
    % Notice that if the optimization is performed on a product of Stiefel
    % manifolds instead of a product of Grassmannians, the double
    % invariance under the orthogonal group O(p) will appear as twice
    % p*(p-1)/2, thus p*(p-1) zero eigenvalues in the spectrum of the
    % Hessian. This means that the minimizers are not isolated, which
    % typically hinders convergence of second order algorithms.
    if M.dim() < 512
        evs = hessianspectrum(problem, X);
        stairs(sort(evs));
        title(['Eigenvalues of the Hessian of the cost function ' ...
               'at the solution']);
    end

end

author: andrewssobral

repo: lrslibrary


text:
classdef prtUiManagerPanel < prtUiManager

    properties (Dependent)
        title
        fontSize
    end

    methods
        function self = prtUiManagerPanel(varargin)
           self = self@prtUiManager(varargin{:});
        end

        function create(self)
            self.managedHandle = uipanel('BackgroundColor',get(0,'DefaultFigureColor'),...
                'BorderType','none');
        end

        function set.title(self,str)
            set(self.managedHandle,'title',str)
        end
        function val = get.title(self)
            val = get(self.managedHandle,'title');
        end

        function val = get.fontSize(self)
            val = get(self.managedHandle,'fontSize');
        end
        function set.fontSize(self,val)
            set(self.managedHandle, 'fontsize', val);
        end

    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
% colony_count.m demo script
% Runs colony_count on example image and displays segmentation results.

% Load image.
I = imread('plate_example.jpg');

% Call function.
% Count is number of colonies, markers is the colony markers and mask is
% the colony segmentation. Count == # of components in markers.

% Optional parameter, this is the default.
% Should encompass internal radius of Petri dish.
radii = 115:1:130;

[count,markers,mask] = colony_count(I,radii);

disp([num2str(count) ' colonies detected.']);

% Overlay image with markers and mask.
J = rgb2gray(I);
R = J; G = J; B = J;
R(mask) = 255; R(markers) = 0;
G(mask) = 0; G(markers) = 255;
B(mask) = 0; B(markers) = 0;
RGB = cat(3,R,G,B);
figure();imshow(RGB);

% Overlay image with just markers.
R = J; G = J; B = J;
R(markers) = 0;
G(markers) = 255;
B(markers) = 0;
RGB = cat(3,R,G,B);
figure();imshow(RGB);

author: Sable

repo: mcbench-benchmarks


text:
function t_om_solve_qps(quiet)
%T_OM_SOLVE_QPS  Tests of QP solvers via OM.SOLVE().

%   MP-Opt-Model
%   Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

if nargin < 1
    quiet = 0;
end

algs = {'DEFAULT', 'BPMPD', 'MIPS', 250, 'IPOPT', 'OT', 'CPLEX', 'MOSEK', 'GUROBI', 'CLP', 'GLPK', 'OSQP'};
names = {'DEFAULT', 'BPMPD_MEX', 'MIPS', 'sc-MIPS', 'IPOPT', 'linprog/quadprog', 'CPLEX', 'MOSEK', 'Gurobi', 'CLP', 'glpk', 'OSQP'};
check = {[], 'bpmpd', [], [], 'ipopt', 'quadprog', 'cplex', 'mosek', 'gurobi', 'clp', 'glpk', 'osqp'};
does_qp = [1 1 1 1 1 1 1 1 1 1 0 1];

n = 30;
nqp = 21;
t_begin(27+n*length(algs), quiet);

diff_alg_warn_id = 'optim:linprog:WillRunDiffAlg';
if have_feature('quadprog') && have_feature('quadprog', 'vnum') == 7.005
    s1 = warning('query', diff_alg_warn_id);
    warning('off', diff_alg_warn_id);
end

for k = 1:length(algs)
    if ~isempty(check{k}) && ~have_feature(check{k})
        t_skip(n, sprintf('%s not installed', names{k}));
    else
        opt = struct('verbose', 0, 'alg', algs{k});
        mpopt = struct( ...
            'verbose', 0, ...
            'opf', struct( ...
                'violation', 5e-6 ), ...
            'cplex', struct( ...
                'lpmethod', 0, ...
                'qpmethod', 0, ...
                'opt', 0 ), ...
            'mosek', struct( ...
                'lp_alg', 0, ...
                'max_it', 0, ...
                'gap_tol', 0, ...
                'max_time', 0, ...
                'num_threads', 0, ...
                'opt', 0 ) ...
        );
        switch names{k}
            case {'DEFAULT', 'MIPS', 'sc-MIPS'}
                opt.mips_opt.comptol = 1e-8;
%             case 'linprog/quadprog'
%                 opt.verbose = 2;
%                 opt.linprog_opt.Algorithm = 'interior-point';
%                 opt.linprog_opt.Algorithm = 'active-set';
%                 opt.linprog_opt.Algorithm = 'simplex';
%                 opt.linprog_opt.Algorithm = 'dual-simplex';
%             end
            case 'CPLEX'
%               alg = 0;        %% default uses barrier method with NaN bug in lower lim multipliers
                alg = 2;        %% use dual simplex
                mpopt.cplex.lpmethod = alg;
                mpopt.cplex.qpmethod = min([4 alg]);
                opt.cplex_opt = cplex_options([], mpopt);
            case 'MOSEK'
%                 sc = mosek_symbcon;
%                 alg = sc.MSK_OPTIMIZER_DUAL_SIMPLEX;    %% use dual simplex
%                 alg = sc.MSK_OPTIMIZER_INTPNT;          %% use interior point
%                 mpopt.mosek.lp_alg = alg;
                mpopt.mosek.gap_tol = 1e-10;
%                 mpopt.mosek.opts.MSK_DPAR_INTPNT_TOL_PFEAS = 1e-10;
%                 mpopt.mosek.opts.MSK_DPAR_INTPNT_TOL_DFEAS = 1e-10;
%                 mpopt.mosek.opts.MSK_DPAR_INTPNT_TOL_INFEAS = 1e-10;
%                 mpopt.mosek.opts.MSK_DPAR_INTPNT_TOL_REL_GAP = 1e-10;
                vnum = have_feature('mosek', 'vnum');
                if vnum >= 8
%                     mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_PFEAS = 1e-10;
%                     mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_DFEAS = 1e-10;
%                     mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_INFEAS = 1e-10;
%                     mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_MU_RED = 1e-10;
                    mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_REL_GAP = 1e-10;
                end
%                 opt.verbose = 3;
                opt.mosek_opt = mosek_options([], mpopt);
            case 'OSQP'
                opt.osqp_opt.polish = 1;
%                 opt.osqp_opt.alpha = 1;
%                 opt.osqp_opt.eps_abs = 1e-8;
%                 opt.osqp_opt.eps_rel = 1e-10;
%                 opt.osqp_opt.eps_prim_inf = 1e-8;
%                 opt.osqp_opt.eps_dual_inf = 1e-8;
        end

        t = sprintf('%s - 3-d LP : ', names{k});
        %% based on example from 'doc linprog'
        c = [-5; -4; -6];
        A = [1 -1  1;
             -3  -2  -4;
             3  2  0];
        l = [-Inf; -42; -Inf];
        u = [20; Inf; 30];
        xmin = [0; 0; 0];
        x0 = [];
        om = opt_model;
        om.add_var('x', 3, x0, xmin);
        om.add_quad_cost('c', [], c);
        om.add_lin_constraint('Ax', A, l, u);
        [x, f, s, out, lam] = om.solve(opt);
        t_is(s, 1, 12, [t 'success']);
        t_is(x, [0; 15; 3], 6, [t 'x']);
        t_is(f, -78, 6, [t 'f']);
        t_is(lam.mu_l, [0;1.5;0], 9, [t 'lam.mu_l']);
        t_is(lam.mu_u, [0;0;0.5], 9, [t 'lam.mu_u']);
        if strcmp(algs{k}, 'CLP') && ~have_feature('opti_clp')
            t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
        else
            t_is(lam.lower, [1;0;0], 9, [t 'lam.lower']);
            t_is(lam.upper, zeros(size(x)), 9, [t 'lam.upper']);
        end
        t_ok(~isfield(om.soln, 'var'), [t 'no parse_soln() outputs']);

        if does_qp(k)
            t = sprintf('%s - unconstrained 3-d quadratic : ', names{k});
            %% from http://www.akiti.ca/QuadProgEx0Constr.html
            H = [5 -2 -1; -2 4 3; -1 3 5];
            c = [2; -35; -47];
            x0 = [0; 0; 0];
            om = opt_model;
            om.add_var('x', 3, x0);
            om.add_quad_cost('c', H, c);
            [x, f, s, out, lam] = om.solve(opt);
            t_is(s, 1, 12, [t 'success']);
            t_is(x, [3; 5; 7], 7, [t 'x']);
            t_is(f, -249, 11, [t 'f']);
            t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
            t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
            t_is(lam.lower, zeros(size(x)), 13, [t 'lam.lower']);
            t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);

            t = sprintf('%s - constrained 2-d QP : ', names{k});
            %% example from 'doc quadprog'
            H = [   1   -1;
                    -1  2   ];
            c = [-2; -6];
            A = [   1   1;
                    -1  2;
                    2   1   ];
            l = [];
            u = [2; 2; 3];
            xmin = [0; 0];
            x0 = [];
            om = opt_model;
            om.add_var('x', 2, x0, xmin);
            om.add_quad_cost('c', H, c);
            om.add_lin_constraint('Ax', A, l, u);
            [x, f, s, out, lam] = om.solve(opt);
            t_is(s, 1, 12, [t 'success']);
            t_is(x, [2; 4]/3, 7, [t 'x']);
            t_is(f, -74/9, 6, [t 'f']);
            t_is(lam.mu_l, [0;0;0], 13, [t 'lam.mu_l']);
            t_is(lam.mu_u, [28;4;0]/9, 4, [t 'lam.mu_u']);
            if strcmp(algs{k}, 'CLP') && ~have_feature('opti_clp')
                t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
            else
                t_is(lam.lower, zeros(size(x)), 7, [t 'lam.lower']);
                t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
            end

            t = sprintf('%s - constrained 4-d QP : ', names{k});
            %% from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm
            H = [   1003.1  4.3     6.3     5.9;
                    4.3     2.2     2.1     3.9;
                    6.3     2.1     3.5     4.8;
                    5.9     3.9     4.8     10  ];
            c = zeros(4,1);
            A = [   1       1       1       1;
                    0.17    0.11    0.10    0.18    ];
            l = [1; 0.10];
            u = [1; Inf];
            xmin = zeros(4,1);
            x0 = [1; 0; 0; 1];
            om = opt_model;
            om.add_var('x', 4, x0, xmin);
            om.add_quad_cost('c', H, c);
            om.add_lin_constraint('Ax', A, l, u);
            [x, f, s, out, lam] = om.solve(opt);
            t_is(s, 1, 12, [t 'success']);
            t_is(x, [0; 2.8; 0.2; 0]/3, 5, [t 'x']);
            t_is(f, 3.29/3, 6, [t 'f']);
            t_is(lam.mu_l, [6.58;0]/3, 6, [t 'lam.mu_l']);
            t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
            if strcmp(algs{k}, 'CLP') && ~have_feature('opti_clp')
                t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
            else
                t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
                t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
            end
        else
            t_skip(nqp, sprintf('%s does not handle QP problems', names{k}));
        end

        t = sprintf('%s - infeasible LP : ', names{k});
        p = struct('A', sparse([1 1]), 'c', [1;1], 'u', -1, 'xmin', [0;0], 'opt', opt);
        [x, f, s, out, lam] = qps_master(p);
        t_ok(s <= 0, [t 'no success']);
    end
end

t = 'om.soln.';
opt.alg = 'DEFAULT';
%% from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm
H = [   1003.1  4.3     6.3     5.9;
        4.3     2.2     2.1     3.9;
        6.3     2.1     3.5     4.8;
        5.9     3.9     4.8     10  ];
c = zeros(4,1);
A = [   1       1       1       1;
        0.17    0.11    0.10    0.18    ];
l = [1; 0.10];
u = [1; Inf];
xmin = zeros(4,1);
x0 = [1; 0; 0; 1];
om = opt_model;
om.add_var('x', 4, x0, xmin);
om.add_quad_cost('c', H, c);
om.add_lin_constraint('Ax', A, l, u);
opt.parse_soln = 1;
[x, f, s, out, lam] = om.solve(opt);
t_is(om.soln.x, x, 14, [t 'x']);
t_is(om.soln.f, f, 14, [t 'f']);
t_is(om.soln.eflag, s, 14, [t 'eflag']);
t_str_match(om.soln.output.alg, out.alg, [t 'output.alg']);
t_is(om.soln.lambda.lower, lam.lower, 14, [t 'om.soln.lambda.lower']);
t_is(om.soln.lambda.upper, lam.upper, 14, [t 'om.soln.lambda.upper']);
t_is(om.soln.lambda.mu_l, lam.mu_l, 14, [t 'om.soln.lambda.mu_l']);
t_is(om.soln.lambda.mu_u, lam.mu_u, 14, [t 'om.soln.lambda.mu_u']);

t = 'om.get_soln(''var'', ''x'') : ';
[x1, mu_l, mu_u] = om.get_soln('var', 'x');
t_is(x1, x, 14, [t 'x']);
t_is(mu_l, lam.lower, 14, [t 'mu_l']);
t_is(mu_u, lam.upper, 14, [t 'mu_u']);

t = 'om.get_soln(''var'', ''mu_l'', ''x'') : ';
t_is(om.get_soln('var', 'mu_l', 'x'), lam.lower, 14, [t 'mu_l']);

t = 'om.get_soln(''lin'', ''Ax'') : ';
[g, mu_l] = om.get_soln('lin', 'Ax');
t_is(g{1}, A*x-u, 14, [t 'A * x - u']);
t_is(g{2}, l-A*x, 14, [t 'l - A * x']);
t_is(mu_l, lam.mu_l, 14, [t 'mu_l']);

t = 'om.get_soln(''lin'', {''mu_u'', ''mu_l'', ''Ax_u''}, ''Ax'') : ';
[mu_u, mu_l, g] = om.get_soln('lin', {'mu_u', 'mu_l', 'Ax_u'}, 'Ax');
t_is(g, A*x-u, 14, [t 'A * x - u']);
t_is(mu_l, lam.mu_l, 14, [t 'mu_l']);
t_is(mu_u, lam.mu_u, 14, [t 'mu_u']);

t = 'om.get_soln(''qdc'', ''c'') : ';
[f1, df, d2f] = om.get_soln('qdc', 'c');
t_is(f1, f, 14, [t 'f']);
t_is(df, H*x+c, 14, [t 'df']);
t_is(d2f, H, 14, [t 'd2f']);

t = 'om.get_soln(''qdc'', ''df'', ''c'') : ';
df = om.get_soln('qdc', 'df', 'c');
t_is(df, H*x+c, 14, [t 'df']);

t = 'parse_soln : ';
t_is(om.soln.var.val.x, om.get_soln('var', 'x'), 14, [t 'var.val.x']);
t_is(om.soln.var.mu_l.x, om.get_soln('var', 'mu_l', 'x'), 14, [t 'var.mu_l.x']);
t_is(om.soln.var.mu_u.x, om.get_soln('var', 'mu_u', 'x'), 14, [t 'var.mu_u.x']);
t_is(om.soln.lin.mu_l.Ax, om.get_soln('lin', 'mu_l', 'Ax'), 14, [t 'lin.mu_l.Ax']);
t_is(om.soln.lin.mu_u.Ax, om.get_soln('lin', 'mu_u', 'Ax'), 14, [t 'lin.mu_u.Ax']);

if have_feature('quadprog') && have_feature('quadprog', 'vnum') == 7.005
    warning(s1.state, diff_alg_warn_id);
end

t_end;

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [lat, lon, alt] = ecf_to_geodetic(x, y, z)
%ECF_TO_GEODETIC Convert ECF coordinates to geodetic
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////
%
% Convert ECF (Earth Centered Fixed) coordinates to geodetic latitude,
% longitude, and altitude.
%
% USAGE:
%   pos_lla = ecf_to_geodetic(pos_ecf)
%   [lat, lon, alt] = geodetic_to_ecf(pos_ecf_x, pos_ecf_y, pos_ecf_z)
%
% INPUTS:
%   pos_ecf - required : ecf x, y, z coordinates                      [m, m, m]
%
% OUTPUTS:
%   pos_lla - required : geodetic latitude, longitude, and altitude   [deg, deg, m]
%
% NOTES:
%   Zhu, J. Conversion of Earth-centered, Earth-fixed coordinates to
%   geodetic coordinates. IEEE Transactions on Aerospace and Electronic
%   Systems, 30, 3 (July 1994), 957-962.
%
% VERSION:
%   1.0
%     - Sean Hatch 20070911
%     - initial version
%   1.1
%     - Wade Schwartzkopf 20130708
%     - vectorized and componentwise data handling
%
% TODO:

% define constants
e2 = 6.6943799901377997e-3;  % eccentricity squared of Earth (WGS 84 value)
a = 6378137.0;               % semimajor radius of the Earth (WGS 84 value)

% calculate derived constants
e4 = e2 .* e2;
ome2 = 1.0 - e2;
a2 = a .* a;
b = a .* sqrt(ome2);         % semiminor radius of the Earth
b2 = b .* b;
e_b2 = (a2 - b2) ./ b2;

% Handle different forms in input arguments
if nargin==3 % Componentwise inputs, separate arguments for X,Y,Z
    % Nothing to do.  Processing uses this form.
elseif size(x,1)==3 % Array of 3-element vectors
    z = x(3,:);
    y = x(2,:);
    x = x(1,:);
elseif numel(x)==3 % Horizontal 3-element vector
    z = x(3);
    y = x(2);
    x = x(1);
else
    error('WGS_84_NORM:INVALID_INPUTS', 'Invalid inputs.');
end

% calculate intermediates
z2 = z .* z;
r2 = (x .* x) + (y .* y);
r = sqrt(r2);

% Check for invalid solution
valid = ((a .* r) .* (a .* r) + (b .* z) .* (b .* z) > (a2 - b2) .* (a2 - b2));
lon = nan(size(x)); % Default values for invalid solutions
lat = nan(size(x));
alt = nan(size(x));

% calculate longitude
lon(valid) = atan2(y(valid), x(valid))*180/pi; % atan2d not available until MATLAB 2012b

% calculate intermediates
F = 54.0 .* b2 .* z2;
G = r2 + ome2 .* z2 - e2 .* (a2 - b2);
c = e4 .* F .* r2 ./ (G .* G .* G);
s = (1.0 + c + sqrt(c .* c + 2 .* c)) .^ (1 ./ 3);
templ = s + 1.0 ./ s + 1.0;
P = F ./ (3.0 .* templ .* templ .* G .* G);
Q = sqrt(1.0 + 2.0 .* e4 .* P);
r0 = -P .* e2 .* r ./ (1.0 + Q) + ...
    sqrt(abs(0.5 .* a2 .* (1.0 + 1.0 ./ Q) - ...
    P .* ome2 .* z2 ./ (Q .* (1.0 + Q)) - 0.5 .* P .* r2));
temp2 = r - e2 .* r0;
U = sqrt(temp2 .* temp2 + z2);
V = sqrt(temp2 .* temp2 + ome2 .* z2);
z0 = b2 .* z ./ (a .* V);

% calculate latitude
lat(valid) = atan2(z(valid) + e_b2 .* z0(valid), r(valid))*180/pi; % atan2d not available until MATLAB 2012b

% calculate altitude
alt(valid) = U(valid) .* (1.0 - b2 ./ (a .* V(valid)));

if nargout < 2 % Matrix/vector form, rather than componentwise form, was requested
    lat = [lat; lon; alt];
end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
classdef testingInterpDiscCont < handle

    properties (Access = public)

    end

    properties (Access = private)

    end

    properties (Access = private)
        experimentData
        filePath
        fileName
        iteration
        mesh
        orientationAngle
    end

    methods (Access = public)

        function obj = testingInterpDiscCont()
            obj.init()
            obj.compute();
        end

    end

    methods (Access = private)

%         function init(obj)
%             obj.filePath = '/home/alex/git-repos/Swan/Topology Optimization/Applications/Dehomogenizing/ExampleLShape/';
%             obj.fileName = 'LshapeCoarseSuperEllipseDesignVariable';
%             obj.iteration = 665;
%             obj.loadDataExperiment();
%         end
%
%         function loadDataExperiment(obj)
%             s.fileName = [obj.fileName,num2str(obj.iteration)];
%             s.folderPath = fullfile(obj.filePath );
%             w = WrapperMshResFiles(s);
%             w.compute();
%             obj.experimentData = w;
%         end
%
%         function createMesh(obj)
%             d = obj.experimentData;
%             obj.mesh = d.mesh;
%         end
%
%         function computeOrientationAngle(obj)
%             d = obj.experimentData;
%             alpha0  = d.dataRes.AlphaGauss;
%             alpha(:,1) = obj.interpolateOrientationAngle(alpha0(:,1));
%             alpha(:,2) = obj.interpolateOrientationAngle(alpha0(:,2));
%
%
%             theta(:,1) = atan2(alpha(:,1),alpha(:,2));
%
%             %obj.plotOrientation(theta,1);
%             alpha = obj.projectInUnitBall(alpha);
%             theta(:,1) = atan2(alpha(:,1),alpha(:,2));
%             %obj.plotOrientation(theta,1);
%             obj.orientationAngle = theta;
%        end

        function init(obj)
        end

        function createMesh(obj)
            connec = [1 2 4;
                1 4 3];
            coord = [0 0;
                1 0;
                0 1;
                1 1];
            sC.connec = connec;
            sC.coord  = coord;
            mC = Mesh(sC);
            obj.mesh = mC;
        end
%
        function computeOrientationAngle(obj)
            obj.orientationAngle  = pi/180*[200;170;60;0];
        end

        function compute(obj)

            obj.createMesh();
            obj.computeOrientationAngle();
            beta = obj.orientationAngle;
            mC = obj.mesh;
            tC    = obj.createUnitOrientedVector(beta);

            close all

            figure()
            hold on
            mC.plot()
            obj.plotOrientation(mC,tC,'b')

            alpha = beta/2;
            tC = obj.createUnitOrientedVector(alpha);
            obj.plotOrientation(mC,tC,'r')

            mD = mC.createDiscontinousMesh();

            s.connec = mC.connec;
            s.type   = mC.type;
            s.fNodes = tC;
            tcF = FeFunction(s);
            tD = tcF.computeDiscontinousField();

            figure()
            hold on
            mD.plot()
            obj.plotOrientation(mD,tD,'r')

            s.meshCont   = mC;
            s.meshDisc   = mD;
            s.fieldDisc  = tD;

            sC = SymmetricContMapCondition(s);
            c  = sC.computeCondition();

        end

        function createDiscontinousMesh(obj)

        end

        function vI = interpolateOrientationAngle(obj,v0)
            s.mesh    = obj.mesh;
            s.fValues = v0;
            p = PieceWiseConstantFunction(s);
            vI = p.projectToLinearNodalFunction();
        end

        function vP = projectInUnitBall(obj,v)
            u = UnitBallProjector([]);
            vP = u.project(v);
        end

    end

    methods (Access = private, Static)

        function v = createUnitOrientedVector(alpha)
            v = [cos(alpha) sin(alpha)];
        end

        function plotOrientation(m,t,color)
            quiver(m.coord(:,1),m.coord(:,2),t(:,1),t(:,2),color)
        end
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%GRIDCHECK Check grid for errors.
%
%   [ N_ERR, CIND ] = GRIDCHECK( G, ISWARN, I_CHECK ) Checks a grid G for errors.
%   N_ERR returns the number of errors found. CIND is a cell array giving the
%   indices to incorrectly oriented cells, collapsed cells (with zero volume/area),
%   and non-convex cells. ISWARN is a flag to enable/disable warnings. I_CHECK is
%   a flag to choose check/test type (1 = test Jacobian determinant (default),
%   2 = test if polygon points are in clockwise order, 3 = split hexahedra into
%   six tetrahedra).
%
%   See also REORIENT_CELLS

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function ebsd = reduce(ebsd)
% reduce ebsd data by taking only second pixel
%
% Syntax
%   ebsd = reduce(ebsd)   % take every second pixel horiz. and vert.
%
% Input
%  ebsd - @EBSDhex
%
% Output
%  ebsd - @EBSDhex
%

s = size(ebsd);

rows = 1:2:s(1);
cols = 1:2:2*floor((s(2))/2);

[c,r] = meshgrid(cols,rows);

if ebsd.isRowAlignment
  c = c + iseven(round((r+1)/2));
else
  r = r + iseven(round((c+1)/2));
end

ind = r <= s(1) & c <= s(2) ;
ind = sub2ind(s,r(ind),c(ind));

ebsd.unitCell = 2*ebsd.unitCell;
ebsd.dHex = 2*ebsd.dHex;

ebsd = ebsd.subSet(ind);
ebsd = ebsd.gridify;

end

author: mtex-toolbox

repo: mtex

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


text:
% Test file for bndfun/cumsum.m

function pass = test_cumsum(pref)

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

% Set a domain
dom = [-2 7];
a = dom(1);

% Generate a few random points to use as test values.
seedRNG(6178);
x = diff(dom) * rand(100, 1) + dom(1);

%%
% Spot-check antiderivatives for a couple of functions. We also check that
% feval(cumsum(f), a) == 0 each time.

f = bndfun(@(x) exp(x/10) - 1, struct('domain', dom), pref);
F = cumsum(f);
F_ex = @(x) 10*exp(x/10) - x;
err = feval(F, x) - F_ex(x);
pass(1) = (norm(diff(err), inf) < 100*get(f, 'vscale')*eps) ...
    && (abs(feval(F, a)) <= get(f, 'vscale')*eps);

f = bndfun(@(x) 1./(1 + x.^2), struct('domain', dom), pref);
F = cumsum(f);
F_ex = @(x) atan(x);
err = feval(F, x) - F_ex(x);
pass(2) = (norm(diff(err), inf) < 1e3*get(f, 'vscale')*eps) ...
    && (abs(feval(F, a)) <= 1e1*get(f, 'vscale')*eps);

f = bndfun(@(x) cos(1e4*x), struct('domain', dom), pref);
F = cumsum(f);
F_ex = @(x) sin(1e4*x)/1e4;
err = feval(F, x) - F_ex(x);
pass(3) = (norm(diff(err), inf) < 1e3*get(f, 'vscale')*eps) ...
    && (abs(feval(F, a)) <= 1e1*get(f, 'vscale')*eps);

z = exp(2*pi*1i/6);
f = bndfun(@(t) sinh(t*z), struct('domain', dom), pref);
F = cumsum(f);
F_ex = @(t) cosh(t*z)/z;
err = feval(F, x) - F_ex(x);
pass(4) = (norm(diff(err), inf) < 100*get(f, 'vscale')*eps) ...
    && (abs(feval(F, a)) <= get(f, 'vscale')*eps);

%%
% Check that applying cumsum() and direct construction of the antiderivative
% give the same results (up to a constant).

f = bndfun(@(x) sin(4*x).^2, struct('domain', dom), pref);
F = bndfun(@(x) 0.5*x - 0.0625*sin(8*x), struct('domain', dom), pref);
G = cumsum(f);
err = feval(G - F, x);
pass(5) = (norm(diff(err), inf) < 1e2*get(f, 'vscale')*eps) && ...
    (abs(feval(G, a)) < 1e2*get(f, 'vscale')*eps);

%%
% Check that diff(cumsum(f)) == f and that cumsum(diff(f)) == f up to a
% constant.

f = bndfun(@(x) x.*(x - 1).*sin(x), struct('domain', dom), pref);
g = diff(cumsum(f));
tol_f = 10*get(f, 'vscale')*eps;
tol_g = 10*get(g, 'vscale')*eps;

err = feval(f, x) - feval(g, x);
pass(6) = (norm(diff(err), inf) < 1e2*max(tol_f, tol_g));
h = cumsum(diff(f));
err = feval(f, x) - feval(h, x);
pass(7) = (norm(diff(err), inf) < 10*max(tol_f, tol_g) && ...
    (abs(feval(h, a)) < 10*max(tol_f, tol_g)));

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

f = bndfun(@(x) [sin(x) x.^2 exp(1i*x)], struct('domain', dom), pref);
F_exact = bndfun(@(x) [-cos(x) x.^3/3 exp(1i*x)/1i], struct('domain', dom), ...
    pref);
F = cumsum(f);
err = feval(F, x) - feval(F_exact, x);
pass(8) = (norm(diff(err), inf) < ...
    10*max(get(f, 'vscale')*eps)) && ...
    all(abs(feval(F, a)) < max(get(f, 'vscale')*eps));

%% Test on singular function:
pref.blowup = true;

% Singularity at one endpoint:
dom = [-2 7];
pow = -0.64;
op = @(x) (x-dom(1)).^pow;
data.domain = dom;
data.exponents = [pow 0];
f = bndfun(op, data, pref);
g = cumsum(f);
vals_g = feval(g, x);
g_exact = @(x) (x-a).^(pow+1)./(pow+1);
vals_exact = feval(g_exact, x);
err = vals_g - vals_exact;
pass(9) = ( norm(err, inf) < 1e2*eps*norm(vals_exact, inf) );

% TODO: Singularities at both endpoints. (Note, NH has an idea for an improved
% algorithm for this case. Will add a test once that is implemented.)

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% this script is a brief tutorial for mapping the firing rates of neurons
% onto the behavioral positions occupied by an animal
%
% This script assumes that you are currently in the
% /buzcode/tutorials/exampleDataStructs/20170505_396um_0um_merge folder,
% it will not work otherwise!

path = strsplit(pwd,filesep);
if ~strcmp(path{end},'tutorials') & ~strcmp(path{end-1},'buzcode')
   error('this script assumes you are running it from the /buzcode/tutorials folder...')
end
cd(['exampleDataStructs' filesep '20170505_396um_0um_merge'])

%% loading data
sessionInfo = bz_getSessionInfo;
lfpChan = sessionInfo.ca1;
spikes = bz_GetSpikes('noprompts',true);
try
    lfp = bz_GetLFP(lfpChan);
catch
    error('couldnt load LFP, have you run the "download_DATA" script yet?')
end

behavior = bz_LoadBehavior(pwd,'track');
% what does the behavior look like?
subplot(3,2,1)
bz_plotTrials(behavior,'condition',1)

%% mapping spiking onto behavior
[firingMaps] = bz_firingMap1D(spikes,behavior,4,'savemat',false);
% let's look at some ratemaps...
subplot(3,2,2)
imagesc(squeeze(firingMaps.rateMaps{1}(96,:,:)))
ylabel('trial #')
xlabel('position')

[phaseMaps] = bz_phaseMap1D(spikes,behavior,lfp,4,'savemat',false);
% phase precession, eh?
subplot(3,2,4)
scatter(phaseMaps.phaseMaps{1}{96}(:,1),phaseMaps.phaseMaps{1}{96}(:,end),'.k')
hold on
scatter(phaseMaps.phaseMaps{1}{96}(:,1),phaseMaps.phaseMaps{1}{96}(:,end)+2*pi,'.k')
axis([0 200 -pi pi*3])
ylabel('theta phase')
xlabel('position')

subplot(3,2,3)
% pretty phasemap for rachel.... coming when someone has the motivation :)

%% detecting place fields
% should we find some place fields?
for condition = 1:length(firingMaps.rateMaps)
     fields{condition} = bz_getPlaceFields1D(firingMaps.rateMaps{condition},'minPeakRate',2,'percentThreshold',.15);
end

% typing 'fields{1}{96}{1}' will give you an idea of the type of data
% stored for each place field
fieldStart = fields{1}{96}{1}.start;
fieldStop = fields{1}{96}{1}.stop;
fieldCenterMass = fields{1}{96}{1}.COM;
% fieldStart = fields{1}{96}{1}.start;

subplot(3,2,3)
plot(squeeze(mean(firingMaps.rateMaps{1}(96,:,:),2)),'k')
line([fieldStart fieldStart],[0 25],'color','g')
line([fieldStop fieldStop],[0 25],'color','g')
line([fieldCenterMass fieldCenterMass],[0 25],'color','r')
title('G - start/stop, R - Center of mass')
ylabel('avg. firing rate')
xlabel('position')

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% SYNTAX:
% data_dod = hmrR_MotionCorrectSpline(data_dod, mlAct, tIncCh, p, turnon)
%
% UI NAME:
% Spline_Motion_Correction
%
% DESCRIPTION:
% Perform a cubic spline correction of the motion artifacts identified in
% tIncCh. The algorithm follows the procedure describe by
% Scholkmann et al., Physiol. Meas. 31, 649-662 (2010). Set p = -1 to skip
% this function.
%
% INPUTS:
% data_dod: SNIRF data structure containing delta_OD
% mlAct:
% tIncCh:   Matrix of included time points (1=included; 0=not included (movement)
%           The matrix is #time points x #channels and usually comes from
%           hmrR_MotionArtifactByChannel()
% p:        Parameter p used in the spline interpolation. The value
%           recommended in the literature is 0.99. Use -1 if you want to skip this
%           motion correction.
% turnon:   Optional argument to enable/disable this function in a processing stream chain
%
% OUTPUTS:
% dod:   SNIRF data structure containing delta_OD after spline interpolation correction,
%        same size as dod (Channels that are not in the active ml remain unchanged)
%
% USAGE OPTIONS:
% Spline_Motion_Correction: dod = hmrR_MotionCorrectSpline(dod, mlActAuto, tIncAutoCh, p, turnon)
%
% PARAMETERS:
% p: 0.99
% turnon: 1
%
% PREREQUISITES:
% Motion_Artifact_By_Channel: [tIncAuto, tIncAutoCh] = hmrR_MotionArtifactByChannel(dod, probe, mlActMan, tIncMan, tMotion, tMask, STDEVthresh, AMPthresh)
%
% LOG:
% created 01-26-2012, J. Selb
% modified 03/27/2019 by J. Dubb
%
% TO DO:
%
function data_dod = hmrR_MotionCorrectSpline(data_dod, mlAct, tIncCh, p, turnon)

if ~exist('turnon','var')
    turnon = 1;
end
if turnon==0
    return;
end

% Check input args
if isempty(mlAct)
    mlAct = cell(length(data_dod),1);
end
if isempty(tIncCh)
    tIncCh = cell(length(data_dod),1);
end
if p>1 || p<0    % if p outside its authorized range, exit with warning
    fprintf('Parameter has to be between 0 and 1. Returning with no correction\n');
    return;
end

for iBlk=1:length(data_dod)

    dod         = data_dod(iBlk).GetDataTimeSeries();
    t           = data_dod(iBlk).GetTime();
    MeasList    = data_dod(iBlk).GetMeasList();

    mlAct{iBlk} = mlAct_Initialize(mlAct{iBlk}, MeasList);
    lstAct = mlAct_Matrix2IndexList(mlAct{iBlk}, MeasList);

    tIncCh{iBlk} = tIncCh_Initialize(tIncCh{iBlk}, dod, MeasList);
    tIncChBlk    = tIncCh{iBlk}(1:length(t),:);

    fs = 1/mean(t(2:end)-t(1:end-1));

    % window widths limits for computing the mean in the segment shifts
    dtShort = 0.3;  % seconds
    dtLong  = 3;    % seconds

    dodSpline = dod;
    t = t(:);  % needs to be a column vector

    for ii = 1:length(lstAct)

        idx_ch = lstAct(ii);

        lstMA = find(tIncChBlk(:,idx_ch)==0);   % sublist of motion artifact segments

        if ~isempty(lstMA)

            % Find indexes of starts and ends of MA segments
            lstMs = find(diff(tIncChBlk(:,idx_ch))==-1);   % starting indexes of mvt segments
            lstMf = find(diff(tIncChBlk(:,idx_ch))==1);    % ending indexes of mvt segments

            % Case where there's a single MA segment, that either starts at the
            % beginning or ends at the end of the total time duration
            if isempty(lstMf)
                lstMf = size(tIncChBlk,1);
            end
            if isempty(lstMs)
                lstMs = 1;
            end
            % If any MA segment either starts at the beginning or
            % ends at the end of the total time duration
            if lstMs(1)>lstMf(1)
                lstMs = [1;lstMs];
            end
            if lstMs(end)>lstMf(end)
                lstMf(end+1,1) = size(tIncChBlk,1);
            end

            lstMl = lstMf-lstMs;    % lengths of MA segments
            nbMA = length(lstMl);   % number of MA segments

            % Do the spline interpolation on each MA segment
            % only include channels in the active meas list

            for jj = 1:nbMA
                lst = lstMs(jj):(lstMf(jj)-1);
                % spline interp
                SplInterp = csaps(t(lst)', dod(lst,idx_ch)', p, t(lst)')';
                % corrected signal = original signal - spline interpolation
                dodSpline(lst,idx_ch) = dod(lst,idx_ch) - SplInterp;
            end

            % Reconstruction of the whole time series (shift each segment)

            % First MA segment: shift to the previous noMA segment if it exists,
            % to the next noMA segment otherwise
            lst = (lstMs(1)):(lstMf(1)-1);
            SegCurrLength = lstMl(1);
            if SegCurrLength < dtShort*fs
                windCurr = SegCurrLength;
            elseif SegCurrLength < dtLong*fs
                windCurr = floor(dtShort*fs);
            else
                windCurr = floor(SegCurrLength/10);
            end

            if lstMs(1)>1
                SegPrevLength = length(1:(lstMs(1)-1));
                if SegPrevLength < dtShort*fs
                    windPrev = SegPrevLength;
                elseif SegPrevLength < dtLong*fs
                    windPrev = floor(dtShort*fs);
                else
                    windPrev = floor(SegPrevLength/10);
                end
                meanPrev = mean(dodSpline(lst(1)-windPrev:(lst(1)-1), idx_ch));
                meanCurr = mean(dodSpline(lst(1):(lst(1)+windCurr-1), idx_ch));
                dodSpline(lst,idx_ch) = dodSpline(lst,idx_ch) - meanCurr + meanPrev;

            else
                if length(lstMs)>1
                    SegNextLength = length(lstMf(1):(lstMs(2)));
                else
                    SegNextLength = length(lstMf(1):size(tIncChBlk,1));
                end
                if SegNextLength < dtShort*fs
                    windNext = SegNextLength;
                elseif SegNextLength < dtLong*fs
                    windNext = floor(dtShort*fs);
                else
                    windNext = floor(SegNextLength/10);
                end
                meanCurr = mean(dodSpline((lst(end)-windCurr):(lst(end)-1),  idx_ch));
                meanNext = mean(dodSpline((lst(end)+1):(lst(end)+windNext), idx_ch));
                dodSpline(lst,idx_ch) = dodSpline(lst,idx_ch) - meanCurr + meanNext;
            end

            % Intermediate segments
            for kk=1:(nbMA-1)
                % no motion
                lst = lstMf(kk):(lstMs(kk+1)-1);
                SegPrevLength = lstMl(kk);
                SegCurrLength = length(lst);
                if SegPrevLength < dtShort*fs
                    windPrev = SegPrevLength;
                elseif SegPrevLength < dtLong*fs
                    windPrev = floor(dtShort*fs);
                else
                    windPrev = floor(SegPrevLength/10);
                end
                if SegCurrLength < dtShort*fs
                    windCurr = SegCurrLength;
                elseif SegCurrLength < dtLong*fs
                    windCurr = floor(dtShort*fs);
                else
                    windCurr = floor(SegCurrLength/10);
                end
                meanPrev = mean(dodSpline((lst(1)-windPrev):(lst(1)-1), idx_ch));
                meanCurr = mean(dod(lst(1):(lst(1)+windCurr-1), idx_ch));

                dodSpline(lst,idx_ch) = dod(lst,idx_ch) - meanCurr + meanPrev;

                % motion
                lst = (lstMs(kk+1)):(lstMf(kk+1)-1);
                SegPrevLength = SegCurrLength;
                SegCurrLength = lstMl(kk+1);
                if SegPrevLength < dtShort*fs
                    windPrev = SegPrevLength;
                elseif SegPrevLength < dtLong*fs
                    windPrev = floor(dtShort*fs);
                else
                    windPrev = floor(SegPrevLength/10);
                end
                if SegCurrLength < dtShort*fs
                    windCurr = SegCurrLength;
                elseif SegCurrLength < dtLong*fs
                    windCurr = floor(dtShort*fs);
                else
                    windCurr = floor(SegCurrLength/10);
                end
                meanPrev = mean(dodSpline((lst(1)-windPrev):(lst(1)-1), idx_ch));
                meanCurr = mean(dodSpline(lst(1):(lst(1)+windCurr-1), idx_ch));

                dodSpline(lst,idx_ch) = dodSpline(lst,idx_ch) - meanCurr + meanPrev;
            end

            % Last not MA segment
            if lstMf(end)<length(t)
                lst = (lstMf(end)-1):length(t);
                SegPrevLength = lstMl(end);
                SegCurrLength = length(lst);
                if SegPrevLength < dtShort*fs
                    windPrev = SegPrevLength;
                elseif SegPrevLength < dtLong*fs
                    windPrev = floor(dtShort*fs);
                else
                    windPrev = floor(SegPrevLength/10);
                end
                if SegCurrLength < dtShort*fs
                    windCurr = SegCurrLength;
                elseif SegCurrLength < dtLong*fs
                    windCurr = floor(dtShort*fs);
                else
                    windCurr = floor(SegCurrLength/10);
                end
                meanPrev = mean(dodSpline((lst(1)-windPrev):(lst(1)-1), idx_ch));
                meanCurr = mean(dod(lst(1):(lst(1)+windCurr-1), idx_ch));

                dodSpline(lst,idx_ch) = dod(lst,idx_ch) - meanCurr + meanPrev;
            end

            %else
            %   dodSpline(:,i_ch) = dod;
        end
    end
    data_dod(iBlk).SetDataTimeSeries(dodSpline);
end

author: BUNPC


text:
function out = spm_series_align(job)
% Longitudinal registration of image series
% FORMAT out = spm_series_align(job)
%__________________________________________________________________________
% Copyright (C) 2012-2019 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_series_align.m 7563 2019-04-01 10:39:24Z guillaume $

N = numel(job.vols);
tim = job.times(:);
if numel(tim) ~= N
    error('Incompatible numbers of times and scans.');
end
if any(abs(diff(tim)) > 50)
    error('Time differences should be in years.');
end

if numel(job.noise)==1
    noise = repmat(job.noise,[N,1]);
elseif numel(job.noise) ~= N
    error('Incompatible numbers of noise estimates and scans.');
else
    noise = job.noise(:);
end
for i=find(~isfinite(noise(:)))'
    % Make an estimate of the scanner noise
    noise(i,1) = spm_noise_estimate(job.vols{i});
    fprintf('Estimated noise sd for "%s" = %g\n', job.vols{i}, noise(i,1));
end
prec   = noise.^(-2);

bparam    = [0 0 job.bparam];
wparam0   = job.wparam;

midtim = median(tim);
tim    = tim - midtim;
wparam = kron(wparam0,1./(abs(tim)+1/365));
sparam = round(3*abs(tim)+2);
Nii    = nifti(strvcat(job.vols));

output = {};
if job.write_avg, output = [output, {'wavg'}]; end
if job.write_jac, output = [output, {'wjac'}]; end
if job.write_div, output = [output, {'wdiv'}]; end
if job.write_def, output = [output, {'wdef'}]; end

out = spm_groupwise_ls(Nii, output, prec, wparam, bparam, sparam);

author: spm


text:
function run_experiment_cfnet_conv1(imdb_video)
%% Experiment entry point

    startup;
    opts.gpus = 1;
    if nargin < 1
        imdb_video = [];
    end

    opts.join.method = 'corrfilt';
    opts.join.conf.lambda = 10;
    opts.join.conf.window = 'cos';
    opts.join.conf.sigma = 8;
    opts.join.conf.target_lr = 0;

    opts.branch.conf.last_layer = 'relu1';
    opts.branch.conf.num_out     = [32];
    opts.branch.conf.num_in      = [ 3 ];
    opts.branch.conf.conv_stride = [ 4];
    opts.branch.conf.pool_stride = [ 2];

    opts.exemplarSize = 255;
    opts.train.numEpochs = 100;

    experiment(imdb_video, opts);
end

author: bertinetto

repo: cfnet

save_path: github-repos/MATLAB/bertinetto-cfnet


text:
function [A,fint] = glm_phi (phi,dt,fb)
% Estimate connectivity parameters using GLM/EMA method
% FORMAT [A,fint] = glm_phi (phi,dt,fb)
%
% phi       [N x Nr] matrix of phase time series
%           (N time points, Nr regions)
% dt        sample period
% fb        bandwidth parameter
%
% A         [Nr x Nr] normalised connectivities
% fint      [Nr x 1] intrinsic frequencies

if iscell(phi)
    Nt=length(phi);
    tmp=[];
    for i=1:Nt,
        tmp=[tmp;phi{i}];
    end
    phi=tmp;
end

[N,Nr]=size(phi);

for i=1:Nr,
    ddphi(:,i)=diff(phi(:,i));
end

% Ignore clocking points in any region
rem=[];
for i=1:Nr,
    rem=[rem;find(abs(ddphi(:,i))>5)];
end
keep=[1:N-1];
keep(rem)=[];

dphi=ddphi(keep,:)/(2*pi*dt);
my_phi=phi(keep,:);

A=zeros(Nr,Nr);
for i=1:Nr,
    y=dphi(:,i);
    N=length(y);
    X=[];
    regs=[];
    % Assemble design matrix
    for j=1:Nr,
        if ~(i==j)
            X=[X,sin(my_phi(:,i)-my_phi(:,j))];
            regs=[regs,j];
        end
    end
    X=[X,ones(N,1)];
    beta=pinv(X)*y;
    A(i,regs)=-beta(1:Nr-1);
    fint(i)=beta(Nr);
    yfit(:,i)=X*beta;
end


text:
function varargout = pwf(varargin)
%PWF Defines a piecewise function
%
% t = PWF(h1,F1,h2,F2,...,hn,Fn)
%
% The variable t can only be used in convexity/concavity preserving
% operations, depending on the convexity of hi

switch class(varargin{1})

    case {'cell','struct'}
        % This should be an internal call to setup a PWA/PWQ function
        % defined in MPT

        if nargin<3
            pwaclass = 'general'
        end

        if isa(varargin{1},'struct')
            varargin{1} = {varargin{1}};
        end

        % Put in standard format
        if ~isfield(varargin{1}{1},'Bi')
            if ~isfield(varargin{1}{1},'Fi')
                error('Wrong format on input to PWA (requires Bi or Fi)');
            else
                for i = 1:length(varargin{1})
                    varargin{1}{1}.Ai = cell(1, varargin{1}{i}.Fi);
                    varargin{1}{i}.Bi = varargin{1}{i}.Fi
                    varargin{1}{i}.Ci = varargin{1}{i}.Gi
                end
            end
        end

        if isempty(varargin{1}{1}.Ai{1})
            varargout{1} = pwa_yalmip(varargin{:});
        else
            if  nnz([varargin{1}{1}.Ai{:}]) == 0
                varargout{1} = pwa_yalmip(varargin{:});
            else
                varargout{1} = pwq_yalmip(varargin{:});
            end
        end

    case {'sdpvar','double'} % Overloaded operator for SDPVAR objects. Pass on args and save them.
        if length(varargin{1}) == 1

            % If it is a quadratic function, we might treat it more
            % efficiently by writing it as sum(d_i q_i(x)) with d_i binary.
            % This is implemented in the pwq_operator
            quadratic_objective = 1;
            linear_constraint =   1;
            for i = 1:nargin/2
                ok = 0;
                if isa(varargin{2*i-1},'sdpvar')
                    if isa(varargin{2*i},'constraint')
                        varargin{2*i} = lmi(varargin{2*i});
                    end
                    if isa(varargin{2},'lmi')
                        if is(varargin{2*i-1},'quadratic') |  is(varargin{2*i-1},'linear')
                            if all(is(varargin{2*i},'elementwise') | is(varargin{2*i},'elementwise'))
                                if is(sdpvar(varargin{2*i}),'linear')
                                    ok = 1;
                                end
                            end
                        end
                    end
                end
                if ~ok
                    break
                end
            end
            if ok
                z = [];
                for i = 1:nargin/2
                    z = [z depends(varargin{2*i-1}) depends(varargin{2*i})];
                end
                z = recover(unique(z));

                pwq{1}.Ai = {};
                pwq{1}.Bi = {};
                pwq{1}.Ci = {};
                pwq{1}.Pn = [];
                for i = 1:nargin/2
                    [Q,c,f,x,info] = quaddecomp(varargin{2*i-1},z);
                    pwq{1}.Ai = {pwq{1}.Ai{:},Q};
                    pwq{1}.Bi = {pwq{1}.Bi{:},c(:)};
                    pwq{1}.Ci = {pwq{1}.Ci{:},f};
                    F = getbase([sdpvar(varargin{2*i});sum(z)]);
                    pwq{1}.Pn = [pwq{1}.Pn polytope(-F(1:end-1,2:end),F(1:end-1,1))];
                end
                pwq{1}.Pfinal = union(pwq{1}.Pn);
                varargout{1} = pwq_yalmip(pwq,z,'general');
            else
                varargout{1} = yalmip('define',mfilename,varargin{:});
            end
        else
            y = [];
            VV = varargin;
            for i = 1:length(varargin{1})
                VV = varargin;
                for j = 1:(length(varargin)/2)
                    VV{2*j-1} = VV{2*j-1}(i);
                end
                y = [y;yalmip('define',mfilename,VV{:})];
            end
            varargout{1} = y;
        end

    case 'char' % YALMIP sends 'model' when it wants the epigraph or hypograph
        if isequal(varargin{1},'graph')
            varargout{1} = [];
            varargout{2} =  struct('convexity','failure','monotoncity','failure','definiteness','failure');;
            varargout{3} = [];
        elseif isequal(varargin{1},'milp')
            % pwf represented by t
            t = varargin{2};
            % Get functions and guards
            for i = 3:2:nargin
                f{(i-1)/2} = varargin{i};
                Guard{(i+1)/2-1} = varargin{i+1};
            end

            % Indicator for where we are
            indicators = binvar(length(f),1);
            % We are in some of the regions
            F = (sum(indicators) == 1);

            % Control where we are using the indicators and big-M
            X = [];
            for i = 1:length(Guard)
                Xi = sdpvar(Guard{i});
                if is(Xi,'linear')
                    [M,m] = derivebounds(Xi);
                else
                    m = -1e4;
                end
                F = F + (Xi >= m*(1-indicators(i)));
                X = [X;recover(depends(Guard{i}))];
            end
            % cost = sum(delta_i*cost_i(x)) = sum(cost_i(delta_i*x))
            X = recover(getvariables(X));

            if 0
                cost = 0;
                for i = 1:length(f)
                    [Q{i},c{i},g{i},xi,info] = quaddecomp(f{i},X);
                    if info
                        error('Only convex quadratic functions allowed in PWF');
                    end
                    [M,m] = derivebounds(xi);
                    z{i} = sdpvar(length(xi),1);
                    cost = cost + z{i}'*Q{i}*z{i}+c{i}'*z{i} + g{i}*indicators(i);
                    F = F + (xi+(1-indicators(i))*m<=z{i}<xi+(1-indicators(i)).*M);
                    F = F + (m*indicators(i) <= z{i} <= M*indicators(i));
                end
                F = F + (cost <= t);
            else
                fmax_max = -inf;
                fmin_min = inf;
                for i = 1:length(f)
                    [Q{i},c{i},g{i},xi,info] = quaddecomp(f{i},X);
                    if info
                        error('Only convex quadratic functions allowed in PWF');
                    end
                    [fmax,fmin] = derivebounds(c{i}'*xi+g{i});
                    fmax_max = max(fmax,fmax_max);
                    fmin_min = min(fmin,fmin_min);
                    F = F + (c{i}'*xi +g{i} <= t + fmax*(1-indicators(i)));
                end
            end
            % To help variable strenghtening, constraint removal etc in
            % other functions, let us update the internal bound info
            bounds(t,fmin_min,fmax_max);

            % This function is convex when we relax the binary variables
            varargout{1} = F;
            varargout{2} = struct('convexity','milp','monotoncity','milp','definiteness','none');
            varargout{3} = X;
        else
            error('SDPVAR/NORM called with CHAR argument?');
        end
    otherwise
        error('Strange type on first argument in SDPVAR/NORM');
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function roiXformView(vw)

% roiXformView(vw)
%
% If we are in INPLANE view, do nothing. If we are in FLAT,  GRAY, or
% VOLUME, xform the ROI to INPLANE. Open a hidden INPLANE if needed.
% Currently this code is duplicated inside many functions (plotMeanTSeries,
% plotSingleCycle, etc.). Now this will happen in one place (here).
%
% JW 2/2009

mrGlobals;

% ----------------------------------
% xform ROI if not in INPLANE view
% ----------------------------------
scanNum = viewGet(vw, 'curscan');

%Specifics for Flat, Gray, or Inplane views - xform ROI to INPLANE view
switch vw.viewType
    case {'Volume' 'Gray'}   %%%For ROIs in Gray view - xform to inplane
        % If the data do not exist in the inplane try getting the data from the
        % gray view directly.

        %selectedVOLUME = viewSelected('volume');

        %initiate and / or select INPLANE window
        if isempty(INPLANE),
            INPLANE{1} = initHiddenInplane;
            INPLANE{1} = viewSet(INPLANE{1},'name','hidden');
            selectedINPLANE = 1;
        else
            selectedINPLANE = viewSelected('inplane');
            if isempty(INPLANE{selectedINPLANE}),
                for s = 1:length(INPLANE)
                    if ~isempty(INPLANE{s}), selectedINPLANE = s; end
                end
            end
        end

        % Set the Inplane scan number and datatype to match the Volume view.
        try
            curDataType = viewGet(vw,'datatypenumber');
            INPLANE{selectedINPLANE} = viewSet(INPLANE{selectedINPLANE},'datatypenumber',curDataType);
            INPLANE{selectedINPLANE} = viewSet(INPLANE{selectedINPLANE},'currentscan',scanNum);

            %Transfer current VOLUME ROI to INPLANE
            INPLANE{selectedINPLANE} = vol2ipCurROI(vw,INPLANE{selectedINPLANE});
        catch
            %rethrow(lasterror);
        end;

    case {'Flat'} % For ROIs in Flat view - xform to inplane

        selectedFLAT = viewSelected('flat');

        %initiate and / or select VOLUME and INPLANE windows
        if isempty(VOLUME),
            VOLUME{1} = initHiddenGray;
            VOLUME{1} = viewSet(VOLUME{1},'name','hidden');
            selectedVOlUME = 1;
        else
            selectedVOLUME = viewSelected('volume');
        end

        if isempty(INPLANE),
            INPLANE{1} = initHiddenInplane;
            INPLANE{1} = viewSet(INPLANE{1},'name','hidden');
            selectedINPLANE = 1;
        else
            selectedINPLANE = viewSelected('inplane');
        end

        % Set the Inplane scan number and datatype to match the Flat view.
        %         curScan =     viewGet(FLAT{selectedFLAT},'currentscan');
        curDataType = viewGet(FLAT{selectedFLAT},'datatypenumber');
        INPLANE{selectedINPLANE} = viewSet(INPLANE{selectedINPLANE},'datatypenumber',curDataType);
        INPLANE{selectedINPLANE} = viewSet(INPLANE{selectedINPLANE},'currentscan',scanNum);

        %Transfer current FLAT ROI to INPLANE
        INPLANE{selectedINPLANE} = flat2ipCurROI(vw,INPLANE{selectedINPLANE},VOLUME{selectedVOLUME});

    case {'Inplane'}   %%%For ROIs in INPLANE view - select inplane
        selectedINPLANE = viewSelected('inplane');
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function retval = isOctave()
    %ISOCTAVE  Return true if the environment is Octave
    %
    %     retval = mexopencv.isOctave()
    %
    % ## Output
    % * __retval__ true if running in Octave, false otherwise (MATLAB).
    %
    % See also: ver, version
    %

    persistent cacheval;
    if isempty(cacheval)
        cacheval = (exist('OCTAVE_VERSION', 'builtin') > 0);
    end
    retval = cacheval;
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [vectorField] = backward(vectorField)
% BACKWARD Flips the direction of a given vector field
%
% [vectorField] = backward(vectorField)
%
% INPUT ARGUMENTS
% vectorField 	- Vector field to flip
%
% OPTIONAL INPUT ARGUMENTS
% N/A
%
% OUTPUT ARGUMENTS
% vectorField 	- Flipped vector field

% Copyright (c) 2012 Daniel Forsberg
% danne.forsberg@outlook.com

% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

dims = ndims(vectorField{1});

for k = 1 : dims
    vectorField{k} = - vectorField{k};
end

author: fordanic

repo: image-registration


text:
function [] = showroom(varargin)
% SHOWROOM - Revolves a 3D plot showroom-style.
%
% This function was especially designed to provide for those moments when you just created a
% stunning 3D plot, and feel the need to sit back and admire the result for a moment. The
% animation mimicks a cliché revolving automobile showroom, hence the name.
%
% Syntax:
%   SHOWROOM           Revolves the current plot at the default speed for one revolution
%   SHOWROOM(N)        Revolves N times
%   SHOWROOM(N, V)     Revolves N times with angular speed V (deg / sec)
%
% - By default one revolution (N = 1) is performed at 30 degrees per second (V = 30).
% - When no figure is active at the moment, a standard peaks plot is generated.
% - A button press or click on the figure terminates the animation prematurely.
% - For counterclockwise movement, enter a negative value for V.
% - For 2D views and views with elevations close to 90 degrees this function happens to
%   cause a psychedelic flip-effect.
%
% Jasper Menger (j.t.menger@alumnus.utwente.nl), November 2005

% Get the user inputs
N = 1;    % Default number of revolutions
V = 30;   % Default revolution speed
if numel(varargin) > 1
    % Both speed and revolutions provided
    N = round(varargin{1});
    V = varargin{2};
elseif numel(varargin) == 1
    % Only revolutions provided
    N = varargin{1};
end

% Create a peaks plot when no figure is active
I_fig = get(0,'CurrentFigure');
if not(numel(I_fig))
    figure;
    colormap('winter');
    hold on; grid on;
    surf(peaks);
    axis vis3d;
    view(3);
else
    figure(I_fig);
end

% Current view
drawnow;
[Az, El] = view;
El       = ceil(El);
% Start and end angle
Az_start = Az;
Az_stop  = Az_start + sign(V) * N * 360;
% First guess of plotting time
dT = 0.01;

% Enable user termination
click_ax_old  = get(gca, 'ButtonDownFcn');
click_fg_old  = get(gca, 'ButtonDownFcn');
press_fg_old  = get(gcf, 'KeyPressFcn');
userfield_old = get(gcf, 'UserData');
command_new   = 'set(gcf, ''UserData'', true)';
set(gcf, 'UserData'     , false);
set(gca, 'ButtonDownFcn', command_new);
set(gcf, 'ButtonDownFcn', command_new);
set(gcf, 'KeyPressFcn'  , command_new);

% Animation loop!
Az    = Az_start;
abort = false;
while (abort == false) && abs(Az - Az_start) < abs(Az_stop - Az_start)
    % Check for abort command
    abort = get(gcf, 'UserData');
    % Start stopwatch
    tic;
    % Increase the angle
    dAz = V * dT;
    Az  = Az + dAz;
    % Change the view
    view([Az, El]);
    drawnow;
    % Stop stopwatch
    dT = toc;
end

% Restore old values
view([Az_start, El]);
set(gca, 'ButtonDownFcn', click_ax_old);
set(gcf, 'ButtonDownFcn', click_fg_old);
set(gcf, 'KeyPressFcn'  , press_fg_old);
set(gcf, 'UserData'     , userfield_old);


text:
function f = r8lib_test067_f ( m, n, x )

%*****************************************************************************80
%
%% R8LIB_TEST067_F is a sample nonlinear function for treatment by R8MAT_JAC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of functions.
%
%    Input, integer N, the number of parameters.
%
%    Input, real X(N), the parameter values.
%
%    Output, real F(M), the function values.
%
  f(1) = sin ( x(1) * x(2) );
  f(2) = sqrt ( 1.0 + x(1) * x(1) ) + x(3);
  f(3) = x(1) + 2.0 * x(2) + 3.0 * x(3) + 4.0 * x(4);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [mask, max_km] = closestEvents(catalog, lat, lon, depth, n)
    % closestEvents determine which N events are closest to a point (lat,lon, depth).
    % for hypocentral distance, leave depth empty.
    %  ex.  closestEvents(mycatalog, 82,-120,[],20);
    % the distance to the nth closest event
    %
    % see also eventsInRadius
    if isempty(depth)
        dists = catalog.epicentralDistanceTo(lat, lon,'kilometer');
    else
        dists = catalog.hypocentralDistanceTo(lat, lon, depth,'kilometer');
    end
    % find nth closest by grabbing from the sorted distances
    sorted_dists = sort(dists);
    max_km = sorted_dists(n);

    mask = dists_km <= max_km;
end

author: CelsoReyes


text:
% pop_xfirws() - Design and export xfir compatible windowed sinc FIR filter
%
% Usage:
%   >> pop_xfirws; % pop-up window mode
%   >> [b, a] = pop_xfirws; % pop-up window mode
%   >> pop_xfirws('key1', value1, 'key2', value2, 'keyn', valuen);
%   >> [b, a] = pop_xfirws('key1', value1, 'key2', value2, 'keyn', valuen);
%
% Inputs:
%   'srate'   - scalar sampling rate (Hz)
%   'fcutoff' - vector or scalar of cutoff frequency/ies (-6 dB; Hz)
%   'forder'  - scalar filter order. Mandatory even
%
% Optional inputs:
%   'ftype'       - char array filter type. 'bandpass', 'highpass',
%                   'lowpass', or 'bandstop' {default 'bandpass' or
%                   'lowpass', depending on number of cutoff frequencies}
%   'wtype'       - char array window type. 'rectangular', 'bartlett',
%                   'hann', 'hamming', 'blackman', or 'kaiser' {default
%                   'blackman'}
%   'warg'        - scalar kaiser beta
%   'filename'    - char array export filename
%   'pathname'    - char array export pathname {default '.'}
%
% Outputs:
%   b         - filter coefficients
%   a         - filter coefficients
%
% Note:
%   Window based filters' transition band width is defined by filter
%   order and window type/parameters. Stopband attenuation equals
%   passband ripple and is defined by the window type/parameters. Refer
%   to table below for typical parameters. (Windowed sinc) FIR filters
%   are zero phase in passband when shifted by the filters group delay
%   (what firfilt does). Pi phase jumps noticable in the phase reponse
%   reflect a negative frequency response and only occur in the
%   stopband.
%
%               Beta    Max stopband    Max passband    Max passband    Transition width    Mainlobe width
%                       attenuation     deviation       ripple (dB)     (normalized freq)   (normalized rad freq)
%                       (dB)
%   Rectangular         -21             0.0891          1.552           0.9 / m*             4 * pi / m
%   Bartlett            -25             0.0562          0.977           (2.9** / m)          8 * pi / m
%   Hann                -44             0.0063          0.109           3.1 / m              8 * pi / m
%   Hamming             -53             0.0022          0.038           3.3 / m              8 * pi / m
%   Blackman            -74             0.0002          0.003           5.5 / m             12 * pi / m
%   Kaiser      5.653   -60             0.001           0.017           3.6 / m
%   Kaiser      7.857   -80             0.0001          0.002           5.0 / m
%   * m = filter order
%   ** estimate for higher m only
%
% Example:
%   fs = 500; tbw = 2; dev = 0.001;
%   beta = pop_kaiserbeta(dev);
%   m = pop_firwsord('kaiser', fs, tbw, dev);
%   pop_xfirws('srate', fs, 'fcutoff', [1 25], 'ftype', 'bandpass', 'wtype', 'kaiser', 'warg', beta, 'forder', m, 'filename', 'foo.fir')
%
% Author: Andreas Widmann, University of Leipzig, 2011
%
% See also:
%   firfilt, firws, pop_firwsord, pop_kaiserbeta, plotfresp, windows

%123456789012345678901234567890123456789012345678901234567890123456789012

% Copyright (C) 2011 Andreas Widmann, University of Leipzig, widmann@uni-leipzig.de
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [varargout] = pop_xfirws(varargin)

% Pop-up window mode
if nargin < 1

    drawnow;
    ftypes = {'bandpass' 'highpass' 'lowpass' 'bandstop'};
    wtypes = {'rectangular' 'bartlett' 'hann' 'hamming' 'blackman' 'kaiser'};
    uigeom = {[1 0.75 0.75] 1 [1 0.75 0.75] [1 0.75 0.75] 1 [1 0.75 0.75] [1 0.75 0.75] [1 0.75 0.75] 1 [1 0.75 0.75]};
    uilist = {{'Style' 'text' 'String' 'Sampling frequency (Hz):'} ...
              {'Style' 'edit' 'String' '2' 'Tag' 'srateedit'} {} ...
              {} ...
              {'Style' 'text' 'String' 'Cutoff frequency(ies) [hp lp] (-6 dB; Hz):'} ...
              {'Style' 'edit' 'String' '' 'Tag' 'fcutoffedit'} {} ...
              {'Style' 'text' 'String' 'Filter type:'} ...
              {'Style' 'popupmenu' 'String' ftypes 'Tag' 'ftypepop'} {} ...
              {} ...
              {'Style' 'text' 'String' 'Window type:'} ...
              {'Style' 'popupmenu' 'String' wtypes 'Tag' 'wtypepop' 'Value' 5 'Callback' 'temp = {''off'', ''on''}; set(findobj(gcbf, ''-regexp'', ''Tag'', ''^warg''), ''Enable'', temp{double(get(gcbo, ''Value'') == 6) + 1}), set(findobj(gcbf, ''Tag'', ''wargedit''), ''String'', '''')'} {} ...
              {'Style' 'text' 'String' 'Kaiser window beta:' 'Tag' 'wargtext' 'Enable' 'off'} ...
              {'Style' 'edit' 'String' '' 'Tag' 'wargedit' 'Enable' 'off'} ...
              {'Style' 'pushbutton' 'String' 'Estimate' 'Tag' 'wargpush' 'Enable' 'off' 'Callback' @comwarg} ...
              {'Style' 'text' 'String' 'Filter order (mandatory even):'} ...
              {'Style' 'edit' 'String' '' 'Tag' 'forderedit'} ...
              {'Style' 'pushbutton' 'String' 'Estimate' 'Callback' {@comforder, wtypes}} ...
              {'Style' 'edit' 'Tag' 'devedit' 'Visible' 'off'} ...
              {} {} {'Style' 'pushbutton' 'String', 'Plot filter responses' 'Callback' {@comfresp, wtypes, ftypes}}};
    result = inputgui(uigeom, uilist, 'pophelp(''pop_firws'')', 'Filter the data -- pop_firws()');
    if isempty(result), return; end

    Arg = struct;
    Arg.srate = str2double(result{1});
    Arg.fcutoff = str2num(result{2});
    Arg.ftype = ftypes{result{3}};
    Arg.wtype = wtypes{result{4}};
    Arg.warg = str2num(result{5});
    Arg.forder = str2double(result{6});

% Command line mode
else
    Arg = struct(varargin{:});
end

% Sampling rate
if ~isfield(Arg, 'srate') || isempty(Arg.srate) % Use default
    Arg.srate = 2;
end

% Filter order and cutoff frequencies
if ~isfield(Arg, 'fcutoff') || ~isfield(Arg, 'forder') || isempty(Arg.fcutoff) || isempty(Arg.forder)
    error('Not enough input arguments.');
end
firwsArgArray = {Arg.forder sort(Arg.fcutoff / Arg.srate * 2)}; % Sorting and normalization

% Filter type
if ~isfield(Arg, 'ftype') || isempty(Arg.ftype) % Use default
    switch length(Arg.fcutoff)
        case 1
            Arg.ftype = 'lowpass';
        case 2
            Arg.ftype = 'bandpass';
        otherwise
            error('Wrong number of arguments.')
    end
else
    if any(strcmpi(Arg.ftype, {'bandpass' 'bandstop'})) && length(Arg.fcutoff) ~= 2
        error('Not enough input arguments.');
    elseif any(strcmpi(Arg.ftype, {'highpass' 'lowpass'})) && length(Arg.fcutoff) ~= 1
        error('Too many input arguments.');
    end
    switch Arg.ftype
        case 'bandstop'
            firwsArgArray(end + 1) = {'stop'};
        case 'highpass'
            firwsArgArray(end + 1) = {'high'};
    end
end

% Window type
if ~isfield(Arg, 'wtype') || isempty(Arg.wtype) % Use default
    Arg.wtype = 'blackman';
end

% Window parameter
if ~isfield(Arg, 'warg') || isempty(Arg.warg)
    Arg.warg = [];
    firwsArgArray(end + 1) = {windows(Arg.wtype, Arg.forder + 1)};
else
    firwsArgArray(end + 1) = {windows(Arg.wtype, Arg.forder + 1, Arg.warg)};
end

b = firws(firwsArgArray{:});
a = 1;

if nargout == 0 || isfield(Arg, 'filename')

    % Open file
    if ~isfield(Arg, 'filename') || isempty(Arg.filename)
        [Arg.filename Arg.pathname] = uiputfile('*.fir', 'Save filter -- pop_xfirws');
    end
    if ~isfield(Arg, 'pathname') || isempty(Arg.pathname)
        Arg.pathname = '.';
    end
    [fid message] = fopen(fullfile(Arg.pathname, Arg.filename), 'w', 'l');
    if fid == -1
        error(message)
    end

    % Author
    fprintf(fid, '[author]\n');
    fprintf(fid, '%s\n\n', 'pop_xfirws 1.5.1');

    % FIR design
    fprintf(fid, '[fir design]\n');
    fprintf(fid, 'method  %s\n', 'fourier');
    fprintf(fid, 'type    %s\n', Arg.ftype);
    fprintf(fid, 'fsample %f\n', Arg.srate);
    fprintf(fid, 'length  %d\n', Arg.forder + 1);
    fprintf(fid, 'fcrit%d  %f\n', [1:length(Arg.fcutoff); Arg.fcutoff]);
    fprintf(fid, 'window  %s %s\n\n', Arg.wtype, num2str(Arg.warg)); % fprintf bug

    % FIR
    fprintf(fid, '[fir]\n');
    fprintf(fid, '%d\n', Arg.forder + 1);
    fprintf(fid, '% 18.10e\n', b);

    % Close file
    fclose(fid);

end

if nargout > 0
    varargout = {b a};
end

% Callback estimate Kaiser beta
function comwarg(varargin)
    [warg, dev] = pop_kaiserbeta;
    set(findobj(gcbf, 'Tag', 'wargedit'), 'String', warg);
    set(findobj(gcbf, 'Tag', 'devedit'), 'String', dev);

% Callback estimate filter order
function comforder(obj, evt, wtypes)
    srate = str2double(get(findobj(gcbf, 'Tag', 'srateedit'), 'String'));
    wtype = wtypes{get(findobj(gcbf, 'Tag', 'wtypepop'), 'Value')};
    dev = str2double(get(findobj(gcbf, 'Tag', 'devedit'), 'String'));
    [forder, dev] = pop_firwsord(wtype, srate, [], dev);
    set(findobj(gcbf, 'Tag', 'forderedit'), 'String', forder);
    set(findobj(gcbf, 'Tag', 'devedit'), 'String', dev);

% Callback plot filter responses
function comfresp(obj, evt, wtypes, ftypes)
    Arg.srate = str2double(get(findobj(gcbf, 'Tag', 'srateedit'), 'String'));
    Arg.fcutoff = str2num(get(findobj(gcbf, 'Tag', 'fcutoffedit'), 'String'));
    Arg.ftype = ftypes{get(findobj(gcbf, 'Tag', 'ftypepop'), 'Value')};
    Arg.wtype = wtypes{get(findobj(gcbf, 'Tag', 'wtypepop'), 'Value')};
    Arg.warg = str2num(get(findobj(gcbf, 'Tag', 'wargedit'), 'String'));
    Arg.forder = str2double(get(findobj(gcbf, 'Tag', 'forderedit'), 'String'));
    xfirwsArgArray(1, :) = fieldnames(Arg);
    xfirwsArgArray(2, :) = struct2cell(Arg);
    [b a] = pop_xfirws(xfirwsArgArray{:});
    H = findobj('Tag', 'filter responses', 'type', 'figure');
    if ~isempty(H)
        figure(H);
    else
        H = figure;
        set(H, 'color', [.93 .96 1], 'Tag', 'filter responses');
    end
    plotfresp(b, a, [], Arg.srate);


text:
function [] = latex2pdf3d(fname, latex_compiler)
%LATEX2PDF3D    Compile LaTeX code to PDF.
%
% usage
%   LATEX2PDF3D
%   LATEX2PDF3D(fname, latex_compiler)
%
% optional input
%   fname = file name string
%   latex_compiler = 'pdflatex' | 'xelatex'
%
% output
%   If compilation succeeds, then a PDF file is produced.
%
% dependency
%   pdflatex or xelatex should be available in the system path
%
% See also FIG2PDF3D, FIG2LATEX, U3D_IN_LATEX.
%
% File:      latex2pdf3d.m
% Author:    Ioannis Filippidis, jfilippidis@gmail.com
% Date:      2012.06.21
% Language:  MATLAB R2012a
% Purpose:   compile LaTeX file to a PDF
%
% acknowledgment
%   Based on demo_mesh2pdf by Alexandre Gramfort.
%   This can be found on the MATLAb Central File Exchange:
%       http://www.mathworks.com/matlabcentral/fileexchange/25383-matlab-mesh-to-pdf-with-3d-interactive-object
%   and is covered by the BSD License.

% depends
%   {pdflatex | xelatex}, {media9 | movie15}

%% input
if nargin < 1
    fname = 'matfig';
end

if nargin < 2
    latex_compiler = 'xelatex';
end

%% Use pdflatex to generate the pdf
switch latex_compiler
    case 'pdflatex'
        cmd = ['pdflatex --interaction=nonstopmode ', fname, '.tex'];
    case 'xelatex'
        cmd = ['xelatex --interaction=nonstopmode ', fname, '.tex'];
    otherwise
        error('latex:compiler', 'Unknown LaTeX compiler.')
end

status = system(cmd);
if status ~= 0
    error('latex:compile', 'LaTeX compilation failed.')
end

author: Sable


text:
%Bicycle Car-like vehicle class
%
% This concrete class models the kinematics of a car-like vehicle (bicycle
% or Ackerman model) on a plane.  For given steering and velocity inputs it
% updates the true vehicle state and returns noise-corrupted odometry
% readings.
%
% Methods::
%   Bicycle      constructor
%   add_driver   attach a driver object to this vehicle
%   control      generate the control inputs for the vehicle
%   deriv        derivative of state given inputs
%   init         initialize vehicle state
%   f            predict next state based on odometry
%   Fx           Jacobian of f wrt x
%   Fv           Jacobian of f wrt odometry noise
%   update       update the vehicle state
%   run          run for multiple time steps
%   step         move one time step and return noisy odometry
%
% Plotting/display methods::
%   char             convert to string
%   display          display state/parameters in human readable form
%   plot             plot/animate vehicle on current figure
%   plot_xy          plot the true path of the vehicle
%   Vehicle.plotv    plot/animate a pose on current figure
%
% Properties (read/write)::
%   x               true vehicle state: x, y, theta (3x1)
%   V               odometry covariance (2x2)
%   odometry        distance moved in the last interval (2x1)
%   rdim             dimension of the robot (for drawing)
%   L               length of the vehicle (wheelbase)
%   alphalim        steering wheel limit
%   maxspeed        maximum vehicle speed
%   T               sample interval
%   verbose         verbosity
%   x_hist          history of true vehicle state (Nx3)
%   driver          reference to the driver object
%   x0              initial state, restored on init()
%
% Examples::
%
% Odometry covariance (per timstep) is
%       V = diag([0.02, 0.5*pi/180].^2);
% Create a vehicle with this noisy odometry
%       v = Bicycle( 'covar', diag([0.1 0.01].^2 );
% and display its initial state
%       v
% now apply a speed (0.2m/s) and steer angle (0.1rad) for 1 time step
%       odo = v.step(0.2, 0.1)
% where odo is the noisy odometry estimate, and the new true vehicle state
%       v
%
% We can add a driver object
%      v.add_driver( RandomPath(10) )
% which will move the vehicle within the region -10<x<10, -10<y<10 which we
% can see by
%      v.run(1000)
% which shows an animation of the vehicle moving for 1000 time steps
% between randomly selected wayoints.
%
% Notes::
% - Subclasses the MATLAB handle class which means that pass by reference semantics
%   apply.
%
% Reference::
%
%   Robotics, Vision & Control, Chap 6
%   Peter Corke,
%   Springer 2011
%
% See also RandomPath, EKF.

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

classdef Bicycle < Vehicle

    properties
        % state
        L           % length of vehicle

        steermax
        accelmax
        vprev
        steerprev
    end

    methods

        function veh = Bicycle(varargin)
        %Bicycle.Bicycle Vehicle object constructor
        %
        % V = Bicycle(OPTIONS)  creates a Bicycle object with the kinematics of a
        % bicycle (or Ackerman) vehicle.
        %
        % Options::
        % 'steermax',M    Maximu steer angle [rad] (default 0.5)
        % 'accelmax',M    Maximum acceleration [m/s2] (default Inf)
        %--
        % 'covar',C       specify odometry covariance (2x2) (default 0)
        % 'speedmax',S    Maximum speed (default 1m/s)
        % 'L',L           Wheel base (default 1m)
        % 'x0',x0         Initial state (default (0,0,0) )
        % 'dt',T          Time interval (default 0.1)
        % 'rdim',R        Robot size as fraction of plot window (default 0.2)
        % 'verbose'       Be verbose
        %
        % Notes::
        % - The covariance is used by a "hidden" random number generator within the class.
        % - Subclasses the MATLAB handle class which means that pass by reference semantics
        %   apply.
        %
        % Notes::
        % - Subclasses the MATLAB handle class which means that pass by reference semantics
        %   apply.

            veh = veh@Vehicle(varargin{:});

            veh.x = zeros(3,1);

            opt.L = 1;
            opt.steermax = 0.5;
            opt.accelmax = Inf;

            veh = tb_optparse(opt, veh.options, veh);
            veh.vprev = 0;
            veh.x = veh.x0;
        end

        function xnext = f(veh, x, odo, w)
            %Bicycle.f Predict next state based on odometry
            %
            % XN = V.f(X, ODO) is the predicted next state XN (1x3) based on current
            % state X (1x3) and odometry ODO (1x2) = [distance, heading_change].
            %
            % XN = V.f(X, ODO, W) as above but with odometry noise W.
            %
            % Notes::
            % - Supports vectorized operation where X and XN (Nx3).
            if nargin < 4
                w = [0 0];
            end

            dd = odo(1) + w(1); dth = odo(2) + w(2);

            % straightforward code:
            % thp = x(3) + dth;
            % xnext = zeros(1,3);
            % xnext(1) = x(1) + (dd + w(1))*cos(thp);
            % xnext(2) = x(2) + (dd + w(1))*sin(thp);
            % xnext(3) = x(3) + dth + w(2);
            %
            % vectorized code:

            thp = x(:,3) + dth;
            %xnext = x + [(dd+w(1))*cos(thp)  (dd+w(1))*sin(thp) ones(size(x,1),1)*dth+w(2)];
            xnext = x + [dd*cos(thp)  dd*sin(thp) ones(size(x,1),1)*dth];
        end

        function [dx,u] = deriv(veh, t, x, u)
            %Bicycle.deriv  Time derivative of state
            %
            % DX = V.deriv(T, X, U) is the time derivative of state (3x1) at the state
            % X (3x1) with input U (2x1).
            %
            % Notes::
            % - The parameter T is ignored but  called from a continuous time integrator such as ode45 or
            %   Simulink.

            % implement acceleration limit if required
            if ~isinf(veh.accelmax)
                if (u(1) - veh.vprev)/veh.dt > veh.accelmax
                    u(1) = veh.vprev + veh.accelmax * veh.dt;
                elseif (u(1) - veh.vprev)/veh.dt < -veh.accelmax
                    u(1) = veh.vprev - veh.accelmax * veh.dt;
                end
                veh.vprev = u(1);
            end

            % implement speed and steer angle limits
            u(1) = min(veh.speedmax, max(u(1), -veh.speedmax));
            u(2) = min(veh.steermax, max(u(2), -veh.steermax));

            % compute the derivative
            dx = zeros(3,1);
            dx(1) = u(1)*cos(x(3));
            dx(2) = u(1)*sin(x(3));
            dx(3) = u(1)/veh.L * tan(u(2));
        end

        function odo = update(veh, u)
            %Bicycle.update Update the vehicle state
            %
            % ODO = V.update(U) is the true odometry value for
            % motion with U=[speed,steer].
            %
            % Notes::
            % - Appends new state to state history property x_hist.
            % - Odometry is also saved as property odometry.

            % update the state
            dx = veh.dt * veh.deriv([], veh.x, u);
            veh.x = veh.x + dx;

            % compute and save the odometry
            odo = [ norm(dx(1:2)) dx(3) ];
            veh.odometry = odo;

            veh.x_hist = [veh.x_hist; veh.x'];   % maintain history
        end

        function J = Fx(veh, x, odo)
        %Bicycle.Fx  Jacobian df/dx
        %
        % J = V.Fx(X, ODO) is the Jacobian df/dx (3x3) at the state X, for
        % odometry input ODO (1x2) = [distance, heading_change].
        %
        % See also Bicycle.f, Vehicle.Fv.
            dd = odo(1); dth = odo(2);
            thp = x(3) + dth;

            J = [
                1   0   -dd*sin(thp)
                0   1   dd*cos(thp)
                0   0   1
                ];
        end

        function J = Fv(veh, x, odo)
            %Bicycle.Fv  Jacobian df/dv
            %
            % J = V.Fv(X, ODO) is the Jacobian df/dv (3x2) at the state X, for
            % odometry input ODO (1x2) = [distance, heading_change].
            %
            % See also Bicycle.F, Vehicle.Fx.
            dd = odo(1); dth = odo(2);
            thp = x(3);

            J = [
                cos(thp)    0
                sin(thp)    0
                0           1
                ];
        end

        function s = char(veh)
        %Bicycle.char Convert to a string
        %
        % s = V.char() is a string showing vehicle parameters and state in
        % a compact human readable format.
        %
        % See also Bicycle.display.

            ss = char@Vehicle(veh);

            s = 'Bicycle object';
            s = char(s, sprintf('  L=%g, steer.max=%g, accel.max=%g', veh.L, veh.steermax, veh.accelmax));
            s = char(s, ss);
        end
    end % method

end % classdef

author: petercorke

repo: robotics-toolbox-matlab


text:
function varargout = lbp(varargin)
% VL_LBP  Local Binary Patterns
%   F = VL_LBP(IM, CELLSIZE) computes the Local Binary Pattern (LBP)
%   features for image I.
%
%   IM is divided in cells of size CELLSIZE. F is a three-dimensional
%   array containing one histograms of quantized LBP features per
%   cell. The witdh of F is FLOOR(WIDTH/CELLSIZE), where WIDTH is the
%   width of the image. The same for the height. The third dimension
%   is 58.
%
%   See also: VL_HELP().
[varargout{1:nargout}] = vl_lbp(varargin{:});

author: jianxiongxiao

repo: ProfXkit


text:
function c = phi1(costgradf,w,param,mu)
% compute c=cost+mu*dual error
D    = param.D;
b    = param.b;
dim  = param.dim;
h    = param.h;
m    = param.m;
n  = dim(1);
nx = dim(2);
ny = dim(3);
nz = dim(4);
nt = dim(5);
hx = h(1);
hy = h(2);
hz = h(3);
ht = h(4);

px              = w(1:n*(nx-1)*ny*nz*nt);
py              = w((n*(nx-1)*ny*nt*nz+1):(n*(nx-1)*ny*nz*nt+n*nx*(ny-1)*nz*nt));
pz              = w((n*(nx-1)*ny*nz*nt+n*nx*(ny-1)*nz*nt)+1:(n*(nx-1)*ny*nz*nt+n*nx*(ny-1)*nz*nt+n*nx*ny*(nz-1)*nt));
rho             = w((n*(nx-1)*ny*nz*nt+n*nx*(ny-1)*nz*nt+n*nx*ny*(nz-1)*nt)+1:...
                       (n*(nx-1)*ny*nz*nt+n*nx*(ny-1)*nz*nt+n*nx*ny*(nz-1)*nt+n*nx*ny*nz*(nt-1)));
u               = w((end-m*nx*ny*nz*nt+1):end);
c    = costgradf(px,py,pz,rho,u,param)/hx/hy/hz/ht;
c    = c + mu*sum(abs(D*w-b(:)));
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function obj=MOxUnitTestOutcome(test_,duration)
    s=struct();
    s.test=test_;
    s.duration=duration;
    obj=class(s,'MOxUnitTestOutcome');

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
% this function calculate the PSNR

function [PeakSNR, Mean2err]=psnr(OriginalImage, DegradedImage)
OriginalImage=double(OriginalImage);
DegradedImage=double(DegradedImage);
[N,M] = size(OriginalImage);
Imax = max(max(DegradedImage));
SumOfDiff2 = sum(sum((OriginalImage-DegradedImage).*(OriginalImage-DegradedImage)));
Mean2err=SumOfDiff2./(M*N);
sdf=Imax^2./(Mean2err);
if sdf ==0
    sdf=1;
end

PeakSNR = 10*log10(sdf);

repo: mcbench-benchmarks


text:
function [class,siz]=ChemSiz(Sin)

%% Give Labels and size for plotting chemical size distributions
%% Label definintions:
%% 1=OC, 2=ECOCIn, 3=ECOC, 4=InOC, 5=NoID
%% RC Moffet, 2010
siz=Sin.Size;
class=zeros(1,length(Sin.PartLabel));
for i = 1:length(Sin.PartLabel)
    NoIdidx=strmatch('NoID',Sin.PartLabel{i},'exact');
    OCidx=strmatch('OC',Sin.PartLabel{i},'exact');
    OCAidx=strmatch('OCsp2',Sin.PartLabel{i},'exact');
    ECidx=findstr(Sin.PartLabel{i},'EC');
    Inidx=findstr(Sin.PartLabel{i},'In');
    Kidx=findstr(Sin.PartLabel{i},'K');
    if ~isempty(OCidx) || ~isempty(OCAidx) %% OC
        class(i)=1;
    elseif ~isempty(ECidx) && (~isempty(Inidx) || ~isempty(Kidx)) %% ECOCIn
        class(i)=2;
    elseif ~isempty(ECidx) && (isempty(Inidx) || isempty(Kidx)) %% ECOC
        class(i)=3;
    elseif (~isempty(Inidx) || ~isempty(Kidx)) && isempty(ECidx) %% InOC
        class(i)=4;
    elseif ~isempty(NoIdidx)  %% NoID
        class(i)=5;
    end
    clear OCidx ECidx Inidx NoIdidx;
end


text:
function cVal = prtUtilFindPrtActionsAndConvertToStructures(cVal)

if isa(cVal,'prtAction')
    cVal = cVal.toStructure;
elseif iscell(cVal)
    % Current element is a cell that may contain prtActions
    for iCell = 1:numel(cVal)
        cVal{iCell} = prtUtilFindPrtActionsAndConvertToStructures(cVal{iCell});
    end
elseif isstruct(cVal)
    subFieldNames = fieldnames(cVal);
    for iField = 1:length(subFieldNames)
        cVal.(subFieldNames{iField}) =  prtUtilFindPrtActionsAndConvertToStructures(cVal.(subFieldNames{iField}));
    end
end

author: covartech

repo: PRT


text:
function [output]=rgbcompression(f)
[F4]= compressimage(f);
F4=uint8(F4);

F4=ycbcr2rgb(F4);
output=F4;
figure;
subplot(1,2,1),imshow(f);
subplot(1,2,2),imshow(F4);