text:
%%***************************************************************************
%% blkeig: compute eigenvalue decomposition of a cell array
%%         whose contents are square matrices or the diagonal
%%         of a diagonal matrix.
%%
%% [d,V] = blkeig(blk,X);
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%***************************************************************************

function [d,V] = blkeig(blk,X)

spdensity = 0.5;

if ~iscell(X);
    if strcmp(blk{1},'s');
        blktmp = blk{2};
        if (length(blktmp) == 1);
            if (nargout == 1);
                d = eig(full(X));
            elseif (nargout == 2);
                [V,d] = eig(full(X));
                d = diag(d);
            end
        else
            if (nargout == 2);
                V = sparse(length(X),length(X));
            end
            d = zeros(sum(blktmp),1);
            xx = mexsvec(blk,X,0);
            blktmp2 = blktmp.*(blktmp+1)/2;
            s2 = [0, cumsum(blktmp2)];
            blksub{1,1} = 's'; blksub{1,2} = 0;
            s = [0, cumsum(blktmp)];
            for i = 1:length(blktmp)
                pos = s(i)+1 : s(i+1);
                blksub{2} = blktmp(i);
                Xsub = mexsmat(blksub,xx(s2(i)+1:s2(i+1)),0);
                if (nargout == 1);
                    lam = eig(Xsub);
                elseif (nargout == 2);
                    [evec,lam] = eig(Xsub);
                    lam = diag(lam);
                    V(pos,pos) = sparse(evec); %#ok
                end
                d(pos,1) = lam;
            end
        end
        n2 = sum(blktmp.*blktmp);
        if (nargout == 2);
            if (nnz(V) <= spdensity*n2);
                V = sparse(V);
            else
                V = full(V);
            end
        end
    elseif strcmp(blk{1},'l');
        if (nargout == 2);
            V = ones(size(X)); d = X;
        elseif (nargout == 1);
            d = X;
        end
    end
else
    if (nargout == 2);
        V = cell(size(X));  d = cell(size(X));
        for p = 1:size(blk,1);
            [d{p},V{p}] = blkeig(blk(p,:),X{p});
        end
    elseif (nargout == 1);
        d = cell(size(X));
        for p = 1:size(blk,1);
            d{p} = blkeig(blk(p,:),X{p});
        end
    end
end
%%***************************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
% z-Transform properties

%	Right shift of x[n]

n=-3:3;
x=0.8.^n;

xminus3=x(1)
xminus2=x(2)
xminus1=x(3)

syms n z
xn1=0.8^(n-1);
Left=ztrans(xn1,z) ;
simplify(Left)

x=0.8^n;
X=ztrans(x,z);
Right=z^-1*X +xminus1;
simplify(Right)

xn2=0.8^(n-2);
Left=ztrans(xn2,z) ;
simplify(Left)

Right=z^-2*X+xminus2+z^-1 *xminus1;
simplify(Right)

author: Sable

repo: mcbench-benchmarks


text:
classdef MixtureTheoryHomogenizer < handle

    properties (Access = public)
        Ch
    end

    properties (Access = private)
        C1
        C0
        angle
        Ex
        Ey
        nu_xy
        nu_yx
        mu
        dir
        ChHorizontal
        Vfrac
    end

    methods (Access = public)

        function obj = MixtureTheoryHomogenizer(C1,C0,Dir,angle,Vfrac)
            obj.init(C1,C0,Dir,angle,Vfrac)
            obj.computeOrthotropicProperties()
            obj.computeHorizontalHomogenizedTensor()
            obj.rotateHorizontalHomogenizedTensor()
        end

    end

    methods (Access = private)

        function init(obj,C1,C0,Dir,angle,Vfrac)
            obj.C1 = C1;
            obj.C0 = C0;
            obj.angle = angle;
            obj.dir = Dir;
            obj.Vfrac = Vfrac;
        end

        function computeOrthotropicProperties(obj)
            E1 = obj.C1.getYoung();
            E0 = obj.C0.getYoung();
            nu1 = obj.C1.getPoisson();
            nu0 = obj.C0.getPoisson();
            mu1 = obj.C1.getMu();
            mu0 = obj.C0.getMu();
            Vol = obj.Vfrac;

            obj.Ex = obj.serialize(E1,E0,Vol);
            obj.Ey = obj.parelalize(E1,E0,Vol);
            obj.nu_xy = obj.serialize(nu1,nu0,Vol);
            obj.nu_yx = obj.nu_xy*obj.Ey/obj.Ex;
            obj.mu = obj.parelalize(mu1,mu0,Vol);
        end

        function computeHorizontalHomogenizedTensor(obj)
            E1    = obj.Ex;
            E2    = obj.Ey;
            nu_12 = obj.nu_xy;
            nu_21 = obj.nu_yx;
            Mu    = obj.mu;

            C = zeros(3,3);
            C(1,1) = E1/(1-nu_12*nu_21);
            C(2,2) = E2/(1-nu_12*nu_21);
            C(1,2) = E1*nu_21/(1-nu_12*nu_21);
            C(2,1) = E2*nu_12/(1-nu_12*nu_21);
            C(3,3) = Mu;
            obj.ChHorizontal  = C;
        end

        function rotateHorizontalHomogenizedTensor(obj)
            d = obj.dir;
            a = obj.angle;
            ChHor = obj.ChHorizontal;
            C = SymmetricFourthOrderPlaneStressVoigtTensor();
            C.setValue(ChHor);
            obj.Ch = Rotator.rotate(C,a,d);
        end

    end

    methods (Access = private, Static)
        function fSerial = serialize(f1,f0,rho)
            fSerial = rho*f1 + (1-rho)*f0;
        end

        function fParalel = parelalize(f1,f0,rho)
            fParalel = 1/(rho/f1 + (1-rho)/f0);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function pass = test_sum3(pref)
% Test file for @chebfun3/sum3.m.

% Obtain preferences.
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = 1e4 * pref.cheb3Prefs.chebfun3eps;

% Check empty case
pass(1) = isempty(sum3(chebfun3()));

% Check constant
f = chebfun3(@(x,y,z) 1);
pass(2) = abs(sum3(f) - 8) < tol;

% Runge function
f = cheb.gallery3('runge');
pass(3) = abs(sum3(f) - 4.28685406230184188268) < tol;

% Different domains
f = chebfun3(@(x,y,z) x, [0, 1, 0, 2, 0, 3]);
pass(4) = abs(sum3(f) - 3) < tol;

f = chebfun3(@(x,y,z) y, [0, 1, 0, 2, 0, 3]);
pass(5) = abs(sum3(f) - 6) < tol;

f = chebfun3(@(x,y,z) z, [0, 1, 0, 2, 0, 3]);
pass(6) = abs(sum3(f) - 9) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function wishart_test04 ( )

%*****************************************************************************80
%
%% WISHART_TEST04 demonstrates the Wishart sampling function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 July 2013
%
%  Author:
%
%    John Burkardt
%

%
%  Access the PDFLIB and RNGLIB libraries.
%
  addpath ( '../pdflib' );
  addpath ( '../rnglib' );
%
%  Initialize the RNGLIB library.
%
  initialize ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'WISHART_TEST04:\n' );
  fprintf ( 1, '  We can compute sample Wishart matrices by:\n' );
  fprintf ( 1, '    W = wishart_sample ( n, df, sigma );\n' );
%
%  Set the parameters and call.
%
  n = 5;
  df = 8;
  sigma = eye ( 5, 5 );
  w = wishart_sample ( n, df, sigma );
  r8mat_print ( n, n, w, '  wishart_sample ( 5, 8, Identity ):' );
%
%  Calling again yields a new matrix.
%
  w = wishart_sample ( n, df, sigma );
  r8mat_print ( n, n, w, '  wishart_sample ( 5, 8, Identity ):' );
%
%  Try a diagonal matrix.
%
  sigma = diag ( [1.0, 2.0, 3.0, 4.0, 5.0 ] );
  w = wishart_sample ( n, df, sigma );
  r8mat_print ( n, n, w, '  wishart_sample ( 5, 8, diag(1,2,3,4,5) ):' );
%
%  Try a smaller matrix.  Sigma must be positive definite symmetric.
%
  n = 3;
  df = 3;
  r = [ 5.0, 1.0, 3.0; ...
        0.0, 4.0, 2.0; ...
        0.0, 0.0, 6.0 ];
  sigma = r' * r;
  r8mat_print ( n, n, sigma, '  Set covariance SIGMA:' );
  w = wishart_sample ( n, df, sigma );
  r8mat_print ( n, n, w, '  wishart_sample ( 3, 3, sigma ):' );
%
%  What is the eigendecomposition of this matrix?
%
  [ v, lambda ] = eigs ( w );
  r8mat_print ( n, n, v, '  Eigenvectors of previous matrix:' );
  r8mat_print ( n, n, lambda, '  Eigenvalues of previous matrix:' );

  rmpath ( '../pdflib' );
  rmpath ( '../rnglib' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ i, j, v ] = find( x )

%Disciplined convex/geometric programming information for FIND:
%   When used in CVX models, FIND cannot deduce the numeric value of a
%   CVX variable. So it returns the indices of elements that are
%   "structurally" nonzero; i.e., that are not *identically* zero.
%   This is similar to the concept of structural nonzeros in sparse
%   matrix factorizations. To illustrate the distinction, consider:
%      variable x(3);
%      x(2) == 0;
%      y = [x(1);0;x(3)];
%   In this case, FIND(x) returns [1;2;3] even though x(2) has been set
%   to zero in an equality constraint. (After all, the overall model may
%   be infeasible, in which case x(2) is arguably not zero.) However,
%   FIND(y) will return [1;3], because y(2) is identically zero.
%
%   When X is a CVX variable, the first two outputs of [I,J,V]=FIND(X)
%   are constant, and the third is a CVX variable containing the
%   structural nonzeros of X.
%
%   FIND(X) places no convexity restrictions on its argument.

ndxs = find( any( x.basis_, 1 ) );
ndxs = ndxs( : );

if nargout > 1,
    i = ndxs - 1;
    j = floor( i / x.size_(1) ) + 1;
    i = rem( i, x.size_(1) ) + 1;
else
	i = ndxs;
end

if nargout > 2,
    v = reshape( cvx_subsref( x, ndxs ), length( ndxs ) );
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function y = evaluate(this, x)
%EVALUATE   Evaluate robust function.
%   EVALUATE(THIS, X) evaluates the robust function THIS at point(s) X.
%
%   This is a member function of the class 'robust_function'.
%
%   Author:  Stefan Roth, Department of Computer Science, TU Darmstadt
%   Contact: sroth@cs.tu-darmstadt.de
%   $Date:  $
%   $Revision: $

% Copyright 2004-2007, Brown University, Providence, RI. USA
% Copyright 2007-2010 TU Darmstadt, Darmstadt, Germany.
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

  y = feval(this.type, x, this.param, 0);

author: kristinbranson

repo: JAABA


text:
function AngleAvg = AverageAnglesRadians(InA)
%       Robert Goldstein
%       Schepens Eye Research Institute
%       robert.goldstein@schepens.harvard.edu
    theNumerator = sum(sin(InA));
    theDenominator = sum(cos(InA));
    if(theDenominator ==0)
        AngleAvg = sign(theNumerator)*pi/2;
    else
        AngleAvg = atan2(theNumerator,theDenominator);
    end

end


function yprime = decay1(t,y)
%DECAY1 Calculates the decay rates of Thorium 227 and Radium 223.
% Function DECAY1 Calculates the rates of change of Thorium 227
% and Radium 223 (yprime) for a given current concentration y.

% Define variables:
%   t         -- Time (in days)
%   y         -- Vector of current concentrations
%
%  Record of revisions:
%      Date       Programmer          Description of change
%      ====       ==========          =====================
%    03/15/07    S. J. Chapman        Original code

% Set decay constants.
lambda_th = 0.03710636;
lambda_ra = 0.0606428;

% Calculate rates of decay
yprime = zeros(2,1);
yprime(1) = -lambda_th * y(1);
yprime(2) = -lambda_ra * y(2) + lambda_th * y(1);


text:
function [image,depth, info]=point2render(coordinate, rgb, image, depth, info)
% render a point cloud to an image
% note that we assume the coordinate is the camera coordinates
% for example, if you have a camera matrix P
% you should input P*X as your coordinate
% coordinate and rgb are 3*n matrix
% rgb is double color from [0,1]

%{
% if you have many point cloud to accumulate. You can use it like this:
for frameID=1:N
    if frameID == 1
        [image,depth, info]=point2render(XYZworld, RGB);
    else
        [image,depth]=point2render(XYZworld, RGB, image, depth, info);
    end
end
%}

if ~exist('info','var')
    info.range.minX = -10;
    info.range.maxX =  10;
    info.range.minY = -10;
    info.range.maxY =  10;
    info.unit = 0.01;
end

if ~exist('image','var') || ~exist('depth','var')
    sizeX = ceil((info.range.maxX - info.range.minX)/info.unit);
    sizeY = ceil((info.range.maxY - info.range.minY)/info.unit);
    image = ones(sizeX,sizeY,3);
    depth = -Inf(sizeX,sizeY);
end
sizeXY = numel(depth);

coordinateX = (coordinate(1,:) - info.range.minX)/info.unit;
coordinateY = (coordinate(2,:) - info.range.minY)/info.unit;

coordinateXi = round(coordinateX);
coordinateYi = round(coordinateY);

valid = find((1<= coordinateXi) & (coordinateXi <= size(depth,1)) & (1<= coordinateYi) & (coordinateYi <= size(depth,2)));

index = sub2ind(size(depth),coordinateXi(valid),coordinateYi(valid));

toColor = depth(index) < coordinate(3,valid);

valid = valid(toColor);
index = index(toColor);

[~, swapOrder]= sort(coordinate(3,valid));
valid = valid(swapOrder);
index = index(swapOrder);

image(index)          = rgb(1,valid);
image(index+sizeXY)   = rgb(2,valid);
image(index+sizeXY*2) = rgb(3,valid);
depth(index)          = coordinate(3,valid);

author: jianxiongxiao

repo: ProfXkit


text:
function x = buncompress( xR, x, sx )
error( nargchk( 2, 3, nargin ) );

if ~isa( xR, 'double' ) && ~isa( xR, 'sparse' ),
    error( 'First argument must be a structure matrix.' );
elseif size( x.basis_, 2 ) ~= size( xR, 1 ),
    error( 'Structure matrix incompatible with vector.' );
elseif nargin < 3 || isempty( sx ),
    sx = size( xR, 2 );
elseif ~cvx_check_dimlist( sx, false ),
    error( 'Third argument must be a size matrix.' );
elseif prod( sx ) ~= isempty( xR ) * prod( x.size_ ) + ~isempty( xR ) * size( xR, 2 ),
    error( 'Incompatible size matrix.' );
end

if isempty( xR ),
    x = cvx( sx, x.basis_ );
else
    x = cvx( sx, x.basis_ * xR );
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.


text:
%% process CW vectors

% Input file
inputFilebase = '../data/mpqa/mpqa';
inputpos = [inputFilebase '.pos'];
inputneg = [inputFilebase '.neg'];
outputname = '../data/mpqa/allSStrMPQA.mat';

%Create hash table for the word list
fid = fopen(inputpos, 'r');
linespos = textscan(fid, '%s', 'delimiter', '\n', 'bufsize', 100000);
fclose(fid);
linespos=linespos{1};

fid = fopen(inputneg, 'r');
linesneg = textscan(fid, '%s', 'delimiter', '\n', 'bufsize', 100000);
fclose(fid);
linesneg=linesneg{1};

allSStrMPQA = {};
for i = 1:length(linespos)
    if mod(i, 1000) == 0;
        fprintf('processing line %d of %d\n', i, length(lines));
    end

    line = regexp(strtrim(linespos{i}),' ','split');

    allSStrMPQA{end+1} = line;

end

for i = 1:length(linesneg)
    if mod(i, 100) == 0;
        fprintf('processing line %d of %d\n', i, length(lines));
    end

    line = regexp(strtrim(linesneg{i}),' ','split');
    allSStrMPQA{end+1} = line;
end
labels = false(1,length(allSStrMPQA));
labels(1:length(linespos)) = true;
save(outputname, 'allSStrMPQA', 'labels');

author: sidaw

repo: nbsvm

save_path: github-repos/MATLAB/sidaw-nbsvm


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

%*****************************************************************************80
%
%% GAMMA_INC_P_VALUES: values of the normalized incomplete Gamma function P(A,X)
%
%  Discussion:
%
%    The (normalized) incomplete Gamma function is defined as:
%
%      P(A,X) = 1/Gamma(A) * Integral ( 0 <= T <= X ) T^(A-1) * exp(-T) dT.
%
%    With this definition, for all A and X,
%
%      0 <= P(A,X) <= 1
%
%    and
%
%      P(A,oo) = 1.0
%
%    In Mathematica, the function can be evaluated by:
%
%      1 - GammaRegularized[A,X]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 April 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real A, the parameter of the function.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 20;

  a_vec = [ ...
     0.10E+00, ...
     0.10E+00, ...
     0.10E+00, ...
     0.50E+00, ...
     0.50E+00, ...
     0.50E+00, ...
     0.10E+01, ...
     0.10E+01, ...
     0.10E+01, ...
     0.11E+01, ...
     0.11E+01, ...
     0.11E+01, ...
     0.20E+01, ...
     0.20E+01, ...
     0.20E+01, ...
     0.60E+01, ...
     0.60E+01, ...
     0.11E+02, ...
     0.26E+02, ...
     0.41E+02  ];

  fx_vec = [ ...
     0.7382350532339351E+00, ...
     0.9083579897300343E+00, ...
     0.9886559833621947E+00, ...
     0.3014646416966613E+00, ...
     0.7793286380801532E+00, ...
     0.9918490284064973E+00, ...
     0.9516258196404043E-01, ...
     0.6321205588285577E+00, ...
     0.9932620530009145E+00, ...
     0.7205974576054322E-01, ...
     0.5891809618706485E+00, ...
     0.9915368159845525E+00, ...
     0.1018582711118352E-01, ...
     0.4421745996289254E+00, ...
     0.9927049442755639E+00, ...
     0.4202103819530612E-01, ...
     0.9796589705830716E+00, ...
     0.9226039842296429E+00, ...
     0.4470785799755852E+00, ...
     0.7444549220718699E+00 ];

  x_vec = [ ...
     0.30E-01, ...
     0.30E+00, ...
     0.15E+01, ...
     0.75E-01, ...
     0.75E+00, ...
     0.35E+01, ...
     0.10E+00, ...
     0.10E+01, ...
     0.50E+01, ...
     0.10E+00, ...
     0.10E+01, ...
     0.50E+01, ...
     0.15E+00, ...
     0.15E+01, ...
     0.70E+01, ...
     0.25E+01, ...
     0.12E+02, ...
     0.16E+02, ...
     0.25E+02, ...
     0.45E+02 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

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

  return
end


text:
classdef OrbitStateEnum < matlab.mixin.SetGet
    %OrbitStateEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        KeplerianOrbit('Keplerian Orbit','KeplerianOrbitStateModel', ...
                       {'SMA','Eccentricity','Inclination','RAAN','Arg Peri','True Aomaly'}, ...
                       {'km','','deg','deg','deg','deg'}, ...
                       'KeplerianOrbitVariable', ...
                       'Central Body', ...
                       'The body around which this orbit is defined.')
        BodyFixed('Body Fixed Orbit','BodyFixedOrbitStateModel', ...
                  {'Latitude','Longitude','Altitude','Velocity Az','Velocity El','Velocity Mag'}, ...
                  {'degN','degE','km','deg','deg','km/s'}, ...
                  'BodyFixedOrbitVariable', ...
                  'Central Body', ...
                  'The body upon which this body-fixed state is defined.  If launching from the surface, the body being launched from.')
        CR3BP('Primary-centered Circular Restricted 3 Body Problem Orbit','CR3BPOrbitStateModel', ...
                  {'Position (X)','Position (Y)','Position (Z)','Velocity (X)','Velocity (Y)','Velocity (Z)'}, ...
                  {'km','km','km','km/s','km/s','km/s'}, ...
                  'BodyFixedOrbitVariable', ...
                  'Secon. Body', ...
                  'The secondary body in the 3 body rotating frame.  The actual orbit will be defined w.r.t. the parent of this body.  E.G. if you want a Kerbin-Mun rotating frame, select "Mun" here.')
    end

    properties
        name char = ''
        class char = ''
        elemNames(6,1) cell
        unitNames(6,1) cell
        varClass char = '';
        cbLabel char = 'Central Body';
        cbTooltip char = '';
    end

    methods
        function obj = OrbitStateEnum(name,class,elemNames,unitNames,varClass,cbLabel,cbTooltip)
            obj.name = name;
            obj.class = class;
            obj.elemNames = elemNames;
            obj.unitNames = unitNames;
            obj.varClass = varClass;
            obj.cbLabel = cbLabel;
            obj.cbTooltip = cbTooltip;
        end
    end

    methods(Static)
        function listBoxStr = getListBoxStr()
            m = enumeration('OrbitStateEnum');
            listBoxStr = {m.name};
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('OrbitStateEnum');
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [ind, enum] = getIndForClass(class)
            m = enumeration('OrbitStateEnum');
            ind = find(ismember({m.class},class),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function determ = poisson_determinant ( nrow, ncol )

%*****************************************************************************80
%
%% POISSON_DETERMINANT returns the determinant of the POISSON matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NROW, NCOL, the number of rows and columns
%    in the grid.
%
%    Output, real DETERM, the determinant.
%
  cr = zeros ( nrow, 1 );
  for i = 1 : nrow
    angle = i * pi / ( nrow + 1 );
    cr(i) = cos ( angle );
  end

  cc = zeros ( ncol, 1 );
  for i = 1 : ncol
    angle = i * pi / ( ncol + 1 );
    cc(i) = cos ( angle );
  end

  determ = 1.0;

  for i = 1 : nrow
    for j = 1 : ncol
      determ = determ * ( 4.0 - 2.0 * cr(i) - 2.0 * cc(j) );
    end
  end

  return
end


text:
function [ Nreg ] = eigRegularizer( N, k )

% -----------------------------------------------------------------
% matrizer regularization by checking eigenvalues
% -----------------------------------------------------------------
%
% N:   input square matrix
% k:   required condition number

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

[V, D] = eig(N); 					 % eigenvectors and eigenvalues
e = diag(D);                         % eigenvalues (ascending order)

% k = 1e6;                           % condition number
% pos = find( e < (e(end)/k) );
% e(pos) = e(end)/k;                 % regularized eigenvalues
% Nreg = V*diag(e)*V';               % regularized matrix

pos = e < (max(e) / k);
e(pos) = max(e) / k;                 % regularized eigenvalues
Nreg = V * diag(e) * V';             % regularized matrix
%Nreg = real(V*diag(e)*V');          % regularized matrix

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

% figure; imagesc(N); colorbar
% figure; imagesc(Nreg); colorbar
% figure; imagesc(N-Nreg); colorbar

author: goGPS-Project


text:
% Minimize sidelobe level of an array with arbitrary 2-D geometry
% "Convex optimization examples" lecture notes (EE364) by S. Boyd
% "Antenna array pattern synthesis via convex optimization"
% by H. Lebret and S. Boyd
% (figures are generated)
%
% Designs an antenna array such that:
% - it minimizes sidelobe level outside the beamwidth of the pattern
% - it has a unit sensitivity at some target direction
% - it has nulls (zero sensitivity) at specified direction(s) (optional)
%
% This is a convex problem (after sampling it can be formulated as an SOCP).
%
%   minimize   max |y(theta)|     for theta outside the beam
%       s.t.   y(theta_tar) = 1
%              y(theta_null) = 0  (optional)
%
% where y is the antenna array gain pattern (complex function) and
% variables are w (antenna array weights or shading coefficients).
% Gain pattern is a linear function of w: y(theta) = w'*a(theta)
% for some a(theta) describing antenna array configuration and specs.
%
% Written for CVX by Almir Mutapcic 02/02/06

% select array geometry
ARRAY_GEOMETRY = '2D_RANDOM';
% ARRAY_GEOMETRY = '1D_UNIFORM_LINE';
% ARRAY_GEOMETRY = '2D_UNIFORM_LATTICE';

% select if the optimal array pattern should enforce nulls or not
HAS_NULLS = 0; % HAS_NULLS = 1;

%********************************************************************
% problem specs
%********************************************************************
lambda = 1;           % wavelength
theta_tar = 60;       % target direction (should be an integer -- discretization)
half_beamwidth = 10;  % half beamwidth around the target direction

% angles where we want nulls (optional)
if HAS_NULLS
    theta_nulls = [95 110 120 140 225];
end

%********************************************************************
% random array of n antenna elements
%********************************************************************
if strcmp( ARRAY_GEOMETRY, '2D_RANDOM' )
  % set random seed to repeat experiments
  rand('state',0);

  % (uniformly distributed on [0,L]-by-[0,L] square)
  n = 40;
  L = 5;
  loc = L*rand(n,2);
  angleRange = 360;

%********************************************************************
% uniform 1D array with n elements with inter-element spacing d
%********************************************************************
elseif strcmp( ARRAY_GEOMETRY, '1D_UNIFORM_LINE' )
  % (unifrom array on a line)
  n = 30;
  d = 0.45*lambda;
  loc = [d*(0:n-1)' zeros(n,1)];
  angleRange = 180;

%********************************************************************
% uniform 2D array with m-by-m element with d spacing
%********************************************************************
elseif strcmp( ARRAY_GEOMETRY, '2D_UNIFORM_LATTICE' )
  m = 6; n = m^2;
  d = 0.45*lambda;

  loc = zeros(n,2);
  for x = 0:m-1
    for y = 0:m-1
      loc(m*y+x+1,:) = [x y];
    end
  end
  loc = loc*d;
  angleRange = 360;

else
  error('Undefined array geometry')
end

%********************************************************************
% construct optimization data
%********************************************************************
% build matrix A that relates w and y(theta), ie, y = A*w
theta = (1:angleRange)';
A = kron(cos(pi*theta/180), loc(:,1)') + kron(sin(pi*theta/180), loc(:,2)');
A = exp(2*pi*1i/lambda*A);

% target constraint matrix
[diff_closest, ind_closest] = min( abs(theta - theta_tar) );
Atar = A(ind_closest,:);

% nulls constraint matrix
if HAS_NULLS
  Anull = []; ind_nulls = [];
  for k = 1:length(theta_nulls)
    [diff_closest, ind_closest] = min( abs(theta - theta_nulls(k)) );
    Anull = [Anull; A(ind_closest,:)]; %#ok
    ind_nulls = [ind_nulls ind_closest]; %#ok
  end
end

% stopband constraint matrix
ind = find(theta <= (theta_tar-half_beamwidth) | ...
           theta >= (theta_tar+half_beamwidth) );
if HAS_NULLS,
    ind = setdiff(ind,ind_nulls);
end;
As = A(ind,:);

%********************************************************************
% optimization problem
%********************************************************************
cvx_begin
  variable w(n) complex
  minimize( max( abs(As*w) ) )
  subject to
    Atar*w == 1;   %#ok target constraint
    if HAS_NULLS   % nulls constraints
      Anull*w == 0; %#ok
    end
cvx_end

% check if problem was successfully solved
disp(['Problem is ' cvx_status])
if ~strfind(cvx_status,'Solved')
  return
end

min_sidelobe_level = 20*log10( max(abs(As*w)) );
fprintf(1,'The minimum sidelobe level is %3.2f dB.\n\n',...
          min_sidelobe_level );

%********************************************************************
% plots
%********************************************************************
figure(1), clf
plot(loc(:,1),loc(:,2),'o')
title('Antenna locations')

% plot array pattern
if angleRange == 180,
    theta = (1:360)';
    A = [ A; -A ];
end
y = A*w;
figure(2), clf
ymin = floor(0.1*min_sidelobe_level)*10-10; ymax = 0;
plot(1:360, 20*log10(abs(y)), ...
     [theta_tar theta_tar],[ymin ymax],'r--',...
     [theta_tar+half_beamwidth theta_tar+half_beamwidth],[ymin ymax],'g--',...
     [theta_tar-half_beamwidth theta_tar-half_beamwidth],[ymin ymax],'g--');
if HAS_NULLS % add lines that represent null positions
  hold on;
  for k = 1:length(theta_nulls)
    plot([theta_nulls(k) theta_nulls(k)],[ymin ymax],'m--');
  end
  hold off;
end
xlabel('look angle'), ylabel('mag y(theta) in dB');
axis([0 360 ymin ymax]);

% polar plot
figure(3), clf
zerodB = -ymin;
dBY = 20*log10(abs(y)) + zerodB;
ind = find( dBY <= 0 ); dBY(ind) = 0;
plot(dBY.*cos(pi*theta/180), dBY.*sin(pi*theta/180), '-');
axis([-zerodB zerodB -zerodB zerodB]), axis('off'), axis('square')
hold on
plot(zerodB*cos(pi*theta/180),zerodB*sin(pi*theta/180),'k:') % 0 dB
plot( (min_sidelobe_level + zerodB)*cos(pi*theta/180), ...
      (min_sidelobe_level + zerodB)*sin(pi*theta/180),'k:')  % min level
text(-zerodB,0,'0 dB')
tt = text(-(min_sidelobe_level + zerodB),0,sprintf('%0.1f dB',min_sidelobe_level));
set(tt,'HorizontalAlignment','right');
theta_1 = theta_tar+half_beamwidth;
theta_2 = theta_tar-half_beamwidth;
plot([0 55*cos(theta_tar*pi/180)], [0 55*sin(theta_tar*pi/180)], 'k:')
plot([0 55*cos(theta_1*pi/180)], [0 55*sin(theta_1*pi/180)], 'k:')
plot([0 55*cos(theta_2*pi/180)], [0 55*sin(theta_2*pi/180)], 'k:')
if HAS_NULLS % add lines that represent null positions
  for k = 1:length(theta_nulls)
    plot([0 55*cos(theta_nulls(k)*pi/180)], ...
         [0 55*sin(theta_nulls(k)*pi/180)], 'k:')
  end
end
hold off

author: yu-jiang

repo: radpbook

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


text:
function r = r(kernel)

% R
%
% This method returns an upper bound on
%
%    R > k(x1, x1) - k(x1, x2)
%
% for arbitrary x1 and x2.  For a Gaussian radial basis kernel R = 1.

%
% File        : @rbf/r.m
%
% Date        : Friday 7th July 2000
%
% Author      : Dr Gavin C. Cawley
%
% Description : Part of an object-oriented implementation of Vapnik's Support
%               Vector Machine, as described in [1].
%
% References  : [1] V.N. Vapnik,
%                   "The Nature of Statistical Learning Theory",
%                   Springer-Verlag, New York, ISBN 0-387-94559-8,
%                   1995.
%
% History     : 07/07/2000 - v1.00
%               12/09/2000 - v1.01 minor improvements to comments and help
%                                  messages
%
% Copyright   : (c) Dr Gavin C. Cawley, July 2000
%
%    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
%

r = 1.0;

% bye bye...

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function varargout = std(varargin)
%STD   Standard deviation of a DISKFUN along one variable.
%   G = STD(F) returns the standard deviation of F in the radial variable
%   (default).If F is defined on the rectangle [-pi,pi] x [0,1] then
%
%                         1
%                        /
%     std(F)^2 = 1/pi    | ( F(THETA,R) - mean(F,1) )^2 dTHETA
%                        /
%                        0
%
%   G = STD(F, FLAG, DIM) takes the standard deviation along the
%   r-variable if DIM = 1 and along the theta-variable (angular) if
%   DIM = 2. The FLAG is ignored and kept in this function so the syntax
%   agrees with the Matlab STD command.
%
% See also CHEBFUN/STD, DISKFUN/MEAN.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

[varargout{1:nargout}] = std@separableApprox(varargin{:});

end


text:
function image = changeclass(class, varargin)
%CHANGECLASS Changes the storage class of an image.
%  I2 = CHANGECLASS(CLASS, I);
%  RGB2 = CHANGECLASS(CLASS, RGB);
%  BW2 = CHANGECLASS(CLASS, BW);
%  X2 = CHANGECLASS(CLASS, X, 'indexed');
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

%  Copyright 1993-2002 The MathWorks, Inc.  Used with permission.

switch class
case 'uint8'
   image = im2uint8(varargin{:});
case 'uint16'
   image = im2uint16(varargin{:});
case 'double'
   image = im2double(varargin{:});
otherwise
   error('Unsupported IPT data class.');
end

author: dipum

repo: dipum-toolbox

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


text:
function cortex = read_nimh_cortex(filename, varargin)

% READ_NIMH_CORTEX
%
% Use as
%  cortex = read_nimh_cortex(filename, ...)
%
% Optional input arguments should come in key-value pairs and may
% include
%   begtrial     = number (default = 1)
%   endtrial     = number (default = inf)
%   epp          = read the EPP data, 'yes' or 'no' (default = 'yes')
%   eog          = read the EOG data, 'yes' or 'no' (default = 'yes')
%   feedback     = display the progress on the screen, 'yes' or 'no' (default = 'no')
%
% The output is a structure array with one structure for every trial that was read.

% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: read_nimh_cortex.m 945 2010-04-21 17:41:20Z roboos $

% get the optional input arguments
feedback = keyval('feedback', varargin); if isempty(feedback), feedback = 'no'; end
begtrial = keyval('begtrial', varargin); if isempty(begtrial), begtrial = 1;   end
endtrial = keyval('endtrial', varargin); if isempty(endtrial), endtrial = inf; end
% reading the epp and eog data is optional
epp = keyval('epp', varargin); if isempty(epp), epp = 'yes'; end
eog = keyval('eog', varargin); if isempty(eog), eog = 'yes'; end

skipepp   = strcmp(epp, 'no');
skipeog   = strcmp(eog, 'no');
clear epp eog

% this will hold the result
cortex  = struct;

% trials are counted with zero-offset
trial = 0;

i_fid = fopen(filename, 'rb', 'ieee-le');

% read until the end of file or until the specified trial number
while ( ~feof (i_fid) && trial<endtrial )

  % the header of all trials must be read, otherwise the size of the subsequent trial data is not known
  % the header of uninteresting trials is not added to the output
  % the epp and eog data may be skipped completely
  skiptrial = (trial+1)<begtrial;

  length = fread(i_fid, 1, 'ushort');

  % allocate space for the header details
  hd = zeros(1,13);

  if ~isempty(length)

    % give some feedback on the screen
    if feedback && (trial+1)>=begtrial
      fprintf('reading trial %d\n', trial+1);
    end

    % Number of bytes for each variable.
    hd(1,1:8)= (fread(i_fid, 8, 'ushort'))';

    % Convert bytes to number of float point values (4 bytes apiece in windows and FreeBSD).
    hd(1,5)=hd(1,5)/4;
    % Convert bytes to number of short values (2 bytes apiece in windows and FreeBSD).
    hd(1,6)=hd(1,6)/2;
    hd(1,7)=hd(1,7)/2;
    hd(1,8)=hd(1,8)/2;

    header.cond_no =   hd(1,1);
    header.repeat_no = hd(1,2);
    header.block_no =  hd(1,3);
    header.trial_no =  hd(1,4);
    header.isi_size =  hd(1,5);
    header.code_size = hd(1,6);
    header.eog_size =  hd(1,7);
    header.epp_size =  hd(1,8);

    hd(1,9:10) = (fread(i_fid, 2, 'uchar'))';

    header.kHz_resolution   = hd(1,9);
    header.eye_storage_rate = hd(1,10);

    hd(1,11:13) = (fread(i_fid, 3, 'ushort'))';

    header.expected_response = hd(1,11);
    header.response          = hd(1,12);
    header.response_error    = hd(1,13);

    if skiptrial
      fseek(i_fid,header.isi_size,  'cof');
      fseek (i_fid,header.code_size, 'cof');
    else
      time  = (fread (i_fid,header.isi_size,  'ulong'));
      event = (fread (i_fid,header.code_size, 'ushort'));
    end % if skiptrial

    if skipepp || skiptrial
      epp1 = [];
      epp2 = [];
      fseek (i_fid,header.epp_size * 2, 'cof');

    else
      epp = fread (i_fid,header.epp_size, 'short');
      epp1 = zeros(1,header.epp_size/2);
      epp2 = zeros(1,header.epp_size/2);

      % fprintf(o_fid, '\nepp(x)\tepp(y)\n\n');

      for i = 1:2:header.epp_size
        % Must extract the data from the raw epp values.
        % The epp value is made up of 12-bits of data, and 4-bits (the
        % low-order 4 bits) of the channel number.
        % To extract the data, must right shift the raw data by 4 bits (to
        % get rid of the channel number and put the data value in the proper
        % location).  After this conversion, you must still add or subtract
        % 2048, since otherwise the value is not right.  (I think that this
        % is because of the way that matlab handles negative values during
        % the bitwise operations.)
        % These calculations cause the results of ctx2txt.m to be the same as
        % for cortview.exe for the EOG and EPP values.

        s = (i+1)/2;
        epp1(s) = bitshift(epp(i),   -4);
        epp2(s) = bitshift(epp(i+1), -4);

        if (epp1(s) < 0)
          epp1(s) = epp1(s) + 2047;
        else
          epp1(s) = epp1(s) - 2048;
        end;

        if (epp2(s) < 0)
          epp2(s) = epp2(s) + 2047;
        else
          epp2(s) = epp2(s) - 2048;
        end;
      end;
    end; % if skipepp

    if skipeog || skiptrial
      eog = [];
      fseek (i_fid,header.eog_size * 2, 'cof');
    else
      eog = fread (i_fid,header.eog_size, 'short');
      eog = reshape(eog, 2, header.eog_size/2);
    end % if skipeog

    if ~skiptrial
      % collect the headers of all trials in a structure array
      cortex(trial+1).header = header;
      cortex(trial+1).length = length;
      cortex(trial+1).time   = time;
      cortex(trial+1).event  = event;
      cortex(trial+1).epp    = [epp1; epp2];
      cortex(trial+1).eog    = eog;
    end

    trial = trial+1;
  end; % if ~isempty(length)
end;

fclose(i_fid);

% remove the initial empty trials, i.e. the ones that were skipped
cortex = cortex(begtrial:end);

author: PatternRecognition

repo: OpenBMI


text:
function x = inv_triu(U)
% INV_TRIU     Invert upper triangular matrix.

% Singularity test:
% inv_triu([1 1; 0 0])

x = solve_triu(U,eye(size(U)));
%x = inv(U);

author: andrewssobral

repo: lrslibrary


text:
function ldldemo
%LDLDEMO demo program for LDL
%
% Example:
%   ldldemo
%
% See also ldlsparse.

% Copyright 2006-2007 by Timothy A. Davis, Univ. of Florida

% compile the LDLSPARSE and LDLSYMBOL mexFunctions
help ldlsparse

fprintf ('\nTesting ldlsparse and ldlsymbol:\n') ;

% create a small random symmetric positive definite sparse matrix
n = 100 ;
d = 0.03 ;
rand ('state', 0) ;
randn ('state', 0) ;
A = sprandn (n, n, d) ;
A = speye (n) + A*A' ;
b = randn (n, 1) ;

figure (1)
clf
subplot (2,2,1) ;
spy (A) ;
title ('original matrix') ;

% permute for sparsity
p = symamd (A) ;
C = A (p,p) ;

subplot (2,2,2) ;
spy (C) ;
title ('permuted matrix') ;
drawnow

% factorize, without using ldlsparse's internal permutation
[L, D, Parent, fl] = ldlsparse (C) ;
L = L + speye (n) ;
err = norm (L*D*L' - C, 1) ;
fprintf ('norm (LDL''-PAP'') = %g\n', err) ;

% solve Ax=b
x = L' \ (D \ (L \ (b (p)))) ;
x (p) = x ;
resid = norm (A*x-b) ;
fprintf ('residual %g for ldlsparse, flops %10.1f\n', resid, fl) ;

% solve Ax=b with one call to ldlsparse
x = ldlsparse (C, [ ], b (p)) ;
x (p) = x ;
resid = norm (A*x-b) ;
fprintf ('residual %g for ldlsparse solve\n', resid) ;

subplot (2,2,3) ;
spy (L + D + L') ;
title ('L+D+L''') ;

subplot (2,2,4) ;
treeplot (Parent)
title ('elimination tree') ;

% try ldlrow (this will be slow)
[L, D] = ldlrow (C) ;
x = L' \ (D \ (L \ (b (p)))) ;
x (p) = x ;
resid = norm (A*x-b) ;
fprintf ('residual %g for ldlrow.m\n', resid) ;

% factorize, using ldlsparse's internal permutation
[L, D, Parent, fl] = ldlsparse (A, p) ;
L = L + speye (n) ;
err = norm (L*D*L' - C, 1) ;
fprintf ('norm (LDL''-PAP'') = %g\n', err) ;

% solve Ax=b
x = L' \ (D \ (L \ (b (p)))) ;
x (p) = x ;
resid = norm (A*x-b) ;
fprintf ('residual %g for ldlsparse, flops %10.1f\n', resid, fl) ;

% solve Ax=b with one call to ldlsparse
x = ldlsparse (A, p, b) ;
resid = norm (A*x-b) ;
fprintf ('residual %g for ldlsparse solve\n\n', resid) ;

% compare ldlsymbol and symbfact
[Lnz, Parent, fl] = ldlsymbol (A) ;
fprintf ('Original matrix: nz in L: %5d  flop count: %g\n', sum (Lnz), fl) ;

Lnz2 = symbfact (A) - 1 ;
Parent2 = etree (A) ;
fl2 = sum (Lnz2 .* (Lnz2 + 2)) ;
if (any (Lnz ~= Lnz2))
    error ('Lnz mismatch') ;
end
if (any (Parent ~= Parent2))
    error ('Parent mismatch') ;
end
if (fl ~= fl2)
    error ('fl mismatch') ;
end

[Lnz, Parent, fl] = ldlsymbol (A, p) ;
fprintf ('Permuted matrix: nz in L: %5d  flop count: %g\n', sum (Lnz), fl) ;

Lnz2 = symbfact (A (p,p)) - 1 ;
Parent2 = etree (A (p,p)) ;
fl2 = sum (Lnz2 .* (Lnz2 + 2)) ;
if (any (Lnz ~= Lnz2))
    error ('Lnz mismatch') ;
end
if (any (Parent ~= Parent2))
    error ('Parent mismatch') ;
end
if (fl ~= fl2)
    error ('fl mismatch') ;
end

fprintf ('\nldldemo: all tests passed\n') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function update_weights_adagrad()
    global config mem;
    for m = 1:length(config.weights)
        config.his_grad{m} = config.his_grad{m} + mem.grads{m} .* mem.grads{m};
        config.weights{m} = config.weights{m} - config.learning_rate * (mem.grads{m} ./ (config.fudge_factor + sqrt(config.his_grad{m})));
    end
end

author: jimmy-ren


text:
clear
clc

ObjectiveFunction = @my_first_SA;   % Function handle to the objective function
X0 = [1 1];   % Starting point
lb = [-2 -2];     % Lower bound
ub = [2 2];       % Upper bound

options = saoptimset('MaxIter',500,'StallIterLim',500,'TolFun',1e-100,'AnnealingFcn',@annealingfast,'InitialTemperature',100,'TemperatureFcn',@temperatureexp,'ReannealInterval',500,'PlotFcns',{@saplotbestx, @saplotbestf, @saplotx, @saplotf,@saplottemperature});

[x,fval] = simulannealbnd(ObjectiveFunction,X0,lb,ub,options);

author: vonsylvia


text:
function [nodes2, edges2] = grMergeNodeClusters(nodes, edges)
%GRMERGENODECLUSTERS Merge cluster of connected nodes in a graph.
%
%   grMergeNodeClusters(nodes, edges)
%   Detects groups of nodes that belongs to the same global node, and
%   replace them by a unique node. Coordinates of reference node is given
%   by the median coordinates of cluster nodes.
%
%   This function is intended to be used as filter after a binary image
%   skeletonization and vectorization.
%
%
%   See also
%   grMergeNodesMedian
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2003-08-13
% Copyright 2003-2022 INRA - TPV URPOI - BIA IMASTE

%% Initialization

% intialize result
nodes2 = nodes;
edges2 = edges;

% compute degree of each node
degrees = grNodeDegree(1:size(nodes, 1), edges)';

% find index of multiple nodes
indMul = find(degrees > 2);

% indices of edges that link several multiple nodes
indEdges = sum(ismember(edges, indMul), 2) == 2;

% associate a label to each cluster
labels = grLabel(nodes, edges(indEdges, :));
clusterLabels = unique(labels(indMul));

%% Replace each cluster by median point

% iterate on clusters
for i = 1:length(clusterLabels)
    % indices of nodes of the current cluster
    inds = find(labels == clusterLabels(i));

    % coordinates of new reference node
    clusterNodes = nodes(inds, :);
    medianNode = median(clusterNodes, 1);

    % replace coordinates of reference node
    refNode = min(inds);
    nodes2(refNode, :) = medianNode;

    % replace node indices in edge array
    edges2(ismember(edges2, inds)) = refNode;
end

%% Clean up

% keep only relevant nodes
inds = unique(edges2(:));
nodes2 = nodes2(inds, :);

% relabeling of edges
for i = 1:length(inds)
    edges2(edges2 == inds(i)) = i;
end

% remove double edges
edges2 = unique(sort(edges2, 2), 'rows');

% remove 'loops'
edges2(edges2(:,1) == edges2(:,2), :) = [];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
%% projection onto T_X(M_r)
%   X=U*diag(S)*V'
%   cf. [Van12]

function P=proj_fr(Z,U,V,sp)

if nargin<4, sp='t'; end

switch sp
    case 't'
        UtZ=U'*Z;
        P=U*UtZ+(Z*V)*V'-(U*(UtZ*V))*V';

    case 'n'
        UtZ=U'*Z;
        P=Z-(U*UtZ+(Z*V)*V'-(U*(UtZ*V))*V');
end

%% old version

% Pu=U*U';
% Pv=V*V';
%
% switch sp
%     case 't'
%         P=Pu*Z+Z*Pv-Pu*Z*Pv;
%
%     case 'n'
%         P=Z-Pu*Z-Z*Pv+Pu*Z*Pv;
% end


text:
%
% SCRIPT TO TEST THE KINEMATICS OF THE 5R robot
%
% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE 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.
%
% ARTE 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 ARTE.  If not, see <http://www.gnu.org/licenses/>.

close all

%load the robot
%the robot is represented (interiorly) by two planar 2dof arms
robot=load_robot('example','5R');

adjust_view(robot)

%Represent a final point in space that will be achieved by the end effector
%P=[Px, Py]=[0.6 0.4]
T=eye(4);
T(1,4)=1.0;%change, for example, to 1.25
%T(2,4)=0.967; %change, for example to 0.967
T(2,4)=0.8; %change, for example to 0.967

%find 4 solutions for the INVERSE KINEMATIC in position
q=inversekinematic(robot, T)

drawrobot3d(robot, q(:,1)), pause(1);
drawrobot3d(robot, q(:,2)), pause(1);
drawrobot3d(robot, q(:,3)), pause(1);
drawrobot3d(robot, q(:,4)), pause(1);
close all
%now solve the DIRECTKINEMATIC problem for each of the above computed
% q
for i=1:4,
    drawrobot3d(robot, q(:,i))
    figure, hold
    % 1 and 3 correspond to the joint variables of the robot
    T=directkinematic(robot, [q(1,i) q(3,i)])
    pause(2);
    close all
end

repo: ARTE


text:
function DIC_3Djoined_results = joinPairs(DIC3DAllPairsResults)
%% function for joining multiple 3D-DIC pair results into one structure (plotting in step3)
%
% Inputs:
% * DIC3DAllPairsResults
%
% Outputs:
% * DIC_3Djoined_results
%

%%
nPairs=numel(DIC3DAllPairsResults);
nFrames=size(DIC3DAllPairsResults{1}.Points3D,1);

% create cells with all the variables for joining
Fi=cell(nPairs,1);
FCi=cell(nPairs,1);
Fisoi=cell(nPairs,1);
Vi=cell(nPairs,nFrames);
% VTi=cell(nPairs,nFrames);
corrCombi=cell(nPairs,nFrames);
FaceCorrCombi=cell(nPairs,nFrames);
DispVi=cell(nPairs,nFrames);
% DispVTi=cell(nPairs,nFrames);
DispMi=cell(nPairs,nFrames);
% DispMTi=cell(nPairs,nFrames);
VFci=cell(nPairs,nFrames);
% VFcTi=cell(nPairs,nFrames);
Ji=cell(nPairs,nFrames);
Lamda1i=cell(nPairs,nFrames);
Lamda2i=cell(nPairs,nFrames);
Emgni=cell(nPairs,nFrames);
emgni=cell(nPairs,nFrames);
Epc1i=cell(nPairs,nFrames);
Epc2i=cell(nPairs,nFrames);
Epc1veci=cell(nPairs,nFrames);
Epc1vecCuri=cell(nPairs,nFrames);
Epc2veci=cell(nPairs,nFrames);
Epc2vecCuri=cell(nPairs,nFrames);
epc1i=cell(nPairs,nFrames);
epc2i=cell(nPairs,nFrames);
epc1veci=cell(nPairs,nFrames);
epc2veci=cell(nPairs,nFrames);

for ip=1:nPairs
    DIC_3Dpair_struct=DIC3DAllPairsResults{ip};

    Fi{ip}=DIC_3Dpair_struct.Faces;
    FCi{ip}=DIC_3Dpair_struct.FaceColors;

    for it=1:nFrames
        Vi{ip,it}=DIC_3Dpair_struct.Points3D{it};
%         VTi{ip,it}=DIC_3Dpair_struct.Points3Dtransformed{it};
        corrCombi{ip,it}=DIC_3Dpair_struct.corrComb{it};
        FaceCorrCombi{ip,it}=DIC_3Dpair_struct.FaceCorrComb{it};
        DispVi{ip,it}=DIC_3Dpair_struct.Disp.DispVec{it};
%         DispVTi{ip,it}=DIC_3Dpair_struct.Disp.DispVecTransformed{it};
        DispMi{ip,it}=DIC_3Dpair_struct.Disp.DispMgn{it};
%         DispMTi{ip,it}=DIC_3Dpair_struct.Disp.DispMgnTransformed{it};
        VFci{ip,it}=DIC_3Dpair_struct.FaceCentroids{it};
%         VFcTi{ip,it}=DIC_3Dpair_struct.FaceCentroidsTransformed{it};
        Ji{ip,it}=DIC_3Dpair_struct.Deform.J{it};
        Lamda1i{ip,it}=DIC_3Dpair_struct.Deform.Lamda1{it};
        Lamda2i{ip,it}=DIC_3Dpair_struct.Deform.Lamda2{it};
        Emgni{ip,it}=DIC_3Dpair_struct.Deform.Emgn{it};
        emgni{ip,it}=DIC_3Dpair_struct.Deform.emgn{it};
        Epc1i{ip,it}=DIC_3Dpair_struct.Deform.Epc1{it};
        Epc2i{ip,it}=DIC_3Dpair_struct.Deform.Epc2{it};
        Epc1veci{ip,it}=DIC_3Dpair_struct.Deform.Epc1vec{it};
        Epc1vecCuri{ip,it}=DIC_3Dpair_struct.Deform.Epc1vecCur{it};
        Epc2veci{ip,it}=DIC_3Dpair_struct.Deform.Epc2vec{it};
        Epc2vecCuri{ip,it}=DIC_3Dpair_struct.Deform.Epc2vecCur{it};
        epc1i{ip,it}=DIC_3Dpair_struct.Deform.epc1{it};
        epc2i{ip,it}=DIC_3Dpair_struct.Deform.epc2{it};
        epc1veci{ip,it}=DIC_3Dpair_struct.Deform.epc1vec{it};
        epc2veci{ip,it}=DIC_3Dpair_struct.Deform.epc2vec{it};
        Fisoi{ip,it}=DIC_3Dpair_struct.FaceIsoInd{it};
    end

end

% join variables that are time invariant
[Fj,~,~]=joinElementSets(Fi,Vi(:,1),[]);
[~,~,FCj]=joinElementSets(Fi,Vi(:,1),FCi);

for it=1:nFrames
    % join variables that are time variant and scalar

    [~,Vj{1,it},DispMj{1,it}]=joinElementSets(Fi,Vi(:,it),DispMi(:,it));
%     [~,VTj{1,it},DispMTj]=joinElementSets(Fi,VTi(:,it),DispMTi(:,it));
    [~,~,corrCombj{1,it}]=joinElementSets(Fi,Vi(:,it),corrCombi(:,it));
    [~,~,FaceCorrCombj{1,it}]=joinElementSets(Fi,Vi(:,it),FaceCorrCombi(:,it));
    [~,~,Jj{1,it}]=joinElementSets(Fi,Vi(:,it),Ji(:,it));
    [~,~,Lamda1j{1,it}]=joinElementSets(Fi,Vi(:,it),Lamda1i(:,it));
    [~,~,Lamda2j{1,it}]=joinElementSets(Fi,Vi(:,it),Lamda2i(:,it));
    [~,~,Emgnj{1,it}]=joinElementSets(Fi,Vi(:,it),Emgni(:,it));
    [~,~,emgnj{1,it}]=joinElementSets(Fi,Vi(:,it),emgni(:,it));
    [~,~,Epc1j{1,it}]=joinElementSets(Fi,Vi(:,it),Epc1i(:,it));
    [~,~,Epc2j{1,it}]=joinElementSets(Fi,Vi(:,it),Epc2i(:,it));
    [~,~,epc1j{1,it}]=joinElementSets(Fi,Vi(:,it),epc1i(:,it));
    [~,~,epc2j{1,it}]=joinElementSets(Fi,Vi(:,it),epc2i(:,it));
    [~,~,Fisoj{1,it}]=joinElementSets(Fi,Vi(:,it),Fisoi(:,it));

    % join variables that are time variant and vectors
    DispVj{1,it}=DispVi{1,it};
%     DispVTj{1,it}=DispVTi{1,it};
    VFcj{1,it}=VFci{1,it};
%     VFcTj{1,it}=VFcTi{1,it};
    Epc1vecj{1,it}=Epc1veci{1,it};
    Epc1vecCurj{1,it}=Epc1vecCuri{1,it};
    Epc2vecj{1,it}=Epc2veci{1,it};
    Epc2vecCurj{1,it}=Epc2vecCuri{1,it};
    epc1vecj{1,it}=epc1veci{1,it};
    epc2vecj{1,it}=epc2veci{1,it};
    for ip=2:nPairs
        DispVj{1,it}=[DispVi{1,it}; DispVi{ip,it}];
%         DispVTj{1,it}=[DispVTi{1,it}; DispVTi{ip,it}];
        VFcj{1,it}=[VFcj{1,it}; VFci{ip,it}];
%         VFcTj{1,it}=[VFcTj{1,it}; VFcTi{ip,it}];
        Epc1vecj{1,it}=[Epc1vecj{1,it}; Epc1veci{ip,it}];
        Epc1vecCurj{1,it}=[Epc1vecCurj{1,it}; Epc1vecCuri{ip,it}];
        Epc2vecj{1,it}=[Epc2vecj{1,it}; Epc2veci{ip,it}];
        Epc2vecCurj{1,it}=[Epc2vecCurj{1,it}; Epc2vecCuri{ip,it}];
        epc1vecj{1,it}=[epc1vecj{1,it}; epc1veci{ip,it}];
        epc2vecj{1,it}=[epc2vecj{1,it}; epc2veci{ip,it}];
    end
end

DIC_3Djoined_results=struct;
DIC_3Djoined_results.Faces=Fj;
DIC_3Djoined_results.FaceColors=FCj;
DIC_3Djoined_results.FaceIsoInd=Fisoj;
DIC_3Djoined_results.FaceCentroids=VFcj;
% DIC_3Djoined_results.FaceCentroidsTransformed=VFcTj;
DIC_3Djoined_results.Points3D=Vj;
% DIC_3Djoined_results.Points3Dtransformed=VTj;
DIC_3Djoined_results.corrComb=corrCombj;
DIC_3Djoined_results.FaceCorrComb=FaceCorrCombj;

DIC_3Djoined_results.Disp.DispVec=DispVj;
DIC_3Djoined_results.Disp.DispMgn=DispMj;
% DIC_3Djoined_results.Disp.DispVecTransformed=DispVTj;
% DIC_3Djoined_results.Disp.DispMgnTransformed=DispMTj;

DIC_3Djoined_results.Deform.J=Jj;
DIC_3Djoined_results.Deform.Lamda1=Lamda1j;
DIC_3Djoined_results.Deform.Lamda2=Lamda2j;
DIC_3Djoined_results.Deform.Emgn=Emgnj;
DIC_3Djoined_results.Deform.emgn=emgnj;
DIC_3Djoined_results.Deform.Epc1=Epc1j;
DIC_3Djoined_results.Deform.Epc2=Epc2j;
DIC_3Djoined_results.Deform.Epc1vec=Epc1vecj;
DIC_3Djoined_results.Deform.Epc1vecCur=Epc1vecCurj;
DIC_3Djoined_results.Deform.Epc2vec=Epc2vecj;
DIC_3Djoined_results.Deform.Epc2vecCur=Epc2vecCurj;
DIC_3Djoined_results.Deform.epc1=epc1j;
DIC_3Djoined_results.Deform.epc2=epc2j;
DIC_3Djoined_results.Deform.epc1vec=epc1vecj;
DIC_3Djoined_results.Deform.epc2vec=epc2vecj;

end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function name = saveName(n, l)
name = num2str(n);
for i=1:l-length(name)
    name = ['0' name];
end;

end

repo: Deep-ADMM-Net


text:
function out = getgenelist(input_dir)
% Reads gene-list text file and returns data.
% Default SIZE for GENELIST.txt is 7 columns! Each column is separated with
% a tab.
%
% USAGE:
%
%    out = getgenelist(input_dir)
%
% INPUT:
%    input_dir:    gene-list text file
%
% OUTPUT:
%    out:          data
%
% .. Author: - Stefan G. Thorleifsson 2011
%
% .. rBioNet is published under GNU GENERAL PUBLIC LICENSE 3.0+
% .. Thorleifsson, S. G., Thiele, I., rBioNet: A COBRA toolbox extension for
% .. reconstructing high-quality biochemical networks, Bioinformatics, Accepted.
% .. rbionet@systemsbiology.is
fid = fopen(input_dir);
line = fgetl(fid);
cnt = 0;
Column = 7; %default
data = cell(1,Column);
while line ~= -1
    cnt = cnt + 1;
    data_line=regexp(line, '\t','split');
    % If mistakes are made when creating genelist, then this removes empty
    % cells.
    if isempty(data_line{1}) || strcmp(line(1),'*') %Lines that start with a star are comments
        %add nothing
        cnt = cnt - 1;
    else
        Sd = size(data_line,2);
        if Sd < Column
            missing = Column-Sd;

            data_line(1,Sd+1:Sd+missing+1) = {'No data'};
        end
        data(cnt,1:Column) = data_line(1:Column);
    end

    line = fgetl(fid);
end

out = data;
fclose(fid);

author: opencobra

repo: cobratoolbox


text:
function stroud_test326 ( )

%*****************************************************************************80
%
%% TEST326 tests SPHERE_CAP_VOLUME_3D, SPHERE_CAP_VOLUME_ND.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;
  ntest = 12;
  r = 1.0;
  center = [ 0.0, 0.0, 0.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST326\n' );
  fprintf ( 1, '  SPHERE_CAP_VOLUME_3D computes the volume of a\n' );
  fprintf ( 1, '    spherical cap, defined by a plane that cuts the\n' );
  fprintf ( 1,'    sphere to a thickness of H units.\n' );
  fprintf ( 1, '  SPHERE_CAP_VOLUME_ND does the same operation,\n' );
  fprintf ( 1, '    but in N dimensions.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Using a radius R = %f\n', r );

  volume1 = sphere_volume_3d ( r );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Volume of the total sphere in 3D = %f\n', volume1 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '        H           Cap        Cap\n' );
  fprintf ( 1, '                    volume_3d  volume_nd\n' );
  fprintf ( 1, '\n' );

  for i = 0 : ntest + 1

    h = 2.0 * r * i / ntest;

    volume1 = sphere_cap_volume_3d ( r, h )

    volume2 = sphere_cap_volume_nd ( dim_num, r, h );

    fprintf ( 1, '  %12.6f  %12.6f  %12.6f\n', h, volume1, volume2 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%RRT Class for rapidly-exploring random tree navigation
%
% A concrete subclass of the abstract Navigation class that implements the rapidly
% exploring random tree (RRT) algorithm.  This is a kinodynamic planner
% that takes into account the motion constraints of the vehicle.
%
% Methods::
%  RRT          Constructor
%  plan         Compute the tree
%  query        Compute a path
%  plot         Display the tree
%  display      Display the parameters in human readable form
%  char         Convert to string
%
% Properties (read only)::
%  graph        A PGraph object describign the tree
%
% Example::
%        goal = [0,0,0];
%        start = [0,2,0];
%        veh = Bicycle('steermax', 1.2);
%        rrt = RRT(veh, 'goal', goal, 'range', 5);
%        rrt.plan()             % create navigation tree
%        rrt.query(start, goal)  % animate path from this start location
%
% References::
% - Randomized kinodynamic planning,
%   S. LaValle and J. Kuffner,
%   International Journal of Robotics Research vol. 20, pp. 378-400, May 2001.
% - Probabilistic roadmaps for path planning in high dimensional configuration spaces,
%   L. Kavraki, P. Svestka, J. Latombe, and M. Overmars,
%   IEEE Transactions on Robotics and Automation, vol. 12, pp. 566-580, Aug 1996.
% - Robotics, Vision & Control, Section 5.2.5,
%   P. Corke, Springer 2011.
%
% See also Navigation, PRM, DXform, Dstar, PGraph.

% 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

% Peter Corke 8/2009.

%TODO
%   more info to the display method
%   distance metric choice or weightings
%   pass time and model options to the simulation

classdef RRT < Navigation

    properties
        npoints         % number of points to find
        graph           % graph Object representing random nodes

        simtime         % path simulation time

        xrange          % range of x coordinates
        yrange          % range of y coordinates

        speed           % speed of vehicle
        vehicle         % Vehicle class object describes kinematics

        revcost         % penalty for going backwards

        root            % coordinate of the root of the tree (3x1)
    end

    methods

        function rrt = RRT(vehicle, varargin)
        %RRT.RRT Create an RRT navigation object
        %
        % R = RRT.RRT(VEH, OPTIONS) is a rapidly exploring tree navigation
        % object  for a vehicle kinematic model given by a Vehicle subclass object VEH.
        %
        % R = RRT.RRT(VEH, MAP, OPTIONS) as above but for a region with obstacles
        % defined by the occupancy grid MAP.
        %
        % Options::
        % 'npoints',N    Number of nodes in the tree (default 500)
        % 'simtime',T    Interval over which to simulate kinematic model toward
        %                random point (default 0.5s)
        % 'goal',P       Goal position (1x2) or pose (1x3) in workspace
        % 'speed',S      Speed of vehicle [m/s] (default 1)
        % 'root',R       Configuration of tree root (3x1) (default [0,0,0])
        % 'revcost',C    Cost penalty for going backwards (default 1)
        % 'range',R      Specify rectangular bounds of robot's workspace:
        %                - R scalar; X: -R to +R, Y: -R to +R
        %                - R (1x2); X: -R(1) to +R(1), Y: -R(2) to +R(2)
        %                - R (1x4); X: R(1) to R(2), Y: R(3) to R(4)
        %
        % Other options are provided by the Navigation superclass.
        %
        % Notes::
        % - 'range' option is ignored if an occupacy grid is provided.
        %
        % Reference::
        % - Robotics, Vision & Control
        %   Peter Corke, Springer 2011.  p102.
        %
        % See also Vehicle, Bicycle, Unicycle.

            % invoke the superclass constructor, it handles some options
            rrt = rrt@Navigation(varargin{:});

            rrt.vehicle = vehicle;

            % handle the options not done by Navigation superclass
            opt.npoints = 500;
            opt.simtime = 0.5;
            opt.speed = vehicle.speedmax;
            opt.revcost = 1;
            opt.root = [0 0 0];

            [rrt,args] = tb_optparse(opt, varargin, rrt);

            if isempty(rrt.occgrid)
                opt = [];
                opt.range = 5;
                [opt,args] = tb_optparse(opt, args);

                % range can be specified as scalar, min/max, different min/max per
                % direction
                switch length(opt.range)
                    case 1
                        rrt.xrange = [-opt.range opt.range];
                        rrt.yrange = [-opt.range opt.range];
                    case 2
                        rrt.xrange = [-opt.range(1) opt.range(1)];
                        rrt.yrange = [-opt.range(2) opt.range(2)];
                    case 4
                        rrt.xrange = [opt.range(1) opt.range(2)];
                        rrt.yrange = [opt.range(3) opt.range(4)];
                    otherwise
                        error('bad range specified');
                end
            else
                rrt.xrange = [1 numcols(rrt.occgrid)];
                rrt.yrange = [1 numrows(rrt.occgrid)];
            end

            rrt.graph = PGraph(3, 'distance', 'SE2', ...
                'dweight', 2*pi/norm(sum([rrt.xrange; rrt.yrange])) );  % graph of points in SE(2)

        end

        function plan(rrt, varargin)
        %RRT.plan Create a rapidly exploring tree
        %
        % R.plan(OPTIONS) creates the tree roadmap by driving the vehicle
        % model toward random goal points.  The resulting graph is kept
        % within the object.
        %
        % Options::
        % 'goal',P        Goal pose (1x3)
        % 'ntrials',N     Number of path trials (default 50)
        % 'noprogress'    Don't show the progress bar
        % 'samples'       Show progress in a plot of the workspace
        %                 - '.' for each random point x_rand
        %                 - 'o' for the nearest point which is added to the tree
        %                 - red line for the best path
        %
        % Notes::
        % - At each iteration we need to find a vehicle path/control that moves it
        %   from a random point towards a point on the graph.  We sample ntrials of
        %   random steer angles and velocities and choose the one that gets us
        %   closest (computationally slow, since each path has to be integrated
        %   over time).

            opt.progress = true;
            opt.samples = false;
            opt.goal = [];
            opt.ntrials = 50;

            opt = tb_optparse(opt, varargin);

            if ~isempty(opt.goal)
                rrt.goal = opt.goal;
            end

            % build a graph over the free space
            rrt.message('create the graph');
            rrt.graph.clear();

            if rrt.verbose
                clf
                %idisp(1-rrt.occgrid, 'ynormal', 'nogui');
                hold on
            end

            % check root node sanity
            if isempty(rrt.root)
                error('no root node specified');
            end
            if ~isvec(rrt.root, 3)
                error('root must be 3-vector');
            end
            assert( ~rrt.isoccupied(rrt.root(1:2)), 'root node cell is occupied')

            % add the goal point as the first node
            vroot = rrt.graph.add_node(rrt.root);
            data.vel = 0;
            data.path = [];
            rrt.graph.setvdata(vroot, data);

            % graphics setup
            if opt.progress
                h = Navigation.progress_init('RRT planning...');
            end
            if opt.samples
                clf
                hold on
                xlabel('x'); ylabel('y');
            end

            npoints = 0;
            while npoints < rrt.npoints       % build the tree

                % Step 3
                % find random state x,y

                % pick a point not in obstacle
                while true
                    xy = rrt.randxy();  % get random coordinate (x,y)

                    if isempty(rrt.occgrid)
                        break
                    else
                        % we have an occgrid
                        xy = round( xy );  % round it to a grid cell coordinate

                        % test if lies in the obstacle map
                        try
                            if ~rrt.isoccupied(xy)
                                break;
                            end
                        catch
                            % index error, point must be off the map
                            continue;
                        end
                    end
                end
                theta = rrt.rand*2*pi;
                xrand = [xy, theta]';
                if opt.samples
                    plot(xy(1), xy(2), '.')
                end

                % Step 4
                % find the existing node closest in state space

                vnear = rrt.graph.closest(xrand);   % nearest vertex
                xnear = rrt.graph.coord(vnear);     % coord of nearest vertex
%                 if rrt.graph.distance_metric(xnear, xrand) < 0.25
%                     continue;
%                 end

                rrt.message('xrand (%g, %g) node %d', xy, vnear);

                % Step 5
                % figure how to drive the robot from xnear to xrand

                best = rrt.bestpath(xnear, xrand, opt.ntrials);

                xnew = best.path(:,best.k);
                if opt.samples
                    plot(xnew(1), xnew(2), 'o');
                    plot2(best.path', 'r');
                    drawnow
                end

%                 % ensure that the path is collision free
%                 if ~rrt.clearpath(y(:,1:2))
%                     disp('path collision');
%                     continue;
%                 end

                % Step 7,8
                % add xnew to the graph, with an edge from xnear
                vnew = rrt.graph.add_node(xnew);

                if rrt.graph.vdata(vnear).vel * best.vel < 0
                    % we changed direction, penalise that
                    cost = rrt.revcost;
                else
                    cost = 1;
                end
                rrt.graph.add_edge(vnear, vnew, cost);

                rrt.graph.setvdata(vnew, best);

                npoints = npoints + 1;
                if opt.progress
                    Navigation.progress(h, npoints / rrt.npoints);
                end

            end

            if opt.progress
                Navigation.progress_delete(h)
            end
            rrt.message('graph create done');
        end

        function p_ = query(rrt, xstart, xgoal)
        %RRT.query Find a path between two points
        %
        % X = R.path(START, GOAL) finds a path (Nx3) from pose START (1x3)
        % to pose GOAL (1x3).  The pose is expressed as [X,Y,THETA].
        %
        % R.path(START, GOAL) as above but plots the path in 3D, where the vertical
        % axis is vehicle heading angle.  The nodes are shown as circles and the
        % line segments are blue for forward motion and red for backward motion.
        %
        % Notes::
        % - The path starts at the vertex closest to the START state, and ends
        %   at the vertex closest to the GOAL state.  If the tree is sparse this
        %   might be a poor approximation to the desired start and end.
        %
        % See also RRT.plot.

            assert(rrt.graph.n > 0, 'RTB:RRT: there is no plan');
            rrt.checkquery(xstart, xgoal);

            g = rrt.graph;
            vstart = g.closest(xstart);
            vgoal = g.closest(xgoal);

            % find path through the graph using A* search
            [path,cost] = g.Astar(vstart, vgoal);

            fprintf('A* path cost %g\n', cost);

            % concatenate the vehicle motion segments
            cpath = [];
            for i = 1:length(path)
                p = path(i);
                data = g.vdata(p);
                if ~isempty(data)
                    if i >= length(path) || g.edgedir(p, path(i+1)) > 0
                        cpath = [cpath data.path];
                    else
                        cpath = [cpath data.path(:,end:-1:1)];
                    end
                end
            end

            if nargout == 0
                % plot the path
                clf; hold on

                plot2(g.coord(path)', 'o');     % plot the node coordinates

                for i = 1:length(path)
                    p = path(i);
                    b = g.vdata(p);            % get path data for segment

                    % draw segment with direction dependent color
                    if ~isempty(b)
                        % if the vertex has a path leading to it

                        if i >= length(path) || g.edgedir(p, path(i+1)) > 0
                            % positive edge
                            %  draw from prev vertex to end of path
                            seg = [g.coord(path(i-1)) b.path]';
                        else
                            % negative edge
                            %  draw reverse path to next next vertex
                            seg = [  b.path(:,end:-1:1)  g.coord(path(i+1))]';
                        end

                        if b.vel > 0
                            plot2(seg, 'b');
                        else
                            plot2(seg, 'r');
                        end
                    end
                end

                xlabel('x'); ylabel('y'); zlabel('\theta');
                grid
            else
                p_ = cpath';
            end
        end

        function plot(rrt, varargin)
        %RRT.plot Visualize navigation environment
        %
        % R.plot() displays the navigation tree in 3D, where the vertical axis is
        % vehicle heading angle.  If an occupancy grid was provided this is also
        % displayed.

            % display the occgrid background
            rrt.plot_bg(varargin{:});

            % display the graph
            %rrt.graph.plot('noedges', 'NodeSize', 3, 'NodeFaceColor', 'm', 'NodeEdgeColor', 'm', 'edges');

            rrt.graph.plot('noedges', 'nocomponentcolor', 'NodeSize', 3, 'NodeFaceColor', 'b', 'NodeEdgeColor', 'b', 'edges');
hold on

            % display the occgrid background
            rrt.plot_fg(varargin{:});
            axis([rrt.xrange rrt.yrange])
            xlabel('x'); ylabel('y'); zlabel('\theta');
            grid on; hold off
            view(0,90);
            axis equal
            rotate3d
        end

        % required by abstract superclass
        function next(rrt)
        end

        function s = char(rrt)
        %RRT.char  Convert to string
        %
        % R.char() is a string representing the state of the RRT
        % object in human-readable form.
        %

            % invoke the superclass char() method
            s = char@Navigation(rrt);

            % add RRT specific stuff information
            s = char(s, sprintf('  region: X %f : %f; Y %f : %f', rrt.xrange, rrt.yrange));
            s = char(s, sprintf('  sim time: %f', rrt.simtime));
            s = char(s, sprintf('  speed: %f', rrt.speed));
            s = char(s, sprintf(' Graph:'));
            s = char(s, char(rrt.graph) );
            if ~isempty(rrt.vehicle)
                s = char(s, char(rrt.vehicle) );
            end
        end

    end % methods

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%    P R I V A T E    M E T H O D S
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    methods (Access='protected')

        function best = bestpath(rrt, x0, xg, N)

            % initial and final state as column vectors
            x0 = x0(:); xg = xg(:);

            best.d = Inf;
            for i=1:N   % for multiple trials

                %choose random direction of motion and random steer angle
                if rand > 0.5
                    vel = rrt.speed;
                else
                    vel = -rrt.speed;
                end
                steer = (2*rrt.rand - 1) * rrt.vehicle.steermax;    % uniformly distributed

                % simulate motion of vehicle for this speed and steer angle which
                % results in a path
                x = rrt.vehicle.run2(rrt.simtime, x0, vel, steer)';

                %% find point on the path closest to xg
                % distance of all path points from goal
                d = colnorm( [bsxfun(@minus, x(1:2,:), xg(1:2)); angdiff(x(3,:), xg(3))] );
                % the closest one
                [dmin,k] = min(d);

                % is it the best so far?
                if dmin < best.d
                    % yes it is!  save it and the inputs that led to it
                    best.d = dmin;
                    best.path = x;
                    best.steer = steer;
                    best.vel = vel;
                    best.k = k;
                end
            end
        end

        % generate a random coordinate within the working region
        function xy = randxy(rrt)
            xy = rrt.rand(1,2) .* [rrt.xrange(2)-rrt.xrange(1) rrt.yrange(2)-rrt.yrange(1)] + ...
                [rrt.xrange(1) rrt.yrange(1)];
        end

        % test if a path is obstacle free
        function c = clearpath(rrt, xy)
            if isempty(rrt.occgrid)
                c = true;
                return;
            end

            xy = round(xy);
            try
                % test that all points along the path do not lie within an obstacle
                for pp=xy'
                    if rrt.isoccupied(pp) > 0
                        c = false;
                        return;
                    end
                end
                c = true;
            catch
                % come here if we index out of bounds
                c = false;
                return;
            end
        end

    end % private methods
end % class

author: petercorke

repo: robotics-toolbox-matlab


text:
function fs_autosegmentToITK(subjID, t1, skipRecon, resample_type)
% Autosegment a t1 weighted anatomical scan using freesurfer
% autosegmentation, and then generate a nifti segmentation file in the
% space and at the resolution of the t1. This function requires that you
% have freesurfer in your shell path.
%
% fs_autosegmentToITK(subjID, t1, [skipRecon], [resample_type])
%
% INPUTS:
%   subjID: directory name in which freesurfer stores outputs
%   t1:     file name (with complete path) of t1 used for segmentation.
%
%   skipRecon: boolean. If true, then we regenerate a t1 class file (nifti)
%       `       from the already complete freesurfer segementation, without
%               re-doing the freesurfer segmentation. [default = false]
%   resample_type: resample type used by freesurfer for resmapling 1x1x1 mm
%               aseg file to the resolution of the t1. Options are
%               'nearest' and 'weighted'. Other options exist (like
%               'trilinear' but do not make sense as the resampling must
%               preserver integer values. [default = 'nearest']
%
% Example:
%   subjID  = 'jw';
%   t1      = fullfile('biac2/wandell2/data/anatomy/winawer/Anatomy20110308', 't1.nii.gz');
%   fs_autosegmentToITK(subjID, t1)
%
% 3/11/2009 Written by JW and HH.
% 6/30/2011 JW: added  'skipRecon' and 'resample_type' as optional input
%               arguments. Change the resample procedure to read the
%               deseried voxel size from the t1 file.
%
%
% see also fs_ribbon2itk.m

%% Check Inputs & Paths

% subjID is the name of the directory that will be created by freesurfer to
% store segmentation and associated files.
if ~exist('subjID', 'var')
    warning('Subject ID is required input'); %#ok<WNTAG>
    eval('help fs_autosegmentToITK');
    return
end

% Get the path to the t1 file if it is not inputed. This can be any
% resolution. Freesurfer will resample to 1x1x1 for autosegmentation. We
% will extract a segmentation at the resolution of the original t1 at the
% end of this function.
if notDefined('t1') || ~exist(t1, 'file'),
    [fname pth] = uigetfile({'t1*.nii.gz', 'T1 files'; '*.nii', '.nii files'; '*.gz', '.gz files'}, 'Cannot locate T1 file. Please find it yourself.', pwd);
    t1 = fullfile(pth, fname);
end
if ~exist(t1, 'file'), error('Cannot locate t1 file'); end

if ~exist('skipRecon', 'var') || isempty(skipRecon)
    skipRecon = false;
end
if ~exist('resample_type', 'var') || isempty(resample_type)
    resample_type = 'nearest';
end

% This is the directory where freesurfer puts subject data. If it is
% defined in the linux shell (e.g., bashrc) then matlab can find it. If it
% is not defined, look for the 'freesurfer_home/subjects', which is the
% default location in freesurfer.
subdir   = getenv('SUBJECTS_DIR');

if isempty(subdir),
    fshome = getenv('FREESURFER_HOME');
    subdir = fullfile(fshome, 'subjects');
end

%% recon all (freesurfer will resample to 1 mm isotropic)
%
if ~skipRecon,
    msg = sprintf('!recon-all -i %s -subjid ''%s'' -all', t1, subjID);
    eval(msg)
end

%
%% Resample aseg into nifti res (if nifit is not [1 1 1] mm)
% -rl flag reslices like the anatomical volume (t1)
% -rt uses nearest neighbour interpolation
%

% check whether t1 is [1 1 1] (freesurfer standard resolution). if not we
% need to resample the segmentation.

ni  = niftiRead(t1);
res = ni.pixdim;
if any(abs(res(1:3) - [1 1 1]) > .000001),
    resampleSeg = true;
else
    resampleSeg = false;
end

if resampleSeg
    msg = sprintf('cd(fullfile(''%s'', ''%s'', ''mri''))', subdir, subjID);
    eval(msg);

    % now make aseg file at resolution of the t1
    % back up 1-mm aseg.mgz and all ribbon files
    if ~exist('1mm', 'dir')
        mkdir('1mm');
        !mv *ribbon* aseg.mgz 1mm/
    end

    msg = sprintf('!mri_convert 1mm/aseg.mgz -o aseg.mgz -rt %s -voxsize %f %f %f', resample_type, res(1), res(2), res(3));
    eval(msg);

    %Recreate ribbon files at t1 resolution
    msg = sprintf('!mris_volmask --label_left_white 2 --label_left_ribbon 3 --label_right_white 41 --label_right_ribbon 42 --save_ribbon --save_distance ''%s'' ', subjID);
    eval(msg);

    % back up new segmentation files
    if ~exist('resampledSeg', 'dir'), mkdir('resampledSeg'); end
    !cp *ribbon* aseg.mgz resampledSeg/
end

%% Convert segmentation from mgz to nifti at resolution of nifti
outfile     = fullfile(fileparts(t1),...
    sprintf('t1_class_fs_%s.nii.gz',  datestr(now, 'yyyy-mm-dd-HH-MM-SS')));
fillWithCSF = true;
alignTo     = t1;
fs_ribbon2itk(subjID, outfile, fillWithCSF, alignTo);

%%
return


text:
%%%%%%%%%%%%%%%%%%%%%%%
%    R E S U L T S
%%%%%%%%%%%%%%%%%%%%%%%

%% Using T3 elements,

% Max y-displacement
% Node           UX           UY
%   8   1.916937e-06    -1.065042e-05
T3_disp = -1.065042e-05;

% Stress in element 5 [MPA]
T3_Sigma_xx_e5 = -650.431082;
T3_Sigma_yy_e5 = -650.431082;
T3_Sigma_xy_e5 = -650.431082;
T3_Vonmises_e5 =  1300.862163;

% Stress in element 6 [MPA]
T3_Sigma_xx_e6 =  650.431082;
T3_Sigma_yy_e6 =  645.504435;
T3_Sigma_xy_e6 = -1349.568918;
T3_Vonmises_e6 =  2425.672941;

% Average stress in Node 8 [MPA]

T3_Sigma_xx = (T3_Sigma_xx_e6+T3_Sigma_xx_e5)/2
T3_Sigma_yy = (T3_Sigma_yy_e6+T3_Sigma_yy_e5)/2
T3_Sigma_xy = (T3_Sigma_xy_e6+T3_Sigma_xy_e5)/2
T3_Vonmises = (T3_Vonmises_e6+T3_Vonmises_e5)/2

% T3_Sigma_xx = 0
% T3_Sigma_yy = -2.4633
% T3_Sigma_xy = -1000
% T3_Vonmises = 1.8633e+03

%% Using Q4 elements,

% Max y-displacement
% Node           UX           UY
%    8   6.1172e-06    -2.6430e-05
Q4_disp = -2.6430e-05;

% Stress in element 3

% Stress in node 3 [MPA]
Q4_Sigma_xx = 1947.025060;
Q4_Sigma_yy = -837.861465;
Q4_Sigma_xy = -543.285579;
Q4_Vonmises = 2647.590101;

%% Using Q8 elements,

% Max y-displacement
% Node           UX           UY
%   18   9.1557e-06    -3.8250e-05
Q8_disp = -3.8250e-05;

% Stress in element 3

% Stress in node 3 [MPA]
Q8_Sigma_xx = 547.919444;
Q8_Sigma_yy = -5070.115307;
Q8_Sigma_xy = -1926.948499;
Q8_Vonmises = 6318.519705;

%% Using Q12 elements,

% Max y-displacement
% Node           UX           UY
%   28   9.3240e-06    -3.8900e-05
Q12_disp = -3.8900e-05;

% Stress in element 3

% Stress in node 3 [MPA]
Q12_Sigma_xx = -209.968713;
Q12_Sigma_yy = -7862.634791;
Q12_Sigma_xy = -2337.689116;
Q12_Vonmises = 8752.632553;

%%%%%%%%%%%%%%%%%%%%%%%
%    C O M P A R E
%%%%%%%%%%%%%%%%%%%%%%%
x = 1:4;
str = ['T3','Q4','Q8','Q12'];
disp_y   = [T3_disp,Q4_disp,Q8_disp,Q12_disp];
Sigma_xx = [T3_Sigma_xx,Q4_Sigma_xx,Q8_Sigma_xx,Q12_Sigma_xx];
Sigma_yy = [T3_Sigma_yy,Q4_Sigma_yy,Q8_Sigma_yy,Q12_Sigma_yy];
Sigma_xy = [T3_Sigma_xy,Q4_Sigma_xy,Q8_Sigma_xy,Q12_Sigma_xy];
VonMises = [T3_Vonmises,Q4_Vonmises,Q8_Vonmises,Q12_Vonmises];
figure(1)
subplot(2,3,1); scatter(x,disp_y,'ob'); title('Max y-displacement'); ylabel('[mm]'); xlabel('T3 - Q4 - Q8 - Q12');
subplot(2,3,2); scatter(x,Sigma_xx,'ok'); title('\sigma_{xx}'); ylabel('[MPa]'); xlabel('T3 - Q4 - Q8 - Q12');
subplot(2,3,3); scatter(x,Sigma_yy,'dr'); title('\sigma_{yy}'); ylabel('[MPa]'); xlabel('T3 - Q4 - Q8 - Q12');
subplot(2,3,4); scatter(x,Sigma_xy,'sb'); title('\sigma_{xy}'); ylabel('[MPa]'); xlabel('T3 - Q4 - Q8 - Q12');
subplot(2,3,5);scatter(x,VonMises,'xk'); title('VonMises'); ylabel('[MPa]'); xlabel('T3 - Q4 - Q8 - Q12');

repo: Aero-matlab


text:
% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB 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.
%
% MVTB 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 MVTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
classdef scamera
    properties
        name
        Tcam
    end

    methods
        function c = scamera
            c.Tcam = eye(4,4);
        end

        function f = plot(c, points, varargin)
            f = c.project(points, varargin{:});
            plot(f(2,:), f(1,:), 'o');
            xlabel('phi (rad)');
            ylabel('theta (rad)');
            grid on
            axis([-pi pi 0 pi]);
        end

        function f = project(c, P, Tcam)
            if nargin < 3
                Tcam = c.Tcam;
            end
            P = transformp(inv(Tcam), P);

            R = sqrt( sum(P.^2) );
            x = P(1,:) ./ R;
            y = P(2,:) ./ R;
            z = P(3,:) ./ R;
            r = sqrt( x.^2 + y.^2);
            theta = atan2(r, z);
            phi = atan2(y, x);
            f = [theta; phi];
        end
    end
end

repo: machinevision-toolbox-matlab


text:
function handles = plotTMPCPrediction(ax, tPred, data, plot_tube, center_points, M_s, be_l, be_u, K_t, ...
    isConstrained, ub, lb, direction, isInput, TargetAvailable, target, LinearizationAvailable, lin)
% function plotStatePrediction
% Authors:      Martin Euler
%               Alexander Wischnewski
% Description:
%   helper function used to plot state predictions
% Inputs/parameters:
%   tPlot: vector containing the time window for plotting [tStart, tEnd]
%   tPred: vector with the time instants corresponding to the prediction
%   data: data of real beahvior
%   plot_tube: flag for tube plots
%   center_points: data for center points
%   M_s: tube shape matrices
%   be_l: lower bound of terminal set constraints
%   be_u: upper bound of terminal set constraints
%   K_t: tube controller
%   isConstrained: flag whether the signal is constrained or not
%   ub: upper limit for the signal
%   lb: lower limit for the signal
%   isInput: set to true in case this is an input (different calculation of tubes)
%   TargetAvailable: target trajectory available and should be plotted
%   target: actual target trajectory
%   LinearizationAvailable: linearization trajectory available and should be plotted
%   lin: actual linearization trajectory

% plot real signal
handles.real = plot(ax, data.Time, data.Data, 'Color', [0 0.3961 0.7412], ...
    'LineWidth', 1, 'DisplayName', 'Real');
% plot predicted signal
handles.pred = plot(ax, tPred, center_points, '*', 'Color', [0.8902, 0.4471, 0.1333], ...
    'DisplayName', 'Prediction');
% plot tube and terminal sets
if plot_tube
    % get tube bounds
    [y_low, y_up] = calcBoundsPlot(center_points, M_s, isInput, K_t, length(tPred), direction);
    % plot tube lower bounds
    handles.tube_ub = plot(ax, tPred, y_low(1,:), '-.', ...
        'Color', [0, 0, 0], 'LineWidth', 1, 'DisplayName', 'Tube');
    % plot tube upper bounds
    handles.tube_lb = plot(ax, tPred, y_up(1,:), '-.', ...
        'Color', [0, 0, 0], 'LineWidth', 1, 'HandleVisibility','off');
    % plot terminal set lower bounds and use a virtual time interval of one second (half a second
    % forth and back) to visualize it properly
    handles.terminal_ub = plot(ax, [tPred(end)-0.5, tPred(end)+0.5], be_l*ones(2, 1), ...
        'Color', [0, 0, 1],'LineWidth', 1.1, 'DisplayName', 'Terminal set');
    % plot terminal set upper bounds
    handles.terminal_lb = plot(ax, [tPred(end)-0.5, tPred(end)+0.5], be_u*ones(2, 1), ...
        'Color', [0, 0, 1], 'LineWidth', 1.1,'HandleVisibility','off');
else
end
if isConstrained
    handles.const_lb = plot(ax, tPred, ub, 'm--', ...
        'Color', [0.6275 0.1255 0.9412], 'LineWidth', 1, 'HandleVisibility', 'off');
    handles.const_ub = plot(ax, tPred, lb, 'm--', ...
        'Color', [0.6275 0.1255 0.9412], 'LineWidth', 1, 'DisplayName', 'Constraint');
end
if TargetAvailable
    handles.target = plot(ax, tPred, target, '-', 'DisplayName', 'Target');
end
if LinearizationAvailable
    handles.lin = plot(ax, tPred, lin, 'c', 'DisplayName', 'Linearization');
end

% adjust x limits appropriately such that it covers a certain range
xlim(ax, [tPred(1)-5, tPred(1)+10]);
legend(ax);

end

author: TUMFTM


text:
%%*****************************************************************************
%% checkdepconst: compute AAt to determine if the
%%             constraint matrices Ak are linearly independent.
%%
%% [At,b,y,idxB,neardepconstr,feasible,AAt] = checkdepconstr(blk,At,b,y,rmdepconstr);
%%
%% rmdepconstr = 1, if want to remove dependent columns in At
%%             = 0, otherwise.
%%
%% idxB = indices of linearly independent columns of original At.
%% neardepconstr = 1 if there is nearly dependent columns in At
%%            = 0, otherwise.
%% Note: the definition of "nearly dependent" is dependent on the
%%       threshold used to determine the small diagonal elements in
%%       the LDLt factorization of A*At.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************************

function  [At,b,y,idxB,neardepconstr,feasible,AAt] = ...
    checkdepconstr(blk,At,b,y,rmdepconstr)

global existlowrank printlevel
global nnzmatold
%%
%% compute AAt
%%
m = length(b);
AAt = sparse(m,m);  numdencol = 0; UU = [];
for p = 1:size(blk,1)
    pblk = blk(p,:);
    if strcmp(pblk{1},'s');
        m1 = size(At{p,1},2); m2 = m - m1;
        if (m2 > 0)
            if (m2 > 0.3*m); AAt = full(AAt); end
            dd = At{p,3};
            lenn = size(dd,1);
            idxD = [0; find(diff(dd(:,1))); lenn];
            ss = [0, cumsum(pblk{3})];
            if (existlowrank)
                AAt(1:m1,1:m1) = AAt(1:m1,1:m1) + At{p,1}'*At{p,1}; %#ok
                for k = 1:m2
                    idx = ss(k)+1 : ss(k+1);
                    idx2 = idxD(k)+1: idxD(k+1);
                    ii = dd(idx2,2)-ss(k); %% undo cumulative indexing
                    jj = dd(idx2,3)-ss(k);
                    len = pblk{3}(k);
                    Dk = spconvert([ii,jj,dd(idx2,4); len,len,0]);
                    tmp = svec(pblk,At{p,2}(:,idx)*Dk*At{p,2}(:,idx)');
                    tmp2 = AAt(1:m1,m1+k) + (tmp'*At{p,1})';
                    AAt(1:m1,m1+k) = tmp2; %#ok
                    AAt(m1+k,1:m1) = tmp2'; %#ok
                end
            end
            DD = spconvert([dd(:,2:4); sum(pblk{3}),sum(pblk{3}),0]);
            VtVD = (At{p,2}'*At{p,2})*DD;
            VtVD2 = VtVD'.*VtVD;
            for k = 1:m2
                idx0 = ss(k)+1 : ss(k+1);
                %%tmp = VtVD(idx0,:)' .* VtVD(:,idx0);
                tmp = VtVD2(:,idx0);
                tmp = tmp*ones(length(idx0),1);
                tmp3 = AAt(m1+1:m1+m2,m1+k) + mexqops(pblk{3},tmp,ones(length(tmp),1),1);
                AAt(m1+1:m1+m2,m1+k) = tmp3; %#ok
            end
        else
            AAt = AAt + abs(At{p,1})'*abs(At{p,1});
        end
    else
        decolidx = checkdense(At{p,1}');
        if ~isempty(decolidx);
            n2 = size(At{p,1},1);
            dd = ones(n2,1);
            len= length(decolidx);
            dd(decolidx) = zeros(len,1);
            AAt = AAt + (abs(At{p,1})' *spdiags(dd,0,n2,n2)) *abs(At{p,1});
            tmp = At{p,1}(decolidx,:)';
            UU = [UU, tmp]; %#ok
            numdencol = numdencol + len;
        else
            AAt = AAt + abs(At{p,1})'*abs(At{p,1});
        end
    end
end
if (numdencol > 0) && (printlevel)
    fprintf('\n number of dense column in A = %d',numdencol);
end
numdencol = size(UU,2);
%%
%%
%%
feasible = 1; neardepconstr = 0;
if ~issparse(AAt); AAt = sparse(AAt); end
nnzmatold = mexnnz(AAt);
rho = 1e-15;
diagAAt = diag(AAt);
mexschurfun(AAt,rho*max(diagAAt,1));
[L.R,indef,L.perm] = chol(AAt,'vector');
L.d = full(diag(L.R)).^2;
if (indef)
    msg = 'AAt is not pos. def.';
    idxB = 1:m;
    neardepconstr = 1;
    if (printlevel); fprintf('\n checkdepconstr: %s',msg); end
    return;
end
%%
%% find independent rows of A
%%
dd = zeros(m,1);
idxB = (1:m)';
dd(L.perm) = abs(L.d);
idxN = find(dd < 1e-13*mean(L.d));
ddB = dd(setdiff(1:m,idxN));
ddN = dd(idxN);
if ~isempty(ddN) && ~isempty(ddB) && (min(ddB)/max(ddN) < 10)
    %% no clear separation of elements in dd
    %% do not label constraints as dependent
    idxN = [];
end
if ~isempty(idxN)
    neardepconstr = 1;
    if (printlevel)
        fprintf('\n number of nearly dependent constraints = %1.0d',length(idxN));
    end
    if (numdencol==0)
        if (rmdepconstr)
            idxB = setdiff((1:m)',idxN);
            if (printlevel)
                fprintf('\n checkdepconstr: removing dependent constraints...');
            end
            [W,resnorm] = findcoeffsub(blk,At,idxB,idxN);
            tol = 1e-8;
            if (resnorm > sqrt(tol))
                idxB = (1:m)';
                neardepconstr = 0;
                if (printlevel)
                    fprintf('\n checkdepconstr: basis rows cannot be reliably identified,');
                    fprintf(' abort removing nearly dependent constraints');
                end
                return;
            end
            tmp = W'*b(idxB) - b(idxN);
            nnorm = norm(tmp)/max(1,norm(b));
            if (nnorm > tol)
                feasible = 0;
                if (printlevel)
                    fprintf('\n checkdepconstr: inconsistent constraints exist,');
                    fprintf(' problem is infeasible.');
                end
            else
                feasible = 1;
                for p = 1:size(blk,1)
                    At{p,1} = At{p,1}(:,idxB);
                end
                b = b(idxB);
                y = y(idxB);
                AAt = AAt(idxB,idxB);
            end
        else
            if (printlevel)
                fprintf('\n To remove these constraints,');
                fprintf(' re-run sqlp.m with OPTIONS.rmdepconstr = 1.');
            end
        end
    else
        if (printlevel)
            fprintf('\n warning: the sparse part of AAt may be nearly singular.');
        end
    end
end
%%*****************************************************************************
%% findcoeffsub:
%%
%% [W,resnorm] = findcoeffsub(blk,At,idXB,idXN);
%%
%% idXB = indices of independent columns of At.
%% idxN = indices of   dependent columns of At.
%%
%% AB = At(:,idxB); AN = At(:,idxN) = AB*W
%%
%% SDPT3: version 3.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 2 Feb 01
%%*****************************************************************************

function [W,resnorm] = findcoeffsub(blk,At,idxB,idxN)

AB = []; AN = [];
for p = 1:size(blk,1)
    AB = [AB; At{p,1}(:,idxB)]; %#ok
    AN = [AN; At{p,1}(:,idxN)]; %#ok
end
n = size(AB,2);
%%
%%-----------------------------------------
%% find W so that AN = AB*W
%%-----------------------------------------
%%
[L,U,P,Q] = lu(sparse(AB));
rhs  = P*AN;
Lhat = L(1:n,:);
W = Q*( U \ (Lhat \ rhs(1:n,:)));
resnorm = norm(AN-AB*W,'fro')/max(1,norm(AN,'fro'));
%%*****************************************************************************

author: zarathustr

repo: LibQPEP


text:
function test_suite = test_vol_grid_convert
% tests for cosmo_vol_grid_convert
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_vol_grid_convert_basics
    aet=@(varargin)assertExceptionThrown(@()...
                            cosmo_vol_grid_convert(varargin{:}),'');

    ds_orig=cosmo_synthetic_dataset('size','normal');
    ds_orig.a.vol=rmfield(ds_orig.a.vol,'xform');
    nfeatures=size(ds_orig.samples,2);
    nfeatures2=round(nfeatures/2);
    while true
        rp=randperm(nfeatures);
        ds=cosmo_slice(ds_orig,repmat(rp(1:nfeatures2),1,2),2);
        ijk=[ds.fa.i; ds.fa.j; ds.fa.k];

        no_missing_features=true;
        for dim=1:3
            no_missing_features=no_missing_features && ...
                                        isequal(unique(ijk(dim,:)),...
                                                    1:ds.a.vol.dim(dim));
        end
        if no_missing_features
            break;
        end
    end

    ds2=cosmo_vol_grid_convert(ds);
    assert(all(~cosmo_isfield(ds2.fa,{'i','j','k'})));
    assertEqual(ds2.a.fdim.labels,{'pos'});
    assertEqual(numel(ds2.a.fdim.values),1);
    assertFalse(isfield(ds2.a,'vol'));

    assert(isfield(ds2.fa,'pos'));
    pos=ds2.a.fdim.values{1}(:,ds2.fa.pos);
    assertEqual(pos, cosmo_vol_coordinates(ds));

    ds3=cosmo_vol_grid_convert(ds,'togrid');
    assertEqual(ds2,ds3);
    ds4=cosmo_vol_grid_convert(ds3,'tovol');
    assertEqual(ds4,ds);
    ds4same=cosmo_vol_grid_convert(ds4,'tovol');
    assertEqual(ds4,ds4same);
    aet(ds4,'foo');

    ds4=cosmo_vol_grid_convert(ds2);
    assertEqual(ds4,ds);

    ds4=cosmo_vol_grid_convert(ds2);
    assertEqual(ds4,ds);

    ds5=cosmo_vol_grid_convert(ds2,'tovol');
    assertEqual(ds5,ds4);

    ds.fa=rmfield(ds.fa,'j');
    aet(ds);

    % irregular grid
    ds3.a.fdim.values{1}(1)=.5;
    aet(ds3);

    % should not work for datasets with missing pos and vol
    ds6=cosmo_synthetic_dataset('type','timelock');
    aet(ds6,'tovol');
    aet(ds6,'togrid');
    aet(ds6);
    aet(ds6,'tovol','togrid');

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function y = horzcat(varargin)
    % Horizontal concatenation for sym arrays.
    %   C = HORZCAT(A, B, ...) horizontally concatenates the sym arrays A,
    %   B, ... .  For matrices, all inputs must have the same number of rows.  For
    %   N-D arrays, all inputs must have the same sizes except in the second
    %   dimension.
    %
    %   C = HORZCAT(A,B) is called for the syntax [A B].
    %
    %   See also VERTCAT.

    X = cellfun(@(x)SymExpression(x),varargin,'UniformOutput',false);
        %     X = SymExpression(varargin);
    str = cellfun(@(x)x.s,X,'UniformOutput',false);
    % evaluate the operation in Mathematica and return the
    % expression string
    sstr = ['Join[Sequence@@ToMatrixForm/@' cell2tensor(str,'ConvertString',false) ',2]'];

    % create a new object with the evaluated string
    y = SymExpression(sstr);

end

author: ayonga

repo: frost-dev

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