text:
function [PBarUVals,dPBarUValsdTheta,d2PBarUValsdTheta2]=NALegendreCosRat(theta,M,scalFactor)
%NALEGENDRECOSRAT Evaluate \bar{P}_{nm}(cos(theta))/u^m for all n from
%                 0 to M and for each n for all m from 0 to n, where
%                 u=abs(sin(theta)). \bar{P}_{nm}(x) is the fully
%                 normalized associated Legendre function of x of degree n
%                 and order m. Also evaluate
%                 D{\bar{P}_{nm}(cos(theta))}/u^m and
%                 D2{\bar{P}_{nm}(cos(theta))}/u^m where D{} is the first
%                 derivative operator with respect to theta and D2{} is the
%                 second derivative operator with respect to theta. All of
%                 the values can be scaled by a factor of scalFac, if
%                 desired, to help prevent overflows with high degrees and
%                 orders.
%
%INPUTS: theta An angle in radians.
%       maxDeg The maximum degree and order of the output. This must be
%              >=0.
%   scalFactor A scale factor to help prevent overflow of the results. In
%              [1], discussed below, a value of 10^(-280) is used.
%
%OUTPUTS: PBarUVals An instance of the CountingClusterSet class such that
%                   PBarUVals(n+1,m+1)=scalFac*\bar{P}_{nm}(cos(theta))/u^m.
%  dPBarUValsdTheta An instance of the CountingClusterSet class such that
%                   dPBarUValsdTheta(n+1,m+1)=scalFac*D{\bar{P}_{nm}(cos(theta))}/u^m
% d2PBarUValsdTheta2 An instance of the CountingClusterSet class such that
%                   d2PBarUValsdTheta2(n+1,m+1)=scalFac*D2{\bar{P}_{nm}(cos(theta))}/u^m
%
%The modified forward row (MFR) algorithm of [1] is used to compute
%PBarUVals and dPBarUValsdTheta. For d2PBarUValsdTheta2, the algorithm of
%[2] is used. However, the paper contains a typo. The first term of the
%first unnumbered equation should be multiplied by an additional (1/u).
%This function uses the correct formulation, which is also described in
%Appendix F of [3].
%
%With the notation P^m_n(x) for an associated Legendre function of degree n
%and order m evaluated at the point x, one generally means
%P^m_n(x)=(-1)^m*(1-x^)^(m/2)*D_m{P_n(x)}
%where P_n(x) is a Legendre polynomial of degree n evaluated at x and
%D_m{P_n(x)} is the mth derivative of the polynomial. These associated
%Legendre functions are the same as those used in the International
%Geomagnetic Reference Field. On the other hand, with the notation
%P_{nm}(x), one generally means (-1)^mP^m_n(x). This notation is also
%called an associated Legendre function. Matlab's built-in function
%legendre(n,x) returns all P^m_n(x) for m=0 to m=n. On the other hand, the
%notation \bar{P}_{nm}(x), refers to a fully normalized associated Legendre
%function. Multiple definitions of what is normalized exist. In this
%function, the normalization is such that
%\bar{P}_{nm}(x)=P_{nm}(x)*sqrt(k*(2*n+1)*factorial(n-m)/factorial(n+m))
%where k=1 if m=0 and k=2 otherwise. This differs from the normalization
%constant that Matlab uses in the normalized version of its built-in
%associated Legendre function. It is the same as the normalization
%constant the NGA uses in its coefficients for the EGM96 and the
%EGM2008 gravitation models.
%
%The fully normalized associated Legendre function ratios that this
%function computes can be used in the synthesis of spherical harmonic
%coefficients.
%
%REFERENCES:
%[1] S. A. Holmes and W. E. Featherstone, "A unified approach to the
%    Clenshaw summation and the recursive computation of very high degree
%    and order normalised associated Legendre functions," Journal of
%    Geodesy, vol. 76, no. 5, pp. 279-299, May 2002.
%[2] S. A. Holmes and W. E. Featherstone, "Short note: Extending simplified
%    high-degree synthesis methods to second latitude derivatives of
%    geopotential," Journal of Geodesy, vol. 76, no. 8, pp. 447-450, Nov.
%    2002.
%[3] D. F. Crouse, "An overview of major terrestrial, celestial, and
%    temporal coordinate systems for target tracking," Naval Research
%    Laboratory, Washington, DC, Tech. Rep. NRL/FR/5344-16-10,279, 10 Aug.
%    2016.
%
%December 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    if(M<0)
       error('M must be >=0.')
    end

    u=sin(theta);
    t=cos(theta);

    %Allocate space for the return variables. n ranges from 0 to M and m
    %ranges from 0 to n. Thus, the number of elements for each n ranges
    %from 1 to M+1. This means that a total of
    %(M+1)*(M+2)/2 elements are necessary.
    numPBarU=(M+1)*(M+2)/2;
    totalP=zeros(numPBarU,1);
    PBarUVals=CountingClusterSet(totalP);

    %The value of PBar_{0,0}(cos(theta)) is independent of theta and is
    %one.
    PBarUVals(0+1,0+1)=1*scalFactor;

    if(M>0)
        %Set the seed value for PBar_{1,1}(cos(theta))/u from which the
        %other values will be computed.
        PBarUVals(1+1,1+1)=sqrt(3)*scalFactor;

        jTerm=1/sqrt(2);
        %First, deal with the case where n=1, m=0;
        n=1;
        m=0;
        %g is given in Equation 18 of the first Holmes and Featherstone paper.
        g=2*(m+1)/sqrt((n-m)*(n+m+1));
        PBarUVals(n+1,m+1)=jTerm*g*t*PBarUVals(n+1,m+1+1);

        %Compute the values along the main diagonal, where m=n
        %starting from m=n=2. This implements equation 28 in the first
        %Holmes and Featherstone paper for the normalized associated
        %Legendre function ratio.
        for m=2:M
            PBarUVals(m+1,m+1)=sqrt((2*m+1)/(2*m))*PBarUVals(m-1+1,m-1+1);
        end

        %Recursively compute the values using Equation 27 from the
        %first Holmes and Featherstone paper, taking into account the
        %fact that the first element of the recursion only has one term.

        %Next, evaluate the values for all other valid n and m.
        for n=2:M
            %Deal with the first element of the recursion,which is
            %where m=n-1.
            m=n-1;
            g=2*(m+1)/sqrt((n-m)*(n+m+1));
            PBarUVals(n+1,m+1)=g*t*PBarUVals(n+1,m+1+1);

            %Recursively compute the values of the rest of the m terms.
            for m=(n-2):-1:1
                g=2*(m+1)/sqrt((n-m)*(n+m+1));
                %h is given in Equation 19 of the first Holmes and
                %Featherstone paper.
                h=sqrt((n+m+2)*(n-m-1)/((n-m)*(n+m+1)));
                PBarUVals(n+1,m+1)=g*t*PBarUVals(n+1,m+1+1)-h*u^2*PBarUVals(n+1,m+2+1);
            end

            %Deal with the special m=0 case.
            m=0;
            g=2*(m+1)/sqrt((n-m)*(n+m+1));
            h=sqrt((n+m+2)*(n-m-1)/((n-m)*(n+m+1)));
            PBarUVals(n+1,m+1)=jTerm*(g*t*PBarUVals(n+1,m+1+1)-h*u^2*PBarUVals(n+1,m+2+1));
        end
    end

    %If the first derivative is desired.
    if(nargout>1)
        %Allocate space.
        dPBarUValsdTheta=CountingClusterSet(totalP);
        %The first derivative of PBar_{0,0}(cos(theta)) is just zero.
        dPBarUValsdTheta(1,1)=0;

        if(M>0)
            m=1;
            n=1;
            %From Equation 30 in the first Holmes and Featherstone paper. This
            %is the seed value from which other values will be computed.
            dPBarUValsdTheta(1+1,1+1)=m*(t/u)*PBarUVals(1+1,1+1);

            n=1;
            m=0;
            %e is given in Equation 22 of the first Holmes and Featherstone
            %paper.
            e=sqrt((n+m+1)*(n-m)/2);
            %This is Equation 30 of the first Holmes and Featherstone paper for
            %m=0.
            dPBarUValsdTheta(n+1,m+1)=-e*u*PBarUVals(n+1,m+1+1);

            %Compute the values along the main diagonal, where m=n starting
            %from m=n=2. This implements Equation 30 in the frist Holmes and
            %Featherstone paper for the ratio of the first derivative.
            for m=2:M
                dPBarUValsdTheta(m+1,m+1)=m*(t/u)*PBarUVals(m+1,m+1);
            end

            %Next, evaluate the values for all other valid n and m.
            for n=2:M
                %Recursively compute the values of the m terms for m>0.
                for m=(n-1):-1:1
                    e=sqrt((n+m+1)*(n-m));
                    dPBarUValsdTheta(n+1,m+1)=m*(t/u)*PBarUVals(n+1,m+1)-e*u*PBarUVals(n+1,m+1+1);
                end
                %Deal with the special m=0 case.
                m=0;
                e=sqrt((n+m+1)*(n-m)/2);
                dPBarUValsdTheta(n+1,m+1)=-e*u*PBarUVals(n+1,m+1+1);
            end
        end
    end

    %If the second derivative is desired
    if(nargout>2)
        %Allocate space
        d2PBarUValsdTheta2=CountingClusterSet(totalP);

        for n=0:M
            for m=0:n
                %From the first (un-numbered) equation in the second Holmes
                %and Featherstone paper AFTER correction.
                d2PBarUValsdTheta2(n+1,m+1)=(m^2/u^2-n*(n+1))*PBarUVals(n+1,m+1)-(t/u)*dPBarUValsdTheta(n+1,m+1);
            end
        end
    end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
classdef MimShowOutputFolder < MimGuiPlugin
    % MimShowOutputFolder. Gui Plugin for opening an explorer/finder window in the
    % cache folder of the currently visualised dataset
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimShowOutputFolder is a Gui Plugin for the MIM Toolkit.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Show output files'
        SelectedText = 'Show output files'
        ToolTip = 'Opens a folder containing the output files for the current dataset'
        Category = 'Results'
        Visibility = 'Dataset'
        Mode = 'Analysis'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        Location = 41
        Icon = 'output_files.png'
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            dataset_cache_path = gui_app.GetOutputPath;
            CoreDiskUtilities.OpenDirectoryWindow(dataset_cache_path);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
classdef SpatialNorm < dagnn.ElementWise
  properties
    param = [2 2 10 2]
  end

  methods
    function outputs = forward(obj, inputs, params)
      outputs{1} = vl_nnspnorm(inputs{1}, obj.param) ;
    end

    function [derInputs, derParams] = backward(obj, inputs, param, derOutputs)
      derInputs{1} = vl_nnspnorm(inputs{1}, obj.param, derOutputs{1}) ;
      derParams = {} ;
    end

    function obj = SpatialNorm(varargin)
      obj.load(varargin) ;
    end
  end
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function varargout = drawEllipse3d(varargin)
%DRAWELLIPSE3D Draw a 3D ellipse
%
%   Possible calls for the function :
%   drawEllipse3d([XC YC ZC A B THETA PHI])
%   drawEllipse3d([XC YC ZC A B THETA PHI PSI])
%   drawEllipse3d([XC YC ZC A B], [THETA PHI])
%   drawEllipse3d([XC YC ZC A B], [THETA PHI PSI])
%   drawEllipse3d([XC YC ZC A B], THETA, PHI)
%   drawEllipse3d([XC YC ZC], A, B, THETA, PHI)
%   drawEllipse3d([XC YC ZC A B], THETA, PHI, PSI)
%   drawEllipse3d([XC YC ZC], A, B, THETA, PHI, PSI)
%   drawEllipse3d(XC, YC, ZC, A, B, THETA, PHI)
%   drawEllipse3d(XC, YC, ZC, A, B, THETA, PHI, PSI)
%
%   where XC, YC, ZY are coordinate of ellipse center, A and B are the
%   half-lengths of the major and minor axes of the ellipse,
%   PHI and THETA are 3D angle (in degrees) of the normal to the plane
%   containing the ellipse (PHI between 0 and 360 corresponding to
%   longitude, and THETA from 0 to 180, corresponding to angle with
%   vertical).
%
%   H = drawEllipse3d(...)
%   return handle on the created LINE object
%
%   ------
%   Author: David Legland
%   e-mail: david.legland@grignon.inra.fr
%   Created: 2008-05-07
%   Copyright 2008 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

%   HISTORY

%   Possible calls for the function, with number of arguments :
%   drawEllipse3d([XC YC ZC A B THETA PHI])             1
%   drawEllipse3d([XC YC ZC A B THETA PHI PSI])         1
%   drawEllipse3d([XC YC ZC A B], [THETA PHI])          2
%   drawEllipse3d([XC YC ZC A B], [THETA PHI PSI])      2
%   drawEllipse3d([XC YC ZC A B], THETA, PHI)           3
%   drawEllipse3d([XC YC ZC A B], THETA, PHI, PSI)      4
%   drawEllipse3d([XC YC ZC], A, B, THETA, PHI)         5
%   drawEllipse3d([XC YC ZC], A, B, THETA, PHI, PSI)    6
%   drawEllipse3d(XC, YC, ZC, A, B, THETA, PHI)         7
%   drawEllipse3d(XC, YC, ZC, A, B, THETA, PHI, PSI)    8

% extract drawing options
ind = find(cellfun(@ischar, varargin), 1, 'first');
options = {};
if ~isempty(ind)
    options = varargin(ind:end);
    varargin(ind:end) = [];
end

if length(varargin)==1
    % get center and radius
    ellipse = varargin{1};
    xc = ellipse(:,1);
    yc = ellipse(:,2);
    zc = ellipse(:,3);
    a  = ellipse(:,4);
    b  = ellipse(:,5);

    % get colatitude of normal
    if size(ellipse, 2)>=6
        theta = ellipse(:,6);
    else
        theta = zeros(size(ellipse, 1), 1);
    end

    % get azimut of normal
    if size(ellipse, 2)>=7
        phi     = ellipse(:,7);
    else
        phi = zeros(size(ellipse, 1), 1);
    end

    % get roll
    if size(ellipse, 2)==8
        psi = ellipse(:,8);
    else
        psi = zeros(size(ellipse, 1), 1);
    end

elseif length(varargin)==2
    % get center and radius
    ellipse = varargin{1};
    xc = ellipse(:,1);
    yc = ellipse(:,2);
    zc = ellipse(:,3);
    a  = ellipse(:,4);
    b  = ellipse(:,5);

    % get angle of normal
    angle = varargin{2};
    theta   = angle(:,1);
    phi     = angle(:,2);

    % get roll
    if size(angle, 2)==3
        psi = angle(:,3);
    else
        psi = zeros(size(angle, 1), 1);
    end

elseif length(varargin)==3
    % get center and radius
    ellipse = varargin{1};
    xc = ellipse(:,1);
    yc = ellipse(:,2);
    zc = ellipse(:,3);
    a  = ellipse(:,4);
    b  = ellipse(:,5);

    % get angle of normal and roll
    theta   = varargin{2};
    phi     = varargin{3};
    psi     = zeros(size(phi, 1), 1);

elseif length(varargin)==4
    % get center and radius
    ellipse = varargin{1};
    xc = ellipse(:,1);
    yc = ellipse(:,2);
    zc = ellipse(:,3);

    if size(ellipse, 2)==5
        a  = ellipse(:,4);
        b  = ellipse(:,5);
    end

    theta   = varargin{2};
    phi     = varargin{3};
    psi     = varargin{4};

elseif length(varargin)==5
    % get center and radius
    ellipse = varargin{1};
    xc      = ellipse(:,1);
    yc      = ellipse(:,2);
    zc      = ellipse(:,3);
    a       = varargin{2};
    b       = varargin{3};
    theta   = varargin{4};
    phi     = varargin{5};
    psi     = zeros(size(phi, 1), 1);

elseif length(varargin)==6
    ellipse = varargin{1};
    xc      = ellipse(:,1);
    yc      = ellipse(:,2);
    zc      = ellipse(:,3);
    a       = varargin{2};
    b       = varargin{3};
    theta   = varargin{4};
    phi     = varargin{5};
    psi     = varargin{6};

elseif length(varargin)==7
    xc      = varargin{1};
    yc      = varargin{2};
    zc      = varargin{3};
    a       = varargin{4};
    b       = varargin{5};
    theta   = varargin{6};
    phi     = varargin{7};
    psi     = zeros(size(phi, 1), 1);

elseif length(varargin)==8
    xc      = varargin{1};
    yc      = varargin{2};
    zc      = varargin{3};
    a       = varargin{4};
    b       = varargin{5};
    theta   = varargin{6};
    phi     = varargin{7};
    psi     = varargin{8};

else
    error('drawEllipse3d: please specify center and radius');
end

% uses 60 intervals
t = linspace(0, 2*pi, 61)';

% polyline approximation of ellipse, centered and parallel to main axes
x       = a * cos(t);
y       = b * sin(t);
z       = zeros(length(t), 1);
base    = [x y z];

% compute transformation from local basis to world basis
trans   = localToGlobal3d(xc, yc, zc, theta, phi, psi);

% transform points composing the ellipse
ellipse = transformPoint3d(base, trans);

% draw the curve
h = drawPolyline3d(ellipse, options{:});

if nargout > 0
    varargout = {h};
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [J, grad] = linearRegCostFunction(X, y, theta, lambda)
%LINEARREGCOSTFUNCTION Compute cost and gradient for regularized linear
%regression with multiple variables
%   [J, grad] = LINEARREGCOSTFUNCTION(X, y, theta, lambda) computes the
%   cost of using theta as the parameter for linear regression to fit the
%   data points in X and y. Returns the cost in J and the gradient in grad

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost and gradient of regularized linear
%               regression for a particular choice of theta.
%
%               You should set J to the cost and grad to the gradient.
%

h = X*theta;
hy = h-y;
J = sum(hy.^2)/(2*m) + lambda*sum(theta(2:end).^2)/(2*m);
hyrep = repmat(hy,[1 size(X,2)]);
grad = sum(hyrep.*X,1)'/m;

grad(2:end) = grad(2:end) + lambda.*theta(2:end)/m;

% =========================================================================

grad = grad(:);

end

repo: machine-learning-coursera


text:
function f = exp(f)
%EXP  Exponential of a CHEBFUN3T object.
%   EXP(F) returns the exponential of a CHEBFUN3T object F.

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

% Empty check:
if ( isempty(f) )
    return
end

op = @(x,y,z) exp(feval(f, x, y, z));    % Resample.
f = chebfun3t(op, f.domain);             % Call constructor.

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function output = callnone(interfacedata)

% Standard interface
output.Primal      = zeros(length(interfacedata.c),1);
output.Dual        = [];
output.Slack       = [];
output.problem     = -2;
output.infostr     = yalmiperror(-2,'YALMIP');
output.solverinput = [];
output.solveroutput= [];
output.solvertime  = 0;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


function [AIC, HQIC, BIC] = IC(S, E, T, K)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 'IC' computes various information criteria

% Inputs:
% - llf, log marginal likelihood
% - T, sample size
% - K, number of regressors

% Filippo Ferroni, 6/1/2015
% Revised, 2/15/2017
% Revised, 3/21/2018
% Revised, 9/11/2019
% Revised, 9/11/2020
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% E   = var.e_ols;
% S   = var.Sigma_ols;
iS  = pinv(S);
N   = size(S,1);
llf = - (T * N / 2) * (1 + log(2 * pi)) - T / 2 * log(det(S));
llf = llf - 1 /2 * trace( iS * E' * E);

AIC = - 2 * llf / T + 2 * K / T;
% SIC = - 2 * llf / T + K * log(T) / T;
HQIC = - 2 * llf / T + 2 * K * log(log(T)) / T;
BIC = - 2 * llf / T + K * log(T) / T;


text:
function y = diag( v, k )

%Disciplined convex/geometric programming information for DIAG:
%   DIAG imposes no convexity restrictions on its arguments.

switch nargin,
    case 0,
        error( 'Not enough arguments.' );
    case 1,
        k = 0;
    case 2,
        if ~isnumeric( k ) || k ~= floor( k ),
            error( 'Second argument must be an integer.' );
        end
end

s = size( v );
if length( s ) ~= 2,
    error( 'First input must be 2D.' );
end

if k < 0,
    absk = -k;
    roff = absk;
    coff = 0;
else
    absk = +k;
    roff = 0;
    coff = absk;
end

if any( s == 1 ),
    nn = prod( s );
    nel = nn + roff + coff;
    y = sparse( roff + 1 : roff + nn, coff + 1 : coff + nn, v, nel, nel );
elseif roff >= s(1) || coff >= s(2),
    y = sparse( 0, 1 );
else
    nel = min( s(1) - roff, s(2) - coff );
    nv = roff + ( coff - 1 ) * s(1) + ( 1 : nel ) * ( s(1) + 1 );
    y = reshape( cvx_subsref( v, nv ), nel, 1 );
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function cs_demo2 (do_pause, matrixpath)
%CS_DEMO2 MATLAB version of the CSparse/Demo/cs_demo2.c program.
%   Solves a linear system using Cholesky, LU, and QR, with various orderings.
%
% Example:
%   cs_demo2
% See also: cs_demo

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

if (nargin < 2)
    matrixpath = [] ;
end

if (isempty (matrixpath))
    try
        % older versions of MATLAB do not have an input argument to mfilename
        p = mfilename ('fullpath') ;
        t = strfind (p, filesep) ;
        matrixpath = [ p(1:t(end)) '../../Matrix' ] ;
    catch
        % assume we are in the C*Sparse/MATLAB/CSparse/Demo directory
        matrixpath = '../../Matrix' ;
    end
end

matrices = { 't1', 'HB/fs_183_1', 'HB/west0067', 'LPnetlib/lp_afiro', ...
'HB/ash219', 'HB/mbeacxc', 'HB/bcsstk01', 'HB/bcsstk16' } ;

if (nargin < 1)
    do_pause = 1 ;
end

for i = 1:length(matrices)
    name = matrices {i} ;
    [C sym] = get_problem (matrixpath, name, 1e-14) ;
    demo2 (C, sym, name) ;
    if (do_pause)
        input ('Hit enter to continue: ') ;
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ signal_filter2 ] = medianfilter_is( signal, fs )
% [signal_filter2] = medianfilter_is( signal, fs )
%   OVERVIEW:   This function estimates the baseline wander signal in the ECG and returns the estimate in the signal_filter2 variable
%
%	INPUT: 	MANDATORY:
%               signal          : a single row of ECG data in samples.
%
%               fs              : sampling frequency for the ecg signal (Hz)
%
%
%   	OUTPUT:
%            	signal_filter2     : baseline wander estimate for the ecg in var signal.
%
%
%	REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Written by Ismail Sadiq
%	COPYRIGHT (C) 2019
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information. The license may be found in the Documents
%       folder of the Physionet-Cardiovascular-Signal-Toolbox.

orderfilter1 = floor(0.2 * fs);
orderfilter2 = floor(0.6 * fs);

signal_filter1 = medfilt1(signal, orderfilter1);
signal_filter2 = medfilt1(signal_filter1, orderfilter2);

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function tutor = svctutor(arg)

% SVCTUTOR
%
% Constructor for a class of tutor objects used to train support vector
% classification networks.  Note this is an abstract base class, you cannot
% train an svc object using an svctutor.

%
% File        : @svctutor/svctutor.m
%
% Date        : Tuesday 12th September 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, September 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
%

if nargin == 0

   % this is the default constructor

   tutor.dummy = [];
   tutor       = class(tutor, 'svctutor');

elseif isa(arg, 'svctutor');

   % this is the copy constructor

   tutor = arg;

else

   % there are no other constructors

   help svctutor

end

% bye bye...

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef MimGuiDatasetState < CoreBaseClass
    % MimGuiDatasetState. Stores information about the currently loaded series
    %
    %     You do not need to modify this file. To add new functionality, create
    %     new plguins in the Plugins and GuiPlugins folders.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (SetAccess = private)
        CurrentSeriesUid
        CurrentPatientId
        CurrentPatientVisibleName
        CurrentSeriesName
        CurrentModality

        CurrentPluginInfo
        CurrentPluginName
        CurrentVisiblePluginName
        CurrentPluginResultIsEdited

        CurrentSegmentationName
    end

    events
        SeriesUidChangedEvent
        PatientIdChangedEvent
        PluginChangedEvent

        ManualSegmentationsChanged % Manual segmentations changed for current dataset
        PreviewImageChanged % Preview images changed for current dataset
        MarkersChanged % Markers changed for current dataset
    end

    methods
        function SetPatientAndSeries(obj, patient_id, series_uid, patient_visible_name, series_name, modality)

            if ~strcmp(patient_id, obj.CurrentPatientId) || ~strcmp(series_uid, obj.CurrentSeriesUid)
                obj.CurrentPatientId = patient_id;
                obj.CurrentSeriesUid = series_uid;
                obj.CurrentPatientVisibleName = patient_visible_name;
                obj.CurrentSeriesName = series_name;
                obj.CurrentModality = modality;
                notify(obj, 'SeriesUidChangedEvent', CoreEventData(series_uid));
            end
        end

        function SetPatientClearSeries(obj, patient_id, patient_visible_name)
            patient_changed = ~strcmp(patient_id, obj.CurrentPatientId);
            series_changed = ~isempty(obj.CurrentSeriesUid);
            obj.CurrentPatientId = patient_id;
            obj.CurrentSeriesUid = [];
            obj.CurrentPatientVisibleName = patient_visible_name;
            obj.CurrentSeriesName = [];
            obj.CurrentModality = [];
            if series_changed
                notify(obj, 'SeriesUidChangedEvent', CoreEventData([]));
            elseif patient_changed
                notify(obj, 'PatientIdChangedEvent', CoreEventData(patient_id));
            end
        end

        function ClearSeries(obj)
            if ~isempty(obj.CurrentPatientId) || ~isempty(obj.CurrentSeriesUid)
                obj.CurrentSeriesUid = [];
                obj.CurrentSeriesName = [];
                obj.CurrentModality = [];
                obj.CurrentSegmentationName = [];
                notify(obj, 'SeriesUidChangedEvent', CoreEventData([]));
            end
        end

        function ClearPatientAndSeries(obj)
            if ~isempty(obj.CurrentPatientId) || ~isempty(obj.CurrentSeriesUid)
                obj.CurrentPatientId = [];
                obj.CurrentSeriesUid = [];
                obj.CurrentPatientVisibleName = [];
                obj.CurrentSeriesName = [];
                obj.CurrentModality = [];
                obj.CurrentSegmentationName = [];
                notify(obj, 'SeriesUidChangedEvent', CoreEventData([]));
            end
        end

        function ClearPlugin(obj)
            obj.CurrentPluginInfo = [];
            obj.CurrentPluginName = [];
            obj.CurrentVisiblePluginName = [];
            obj.CurrentPluginResultIsEdited = false;
            obj.CurrentSegmentationName = [];
            notify(obj, 'PluginChangedEvent', CoreEventData([]));
        end

        function SetPlugin(obj, plugin_info, plugin_name, plugin_visible_name, is_edited)
            obj.CurrentPluginInfo = plugin_info;
            obj.CurrentPluginName = plugin_name;
            obj.CurrentVisiblePluginName = plugin_visible_name;
            obj.CurrentPluginResultIsEdited = is_edited;
            obj.CurrentSegmentationName = [];
            notify(obj, 'PluginChangedEvent', CoreEventData(plugin_name));
        end

        function SetSegmentation(obj, segmentation_name)
            obj.ClearPlugin;
            obj.CurrentSegmentationName = segmentation_name;
            notify(obj, 'PluginChangedEvent', CoreEventData(segmentation_name));
        end

        function UpdateEditStatus(obj, is_edited)
            if ~isequal(obj.CurrentPluginResultIsEdited, is_edited)
                obj.CurrentPluginResultIsEdited = is_edited;
                notify(obj, 'PluginChangedEvent', CoreEventData(obj.CurrentPluginName));
            end
        end
    end
end


text:
function grad = B_logistic(input_layers, CostLayer)

[~, output, target] = prepareCostEvaluation(input_layers, CostLayer);
m = length(output);
output = sigmoid(output);

grad = -(target - output)/m;

end

author: singaxiong

repo: SignalGraph


text:
function W = compute_mesh_weight(vertex,face,type,options)

% compute_mesh_weight - compute a weight matrix
%
%   W = compute_mesh_weight(vertex,face,type,options);
%
%   W is sparse weight matrix and W(i,j)=0 is vertex i and vertex j are not
%   connected in the mesh.
%
%   type is either
%       'combinatorial': W(i,j)=1 is vertex i is conntected to vertex j.
%       'distance': W(i,j) = 1/d_ij^2 where d_ij is distance between vertex
%           i and j.
%       'conformal': W(i,j) = cot(alpha_ij)+cot(beta_ij) where alpha_ij and
%           beta_ij are the adjacent angle to edge (i,j)
%
%   If options.normalize=1, the the rows of W are normalize to sum to 1.
%
%   Copyright (c) 2007 Gabriel Peyre

options.null = 0;
[vertex,face] = check_face_vertex(vertex,face);

nface = size(face,1);
n = max(max(face));

if isfield(options, 'verb')
    verb = options.verb;
else
    verb = n>5000;
end

if nargin<3
    type = 'conformal';
end

switch lower(type)
    case 'combinatorial'
        W = triangulation2adjacency(face);
    case 'distance'
        W = my_euclidean_distance(triangulation2adjacency(face),vertex);
        W(W>0) = 1./W(W>0);
        W = (W+W')/2;
    case 'conformal'
        % conformal laplacian
        W = sparse(n,n);
        ring = compute_vertex_face_ring(face);
        for i = 1:n
            if verb
                progressbar(i,n);
            end
            for b = ring{i}
                % b is a face adjacent to a
                bf = face(:,b);
                % compute complementary vertices
                if bf(1)==i
                    v = bf(2:3);
                elseif bf(2)==i
                    v = bf([1 3]);
                elseif bf(3)==i
                    v = bf(1:2);
                else
                    error('Problem in face ring.');
                end
                j = v(1); k = v(2);
                vi = vertex(:,i);
                vj = vertex(:,j);
                vk = vertex(:,k);
                % angles
                alpha = myangle(vk-vi,vk-vj);
                beta = myangle(vj-vi,vj-vk);
                % add weight
                W(i,j) = W(i,j) + cot( alpha );
                W(i,k) = W(i,k) + cot( beta );
            end
        end
    otherwise
        error('Unknown type.')
end

if isfield(options, 'normalize') && options.normalize==1
    W = diag(sum(W,2).^(-1)) * W;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function beta = myangle(u,v);

du = sqrt( sum(u.^2) );
dv = sqrt( sum(v.^2) );
du = max(du,eps); dv = max(dv,eps);
beta = acos( sum(u.*v) / (du*dv) );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function W = my_euclidean_distance(A,vertex)

if size(vertex,1)<size(vertex,2)
    vertex = vertex';
end

[i,j,s] = find(sparse(A));
d = sum( (vertex(i,:) - vertex(j,:)).^2, 2);
W = sparse(i,j,d);

author: gpeyre

repo: matlab-toolboxes


text:
function value = r4_asin ( x )

%*****************************************************************************80
%
%% R4_ASIN evaluates the arc-sine of an R4 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the arc-sine of X.
%
  persistent asincs
  persistent nterms
  persistent sqeps

  if ( isempty ( nterms ) )
    asincs = [ ...
      0.10246391753227159, ...
      0.054946487221245833, ...
      0.004080630392544969, ...
      0.000407890068546044, ...
      0.000046985367432203, ...
      0.000005880975813970, ...
      0.000000777323124627, ...
      0.000000106774233400, ...
      0.000000015092399536, ...
      0.000000002180972408, ...
      0.000000000320759842, ...
      0.000000000047855369, ...
      0.000000000007225128, ...
      0.000000000001101833, ...
      0.000000000000169476, ...
      0.000000000000026261, ...
      0.000000000000004095, ...
      0.000000000000000642, ...
      0.000000000000000101, ...
      0.000000000000000016 ]';
    nterms = r4_inits ( asincs, 20, 0.1 * r4_mach ( 3 ) );
    sqeps = sqrt ( 6.0 * r4_mach ( 3 ) );
  end

  y = abs ( x );

  if ( x < - 1.0 - sqeps )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_ASIN - Fatal error!\n' );
    fprintf ( 1, '  X < - 1.0\n' );
    error ( 'R4_ASIN - Fatal error!' )

  elseif ( x < - 1.0 )

    value = - 0.5 * pi;

  elseif ( x < 1.0 )

    z = 0.0;
    if ( sqeps < y )
      z = y * y;
    end

    if ( z <= 0.5 )
      value = x * ( 1.0 + r4_csevl ( 4.0 * z - 1.0, asincs, nterms ) );
    else
      value = 0.5 * pi - sqrt ( 1.0 - z ) * ( 1.0 + ...
        r4_csevl ( 3.0 - 4.0 * z, asincs, nterms ) );
    end

    if ( x < 0.0 )
      value = - abs ( value );
    elseif ( 0.0 < x )
      value = + abs ( value );
    end

  elseif ( x < 1.0 + sqeps )

    value = 0.5 * pi;

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_ASIN - Fatal error!\n' );
    fprintf ( 1, '  1.0 < X\n' );
    error ( 'R4_ASIN - Fatal error!' )

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%filename='Cantilever_tetrahedra';
%filename='Cantileverbeam_Tetrahedra_Linear_Structured_Fine';
%filename = 'Cantilever3D';
%filename = 'Cantilever3DLarge';
filename = 'Cantilever3DLargeFineFine';
%filename='Cantilever_quad_fine';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance','perimeter'};
%cost = {'compliance'};
weights = [1, 1];
%weights = 1;
constraint = {'volumeConstraint'};
optimizerUnconstrained = 'SLERP';
incrementFactor = 1;
designVariable = 'LevelSet';

%optimizerUnconstrained = 'PROJECTED GRADIENT';
%designVariable = 'Density';

filterType = 'PDE';
optimizer = 'DualNestedInPrimal';
%optimizer = 'AlternatingPrimalDual';

nsteps = 80;
Vfrac_final = 0.15;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final = 1e-2;

Vfrac_initial      = 0.15;
optimality_initial = 0.5*1e-2;
constr_initial     = 1e-2;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

plotting = true;
printing = false;
monitoring = true;
monitoring_interval = 1;

maxiter = 80;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


% Convert the samples in the input file to a different format (ex: float complex to int16 complex)
%
% The input and output types need to be string names like 'single', 'int16', or 'int8'.
% To convert a file containing complex 32-bit floats to complex 16-bit ints for a radio with a 12-bit DAC:
%     convert_complex_file_type('/tmp/foo.fc32', 'single', '/tmp/foo.sc16', 'int16', 2^12);
%
% @param input_file_path Path to the input file containing complex samples of `input_type`
% @param input_type Sample type of each I and Q value in `input_file_path`
% @param output_file_path Path to the new file containing samples from the input file converted to `output_type`
% @param output_type New sample type of each I and Q value (converted from `input_type`)
% @param scalar Value to scale the input samples by before writing to disk
% @return sample_count Number of complex samples written to `output_file_path`
function [sample_count] = convert_complex_file_type(input_file_path, input_type, output_file_path, output_type, scalar)
    assert(isstring(input_file_path) || ischar(input_file_path) || iscellstr(input_file_path), ...
        'Input file path must be a string, char array, or cell string');
    assert(isstring(input_type) || ischar(input_type), 'Input type must be a string or char array');
    assert(isstring(output_file_path) || ischar(output_file_path) || iscellstr(output_file_path), ...
        'Output file path must be a string, char array, or cell string');
    assert(isstring(output_file_path) || ischar(output_file_path), 'Output type must be a string or char array');
    assert(isnumeric(scalar), 'Scalar must be a number');

    input_file_handle = fopen(input_file_path, 'r');
    assert(input_file_handle ~= -1, 'Could not open input file "%s"', input_file_path);

    output_file_handle = fopen(output_file_path, 'w');
    assert(output_file_handle ~= -1, 'Could not open output file "%s" for writing', output_file_path);

    input_file_sample_bytes = get_bytes_per_sample(input_type);

    % How many samples to work on at a time.  Keep this value as high as possible to reduce runtime.  But, don't go
    % crazy with it as it will plow through RAM
    chunk_size = 1e6;

    % Track the number of samples written
    sample_count = 0;

    % Read until there are no samples left in the input file
    while (~ feof(input_file_handle))
        % Read up to one chunk of samples (could be less if there aren't enough samples left in the file)
        input_samples = fread(input_file_handle, chunk_size * input_file_sample_bytes, input_type);

        % Scale and cast the input samples to the output
        output_samples = cast(input_samples * scalar, output_type);

        % Write the scaled samples to disk and update how many samples have been written
        sample_count = sample_count + fwrite(output_file_handle, output_samples, output_type);
    end

    % Convert to complex sample count
    sample_count = sample_count / 2;
end


text:
function [aveSAR,avePower,aveKGram]=DoSARAverage(SAR,Power,VMass,N_KGram)

global VObj;

%% average SAR
try
    % create 3D searching sphere voxel index list
    Mxdims=size(VObj.MassDen);
    [xgrid,ygrid,zgrid]=meshgrid((-(Mxdims(2)-1)/2)*VObj.XDimRes:VObj.XDimRes:((Mxdims(2)-1)/2)*VObj.XDimRes,...
                                 (-(Mxdims(1)-1)/2)*VObj.YDimRes:VObj.YDimRes:((Mxdims(1)-1)/2)*VObj.YDimRes,...
                                 (-(Mxdims(3)-1)/2)*VObj.ZDimRes:VObj.ZDimRes:((Mxdims(3)-1)/2)*VObj.ZDimRes);
    dist=sqrt(xgrid.^2+ygrid.^2+zgrid.^2);
    [B,idx]=sort(dist(:));
    idx=idx-idx(1);

    aveSAR=zeros(size(VObj.MassDen));
    avePower=zeros(size(VObj.MassDen));
    aveKGram=zeros(size(VObj.MassDen));
    totnum=numel(VObj.MassDen);
    num=numel(find(VObj.MassDen~=0));
    vrange=1;
    for i=(find(VObj.MassDen~=0))'
        tidx=idx+i;
        tidx(tidx<1 | tidx>totnum)=[];
        signal=1;
        direct=0;
        vrange=max(1,min(length(tidx),vrange));
        while vrange>=1 & vrange<=length(tidx)
            mass=sum(VMass(tidx(1:vrange)));
            if mass<N_KGram
                if direct==-1
                    break;
                end
                if signal==1
                    direct=1;
                end
                flag=1;
            else
                if direct==1
                    break;
                end
                if signal==1
                    direct=-1;
                end
                flag=-1;
            end
            vrange=vrange+flag;
            signal=0;
        end
        aveKGram(i)=mass;
        tmpSAR=SAR(tidx(1:max(1,min(length(tidx),vrange))));
        aveSAR(i)=mean(tmpSAR(tmpSAR~=0));
        tmpPower=Power(tidx(1:max(1,min(length(tidx),vrange))));
        avePower(i)=mean(tmpPower(tmpPower~=0));
    end

catch me
    error_msg{1,1}='ERROR!!! SAR averaging process aborted.';
    error_msg{2,1}=me.message;
    errordlg(error_msg);
    return;
end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function linplus_test29 ( )

%*****************************************************************************80
%
%% TEST29 tests R8GE_DET.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 4;
  x = 2.0E+00;
  y = 3.0E+00;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST29\n' );
  fprintf ( 1, '  For a matrix in general storage,\n' );
  fprintf ( 1, '  R8GE_DET computes the determinant.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
%
%  Set the matrix.
%
  for i = 1 : n
    for j = 1 : n
      if ( i == j )
        a(i,j) = x + y;
      else
        a(i,j) = y;
      end
    end
  end
%
%  Factor the matrix.
%
  [ a_lu, pivot, info ] = r8ge_fa ( n, a );
%
%  Compute the determinant.
%
  det = r8ge_det ( n, a_lu, pivot );

  exact = x^(n-1) * ( x + n * y );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  R8GE_DET computes the determinant = %f\n', det );
  fprintf ( 1, '  Exact determinant =                %f\n', exact );

  return
end


text:
function test19
%TEST19 look for NaN's from lchol (caused by Intel MKL 7.x bug)
% Example:
%   test19
% See also cholmod_test

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

fprintf ('=================================================================\n');
fprintf ('test19: look for NaN''s from lchol (caused by Intel MKL 7.x bug)\n') ;

Prob = UFget (936)							    %#ok
A = Prob.A ;
[p count] = analyze (A) ;
A = A (p,p) ;
tic
L = lchol (A) ;
t = toc ;
fl = sum (count.^2) ;
fprintf ('mflop rate: %8.2f\n', 1e-6*fl/t) ;
n = size (L,1) ;
for k = 1:n
    if (any (isnan (L (:,k))))
	k								    %#ok
	error ('!') ;
    end
end

fprintf ('test19 passed; you have a NaN-free BLAS (must not be MKL 7.x...)\n') ;


text:
%% Calculate geodesic error curves
clear all; close all; clc

addpath(genpath('./'))
addpath(genpath('./../Tools/'))

mesh_0 = load('./faust_synthetic/shapes/tr_reg_080'); %Choose the indices of the test pair
mesh_1 = load('./faust_synthetic/shapes/tr_reg_087'); %Choose the indices of the test pair

X = load('./Results/test_faust_synthetic/080_087.mat'); %Choose the indices of the test pair
[~, matches] = max(squeeze(X.softCorr),[],1);

D_model = load('.\faust_synthetic\distance_matrix\tr_reg_087.mat'); %Choose the indices of the test pair
D_model = D_model.D;

gt_matches = 1:6890;
errs = calc_geo_err(matches, gt_matches, D_model);
curve = calc_err_curve(errs, 0:0.001:1.0)/100;
plot(0:0.001:1.0, curve); set(gca, 'xlim', [0 0.1]); set(gca, 'ylim', [0 1])

xlabel('Geodeisc error')
ylabel('Correspondence Accuracy %')

author: OshriHalimi


text:
function [Xi,ODEs] = sparsifyDynamics(Theta,dXdt,LHS_Sym,lambda,N,Sym_Struct,disp,NormalizeLib)
% Copyright 2015, All Rights Reserved
% Code by Steven L. Brunton
% For Paper, "Discovering Governing Equations from Data:
%        Sparse Identification of Nonlinear Dynamical Systems"
% by S. L. Brunton, J. L. Proctor, and J. N. Kutz

% compute Sparse regression: sequential least squares
%%
% Modified By: K
% Last Updated
%% Normalize the library data
if NormalizeLib==1
    % Change thisinto parfor to increase the speed if you have a large
    % library

    %parfor norm_k=1:size(Theta,2)

    % Else use the normal for loop
     for norm_k=1:size(Theta,2)
        normLib(norm_k) = norm(Theta(:,norm_k));
        Theta(:,norm_k) = Theta(:,norm_k)/normLib(norm_k);
    end
end

%% Peform sparse regression
Xi = Theta\dXdt;  % initial guess: Least-squares
[n,m]=size(dXdt);

% lambda is our sparsification knob.
for k=1:N
    smallinds = (abs(Xi)<lambda);   % find small coefficients
    Xi(smallinds)=0;                % and threshold
    for ind = 1:m                   % n is state dimension
        biginds = ~smallinds(:,ind);
        % Regress dynamics onto remaining terms to find sparse Xi
        Xi(biginds,ind) = Theta(:,biginds)\dXdt(:,ind);
    end
end

%% Now output the SINDy Identified ODEs

% Now retrive the parameters
if NormalizeLib==1
    for norm_k=1:length(Xi)
        Xi(norm_k,:) = Xi(norm_k,:)/normLib(norm_k);
    end
end

for i=1:m
     ODEs(i,1)=vpa(cell2sym(Sym_Struct)*Xi(:,i));
end

%% Choose whether you want to display the final discovered equation

if disp==1
     fprintf('The SINDy-PI discovered PDE is:\n')
     digits(6)
     for i=1:m
          fprintf(strcat('\t',char(LHS_Sym),'=',char(simplify(ODEs(i,1))),'\n'));
     end
end

author: dynamicslab

repo: SINDy-PI


text:
function L = obslike_ehmm(ehmm,Gamma,residuals,XX,k)
%
% Evaluate likelihood of data given observation model
% for chain k, for one continuous trial
% (It's pretty much the same than obslike but simplified)
%
% INPUT
% X          N by ndim data matrix
% ehmm        ehmm data structure
% residuals  in case we train on residuals, the value of those.
% XX        alternatively to X (which in this case can be specified as []),
%               XX can be provided as computed by setxx.m
% OUTPUT
% B          Likelihood of N data points
%
% Author: Diego Vidaurre, University of Oxford / Aarhus University (2022)

K = ehmm.K;

[T,ndim] = size(residuals);
setstateoptions;
ltpi = sum(regressed)/2 * log(2*pi);
L = zeros(T+ehmm.train.maxorder,2);

ldetWishB = 0;
PsiWish_alphasum = 0;
for n = 1:ndim % only diagonal?
    if ~regressed(n), continue; end
    ldetWishB = ldetWishB+0.5*log(ehmm.Omega.Gam_rate(n));
    PsiWish_alphasum = PsiWish_alphasum+0.5*psi(ehmm.Omega.Gam_shape);
end
C = ehmm.Omega.Gam_shape ./ ehmm.Omega.Gam_rate;

for l = 1:2

    if l == 1, Gamma(:,k) = 1;
    else, Gamma(:,k) = 0;
    end

    [Xhat,XXstar] = computeStateResponses(XX,ehmm,Gamma,1:K+1);

    d = residuals(:,regressed) - Xhat(:,regressed);
    Cd = bsxfun(@times,C(regressed),d)';
    dist = zeros(T,1);
    for n = 1:sum(regressed)
        dist = dist - 0.5 * (d(:,n).*Cd(n,:)');
    end

    NormWishtrace = zeros(T,1);
%     if ndim == 1
%           NormWishtrace = NormWishtrace + 0.5 * C * ...
%                 sum( (XXstar * ehmm.state_shared(1).W.S_W) .* XXstar, 2);
%     else
%         for n = 1:ndim
%             if ~regressed(n), continue; end
%             Sind_all = repmat(Sind(:,n),K+1,1) == 1;
%             NormWishtrace = NormWishtrace + 0.5 * C(n) * ...
%                 sum( (XXstar(:,Sind_all) * ehmm.state_shared(n).S_W(Sind_all,Sind_all)) ...
%                 .* XXstar(:,Sind_all), 2);
%         end
%     end

    L(ehmm.train.maxorder+1:end,l) = - ltpi - ldetWishB + ...
        PsiWish_alphasum + dist - NormWishtrace;

end

Lneg = all(L<0,2);
if any(Lneg)
    L(Lneg,:) = L(Lneg,:) - repmat(max(L(Lneg,:),[],2),1,2);
end
%if any(sum(exp(L),2)==0), keyboard; end

L = exp(L);

end

author: OHBA-analysis

repo: HMM-MAR


text:
function [states] = tapas_vlinear_init_states(data, model, inference)
%%
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

states = cell(inference.niter, 1);

end

author: translationalneuromodeling

repo: tapas


text:
function model=rmInterpolate(view,model,params)
% rmInterpolate - interpolate missing values across the Gray surface
% that were not sampled with the coarse sampling
%
% model = rmInterpolate(view,model,params);
%
% 2007/03 SOD: wrote it.

if notDefined('view'),   error('Need view');   end;
if notDefined('model'),  error('Need model');  end;
if notDefined('params'), error('Need params'); end;

numNeighbors = double(view.nodes(4,:));
myzeros      = zeros(size(numNeighbors));
grayConMat   = [];

roiCoords = rmGet(params,'roiCoords');
if ~isempty(roiCoords),
    % find the index' that cover the roi
    doROI = true;
    inROI = rmGet(params,'roiIndex');
    % coords with data within roi
    coarseIndex     = rmCoarseSamples(roiCoords,params.analysis.coarseSample);
    withData        = logical(myzeros);
    withData(inROI(coarseIndex)) = true;

    % coords without data within roi
    toInterp        = logical(myzeros);
    toInterp(inROI(~coarseIndex)) = true;
else
    doROI    = false;
    coarseIndex     = rmCoarseSamples(viewGet(view,'coords'),params.analysis.coarseSample);
    % coords with data
    withData = coarseIndex;
    % coords without data
    toInterp = ~coarseIndex;
end;

% Sanity check: we only need interpolation if there are points to
% interpolate, e.g. very small ROIs will not be coarsely sampled.
if all(coarseIndex), % if all data points have estimates,
    return;
end;

% now we have to loop over models and parameters
%rp = {'x','y','s','rss','rawrss','x02','y02','s2','rss2','rawrss2', 'exponent','sigmamajor','sigmaminor'};
for n=1:numel(model),
    if strcmp(model{n}.description,'radial oval 2D pRF fit (x,y,sigma_major,sigma_minor)')
        rp = {'x','y','rss','rawrss','x02','y02','rss2','rawrss2', 'exponent','sigmamajor','sigmaminor','theta'};
    else
        rp = {'x','y','s','rss','rawrss','x02','y02','s2','rss2','rawrss2', 'exponent'};
    end

    % reset for different models (should be the same i think)
    fulval = myzeros;
    for p=1:numel(rp),
        val = rmGet(model{n},rp{p});
        if ~isempty(val) && numel(val)==sum(coarseIndex),
            % put data in full datastructure and fill the rest with zeros
            if doROI,
                fulval(inROI(coarseIndex)) = double(val);
            else
                fulval(coarseIndex)        = double(val);
            end;
            % now interpolate the data to some voxels without data
            %newval = myinterp(fulval,withData,toInterp,edges,numNeighbors,edgeOffsets);
            [newval grayConMat] = myinterp(fulval,withData,toInterp,view,grayConMat);
            % if roi we should store the roi voxels only otherwise store
            % the whole thing
            if doROI,
                model{n} = rmSet(model{n},rp{p},newval(inROI));
            else
                model{n} = rmSet(model{n},rp{p},newval);
            end;
        elseif ~isempty(val) && numel(val)==numel(fulval),
            model{n} = rmSet(model{n},rp{p},val);
        end;
    end;
    if isfield(model{n}, 'exponent') && ~isempty(model{n}.exponent) && ~strcmp(model{n}.description,'radial oval 2D pRF fit (x,y,sigma_major,sigma_minor)')
        model{n} = rmSet(model{n}, 's', model{n}.sigma.major .* sqrt(model{n}.exponent));
    end
end;

% separate loop for the betas, because they are saved in a different way
for n=1:numel(model),
    % reset for different models (should be the same i think)
    fulval = myzeros;
    val = double(rmGet(model{n},'b'));
    nBetas  = size(val,3);
    newBeta = zeros(size(val,1),numel(coarseIndex),nBetas);

    for ii=1:nBetas,
        if doROI,
            fulval(inROI(coarseIndex)) = double(val(:,:,ii));
        else
            fulval(coarseIndex) = double(val(:,:,ii));
        end
        newval = myinterp(fulval,withData,toInterp,view,grayConMat);
        % if roi we should store the roi voxels only otherwise store
        % the whole thing
        if doROI,
            newBeta(:,:,ii) = newval(inROI);
        else
            newBeta(:,:,ii) = newval;
        end;
    end;
    model{n} = rmSet(model{n},'b',newBeta);
end;

return;
%---------------------------------------

%---------------------------------------
function [val grayConMat] = myinterp(val,withData,toInterp,view,grayConMat)
% actual interpolation function

lzeros     = false(size(withData));
doTouch    = toInterp;  % keep refining these estimates
oldtoInterp = toInterp;
while (1),
    % 'smooth' data. Points without data are zero and do contribute to this
    % initial estimate:
    [tmp grayConMat] = dhkGraySmooth(view,val,[1 1],grayConMat);

    % Now check which surrounding points had data, 0=no data in
    % neighborhood, 1=all neighbors have data, intermediate is the weight
    % that should have been used to create the new estimate and will be
    % used to remove the contributions of no-data neighbors (0).
    c2               = dhkGraySmooth(view,double(withData),[1 1],grayConMat);

    % All point that we want to interpolate (doTouch) that have neighbors
    % with data and that thus will have a new estimate.
    toInterp          = lzeros;
    toInterp(doTouch) = c2(doTouch)>0.001;

    % Compute new final values taking into account the number of points
    % that had data (c2)
    val(toInterp)     = tmp(toInterp)./c2(toInterp);

    % update toInterp: find voxels that still need interpolating
    toInterp       = (oldtoInterp-toInterp)>0;

    % sanity check: if no changes in toInterp than quit. Ideally
    % toInterp should go to zero but this does not have to be the
    % case if some voxels are not connected to voxels with data.
    if any(oldtoInterp-toInterp) || sum(toInterp)==0,
        break;
    else
        % 'grow' logical matrix that keeps track of which voxels have data
        withData       = toInterp | withData;
        % update
        oldtoInterp = toInterp;
    end;
end;

return;
%---------------------------------------


text:
clear;
options.J = 1;
options.M = 2;
options.L = 3;

field_list = {'J', 'M', 'aa'};

%%
options2 = sub_options(options, field_list);

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function varargout=abs(varargin)
%ABS (overloaded)
%
% t = abs(x)
%
% The variable t can only be used in convexity preserving
% operations such as t<=0, minimize t etc.

%% ***************************************************
% This file defines a nonlinear operator for YALMIP
%
% It can take three different inputs
% For double inputs, it returns standard double values
% For sdpvar inputs, it genreates a an internal variable
% When first input is 'model' it generates the epigraph

%% ***************************************************
switch class(varargin{1})

    case 'double'
        error('Overloaded SDPVAR/ABS CALLED WITH DOUBLE. Report error')

    case 'sdpvar' % Overloaded operator for SDPVAR objects. Pass on args and save them.
        if isreal(varargin{1})
            varargout{1} = yalmip('addextendedvariable',mfilename,varargin{1});
        else
            % For complex args, abs(X) is defined [norm(X(i,j),2)] in MATLAB
            y = [];
            x = varargin{1};
            for i = 1:size(x,1)
                temp = [];
                for j = 1:size(x,2)
                    temp = [temp norm(extsubsref(x,i,j))];
                end
                y = [y;temp];
            end
            varargout{1} = y;
        end

    case 'char' % YALMIP send 'graph' when it wants the epigraph or hypograph
        switch varargin{1}
            case 'graph'
                % Description using epigraphs
                t = varargin{2};
                X = varargin{3};
                varargout{1} = (-t <= X <= t);
                varargout{2} = struct('convexity','convex','monotonicity','none','definiteness','positive');
                varargout{3} = X;

            case 'milp'
                % Exact description using binary variables
                t = varargin{2};
                X = varargin{3};
                F = ([]);
                [M,m]=derivebounds(X);
                if m>=0
                    F = F + (t == X);
                elseif M<0
                    F = F + (t == -X);
                else
                    d = binvar(1,1);
                    F = F + (X <= M*d)     + (-2*(M-m)*d     <= t+X <= 2*(M-m)*d);
                    F = F + (X >= m*(1-d)) + (-2*(M-m)*(1-d) <= t-X <= 2*(M-m)*(1-d));
                end
                varargout{1} = F;
                varargout{2} = struct('convexity','milp','monotonicity','milp','definiteness','positive');
                varargout{3} = X;
            otherwise
                error('SDPVAR/ABS called with CHAR argument?');
        end
    otherwise
        error('Strange type on first argument in SDPVAR/ABS');
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function svmPlotConfusionMatrix(confusionImage,grouplabels,labelsoff)
% Function to plot a confusion image that was created in
% svmConfusionMatrix.

if notDefined('labelsoff'), labelsoff = 0; end
nConditions = length(grouplabels);

figure;
% Clean up the figure junk matlab generates
set(gcf, 'DockControls', 'off');
set(gcf, 'ToolBar', 'none');

% Plot data and begin tidying it up
imagesc(confusionImage);
axis square; % Correct for stretching
set(gca, 'TickLength', [0 0]); % No ticks

if ~labelsoff
    set(gca, 'YTickLabel', grouplabels);
else
    set(gca, 'YTickLabel',[]);
end

% Create rotated labels on the confusion matrix
% Inspired by a script from Andrew Bliss on MATLAB Central, nice!
% http://www.mathworks.com/matlabcentral/fileexchange/8722-rotate-tick-label
set(gca, 'XTickLabel', []); % Clear original x ticks, replace with text command
if ~labelsoff
    xVals = 1:nConditions;
    yVals = repmat(nConditions + (.1 * nConditions), nConditions, 1);
    text(xVals, yVals, grouplabels, 'HorizontalAlignment', 'right', 'rotation', 45, 'FontWeight', 'bold');
    set(gca, 'FontWeight', 'bold');
end

% Shift over and up a bit, to make room for axis labels and legend
set(gca, 'OuterPosition', [0 .05 1 1]);

xLabelPosition = get(get(gca, 'XLabel'), 'Position');
xLabelPosition(2) = nConditions + (.25 * nConditions); % Set position of x label lower than normal
xlabel('Prediction', 'FontSize', 15, 'FontWeight', 'bold', 'Position', xLabelPosition);
ylabel('Query', 'FontSize', 15, 'FontWeight', 'bold');
colorbar;
%colormap gray
colormap autumn


text:
function indx = r82vec_sort_heap_index_a ( n, a )

%*****************************************************************************80
%
%% R82VEC_SORT_HEAP_INDEX_A does an indexed heap ascending sort of an R82VEC.
%
%  Discussion:
%
%    The sorting is not actually carried out.  Rather an index array is
%    created which defines the sorting.  This array may be used to sort
%    or index the array, or to sort or index related arrays keyed on the
%    original array.
%
%    Once the index array is computed, the sorting can be carried out
%    "implicitly:
%
%      A(1:2,INDX(I)), I = 1 to N is sorted,
%
%    or explicitly, by the call
%
%      A = R82VEC_PERMUTE ( N, A, INDX )
%
%    after which A(1:2,I), I = 1 to N is sorted.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the array.
%
%    Input, real A(2,N), an array to be index-sorted.
%
%    Output, integer INDX(N), the sort index.  The
%    I-th element of the sorted array is A(1:2,INDX(I)).
%
  if ( n < 1 )
    return
  end

  if ( n == 1 )
    indx(1) = 1;
    return
  end

  indx = i4vec_indicator ( n );

  l = floor ( n / 2 ) + 1;
  ir = n;

  while ( 1 )

    if ( 1 < l )

      l = l - 1;
      indxt = indx(l);
      aval(1:2) = a(1:2,indxt);

    else

      indxt = indx(ir);
      aval(1:2) = a(1:2,indxt);
      indx(ir) = indx(1);
      ir = ir - 1;

      if ( ir == 1 )
        indx(1) = indxt;
        break
      end

    end

    i = l;
    j = l + l;

    while ( j <= ir )

      if ( j < ir )
        if (   a(1,indx(j)) <  a(1,indx(j+1)) | ...
             ( a(1,indx(j)) == a(1,indx(j+1)) & ...
               a(2,indx(j)) <  a(2,indx(j+1)) ) )
          j = j + 1;
        end
      end

      if (   aval(1) <  a(1,indx(j)) | ...
           ( aval(1) == a(1,indx(j)) & ...
             aval(2) <  a(2,indx(j)) ) )
        indx(i) = indx(j);
        i = j;
        j = j + j;
      else
        j = ir + 1;
      end

    end

    indx(i) = indxt;

  end

  return
end


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function MU_funcACos(Temp,Event,handles)
handles = guidata(handles.MU_matrix_display);

MergeM=get(handles.Matrix_name_edit,'String');
if isfloat(handles.TMatrix(1))
    set(handles.Matrix_name_edit,'String',['acos([' MergeM '])']);
else
    set(handles.Matrix_name_edit,'String',['acos(double([' MergeM ']))']);
end
MU_calc_matrix(handles);

end


text:
function initialize_native()
% initialize_native Initialize all native components
%
% A script that compiles all native components (MEX functions) and places
% them in a given output directory.
%
% Input:
% - output_path (string): Path to output directory.
%

toolkit_path = get_global_variable('toolkit_path');
output_path = get_global_variable('native_path');

trax_path = get_global_variable('trax_source', fullfile(output_path, 'trax'));
if ~download_trax_source(trax_path)
    error('Unable to compile native resources.');
end;

print_text('Verifying native components ...');

success = true;

include_paths = {fullfile(trax_path, 'src'), fullfile(trax_path, 'include')};

success = success && compile_mex('region_overlap', {fullfile(toolkit_path, 'sequence', 'region_overlap.cpp'), ...
    fullfile(trax_path, 'src', 'region.c')}, include_paths, output_path, '-DTRAX_STATIC_DEFINE');

success = success && compile_mex('region_mask', {fullfile(toolkit_path, 'sequence', 'region_mask.cpp'), ...
    fullfile(trax_path, 'src', 'region.c')}, include_paths, output_path, '-DTRAX_STATIC_DEFINE');

success = success && compile_mex('region_convert', {fullfile(toolkit_path, 'sequence', 'region_convert.cpp'), ...
    fullfile(trax_path, 'src', 'region.c')}, include_paths, output_path, '-DTRAX_STATIC_DEFINE');

success = success && compile_mex('read_trajectory', {fullfile(toolkit_path, 'sequence', 'read_trajectory.cpp'), ...
    fullfile(trax_path, 'src', 'region.c')}, include_paths, output_path, '-DTRAX_STATIC_DEFINE');

success = success && compile_mex('write_trajectory', {fullfile(toolkit_path, 'sequence', 'write_trajectory.cpp'), ...
    fullfile(trax_path, 'src', 'region.c')}, include_paths, output_path, '-DTRAX_STATIC_DEFINE');

success = success && compile_mex('benchmark_native', {fullfile(toolkit_path, 'tracker', 'benchmark_native.cpp')}, ...
    {}, output_path);

success = success && compile_mex('md5hash', {fullfile(toolkit_path, 'utilities', 'md5hash.cpp')}, ...
    {}, output_path);

trax_mex_path = fullfile(output_path, 'mex');
mkpath(trax_mex_path);

success = success && compile_mex('traxserver', ...
    {fullfile(trax_path, 'support', 'matlab', 'traxserver.cpp'), ...
     fullfile(trax_path, 'support', 'matlab', 'helpers.cpp'), ...
     fullfile(trax_path, 'src', 'trax.c'), ...
     fullfile(trax_path, 'src', 'region.c'), ...
     fullfile(trax_path, 'src', 'strmap.c'), ...
     fullfile(trax_path, 'src', 'message.c'), ...
     fullfile(trax_path, 'src', 'traxpp.cpp'), ...
     fullfile(trax_path, 'src', 'base64.c')}, ...
    {fullfile(trax_path, 'src'), fullfile(trax_path, 'include')}, trax_mex_path, '-DTRAX_STATIC_DEFINE');

% Additional OS-specific flags for traxclient MEX
os_specific = {};
if isunix() && ~ismac()
    % clock_gettime() in trax/support/client/timer.cpp requires librt on
    % linux systems with glibc < 2.17 (so to be safe, we always add it)
    os_specific{end+1} = '-lrt';
end

success = success && compile_mex('traxclient', ...
    {fullfile(trax_path, 'support', 'matlab', 'traxclient.cpp'), ...
     fullfile(trax_path, 'support', 'matlab', 'helpers.cpp'), ...
     fullfile(trax_path, 'support', 'client', 'client.cpp'), ...
     fullfile(trax_path, 'support', 'client', 'process.cpp'), ...
     fullfile(trax_path, 'support', 'client', 'threads.cpp'), ...
     fullfile(trax_path, 'support', 'client', 'timer.cpp'), ...
     fullfile(trax_path, 'src', 'trax.c'), ...
     fullfile(trax_path, 'src', 'traxpp.cpp'), ...
     fullfile(trax_path, 'src', 'region.c'), ...
     fullfile(trax_path, 'src', 'strmap.c'), ...
     fullfile(trax_path, 'src', 'message.c'), ...
     fullfile(trax_path, 'src', 'base64.c')}, ...
    {fullfile(trax_path, 'src'), ...
     fullfile(trax_path, 'include'), ...
     fullfile(trax_path, 'support', 'client', 'include'), ...
     fullfile(trax_path, 'support', 'matlab')}, output_path, '-DTRAX_STATIC_DEFINE', os_specific{:});

if ~success
    error('Unable to compile all native resources.');
end;

set_global_variable('trax_mex', trax_mex_path);

end

function success = download_trax_source(trax_dir)
% download_trax_source Download external components from TraX repository.
%
% To reduce redundant code, a part of the source for MEX files is provided
% by the TraX library. This function downloads and unpacks the source of
% the library and places it in a desired directory.
%
% Input:
% - trax_dir (string): Path to the destination directory.
%
% Output:
% - success (boolean): True on success.
%

trax_branch = get_global_variable('trax_source_branch', 'master');

trax_url = get_global_variable('trax_url', sprintf('https://codeload.github.com/votchallenge/trax/zip/%s', trax_branch));

if isempty(trax_url)
    success = false;
    return;
end;

trax_header = fullfile(trax_dir, 'include', 'trax.h');

if ~exist(trax_header, 'file')
    print_text('Downloading TraX source from "%s". Please wait ...', trax_url);
    working_directory = tempname;
    mkdir(working_directory);
    mkdir(trax_dir);
    bundle = [tempname, '.zip'];
    try
        urlwrite(trax_url, bundle);
        unzip(bundle, working_directory);
        delete(bundle);
        copyfile(fullfile(working_directory, sprintf('trax-%s', trax_branch), '*'), trax_dir);
        success = true;
		delpath(working_directory);
    catch
        print_text('Unable to unpack TraX source code.');
        success = false;
    end;
    delpath(working_directory);
else
    print_debug('TraX source code already present.');
    success = true;
end;

if success
    if exist(fullfile(trax_dir, 'python'), 'dir') == 7
        set_global_variable('trax_python', fullfile(trax_dir, 'support', 'python'));
    end
end

end

author: votchallenge

repo: toolkit-legacy


text:
function out = optPriceVal(type, Price, Strike, Rate, Time, Vol, Yield, Strike2, Spread)
% OPTPRICE: Pricing function for basic and composite European options.
%
% OPTPRICE provides the means to calculate the price of certian options for
% which a closed-form solution via the Black-Scholes model is known.  It is
% meant to be called using the GUI in one of the following syntaxes:
%

% Sanitize the inputs
if any(isnan(Price(:))) || ~isnumeric(Price) || isempty(Price) || ~isreal(Price) || any(Price(:) < 0)
    error('optPriceVal:InvalidPrice', ...
        'Spot price must be a positive number')
elseif any(isnan(Strike(:))) || ~isnumeric(Strike) || isempty(Strike) || ~isreal(Strike) || any(Strike(:) < 0)
    error('optPriceVal:InvalidStrike', ...
        'Strike price must be a positive number')
elseif any(isnan(Rate(:))) ||~isnumeric(Rate) || isempty(Rate) || ~isreal(Rate) || any(Rate(:) < 0)
    error('optPriceVal:InvalidRate', ...
        'Risk free rate must be a positive number')
elseif any(isnan(Time(:))) ||~isnumeric(Time) || isempty(Time) || ~isreal(Time) || any(Time(:) < 0)
    error('optPriceVal:InvalidTime', ...
        'Time to expiration must be a positive number')
elseif any(isnan(Vol(:))) ||~isnumeric(Vol) || isempty(Vol) || ~isreal(Vol) || any(Vol(:) < 0)
    error('optPriceVal:InvalidVol', ...
        'Volatility must be a positive number')
elseif any(isnan(Yield(:))) ||~isnumeric(Yield) || isempty(Yield) || ~isreal(Yield) || any(Yield(:) < 0)
    error('optPriceVal:InvalidYield', ...
        'Yield must be a positive number')
end
% Optional sanitation of inputs, depending on the option:
switch type
    case {'Strangle', 'Bull Spread', 'Bear Spread'}
        if isnan(Strike2) || ~isnumeric(Strike2) || isempty(Strike2) || ~isreal(Strike2) || Strike2 < 0
            error('optPriceVal:InvalidStrike', ...
                'Strike price must be a positive number')
        end
    case 'Butterfly'
        if isnan(Spread) || ~isnumeric(Spread) || isempty(Spread) || ~isreal(Spread) || Spread < 0 || Spread > 1
            error('optPriceVal:InvalidSpread', ...
                'Spread must be a positive number between 0 and 1')
        end
end

switch type
    case 'Call'
        out = blsprice(Price, Strike, Rate, Time, Vol, Yield);
    case 'Put'
        [dummy, out] = blsprice(Price, Strike, Rate, Time, Vol, Yield);
    case 'Straddle'
        % A Straddle is just a Call and Put option bought at the same
        % strike price and with the same dates involved.
        [Call, Put] = blsprice(Price, Strike, Rate, Time, Vol, Yield);
        out = Call + Put;
    case 'Strangle'
        % A strangle is similar to the straddle, except that there are two
        % strike prices involved: one for the put, and the other for the
        % call.
        Call = blsprice(Price, Strike2, Rate, Time, Vol, Yield);
        [dummy, Put] = blsprice(Price, Strike, Rate, Time, Vol, Yield);
        out = Call + Put;
    case 'Bull Spread'
        % A Bull Spread is formed by buying a call option with a given
        % strike price and then selling another call option with a higher
        % strike price.
        Call1 = blsprice(Price, Strike, Rate, Time, Vol, Yield);
        Call2 = blsprice(Price, Strike2, Rate, Time, Vol, Yield);
        out = Call1 - Call2;
        % Strike2 should be larger than Strike; if not, we'll just
        % tacitly invert them rather than giving an error message.
        if Strike2 < Strike
            out = -out;
        end
    case 'Bear Spread'
        % A Bear Spread is formed by buying a call option with a given
        % strike price and then selling another call option with a lower
        % strike price.
        Call1 = blsprice(Price, Strike, Rate, Time, Vol, Yield);
        Call2 = blsprice(Price, Strike2, Rate, Time, Vol, Yield);
        out = Call2 - Call1;
        % Strike2 should be larger than Strike; if not, we'll just
        % tacitly invert them rather than giving an error message.
        if Strike2 < Strike
            out = -out;
        end
    case 'Butterfly'
        % A Butterfly option involves 4 call options at differing strike
        % prices:
        % Long 1 call at (1 - Spread)*Strike
        % Short 2 calls at Strike
        % Long 1 call at (1 + Spread)*Strike
        Call1 = blsprice(Price, (1-Spread)*Strike, Rate, Time, Vol, Yield);
        Call2 = blsprice(Price, Strike, Rate, Time, Vol, Yield);
        Call3 = blsprice(Price, (1+Spread)*Strike, Rate, Time, Vol, Yield);
        out = Call1 - 2*Call2 + Call3;
    otherwise
        error('optPriceVal:UnknownOption', ...
        ['Invalid option type.  Valid choices are Call, Put, Straddle, '...
        'Strangle, Bull Spread, Bear Spread, or Butterfly'])
end

author: Sable

repo: mcbench-benchmarks


text:
function [ int, count, msg ] = pnmgeti( fid, n )
%PNMGETI Get integers from an ascii encoded PBM/PGM/PPM file.
%
%   [ INT, COUNT, MSG ] = PNMGETI( FID, N ) tries to read N integers
%   from the ascii encoded PBM/PGM/PPM file with file identifier FID and
%   returns the integers in the vector INT. COUNT is an optional output
%   argument that returns the number of elements successfully read. MSG
%   is an optional output argument that returns an error message string
%   if an error occurred or an empty matrix if an error did not occur.
%
%   If N is omitted, the whole remaining of the file is read.
%
%   The main difference between PNMGETI( FID ) and FSCANF( FID, '%d' )
%   is that PNMGETI ignores comments (from # to end of line). PNMGETI
%   also ignores garbage, which is anything that is neither whitespace,
%   digit nor comment.

%   Author:      Peter J. Acklam
%   Time-stamp:  1998-04-06 02:38:34
%   E-mail:      jacklam@math.uio.no (Internet)
%   URL:         http://www.math.uio.no/~jacklam

%
% Check number of input arguments and assign default value to omitted
% argument.
%
error( nargchk( 1, 2, nargin ) );
if nargin < 2
   n = Inf;
end

% Initialize output arguments.
int   = [];             % Image data vector.
count = 0;              % Number of elements read. Same as length(int).
msg   = '';             % Error message string.

while 1

   % Calculate number of integers missing and try to read that many.
   ints_missing = n - count;
   [ x, this_count ] = fscanf( fid, '%d', ints_missing );

   % Append new data to main data vector and increment counter.
   int = [ int ; x ];
   count = count + this_count;

   % Return if we have got the desired number of elements.
   if count == n
      return
   end

   % Return if we have reached EOF.
   if feof( fid )
      msg = 'End of file reached too early.';
      return
   end

   %
   % If we get here we have reached a comment or some garbage. Garbage
   % is anything that is neither whitespace, digit nor comment.
   %

   char = fscanf( fid, '%c', 1 );
   if ( char == '#' )

      % Found a comment, so skip the rest of the line and redo the loop.
      fgetl( fid );

   else

      % We found some garbage, so give a message.
      msg = 'Garbage found where image data was expected.';

      %
      % Read past the garbage and following whitespace (i.e., until
      % first number character). This would probably be faster if we
      % read and checked a whole vector of data, but it probably won't
      % be executed often. char is empty at EOF.
      %
      while ~isempty( char ) & ( char < '0' | char > '9' )
         char = fscanf( fid, '%c', 1 );
      end

      %
      % Return if we have reached EOF. This error message may overwrite
      % any error message telling about garbage, but that doesn't matter
      % since reaching EOF too early is a more serious error.
      %
      if feof( fid )
         msg = 'End of file reached too early.';
         return
      end

      %
      % We found a number character, so jump back one byte so fscanf
      % catches it.
      %
      fseek( fid, -1, 0 );

   end

end


text:
function outname = writeww3(outfiname,EToV,VX,B,opedat, title)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function writes WW3 grid including nodes             %
% (longitude,latitude,depth), open bounday nodes            %
% and element connections (triangles)                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT:                                                    %
% EToV: triangle (nelem,3);                                 %
% VX(:,1): longitude (nnode,1);                             %
% VX(:,2): laritude (nnode,1);                              %
% B: depth (nnode,1);                                       %
% opedat.nbdv: Open boundary nodes ID                       %
% outfiname: mesh name                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%      Ali Abdolali August 2018 ali.abdolali@noaa.gov       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('write WW3 msh') ;

node(:,1)=(1:length(VX(:,1)));
node(:,2)=VX(:,1);
node(:,3)=VX(:,2);
node(:,4)=B;

fid = fopen(outfiname,'w');
outname = outfiname ;
disp( title )  ;
fprintf(fid,'%s\n', '$MeshFormat');
fprintf(fid,'%s\n', '2 0 8');
fprintf(fid,'%s\n', '$EndMeshFormat');
fprintf(fid,'%s\n', '$Nodes');
fprintf(fid,'%d\n', length(node(:,1)));
   for i=1:length(node(:,1))
        fprintf(fid,['%d %s %5.5f %s %5.5f %s %5.5f\n'], node(i,1),'', VX(i,1),'',VX(i,2),'',B(i,1));
    end
fprintf(fid,'%s\n', '$EndNodes');
fprintf(fid,'%s\n', '$Elements');
fprintf(fid,'%d\n', length(EToV(:,1))+opedat.neta);
m=0;
for i=1:opedat.nope
    for j=1:opedat.nvdll(i)
        m=m+1;
        fprintf(fid,['%d %s %d %s %d %s %d %s %d %s %d\n'], m,'',15,'',2,'',0,'',0,'',opedat.nbdv(j,i));
    end
end
for i=1:length(EToV(:,1))
   m=m+1;
   fprintf(fid,['%d %s %d %s %d %s %d %s %d %s %d %s %d %s %d %s %d\n'], m,'',2,'',3,'',0,'',i,'',0,'',EToV(i,1),'',EToV(i,2),'',EToV(i,3));
end
fprintf(fid,'%s', '$EndElements');
fclose(fid) ;
return

author: CHLNDDEV


text:
% Profiler extension template.
%
% Profiler extension for [algorithm]
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab.
% https://github.com/steven2358/kafbox/

classdef kafbox_template_profiler < kafbox_template

    properties (GetAccess = 'public', SetAccess = 'private')
        elapsed = 0; % elapsed time
    end

    methods

        function kaf = kafbox_template_profiler(parameters) % constructor
            if nargin<1, parameters = struct(); end
            kaf = kaf@kafbox_template(parameters);
        end

        function flops = lastflops(kaf) % flops for last iteration
            % numbers of operations
            m1 = 1;
            m2 = 1;
            m3 = 1;
            m4 = 1;

            floptions = struct(...
                'sum', m1, ...
                'mult', m2, ...
                'div', m3, ...
                sprintf('%s_kernel',kaf.kerneltype), [m4,1,size(kaf.dict,2)]);

            flops = kflops(floptions);
        end

        %% flops breakdown

        % [space for remarks on number of operations used above]

        %%

        function train_profiled(kaf,x,y)
            t1 = tic;
            kaf.train(x,y);
            t2 = toc(t1);
            kaf.elapsed = kaf.elapsed + t2;
        end

        function bytes = lastbytes(kaf) % bytes used in last iteration
            m = size(kaf.dict,1);
            m2 = 1;
            bytes = (m2 + m*size(kaf.dict,2)); % 8 bytes for double precision
            % [list variables]
        end

    end
end

repo: kafbox


text:
function y=obj1(x)

global counterf
counterf=counterf+1;

y=x^2+2*x;

repo: Optimization-Theory-and-Algorithm


text:
classdef testPrintingInputFileForFem < testNotShowingError

    properties (Access = private)
        fileName
        gmsFile
        resultsDir
        fullFileName
    end

    methods (Access = public)

        function obj = testPrintingInputFileForFem()
            obj.init();
            g  = obj.gmsFile;
            oD = obj.resultsDir;
            oF = obj.fullFileName;
            c = GmsFile2SwanFileConverter(g,oD,oF);
            c.convert();
        end

    end

    methods (Access = protected)

        function hasPassed = hasPassed(obj)
            sF = fullfile('tests','PrintingTests',['test',obj.fileName],['test',obj.fileName,'.m']);
            oF = obj.fullFileName;
            hasChanged = FileComparator().areFilesDifferent(sF,oF);
            hasPassed = ~hasChanged;
        end

    end

    methods (Access = private)

        function init(obj)
            obj.fileName = 'InputFileForFem';
            obj.gmsFile  = 'testReadingGmsh.msh';
            obj.resultsDir   = fullfile('Output',obj.fileName);
            obj.fullFileName = fullfile(obj.resultsDir,[obj.fileName,'.m']);
        end

    end

end


text:
function model = svmTrain(X, Y, opts)
    % trains an SVM using minFunc
    % X is NxD features
    % Y is discrete array of labels that must (for now) be 1...K

    % opts.type specifies which method to use.
    % 1= linear SVM
    % 2= polynomial Kernel SVM of order opts.polyOrder
    % 3= rbf Kernel SVM with scale parameter opts.rbfScale
    %    the scale is used in equation Z = 1/sqrt(2*pi*sigma^2);
    %    to scale the output of the exponential, where sigma is the scale
    % 4= l2svm, which unlike normal svm uses squared hinge loss

    % model is to be plugged in to yhat = svmTest(model, X)

    % TODO: Support labels that are not in 1..K
    % TODO: n-fold Cross validation support

    addOnes= false;
    C= 1;
    type= 1;
    rbfScale= 1;
    polyOrder= 2;
    if nargin < 3, opts= struct; end
    if isfield(opts, 'addOnes'), addOnes= opts.addOnes; end
    if isfield(opts, 'C'), C= opts.C; end
    if isfield(opts, 'type'), type= opts.type; end
    if isfield(opts, 'rbfScale'), rbfScale= opts.rbfScale; end
    if isfield(opts, 'polyOrder'), polyOrder= opts.polyOrder; end

    if(addOnes), X= [X, ones(size(X,1), 1)]; end

    [N, D]= size(X);
    u= unique(Y);
    K = length(u);

    minFuncOpts= struct;
    minFuncOpts.display= 0;

    if type == 1

        % Linear SVM
        funObj = @(w) SSVMMultiLoss(w, X, Y, K);
        wLinear = minFunc(@penalizedL2, zeros(D*K,1), minFuncOpts, funObj, C);
        model.w = reshape(wLinear,[D K]);

    elseif type == 2

        % Polynomial SVM
        Kpoly = kernelPoly(X, X, polyOrder);
        funObj = @(v) SSVMMultiLoss(v, Kpoly, Y, K);
        uPoly = minFunc(@penalizedKernelL2_matrix, randn(N*K,1), minFuncOpts, Kpoly, K, funObj, C);
        model.X = X; % must save all the training data... (ok only support vectors, todo)
        model.uPoly= reshape(uPoly, [N, K]);
        model.polyOrder= polyOrder;

    elseif type == 3

        % RBF SVM
        Krbf = kernelRBF(X, X, rbfScale);
        funObj = @(v) SSVMMultiLoss(v, Krbf, Y, K);
        uRBF = minFunc(@penalizedKernelL2_matrix, randn(N*K,1), minFuncOpts, Krbf, K, funObj, C);
        model.X = X; % must save all the training data... (ok only support vectors, todo)
        model.urbf= reshape(uRBF,[N, K]);
        model.rbfScale= rbfScale;

    elseif type == 4

        % L2 SVM (squares the slack variables, i.e. squared hinge loss)
        funObj = @(v) l2svmloss(v, X, Y, K, C);
        wLinear = minFunc(@penalizedL2, zeros(D*K,1), minFuncOpts, funObj, C);
        model.w = reshape(wLinear,[D K]);

    else

        fprintf('Unrecognized type %d, exitting...\n', type);
        model= struct;
        return;

    end

    model.type= type;
    model.addOnes= addOnes;
    model.C= C;
    model.u= u;

    % 1-vs-all L2-svm loss function;  similar to LibLinear.
    % Originally taken from Adam Coates' code, slightly adapted
    function [loss, g] = l2svmloss(w, X, y, K, C)

      [M,N] = size(X);
      theta = reshape(w, N,K);
      Y = bsxfun(@(y,ypos) 2*(y==ypos)-1, y, 1:K);
      margin = max(0, 1 - Y .* (X*theta));
      loss = (0.5 * sum(theta.^2)) + C*sum(margin.^2);
      loss = sum(loss);
      g = theta - 2 * C * (X' * (margin .* Y));
      g = g(:);

author: karpathy

repo: Random-Forest-Matlab


text:
function [x, y, z]=splat(r, theta, phi)
% % mapping, delaunay, flat projection
% %
% % Syntax;
% %
% % [x, y, z]=splat(r, theta, phi);
% %
% % ***********************************************************
% %
% % Description
% %
% % Program maps the spherical coordinates r, theta, and phi from a
% % sphere to a flat plane to aid in triangularization of hemispherical
% % coordinates.
% %
% % ***********************************************************
% %
% % Input Variables
% %
% % r, theta, phi, are matrices of spherical coordinates.
% %
% % ***********************************************************
% %
% % Output Variables
% %
% % x, y, and z, are matrices of rectangular coordinates.
% %
% % ***********************************************************
% %
% Example
%
% phi=pi/4:pi/320:pi/2;
% theta=0:(pi/16):(5*pi);
% r=ones(size(phi));
% [x, y, z]=spherical_to_rectangular(r, theta, phi);
% [x10, y10, z]=splat(r, theta, phi);
% TRI = delaunay(x10, y10);
%
% figure(1);
% h=trisurf(TRI,x,y,z, 'FaceColor', [1 1 1], 'EdgeColor', 0.5*[1 1 1],...
% 'LineWidth', 1 );
% hold on;
% plot3(x,y,z, '-r', 'Linewidth', 3 );
%
% TRI2 = delaunay(x, y);
%
% figure(2);
% h=trisurf(TRI2,x,y,z, 'FaceColor', [1 1 1], 'EdgeColor', 0.5*[1 1 1],...
% 'LineWidth', 1 );
% hold on;
% plot3(x,y,z, '-r', 'Linewidth', 3 );
%
% % ***********************************************************
% %
% % This program was written by Edward L. Zechmann
% %
% %     date     January 2007
% %
% % modified 11   March   2008  added comments
% %
% % ***********************************************************
% %
% % Feel free to modify this code.
% %

x=r.*phi.*cos(theta);
y=r.*phi.*sin(theta);
z=r.*cos(phi);


text:
%% Example of using KST class for interfacing with KUKA iiwa robots

% Using the gamepad for controlling the joints positions of KUKA
% iiwa 7 R 800 robot.
% This script implements a graphical user interface, which gives a feedback
% about the position of the joints' angles.

% This example works with Sunrise application version KST_1.7  and higher.
% Copyright: Mohammad SAFEEA, 29th of June 2018

close all;clear;clc;
warning('off')
% Start the joystick, make sure that the joystick is connected
instrreset;
ID=1;
joy=vrjoystick(ID);
%% Add path of KST class to work space
cDir = pwd;
cDir=getTheKSTDirectory(cDir);
addpath(cDir);

%% Instantiate the KST object
ip='172.31.1.147'; % The IP of the controller
arg1=KST.LBR7R800; % choose the robot iiwa7R800 or iiwa14R820
arg2=KST.Medien_Flansch_elektrisch; % choose the type of flange
Tef_flange=eye(4); % transofrm matrix of EEF with respect to flange
iiwa=KST(ip,arg1,arg2,Tef_flange); % create the object

%% Start a connection with the server
flag=iiwa.net_establishConnection();
if flag==0
	fprintf('Can not connect to KST \n');
    fprintf('Program terminated \n');
    return;
end

%% Move robot to some initial position
jPos={0,0,0,0,0,0,0};
relVel=0.15;
iiwa.movePTPJointSpace(jPos, relVel);

%% Draw the interface
[figureHandle,anglesTextHandlesCellArray,labelTextHandlesCellArray]=constructInterface();

%% max angular velocity
w=10*pi/180; % rad/sec

%% Joint space control
firstExecution=0;
jointIndex=1; % the index of the selcted joint for realtime control.

setBackGroundColor(labelTextHandlesCellArray, 1 );
pause(0.1);
frameUpdateCount=0;
maxDisp=1*pi/180;
while true
    joyStatus=read(joy);
    %% Remove the bias in the analog signal
    analogPrecission=1/50;
    for i=1:4
        if abs( joyStatus(i))<analogPrecission
            joyStatus(i)=0;
        end
    end

    % About the variable (joyStatus)
    % joyStatus: is a 4x1 vector, the first and the second elements of this
    % vector correspond to the analog values of the left analog-stick
    % psoition. the third and the fourth elements of this
    % vector correspond to the analog values of the right analog-stick
    % psoition.
    %
    % Only the seconds and the third elements of this vector are used for
    % controlling the joints motion of the robot one at a time.
    %
    % The 2d element is used to select the robot axes, this element corresponds
    % to the (up,down) directions of the left analog stick of the gamepad.
    %
    % The 3rd element is used to control the motion of the selcted axes,
    % this element corresponds to the (right, left) directions of the right
    % analog stick of the gamepad.

    %% Calculate the elapsed time
    if(firstExecution==0)
        firstExecution=1;
        a=datevec(now);
        timeNow=a(6)+a(5)*60+a(4)*60*60; % calculate time at this instant
        dt=0; % elapsed is zero at first excution
        time0=timeNow;
    else
        a=datevec(now);
        timeNow=a(6)+a(5)*60+a(4)*60*60; % calculate time at this instant
        dt=timeNow-time0; % calculate elapsed time interval
        time0=timeNow;
    end

%% Change index  of controlled axes,
% This happens when the left joiystic of the gamepad is moved in up/down
% direction.
%
    if(abs(joyStatus(2))>0.9)
        step=sign(joyStatus(2))*1;
        jointIndex=jointIndex+step;

        if(jointIndex==0)
            jointIndex=7;
        end

        if(jointIndex==8)
            jointIndex=1;
        end

        setBackGroundColor(labelTextHandlesCellArray, jointIndex );
        pause(0.5);
        continue;
    end

 %% Control axes motion
     if(joyStatus(2)==0)
        % move the axes only if the analog signal corresponding to change
        % axes command is zero.
        dq=w*dt*joyStatus(3);
        %% Check for joint limit avoidance
        dq=dq*jointLimitAvoidance(jointIndex,jPos,dq);
        if dq>maxDisp
            dq=maxDisp;
        end
        jPos{jointIndex}=jPos{jointIndex}+dq;
        frameUpdateCount=frameUpdateCount+1;
        if(frameUpdateCount==50)
            setJointAngleToText(anglesTextHandlesCellArray, jointIndex,jPos{jointIndex} )
            frameUpdateCount=0;
            pause(0.001);
        end
        iiwa.sendJointsPositions(jPos);
     end

     %% A condition to break the loop, when x button of the gamepad is pressed
     b=button(joy,1);
     if(b==1)
         break;
     end

end
close(figureHandle);
%% turn off the server
iiwa.net_turnOffServer();

repo: KST-Kuka-Sunrise-Toolbox


text:
function grad = B_MVDR_spatialCov(X, curr_layer, beamform_layer, after_power_layer)

% X is the multichannel complex spectrum inputs
[D,C,T,N] = size(X);
% weight is the beamforming weight
weight = reshape(curr_layer.a, D,C,N);
lambda = curr_layer.lambda;
phi_s = curr_layer.phi_s;
phi_n = curr_layer.phi_n;

if isfield(curr_layer, 'noiseCovL2')
    noiseCovL2 = curr_layer.noiseCovL2;
else
    noiseCovL2 = 0;  % add noiseCovRegularization*\lambda*I to noise covariance, where \lambda is the maximum eigenvalue
end

% Y is the beamforming's output
Y = beamform_layer.a;
% future_grad is the gradient of the power specturm of beamformed signal
future_grad = after_power_layer.grad;

if N>1
    [validMask, variableLength] = getValidFrameMask(curr_layer);
    future_gradUtt = ExtractVariableLengthTrajectory(future_grad, validMask);
    YUtt = ExtractVariableLengthTrajectory(Y, validMask);
    X = reshape(X, D*C, T, N);
    XUtt = ExtractVariableLengthTrajectory(X, validMask);
    for i=1:N
        grad{i} = GetGradUtt(reshape(XUtt{i},D,C,size(XUtt{i},2)),YUtt{i},phi_s(:,:,:,:,i), phi_n(:,:,:,:,i), ...
            weight(:,:,i), lambda(:,:,:,:,i), future_gradUtt{i}, noiseCovL2);
    end
    grad = cell2mat_gpu(grad);
    grad = permute(grad, [1 3 2]);
else
    grad = GetGradUtt(X,Y,phi_s, phi_n, weight, lambda, future_grad, noiseCovL2);
end

end

%%
function grad = GetGradUtt(X,Y,phi_s, phi_n, weight, lambda, future_grad, noiseCovL2)
[D,C,T,N] = size(X);
u = zeros(C,1);
u(1) = 1;
for f=1:D
    x = squeeze(X(f,:,:));
    xx = x*x';
    ww = weight(f,:).' * conj(weight(f,:));
    phi_n_inv = inv(phi_n{1,1,f});
    yy = abs(Y(f,:).*conj(Y(f,:)));
    dyyy = sum(future_grad(f,:).*yy);
    dyxx = bsxfun(@times, x, future_grad(f,:)) * x';

    if 1
        grad_phi_s(:,:,f) = 2 * phi_n_inv * dyxx * weight(f,:).' * u' / lambda{1,1,f} -2 * dyyy * phi_n_inv / lambda{1,1,f};
        grad_phi_n(:,:,f) = - 2* phi_n_inv * dyxx * ww + 2 * dyyy * phi_n_inv * phi_s{1,1,f} * phi_n_inv / lambda{1,1,f};
    elseif 0    % no lambda
        grad_phi_s(:,:,f) = 2 * phi_n_inv * dyxx * weight(f,:).' * u';
        grad_phi_n(:,:,f) = - 2* phi_n_inv * dyxx * ww;
    else    % only lambda
        grad_phi_s(:,:,f) = -2 * dyyy * phi_n_inv / lambda{1,1,f};
        grad_phi_n(:,:,f) = 2 * dyyy * phi_n_inv * phi_s{1,1,f} * phi_n_inv / lambda{1,1,f};
    end
end

grad = conj([reshape(grad_phi_s,C*C*D,1); reshape(grad_phi_n,C*C*D,1)]);
end


text:
function tests = PGraphTest
  tests = functiontests(localfunctions);
end
function PGraph2_test(testCase)

    stream = RandStream.getGlobalStream;
    stream.reset()

    g = PGraph(2);

    % add first component, nodes 1-5
    g.add_node( rand(2,1));
    g.add_node( rand(2,1));
    g.add_node( rand(2,1));
    g.add_node( rand(2,1));
    g.add_node( rand(2,1));

    g.add_edge(1,2);
    g.add_edge(1,3);
    g.add_edge(1,4);
    g.add_edge(2,3);
    g.add_edge(2,4);
    g.add_edge(4,5);

    verifyEqual(testCase, g.n, 5);
    verifyEqual(testCase, g.ne, 6);
    verifyEqual(testCase, g.nc, 1);
    verifyEqual(testCase, g.neighbours(2), [3 4 1]);
    z = g.coord(1);
    %verifyEqual(testCase,  g.distance(1,2), 0.6878, 'absTol',1e-4)
%     verifyEqual(testCase,  g.distances([0 0]), ...
%         [0.6137    0.6398    0.9222    1.2183    1.3593], 'absTol',1e-4)

    D = g.degree();
    I = g.incidence();
    A = g.adjacency();
    L = g.laplacian();

    verifyEqual(testCase, D, diag([3 3 2 3 1]));
    verifyEqual(testCase, sum(I), [2 2 2 2 2 2]);
    verifyEqual(testCase, sum(I'), [3 3 2 3 1]);
    verifyEqual(testCase, max(max(A-A')), 0);

    verifyEqual(testCase,  g.edgedir(1,3), 1);
    verifyEqual(testCase,  g.edgedir(3,1), -1);
    verifyEqual(testCase,  g.edgedir(2,5), 0);
    verifyEqual(testCase, g.neighbours_d(2), [-1 3 4]);

    s = g.char();

    % test node data
    g.setvdata(1, [1 2 3]);
    d = g.vdata(1);
    verifyEqual(testCase, d, [1 2 3]);

    % test edge data
    g.setedata(1, [4 5 6 7]);
    d = g.edata(1);
    verifyEqual(testCase, d, [4 5 6 7]);

    % add second component, nodes 6-8

    g.add_node( rand(2,1));
    g.add_node( rand(2,1), 6, 1);
    g.add_node( rand(2,1), 7);
    g.add_edge(6,7);
    g.add_edge(6,8);

    verifyEqual(testCase, g.nc, 2);

    clf
    g.plot('labels')
    g.plot('edgelabels')
    g.plot('componentcolor')

    [path,cost]=g.Astar(3,5);
    g.highlight_path(path);
    %verifyEqual(testCase, path, [3 2 4 5]);
    %verifyEqual(testCase,  cost, 2.1536, 'absTol',1e-4)

    %verifyEqual(testCase, g.closest([0 0]), 4);

    g.setcost(1, 99);
    verifyEqual(testCase,  g.cost(1), 99);

    g.Astar(3, 5);
    %verifyEqual(testCase, path, [3 2 4 5]);
end

function PGraph3_test(testCase)

    stream = RandStream.getGlobalStream;
    stream.reset()

    g = PGraph(3);

    g.add_node( rand(3,1));
    g.add_node( rand(3,1));
    g.add_node( rand(3,1));
    g.add_node( rand(3,1));
    g.add_node( rand(3,1));

    g.add_edge(1,2);
    g.add_edge(1,3);
    g.add_edge(1,4);
    g.add_edge(2,3);
    g.add_edge(2,4);
    g.add_edge(4,5);

    clf
    g.plot('labels')
    g.plot('edgelabels')
end

author: petercorke

repo: spatialmath-matlab


text:
function [qualMeas]=Measure_Quality(res_prev,res,QualMeasOpts)
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Manasavee Lohvithee
%--------------------------------------------------------------------------

%for loop over parameters to get the name
for ii=1:length(QualMeasOpts)
    opt=QualMeasOpts{ii};

    switch opt
        case 'RMSE'
         q=RMSE(res_prev,res);

        case 'CC'
         q=CC(res_prev,res);

        case 'MSSIM'
         q=MSSIM(res_prev,res);

        case 'UQI'
         q=UQI(res_prev,res);
        case 'error_norm'
         q=im3Dnorm(res_prev-res,'L2');

    end

    qualMeas(ii)=q;
end

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE