text:
function header = text2header(requestdata,config)
request_text=char(requestdata);
request_lines = regexp(request_text, '\r\n+', 'split');
request_words = regexp(request_lines, '\s+', 'split');
for i=1:length(request_lines)
    line=request_lines{i};
    if(isempty(line)), break; end
    type=request_words{i}{1};
    switch(lower(type))
        case 'get'
            header.Get.Filename=request_words{i}{2};
            header.Get.Protocol=request_words{i}{3};
        case 'post'
            header.Post.Filename=request_words{i}{2};
            header.Post.Protocol=request_words{i}{3};
        case 'host:'
            header.Host=rmvp(line(7:end));
        case 'user-agent:'
            header.UserAgent=rmvp(line(13:end));
        case 'accept:'
            header.Accept=rmvp(line(9:end));
        case 'accept-language:'
            header.AcceptLanguage=rmvp(line(18:end));
        case 'accept-encoding:'
            header.AcceptEncoding=rmvp(line(18:end));
        case 'accept-charset:'
            header.AcceptCharset=rmvp(line(17:end));
        case 'keep-alive:'
            header.KeepAlive=rmvp(line(13:end));
        case 'connection:'
            header.Connection=rmvp(line(13:end));
        case 'content-length:'
            header.ContentLength=rmvp(line(17:end));
        case 'content-type:'
            %lines=rmvp(line(15:end));
            switch rmvp(request_words{i}{2})
                case {'application/x-www-form-urlencoded','application/x-www-form-urlencoded;'}
                    header.ContentType.Type='application/x-www-form-urlencoded';
                    header.ContentType.Boundary='&';
                case {'multipart/form-data','multipart/form-data;'}
                    header.ContentType.Type='multipart/form-data';
                    str=request_words{i}{3};
                    header.ContentType.Boundary=str(10:end);
                otherwise
                    disp('unknown')
            end
        otherwise
    end
end
header.Content=struct;
if(isfield(header,'ContentLength'))
    cl=str2double(header.ContentLength);
    str=request_text(end-cl+1:end);
    data=requestdata(end-cl+1:end);
    if(~isfield(header,'ContentType'))
        header.ContentType.Type=''; header.ContentType.Boundary='&';
    end
    switch (header.ContentType.Type)
        case {'application/x-www-form-urlencoded',''}
            str=rmvp(str);
            words = regexp(str, '&', 'split');
            for i=1:length(words)
                words2 = regexp(words{i}, '=', 'split');
                header.Content.(words2{1})=words2{2};
            end
        case 'multipart/form-data'
            pos=strfind(str,header.ContentType.Boundary);
            while((pos(1)>1)&&(str(pos(1)-1)=='-'))
                header.ContentType.Boundary=['-' header.ContentType.Boundary];
                pos=strfind(str,header.ContentType.Boundary);
            end

            for i=1:(length(pos)-1)
                pstart=pos(i)+length(header.ContentType.Boundary);
                pend=pos(i+1)-3; % Remove "13 10" End-line characters
                subrequestdata=data(pstart:pend);
                subdata= multipart2struct(subrequestdata,config);
                header.Content.(subdata.Name).Filename=subdata.Filename;
                header.Content.(subdata.Name).ContentType=subdata.ContentType;
                header.Content.(subdata.Name).ContentData=subdata.ContentData;
            end
        otherwise
            disp('unknown')
    end
end

author: jaderberg

repo: python-matlab-bridge


text:
function X = positive_definite_karcher_mean(A)
% Computes a Karcher mean of a collection of positive definite matrices.
%
% function X = positive_definite_karcher_mean(A)
%
% Input:  A 3D matrix A of size nxnxm such that each slice A(:,:,k) is a
%         positive definite matrix of size nxn.
%
% Output: A positive definite matrix X of size nxn which is a Karcher mean
%         of the m matrices in A, that is, X minimizes the sum of squared
%         Riemannian distances to the matrices in A:
%            f(X) = sum_k=1^m .5*dist^2(X, A(:, :, k))
%         The distance is defined by the natural metric on the set of
%         positive definite matrices: dist(X,Y) = norm(logm(X\Y), 'fro').
%
% This simple example is not the best way to compute Karcher means. Its
% purpose it to serve as base code to explore other algorithms. In
% particular, in the presence of large noise, this algorithm seems to not
% be able to reach points with a very small gradient norm. This may be
% caused by insufficient accuracy in the gradient computation.

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

    % Generate some random data to test the function if none is given.
    if ~exist('A', 'var') || isempty(A)
        n = 5;
        m = 10;
        A = zeros(n, n, m);
        ref = diag(max(.1, 1+.1*randn(n, 1)));
        for i = 1 : m
            noise = 0.01*randn(n);
            noise = (noise + noise')/2;
            [V D] = eig(ref + noise);
            A(:, :, i) = V*diag(max(.01, diag(D)))*V';
        end
    end

    % Retrieve the size of the problem:
    % There are m matrices of size nxn to average.
    n = size(A, 1);
    m = size(A, 3);
    assert(n == size(A, 2), ...
           ['The slices of A must be square, i.e., the ' ...
	        'first and second dimensions of A must be equal.']);

    % Our search space is the set of positive definite matrices of size n.
    % Notice that this is the only place we specify on which manifold we
    % wish to compute Karcher means. Replacing this factory for another
    % geometry will yield code to compute Karcher means on that other
    % manifold, provided that manifold is equipped with a dist function and
    % a logarithmic map log.
    M = sympositivedefinitefactory(n);

    % Define a problem structure, specifying the manifold M, the cost
    % function and its gradient.
    problem.M = M;
    problem.cost = @cost;
    problem.grad = @grad;

    % The functions below make many redundant computations. This
    % performance hit can be alleviated by using the caching system. We go
    % for a simple implementation here, as a tutorial example.

    % Cost function
    function f = cost(X)
        f = 0;
        for k = 1 : m
            f = f + M.dist(X, A(:, :, k))^2;
        end
        f = f/(2*m);
    end

    % Riemannian gradient of the cost function
    function g = grad(X)
        g = M.zerovec(X);
        for k = 1 : m
            % Update g in a linear combination of the form
            % g = g - [something]/m.
            g = M.lincomb(X, 1, g, -1/m, M.log(X, A(:, :, k)));
        end
    end

    % Execute some checks on the derivatives for early debugging.
    % These things can be commented out of course.
    % The slopes should agree on part of the plot at least. In this case,
    % it is sometimes necessary to inspect the plot visually to make the
    % call, but it is indeed correct.
    % checkgradient(problem);
    % pause;

    % Execute this if you want to force using a proper parallel vector
    % transport. This is not necessary. If you omit this, the default
    % vector transport is the identity map, which is (of course) cheaper
    % and seems to perform well in practice.
    % M.transp = M.paralleltransp;

    % Issue a call to a solver. Default options are selected.
    % Our initial guess is the first data point.
    X = trustregions(problem, A(:, :, 1));

end

author: andrewssobral

repo: lrslibrary


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: midpoint quadrature rule for 2D spline
%
%==============================================================================

clear, close all, help(mfilename)

omega = [0,6,0,8]; I = 36; T = zeros(6,8); T(3,4) = 1; h = []; Q = [];
psi = @(xc) splineInter(T,omega,xc);
for j=1:10,
  m    = 2^j*[1,1];
  h(j) = prod((omega(2:2:end)-omega(1:2:end))./m);
  xc   = getCellCenteredGrid(omega,m);
  Q(j) = h(j)*sum(psi(xc));
end;
figure(1); clf; p1=semilogx(h/h(1),Q+eps,'kx',h/h(1),Q,'k-');
figure(2); clf; p2=loglog(h/h(1),abs(I-Q)+eps,'kx',h/h(1),abs(I-Q),'k-');
%==============================================================================


text:
function [M11,M12,M21,M22,f1,f2] = autoGen_acrobotDynamics(q1,q2,dq1,dq2,u,m1,m2,g,l1,l2)
%AUTOGEN_ACROBOTDYNAMICS
%    [M11,M12,M21,M22,F1,F2] = AUTOGEN_ACROBOTDYNAMICS(Q1,Q2,DQ1,DQ2,U,M1,M2,G,L1,L2)

%    This function was generated by the Symbolic Math Toolbox version 6.2.
%    11-Jul-2015 20:41:44

t2 = cos(q1);
t3 = l1.^2;
t4 = sin(q1);
t5 = cos(q2);
t6 = l1.*t2;
t7 = l2.*t5;
t8 = t6+t7;
t9 = sin(q2);
t10 = l1.*t4;
t11 = l2.*t9;
t12 = t10+t11;
M11 = -m1.*t2.^2.*t3-m1.*t3.*t4.^2-l1.*m2.*t2.*t8-l1.*m2.*t4.*t12;
if nargout > 1
    M12 = -l2.*m2.*t5.*t8-l2.*m2.*t9.*t12;
end
if nargout > 2
    M21 = -l1.*l2.*m2.*t2.*t5-l1.*l2.*m2.*t4.*t9;
end
if nargout > 3
    t13 = l2.^2;
    M22 = -m2.*t5.^2.*t13-m2.*t9.^2.*t13;
end
if nargout > 4
    t14 = dq1.^2;
    t15 = dq2.^2;
    t16 = l1.*t2.*t14;
    t17 = l2.*t5.*t15;
    t18 = t16+t17;
    t19 = l1.*t4.*t14;
    t20 = l2.*t9.*t15;
    t21 = t19+t20;
    t22 = m2.*t12.*t18;
    t23 = g.*m2.*t12;
    t24 = g.*l1.*m1.*t4;
    f1 = t22+t23+t24-m2.*t8.*t21;
end
if nargout > 5
    t25 = l2.*m2.*t9.*t18;
    t26 = g.*l2.*m2.*t9;
    f2 = t25+t26-u-l2.*m2.*t5.*t21;
end

author: MatthewPeterKelly

repo: OptimTraj


text:
% [VEC] = columnize(MTX)
%
% Pack elements of MTX into a column vector.  Just provides a
% function-call notatoin for the operation MTX(:)

function vec = columnize(mtx)

vec = mtx(:);

repo: SelfExSR


text:
function pass = test_real(pref)
% Test REAL

if ( nargin == 0 )
    pref = chebfunpref;
end
tol = 100*pref.cheb3Prefs.chebfun3eps;

f = chebfun3v(@(x,y,z) cos(x.*y.*z), @(x,y,z) cos(x.*y.*z));
g = real(f);
pass(1) = norm(g - f) < tol;

f = chebfun3v(@(x,y,z) cos(x.*y.*z), @(x,y,z) cos(x.*y.*z));
g = real(1i*f);
pass(2) = norm(g) < tol;

f1 = chebfun3v(@(x,y,z) cos(x.*y.*z), @(x,y,z) cos(x.*y.*z));
f2 = chebfun3v(@(x,y,z) sin(x + y.^2+z.^3), @(x,y,z) sin(x + y.^2+z.^3));
g = real(f1 + 1i*f2);
pass(3) = norm(g - f1) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [nCoords, nDim, nVTMode, c] = plx_vt_interpret(ts, sv);
% plx_vt_interpret - interpret CinePlex video tracking data
%
% [nCoords, nDim, nVTMode, c] = plx_vt_interpret(ts, sv);
%
% INPUT:
%   ts - array of timestamps (in seconds) (see plx_event_ts.m)
%   sv - array of strobed event values (see plx_event_ts.m)
%
% OUTPUT:
%   nCoords - number of produced coordinates
%   nDim    - number of elemnts in produced coordinates
%             nDim = 3 for CENTROID, LED_1, LED_2, LED3
%             nDim = 4 for CENTROID_WITH_MOTION
%             nDim = 5 for LED_12, LED_23, LED_13
%             nDim = 7 for LED_123
%   nVTMode - VT mode:
%              0 = UNKNOWN
%              1 = CENTROID                // 1 set of coordinates, no motion
%              2 = CENTROID_WITH_MOTION    // 1 set of coordinates, with motion
%              3 = LED_1                   // 1 set of coordinates
%              4 = LED_2
%              5 = LED_3
%              6 = LED_12                  // 2 sets of coordinates
%              7 = LED_13
%              8 = LED_23
%              9 = LED_123                 // 3 sets of coordinates
%   c       - nCoords by nDim matrix of produced coordinates
%             c(:, 1) - timestamp
%             c(:, 2) - x1
%             c(:, 3) - y1
%             c(:, 4) - x2 or motion (if present)
%             c(:, 5) - y2 (if present)
%             c(:, 6) - x3 (if present)
%             c(:, 7) - y3 (if present)
%

if nargin < 2
    error 'Expected 2 input arguments';
end

[nCoords, nDim, nVTMode, c] = mexPlex(21, '', ts, sv);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ y, m, d, f ] = jed_to_ymdf_julian ( jed )

%*****************************************************************************80
%
%% JED_TO_YMDF_JULIAN converts a JED to a Julian YMDF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Richards,
%    Algorithm F,
%    Mapping Time, The Calendar and Its History,
%    Oxford, 1999, pages 324-325.
%
%  Parameters:
%
%    Input, real JED, the Julian Ephemeris Date.
%
%    Output, integer Y, integer M, integer D, real F,
%    the YMDF date.
%

%
%  Determine the computational date (Y'/M'/D').
%
  j = floor ( jed + 0.5 );
  f = ( jed + 0.5 ) - j;

  j_prime = j + 1401;

  y_prime = floor (     ( 4 * j_prime + 3 ) / 1461 );
  t_prime = floor ( mod ( 4 * j_prime + 3, 1461 ) / 4 );
  m_prime = floor (     ( 5 * t_prime + 2 ) / 153 );
  d_prime = floor ( mod ( 5 * t_prime + 2, 153 ) / 5 );
%
%  Convert the computational date to a calendar date.
%
  d = d_prime + 1;
  m = mod ( m_prime + 2, 12 ) + 1;
  y = y_prime - 4716 + floor ( ( 14 - m ) / 12 );
%
%  Any year before 1 AD must be moved one year further back, since
%  this calendar does not include a year 0.
%
  y = y_astronomical_to_common ( y );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% This file should be run before running the HOWTO scripts
function startup

current_directory = fileparts(mfilename('fullpath'));
entries = strsplit(genpath(current_directory), pathsep);

% helper function to find matching entries
entry_matcher = @(entries, pattern) ~cellfun(@isempty, strfind(entries, pattern));

% Remove all .git and .svn entries from the path
entries(entry_matcher(entries, '.git') | entry_matcher(entries, '.svn')) = [];

% as well as the documentation folder
entries(entry_matcher(entries, strcat(filesep, 'documentation'))) = [];

if ~exist('OCTAVE_VERSION', 'builtin')
    % Matlab is running
    % Thus remove also the octave helper functions
    octave_path = strcat(filesep, 'external', filesep, 'octave_only');
    entries(entry_matcher(entries, octave_path)) = [];
else
    % Octave is running
    more off;
    pkg load tsa;
    pkg load signal;
end

% remove backward compatible functions
MATLABVERSION=version('-release');
MATLABVERSION=MATLABVERSION(1:end-1);
if str2double(MATLABVERSION)>=2015
    entries(entry_matcher(entries, 'backcompatibility_2015')) = [];
end

% addpaths
addpath(strjoin(entries, pathsep));
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
clear all
close all
clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demo file for deconvtv
% Image 'salt and pepper' noise removal
%
% Stanley Chan
% University of California, San Diego
% 20 Jan, 2011
%
% Copyright 2011
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Prepare images
f_orig  = im2double(imread('C:\Users\Stanley Chan\Dropbox\TestImages\Lena.bmp'));
[rows cols frames] = size(f_orig);
H       = fspecial('gaussian', [9 9], 2);
g       = imfilter(f_orig, H, 'circular');
g       = imnoise(g, 'salt & pepper', 0.05);

% Setup parameters (for example)
opts.rho_r   = 5;
opts.rho_o   = 100;
opts.beta    = [1 1 0];
opts.print   = true;
opts.alpha   = 0.7;
opts.method  = 'l1';

% Setup mu
mu           = 20;

% Main routine
tic
out = deconvtv(g, H, mu, opts);
toc

% Display results
figure(1);
imshow(g);
title('input');

figure(2);
imshow(out.f);
title('output');

author: Sable

repo: mcbench-benchmarks


text:
classdef ConstantHeadingModelX <  DynamicModelX
% ConstantHeadingModelX class
%
% Summary of ConstantHeadingModelX
% This is a class implementation of a time-varying 2D Linear-Gaussian
% Constant Heading Dynamic Model [1].
%
% The model is described by the following SDEs:
%   dx = v*cos(h)*dt                      | Position on X-axis (m)
%   dy = v*sin(h)*dt                      | Position on Y axis (m)
%   dv = q_vel*dW_t,  W_t~N(0,t)          | Absolute velocity  (m/s)
%   dh = q_head*dB_t, B_t~N(0,t)          | Heading            (rad)
%
% Or equivalently:
%   x(t) = f(x(t-1),Dt) + w(t),  w(t)~ N(0,Q)
%
% where: (using Euler discretisation of the SDEs)
%   x = [x; y; v; h]
%   F = [x + v*cos(h)*Dt; y+ v*sin(h)*Dt; v; h]
%   Q = [0, 0, 0, 0;
%        0, 0, 0, 0;
%        0, 0, q_vel*sqrt(Dt), 0;
%        0, 0, 0, q_head*sqrt(Dt)];
%
% ConstantVelocityModelX_2D Properties:
%   - VelocityErrVariance  Value of the velocity noise diffusion coefficient q_vel
%   - HeadingErrVariance   Value of the heading noise diffusion coefficient q_head
%   - TimeVariant          Value of the time variant Dt
%
% ConstantVelocityModelX_2D Methods:
%   - feval(~)         Equivalent to applying the model transition equations
%   - random(~)        Process noise sample generator function
%   - pdf(~)           Function to evaluate the probability p(x_k|x_{k-1}) of
%                       a set of new states, given a set of (particle) state vectors
%                       e.g. eval = @(xk,xkm1) mvnpdf(xk,xkm1,Q);
%   - covariance(~)    Returns the state covariance process Q_k
%
%  February 2018 Lyudmil Vladimirov, University of Liverpool

    properties
        VelocityErrVariance
        HeadingErrVariance
        TimeVariant
    end

    properties (Access = private)
        f = @(xkm1,Dt) [xkm1(1,:)+Dt*xkm1(3,:).*cos(xkm1(4,:));
                        xkm1(2,:)+Dt*xkm1(3,:).*sin(xkm1(4,:));
                        xkm1(3,:);
                        xkm1(4,:)];
        Q = @(Dt, q_vel, q_head) Dt*blkdiag(q_vel^2, q_vel^2, q_vel, q_head);
    end

    methods
        function this = ConstantHeadingModelX(varargin)
        % CONSTANTHEADINGMODELX Constructor method
        %
        % DESCRIPTION:
        % * ConstantHeadingModelX(q_vel,q_head) instantiates aan object handle
        %   configured with the provided velocity and heading process noise
        %   diffusion coefficients q_vel and q_head (scalar).
        % * ConstantVelocityModelX_2D(config) instantiates an object handle
        %   configured with the provided velocity and heading process noise
        %   diffusion coefficients config.VelocityErrVariance and
        %   config.HeadingErrVariance (scalar).
        % * ConstantVelocityModelX_2D(___,Name,Value) instantiates an object
        %   handle, configured with additional options specified by one or
        %   more Name,Value pair arguments.
        %
        % PARAMETERS
        % * VelocityErrVariance - (Required) The VelocityErrVariance is a scalar
        %   value describing the velocity noise diffusion coefficient.
        % * HeadingErrVariance - (Required) The VelocityErrVariance is a scalar
        %   value describing the heading noise diffusion coefficient (in rads/sec).
        %
        %  See also apply, rnd, pdf, covariance.

            % Call SuperClass method
            this@DynamicModelX;

            % Return quickly if no arguments are passed
            if(nargin==0)
                this.TimeVariant = 1;
                this.VelocityErrVariance = 1;
                this.HeadingErrVariance = 1;
            end

            % First check to see if a structure was received
            if(nargin==1)
                if(isstruct(varargin{1}))
                    this.VelocityErrVariance = varargin{1}.VelocityErrVariance;
                    this.HeadingErrVariance = varargin{1}.HeadingErrVariance;
                    this.TimeVariant = varargin{1}.TimeVariant;
                    return;
                end
            end

            % Otherwise, fall back to input parser
            parser = inputParser;
            parser.KeepUnmatched = true;
            parser.addOptional('VelocityErrVariance',[]);
            parser.addOptional('HeadingErrVariance',[]);
            parser.addOptional('TimeVariant',1);
            parser.parse(varargin{:});

            this.VelocityErrVariance = parser.Results.VelocityErrVariance;
            this.HeadingErrVariance = parser.Results.HeadingErrVariance;
            this.TimeVariant = parser.Results.TimeVariant;

            this.NumStateDims = 4;

        end

        function xk = feval(this, xkm1, wk, Dt)
        % FEVAL Propagate a given state through the dynamic model
        %
        % Parameters
        % ----------
        % xkm1: (NumStateDims x Ns) matrix, optional
        %   - A matrix, whose columns correspond to individual state vectors.
        %   - If not provided, then the state transition matrix (F) will be
        %     returned.
        % wk: (NumStateDims x Ns) matrix or boolean, optional
        %   - If wk is a boolean and set True, then the default model noise
        %     generation function (this.random()) will be used to generate the
        %     noise.
        %   - Otherwise, wk should be a matrix of equal dimensions to xk,
        %     where each column corresponds to the random noise which will
        %     be added to each state vector.
        %   - If not provided, then no noise will be added to the state vectors.
        % Dt: scalar, optional
        %   A time variant. (default=this.TimeVariant)
        %
        % Returns
        % -------
        % xk: (NumStateDims x Ns) matrix or function handle
        %   - If no parameters are passed to the function, then xk will be the
        %     state transition function f.
        %   - Else, xk will be a (NumStateDims x Ns) matrix, whose columns
        %     correspond to the columns of xkm1, each propagated through the
        %     dynamic model
        %
        % Usage
        % -----
        % * xk = FEVAL(this,xkm1,wk,Dt) returns the new state xk, produced by
        %   propagating the given state xkm1 through the dynamic model for
        %   time Dt and adding random noise wk. xk, xkm1 and wk must have
        %   the same dimensions, i.e (NumStateDim x Ns), where Ns is the
        %   number of states/columns in xkm1.
        % * xk = FEVAL(this,xkm1,wk) returns the new state xk, produced by
        %   propagating the given state xkm1 through the dynamic model for
        %   time this.TimeVariant and adding random noise wk.
        %   (i.e. Default Dt = this.TimeVariant)
        % * xk = FEVAL(this,xkm1) returns the new state xk, produced by
        %   propagating the given state xkm1 through the dynamic model for
        %   time this.TimeVariant without the addition of noise.
        %   (i.e. Default wk = 0)
        % * xk = FEVAL(this) returns the model's transition matrix resulting
        %   through the application of this.TimeVariant.
        %   (i.e. Default xkm1 = 1)
        %
        % See also PDF, RANDOM, COVARIANCE.

            switch(nargin)
                case 1
                    xk = this.f;
                    return;
%                     xkm1 = 1;
%                     wk   = 0;
%                     Dt = this.TimeVariant;
                case 2
                    wk   = 0;
                    Dt = this.TimeVariant;
                case 3
                    Dt = this.TimeVariant;
                    if(islogical(wk) && wk)
                        wk = this.random(size(xkm1,2),Dt);
                    end
            end

            % Compute result
            xk = this.f(xkm1,Dt) + wk;
        end

        function cov = covariance(this,Dt)
        % COVARIANCE Returns process covariance matrix.
        %
        % Parameters
        % ----------
        % Dt: scalar, optional
        %   A time variant. (default=this.TimeVariant)
        %
        % Returns
        % -------
        % cov: (NumStateDims x NumStateDims) matrix
        %   The process noise covariance matrix
        %
        % Usage
        % -----
        % * Qk = covariance(this,Dt) returns process covariance matrix, upon
        %   application of Dt.
        % * Qk = covariance(this) returns process covariance matrix, upon
        %   application of this.TimeVariant.
        % (i.e. Default Dt = this.TimeVariant)
        %
        % See also FEVAL, RANDOM, PDF.
            switch(nargin)
                case 1
                    Dt = this.TimeVariant;
            end

            % Return process covariance
            cov = this.Q(Dt,this.VelocityErrVariance, this.HeadingErrVariance); % Time variant
        end

        function wk = random(this, Ns, Dt)
        % RANDOM Generates random samples from the dynamic model's noise
        % distribution.
        %
        % Parameters
        % ----------
        % Ns: scalar, optional
        %   The number of samples to be generated.
        %   (default = 1)
        % Dt: scalar, optional
        %   A time variant.
        %   (default=this.TimeVariant)
        %
        % Returns
        % -------
        % wk: (NumStateDims x Ns)
        %   A matrix, whose columns correspond to indivual/independent
        %   noise samples.
        %
        % Usage
        % -----
        % * wk = random(this,Ns,Dt) generates and returns a set of Ns samples
        %   generated from the noise distribution of the dynamic model, i.e.
        %   wk ~ N(0,Q(Dt)), where Q is the noise covariance upon application
        %   of the time variant Dt.
        % * wk = random(this,Ns) generates and returns a set of Ns samples
        %   generated from the noise distribution of the dynamic model, i.e.
        %   wk ~ N(0,Q), where Q is the noise covariance upon application
        %   of the time variant this.TimeIndex.
        % (i.e. Default Dt = this.TimeVariant)
        %
        % See also FEVAL, PDF, COVARIANCE

            switch(nargin)
                case 1
                    Ns = 1;
                    Dt = this.TimeVariant;
                case 2
                    Dt = this.TimeVariant;
            end

            wk = mvnrnd(zeros(this.NumStateDims,1),this.Q(Dt,this.VelocityErrVariance, this.HeadingErrVariance),Ns)';
        end

        function prob = pdf(this, xk, xkm1, Dt)
        % PDF Evaluates the probability/likelihood p(x_k|x_{k-1}) of a
        % (set of) new state vector(s), given a (set of) old state vector(s)
        %
        % Parameters
        % ----------
        % xk: (NumStateDims x Ns) matrix
        %   A matrix, whose columns correspond to individual new state vectors.
        % xkm1: (NumStateDims x Np) matrix
        %   A matrix, whose columns correspond to individual old state vectors
        % Dt: scalar, optional
        %   A time variant.
        %   (default=this.TimeVariant)
        %
        % Returns
        % -------
        % prob: (Np x Ns) matrix
        %   A matrix, where each element (i,j) corresponds to the evaluated
        %   probability p(xk(:,j)|xkm1(:,i))
        %
        % Usage
        % -----
        % * prob = pdf(x_k, x_km1, Dt) evaluates and returns a (a x b)
        %   probability/likelihood matrix given the (NumStatesDim x a) x_k
        %   and (NumStatesDim x b) x_km1 state matrices. Dt is an optional
        %   argument (Default = this.TimeVariant) time variable, which is
        %   used for computing the model's covariance.
        %
        % See also FEVAL, RANDOM, COVARIANCE

            if(nargin<4)
                Dt = this.TimeVariant;
            end

            xk_km1 = this.feval(xkm1);
            prob = zeros(size(xk,2), size(xkm1,2));
            if(size(xkm1,2)>size(xk,2))
                for i=1:size(xk,2)
                    prob(i,:) = gauss_pdf(xk(:,i), xk_km1, this.Q(Dt,this.VelocityErrVariance,this.HeadingErrVariance));
                end
            else
                for i=1:size(xkm1,2)
                    prob(:,i) = gauss_pdf(xk, xk_km1(:,i), this.Q(Dt,this.VelocityErrVariance,this.HeadingErrVariance))';
                end
            end

        end
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function ppma_write ( file_name, r, g, b )

%*****************************************************************************80
%
%% PPMA_WRITE writes RGB data into an ASCII PPM file.
%
%  The PPMA file will have the format:
%
%    P3
%    # file_name created by Matlab PPMA_WRITE.
%    ncol nrow
%    maxrgb, the maximum legal value for the data.
%    row 1, a list of ncol triples between 0 and maxrgb.
%    row 2, a list of ncol triples between 0 and maxrgb.
%    ...
%    row nrow, a list of ncol triples between 0 and maxrgb.
%
%    Lines beginning with '#' are comments.
%
%  Example:
%
%    P3
%    # feep.ppma
%    4 4
%    15
%     0  0  0  0  0  0    0  0  0   15  0 15
%     0  0  0  0 15  7    0  0  0    0  0  0
%     0  0  0  0  0  0    0 15  7    0  0  0
%    15  0 15  0  0  0    0  0  0    0  0  0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILENAME, the name of the PPMA file to create.
%
%    Input, integer R(NROW,NCOL), G(NROW,NCOL), B(NROW,NCOL), the RGB data.
%
  fid = fopen ( file_name, 'wt' );

  if ( fid < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PPMA_WRITE - Fatal error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'PPMA_WRITE - Fatal error!' );
    return;
  end
%
%  First line identifies the file as a PPMA file.
%
  fprintf ( fid, 'P3\n' );
%
%  Second line is a comment with the filename and creator.
%
  fprintf ( fid, '# %s created by PPMA_IO::PPMA_WRITE.\n', file_name );
%
%  Third line holds NCOL and NROW.
%
  [nrow, ncol] = size ( r );
  fprintf ( fid, '%d %d\n', ncol, nrow );
%
%  Fourth line contains the maximum RGB value.
%
  maxr = max ( max ( r ) );
  maxg = max ( max ( g ) );
  maxb = max ( max ( b ) );
  maxrgb = max ( maxr, max ( maxg, maxb ) );
  fprintf ( fid, '%d\n', maxrgb );
%
%  Subsequent lines contain the RGB values, no more than 3 triples to a line.
%
  for i = 1 : nrow
    for j = 1 : ncol
      fprintf ( fid, '%d ', r(i,j), g(i,j), b(i,j) );
      if ( j == ncol || rem ( j, 3 ) == 0 )
        fprintf ( fid, '\n' );
      end
    end
  end

  fclose ( fid );

  return
end


text:
function u = rgbtogrey(v)

% rgbtogrey -- Create greyscale image from RGB image
%
% Usage:
%       u = rgbtogrey(v)
%
% Input:
%       v           Input RGB image
%
% Output:
%       u           Output greyscale image
%
%
% Author: Brendt Wohlberg <brendt@lanl.gov>  Modified: 2014-05-29
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'License' file distributed with
% the library.

if size(v,3) == 1,
  u = v;
elseif size(v,3) == 3,
  C = [0.30 0.59 0.11];
  u = C(1)*v(:,:,1) +  C(2)*v(:,:,2) + C(3)*v(:,:,3);
else
  warning('rgbtogrey: input image must have 1 or 3 bands');
  u = [];
end

return

author: thfylsty


text:
function [degree, node] = grNodeDegree(node, edges)
%GRNODEDEGREE Degree of a node in a (undirected) graph.
%
%   DEGREE = grNodeDegree(NODE_INDEX, EDGES);
%   return the degree of a node in the given edge list, that is the number
%   of edges connected to it.
%   NODE_INDEX is the index of the node, and EDGES is a liste of couples of
%   indices (origin and destination node).
%   This degree is the sum of inner degree (number of edges arriving on the
%   node) and the outer degree (number of emanating edges).
%
%   Note: Also works when NODE_INDEX is a vector of indices
%
%   DEGREE = grNodeDegree(EDGES);
%   Return the degree of each node references by the array EDGES. DEGREE is
%   a column vector with as many rows as the number of nodes referenced by
%   edges.
%
%   [DEG, INDS] = grNodeDegree(EDGES);
%   Also returns the indices of the nodes that were referenced.
%
%   Example
%     edges = [1 2;1 3;2 3;2 4;3 4];
%     grNodeDegree(2, edges)
%     ans =
%          3
%     grNodeDegree(edges)'
%     ans =
%          2     3     3     2
%
%   See also grNodeInnerDegree, grNodeOuterDegree
%

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2003-08-13
% Copyright 2003-2022 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas)

% If only edge array is given, assume we want the degree of each node
if nargin == 1
    edges = node;
    node = unique(edges(:));
end

% allocate array for result
degree = zeros(size(node));

% for each node ID, count the number of edges containing it
for i = 1:length(node(:))
    degree(i) = sum(edges(:,1) == node(i)) + sum(edges(:,2) == node(i));
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function varargout = rotation3dToEulerAngles(mat, varargin)
%ROTATION3DTOEULERANGLES Extract Euler angles from a rotation matrix.
%
%   [PHI, THETA, PSI] = rotation3dToEulerAngles(MAT)
%   Computes Euler angles PHI, THETA and PSI (in degrees) from a 3D 4-by-4
%   or 3-by-3 rotation matrix.
%
%   ANGLES = rotation3dToEulerAngles(MAT)
%   Concatenates results in a single 1-by-3 row vector. This format is used
%   for representing some 3D shapes like ellipsoids.
%
%   ... = rotation3dToEulerAngles(MAT, CONVENTION)
%   CONVENTION specifies the axis rotation sequence.
%   Supported conventions are: 'ZYX', 'ZYZ'. Default is 'ZYX'
%
%   Example
%   rotation3dToEulerAngles
%
%   References
%   Code from Graphics Gems IV on euler angles
%   http://tog.acm.org/resources/GraphicsGems/gemsiv/euler_angle/EulerAngles.c
%   Modified using explanations in:
%   http://www.gregslabaugh.name/publications/euler.pdf
%
%   See also
%   transforms3d, rotation3dAxisAndAngle, createRotation3dLineAngle,
%   eulerAnglesToRotation3d
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-08-11,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

p = inputParser;
validStrings = {'ZYX','ZYZ'};
addOptional(p,'convention','ZYX',@(x) any(validatestring(x,validStrings)));
parse(p,varargin{:});
convention=p.Results.convention;

% conversion from radians to degrees
k = 180 / pi;

switch convention
    case 'ZYX'
        % extract |cos(theta)|
        cy = hypot(mat(1,1), mat(2,1));
        % avoid dividing by 0
        if cy > 16*eps
            % normal case: theta <> 0
            phi   = k * atan2( mat(2,1), mat(1,1));
            theta = k * atan2(-mat(3,1), cy);
            psi   = k * atan2( mat(3,2), mat(3,3));
        else
            %
            phi   = 0;
            theta = k * atan2(-mat(3,1), cy);
            psi   = k * atan2(-mat(2,3), mat(2,2));
        end
    case 'ZYZ'
        cy = hypot(mat(3,2), mat(3,1));
        if cy > 16*eps
            phi   = k * -atan2(mat(2,3), -mat(1,3));
            theta = k * -atan2(cy, mat(3,3));
            psi   = k * -atan2(mat(3,2), mat(3,1));
        else
            phi   = 0;
            theta = k * atan2(cy, mat(3,3));
            psi   = k * atan2(mat(2,1), mat(2,2));
        end
end

% format output arguments
if nargout <= 1
    % one array
    varargout{1} = [phi theta psi];
else
    % three separate arrays
    varargout = {phi, theta, psi};
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function v = so3_hatinv(xi)
v = [xi(3,2); xi(1,3); xi(2,1)];

author: RomaTeng

repo: EKF-SLAM-on-Manifold


text:
classdef BMatrixComputer < handle

    properties (Access = private)
        fun
        dNdx
        nVoigt
    end

    methods (Access = public)

        function obj = BMatrixComputer(cParams)
            obj.init(cParams);
        end

        function B = compute(obj,igaus)
            ndimf = obj.fun.ndimf;
            switch ndimf
                case 1
                    obj.nVoigt = 2;
                    B = obj.computeBin1D(igaus);
                case 2
                    obj.nVoigt = 3;
                    B = obj.computeBin2D(igaus);
                case 3
                    obj.nVoigt = 6;
                    B = obj.computeBin3D(igaus);
            end
        end

    end

    methods (Access = private)

        function init(obj, cParams)
            obj.fun  = cParams.fun;
            obj.dNdx = cParams.dNdx;
        end

        function B = computeBin2D(obj,iGaus)
            deriv = obj.dNdx;
            nStre = obj.nVoigt;
            nDimf = obj.fun.ndimf;
            nNodE = size(deriv,2);
            nDofE = nNodE*nDimf;
            nElem = size(deriv,3);
            B = zeros(nStre,nDofE,nElem);
            for iNode = 1:nNodE
                j = nDimf*(iNode-1)+1;
                B(1,j,:)   = deriv(1,iNode,:,iGaus);
                B(2,j+1,:) = deriv(2,iNode,:,iGaus);
                B(3,j,:)   = deriv(2,iNode,:,iGaus);
                B(3,j+1,:) = deriv(1,iNode,:,iGaus);
            end
        end

        function B = computeBin3D(obj,iGaus)
            deriv = obj.dNdx;
            nNode = size(deriv,2);
            nElem = size(deriv,3);
            B = zeros(obj.nVoigt,nNode,nElem);
            for inode = 1:nNode
                j = obj.fun.ndimf*(inode-1)+1;
                % associated to normal strains
                B(1,j,:)   = deriv(1,inode,:,iGaus);
                B(2,j+1,:) = deriv(2,inode,:,iGaus);
                B(3,j+2,:) = deriv(3,inode,:,iGaus);
                % associated to shear strain, gamma12
                B(4,j,:)   = deriv(2,inode,:,iGaus);
                B(4,j+1,:) = deriv(1,inode,:,iGaus);
                % associated to shear strain, gamma13
                B(5,j,:)   = deriv(3,inode,:,iGaus);
                B(5,j+2,:) = deriv(1,inode,:,iGaus);
                % associated to shear strain, gamma23
                B(6,j+1,:) = deriv(3,inode,:,iGaus);
                B(6,j+2,:) = deriv(2,inode,:,iGaus);
            end
        end

        function [B] = computeBin1D(obj, igaus)
            deriv  = obj.dNdx(:,:,:,igaus);
            nDimf = obj.fun.ndimf;
            nNode = size(deriv,2);
            nElem = size(obj.dNdx,3);
            nDofs = nDimf*nNode;
            B = zeros(obj.nVoigt,nDofs,nElem);
            for inode = 1:nNode
                j = nDimf*(inode-1) + 1;
                B(1,j,:) = deriv(1,inode,:);
                B(2,j,:) = deriv(2,inode,:);
            end
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function psi = psiEval(l, C, LR, N, dim, nVars)
%PSIEVAL   Evaluate a psi-function.
%   PSI = PSIEVAL(L, C, LR, N, DIM, NVARS) evaluates the psi-function of index L
%   and coefficient C with the contour LR, N grid points, in dimension DIM and
%   with NVARS variables.
%
% See also EXPINTEG/PHIEVAL, EXPINTEG/PHIFUN.

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

% Get a function handle to the phi-function of index L:
phi = expinteg.phiFun(l);

% Evaluate the psi-function with a contour integral:
psi = mean(C^l*feval(phi, C*LR), 2);

% Reshape it when nVars>1 or/and dim>1:
psi = reshape(psi, nVars*N, N^(dim>1), N^(dim>2));

end


text:
function pde = constant_fun(bm,bp,r,x0,y0,z0,coef)
%% USAGE: polynomial solution for Poisson equation
%  Last Modified: 02/21/2020 by Xu Zhang

%% PDE Structure
pde = struct('intf',@intf,...
    'exactu1',@exactu1,'exactu2',@exactu2,'exactu3',@exactu3,...
    'um1',@um1,'um2',@um2,'um3',@um3,'up1',@up1,'up2',@up2,'up3',@up3,...
    'Dxu',@Dxu,'Dxum',@Dxum,'Dxup',@Dxup,'Dyu',@Dyu,...
    'Dyum',@Dyum,'Dyup',@Dyup,'Dzu',@Dzu,'Dzum',@Dzum,'Dzup',@Dzup,...
    'f1',@f1,'f2',@f2,'f3',@f3,...
    'fm1',@fm1,'fm2',@fm2,'fm3',@fm3,...
    'fp1',@fp1,'fp2',@fp2,'fp3',@fp3,...
    'A',@A,'Am',@Am,'Ap',@Ap,'one',@one);

pde.bm = bm;
pde.bp = bp;
%% interface function
    function u = intf(x,y,z)
        u = ((x-x0).^2 + (y-y0).^2 + (z-z0).^2).^(1/2)/r-1;
    end

%% exact solution
    function u = exactu1(x,y,z)
        u = um1(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up1(x(id),y(id),z(id));
    end
    function u = exactu2(x,y,z)
        u = um2(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up2(x(id),y(id),z(id));
    end
    function u = exactu3(x,y,z)
        u = um3(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up3(x(id),y(id),z(id));
    end
    function u = um1(x,y,z)
        u = coef*ones(size(x));
    end
    function u = um2(x,y,z)
        u = coef*ones(size(x));
    end
    function u = um3(x,y,z)
        u = coef*ones(size(x));
    end
    function u = up1(x,y,z)
        u = coef*ones(size(x));
    end
    function u = up2(x,y,z)
        u = coef*ones(size(x));
    end
    function u = up3(x,y,z)
        u = coef*ones(size(x));
    end

%% Boundary Function
    function u = gD1(x,y,z)
        u = exactu1(x,y,z);
    end
    function u = gD2(x,y,z)
        u = exactu2(x,y,z);
    end
    function u = gD3(x,y,z)
        u = exactu3(x,y,z);
    end
%% Derivative of the exact solution
    function u = Dxu(x,y,z)
        u = Dxum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dxup(x(id),y(id),z(id));
    end
    function u = Dyu(x,y,z)
        u = Dyum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dyup(x(id),y(id),z(id));
    end
    function u = Dzu(x,y,z)
        u = Dzum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dzup(x(id),y(id),z(id));
    end
    function u = Dxum(x,y,z)
        u = zeros(size(x));
    end
    function u = Dyum(x,y,z)
        u = zeros(size(x));
    end
    function u = Dzum(x,y,z)
        u = zeros(size(x));
    end
    function u = Dxup(x,y,z)
        u = zeros(size(x));
    end
    function u = Dyup(x,y,z)
        u = zeros(size(x));
    end
    function u = Dzup(x,y,z)
        u = zeros(size(x));
    end

    function u = Duker(x,y,z)
        u = zeros(size(x));
    end

%% right hand side function
    function u = f1(x,y,z)
        u = fm1(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp1(x(id),y(id),z(id));
    end
    function u = f2(x,y,z)
        u = fm2(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp2(x(id),y(id),z(id));
    end
    function u = f3(x,y,z)
        u = fm3(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp3(x(id),y(id),z(id));
    end

    function u = fm1(x,y,z)
        u = coef*ones(size(x));
    end
    function u = fm2(x,y,z)
        u = coef*ones(size(x));
    end
    function u = fm3(x,y,z)
        u = coef*ones(size(x));
    end
    function u = fp1(x,y,z)
        u = coef*ones(size(x));
    end
    function u = fp2(x,y,z)
        u = coef*ones(size(x));
    end
    function u = fp3(x,y,z)
        u = coef*ones(size(x));
    end

%% Diffusion coefficient function
    function u = A(x,y,z)
        u = Am(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Ap(x(id),y(id),z(id));
    end
    function u = Am(x,y,z)
        u = bm*ones(size(x));
    end
    function u = Ap(x,y,z)
        u = bp*ones(size(x));
    end

%% Other function
    function u = one(x,y,z)
        u = ones(size(x));
    end
end

repo: ifem


text:
% @FACTORIZE: an object-oriented method for solving linear system, solving
% least-square problems, and for efficient computation of mathematical
% expressions that use inv(A).
%
%   disp      - displays the factorization F
%   double    - returns the factorization as a single matrix, A or inv(A)
%   end       - returns index of last item for use in subsref
%   factorize - an object-oriented method for solving linear systems
%   inverse   - "inverts" F by flagging it as the factorization of inv(A).
%   mldivide  - x = A\b using the factorization F = factorize(A)
%   mrdivide  - x = b/A using the factorization F = factorize(A)
%   mtimes    - A*b, inv(A)*b, b*A, or b*inv(A), without computing inv(A)
%   size      - returns the size of the matrix F.A in the factorization F
%   subsref   - A(i,j) or (i,j)th entry of inv(A) if F is inverted.
%   plus      - update a dense Cholesky factorization
%   minus     - downdate a dense Cholesky factorization
%
% Example
%
%   F = factorize(A) ;
%   x = F\b ;           % same as x=A\b

% Copyright 2009, Timothy A. Davis, University of Florida

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% The COBRAToolbox: testGeometricFBA.m
%
% Purpose:
%     - Tests whether geometricFBA works as intended.
%
% Authors:
%     - Thomas Pfau, Oct 2018

% save the current path
currentDir = pwd;

% The testmodel used is structured as follows:
%
%   <-> A -> B ---> C --> E <->
%        \               ^
%         \             /
%           -> F -->  G
%
% Thus the flux can be split between A -> B -> C -> E and A -> F -> G -> E

% load the test models
model = createToyModelForGeoFBA();
% model with Constraints
modelWConst = addCOBRAConstraints(model,{'R1'},400); % R1 restricted to a maximum flux of 400

%When detectDeadEnds is changed according to Ronans suggestion, we need to test
%multiple solvers.
solverPkgs = prepareTest('needsLP', true, 'minimalMatlabSolverVersion',8.0, 'excludeSolvers',{'pdco','ddqMinos','quadMinos','matlab'});
tol = 1e-4;
for k = 1:length(solverPkgs.LP)
    % set the solver
    changeCobraSolver(solverPkgs.LP{k}, 'LP', 0);
    fprintf('Testing geometricFBA with solver %s ...\n',solverPkgs.LP{k});
    sol = geometricFBA(model);
    assert(all(abs(sol(1:6) - 500) < tol))
    sol = geometricFBA(modelWConst);
    % flux is now centered between 0 and 400 for R11.. R3 and between 600
    % and 1000 for R4..6
    assert(all(abs(sol(1:3) - 200) < tol))
    assert(all(abs(sol(4:6) - 800) < tol))
end

fprintf('Done...\n');
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function [xi,yi,curvatures] = maxr(a)
%MAXR   Find interpolated maximizer(s) and max value(s)
%       for (each column of) a.
%
%               [xi,yi,curvatures] = maxr(a)
%
%       Calls max() followed by qint() for quadratic interpolation.
%
   [m,n] = size(a);
   if m==1, a=a'; t=m; m=n; n=t; end;
   [y,x] = max(a);
   xi=x;    % vector of maximizer locations, one per col of a
   yi=y;    % vector of maximum values, one per column of a
   if nargout>2, curvatures = zeros(1,n); end
   for j=1:n,   % loop over columns
     if x(j)>1, % only support interior maxima
       if x(j)<m,
         [xdelta,yij,cj] = qint(a(x(j)-1,j),y(j),a(x(j)+1,j));
         xi(j) = x(j) + xdelta;
         if nargout>2, curvatures(j) = cj; end
         if (nargout>1), yi(j) = yij; end
       end;
     end;
   end;

repo: REST


text:
function visualizetiles(board, orientation, tiles, visualizationMethod)
% VISUALIZETILES Draws the board and the tiles

% The MATLAB Contest Team
% Copyright 2012 The MathWorks, Inc.

clf;
if visualizationMethod == 1
    BoardVisualizer(tiles, board, orientation);
elseif visualizationMethod == 2
    ScoreVisualizer(tiles, board ,orientation);
end

end

function h = BoardVisualizer(tiles, pos, rot)

[row,col] = size(pos);
boardcapacity = nnz(pos);
y_tile_loc = (pos.') ~= 0;
x_tile_loc = (pos.') ~= 0;

%Rotate tiles matrix into correct position
tiles = rotateTiles(tiles,rot);

%Rearrange tiles and add zeros
inds = pos';
inds = inds(:);
inds(inds == 0) = [];
tiles = tiles(inds,:);

%Set up figure
axis equal
axis off
hold on
colormap(rot90(pink,2));
set(gcf,'color','w')

%Plot background mesh
h.mesh = mesh(0:col,0:row,zeros(row+1,col+1));
set(h.mesh,'facealpha',0,'linewidth',1,'edgecolor','k');

if isempty(tiles)
    return;  % empty board -- nothing else to do
end

%Plot patch for north, east, south, and west wedges
X = repmat([0:col-1; .5:col; 1:col],1,row);
X = X(:,x_tile_loc);
y1 = repmat(row:-1:1,col,1);
y1 = y1(y_tile_loc);
Y = [reshape(y1,1,boardcapacity);
    reshape(y1,1,boardcapacity)-.5;
    reshape(y1,1,boardcapacity)];
h.north = patch('XData',X, 'YData',Y, 'CDataMapping','scaled', 'CData',tiles(:,1)','linestyle','none', 'FaceColor','flat');
hold on
Xt = X(2,:);
Yt = Y(1,:) - 1/6;
h.northText = text(Xt,Yt,num2str(tiles(:,1)),'fontweight','bold');

X = repmat([1:col; .5:col; 1:col],1,row);
X = X(:,x_tile_loc);
y1 = repmat(row:-1:1,col,1);
y1 = y1(y_tile_loc);
Y = [reshape(y1,1,boardcapacity);
    reshape(y1,1,boardcapacity)-.5;
    reshape(y1,1,boardcapacity)-1];
h.east = patch('XData',X, 'YData',Y, 'CDataMapping','scaled', 'CData',tiles(:,2)','linestyle','none', 'FaceColor','flat');
Xt = X(1,:) - 1/6;
Yt = Y(2,:);
h.eastText = text(Xt,Yt,num2str(tiles(:,2)),'fontweight','bold');

X = repmat([0:col-1; .5:col; 1:col],1,row);
X = X(:,x_tile_loc);
y1 = repmat(row-1:-1:0,col,1);
y1 = y1(y_tile_loc);
Y = [reshape(y1,1,boardcapacity);
    reshape(y1,1,boardcapacity)+.5;
    reshape(y1,1,boardcapacity)];
h.south = patch('XData',X, 'YData',Y, 'CDataMapping','scaled', 'CData',tiles(:,3)','linestyle','none', 'FaceColor','flat');
Xt = X(2,:);
Yt = Y(1,:) + 1/6;
h.southText = text(Xt,Yt,num2str(tiles(:,3)),'fontweight','bold');

X = repmat([0:col-1; .5:col; 0:col-1],1,row);
X = X(:,x_tile_loc);
y1 = repmat(row:-1:1,col,1);
y1 = y1(y_tile_loc);
Y = [reshape(y1,1,boardcapacity);
    reshape(y1,1,boardcapacity)-.5;
    reshape(y1,1,boardcapacity)-1];
h.west = patch('XData',X, 'YData',Y, 'CDataMapping','scaled', 'CData',tiles(:,4)','linestyle','none', 'FaceColor','flat');
Xt = X(1,:) + 1/6;
Yt = Y(2,:);
h.westText = text(Xt,Yt,num2str(tiles(:,4)),'fontweight','bold');

clim = [min(tiles(:)), max(tiles(:))];
range = diff(clim);
clim(1) = clim(1) - 0.1*range;
clim(2) = clim(2) + 0.1*range;
set(gca,'CLim',clim);
end

function h = ScoreVisualizer(tiles,pos,rot)

%Get size information and initialize values matrix
brow = size(pos,1);
bcol = size(pos,2);
bsize = [brow, bcol];
vsize = [2*brow+1,bcol+1];

%Rotate tiles matrix into correct position
tiles = rotateTiles(tiles,rot);
%Determine the score of each element in values matrix
values = evaluateBoard(tiles,pos,vsize,bsize);

h = drawResults(values,bsize);
end

function tout = rotateTiles(tin,rot)
tout = zeros(size(tin));
for i = 1:length(rot)
    tout(i,:) = circshift(tin(i,:), [0 1-rot(i)]);
end
end

function values = evaluateBoard(tiles,pos,vsize,bsize)
values = zeros(vsize(1),vsize(2));
values(logical(mod(1:vsize(1),2)),end) = NaN(bsize(1)+1,1);
for r = 1:2:vsize(1)
    for c = 1:vsize(2)-1
        above = (r-1)/2;
        below = (r+1)/2;
        if r == 1 %first row
            if pos(below,c) ~= 0
                values(r,c) = abs(tiles(pos(below,c),1));
            end
        elseif r == vsize(1) %last row
            if pos(above,c) ~= 0
                values(r,c) = abs(tiles(pos(above,c),3));
            end
        else %middle rows
            if pos(below,c) ~= 0 && pos(above,c) ~= 0
                values(r,c) = abs(tiles(pos(above,c),3)-...
                    tiles(pos(below,c),1));
            elseif pos(below,c) ~= 0 && pos(above,c) == 0
                values(r,c) = abs(tiles(pos(below,c),1));
            elseif pos(below,c) == 0 && pos(above,c) ~= 0
                values(r,c) = abs(tiles(pos(above,c),3));
            end
        end
    end
end
for r = 2:2:vsize(1)
    for c = 1:vsize(2)
        left = c-1;
        right = c;
        if c == 1 % first column
            if pos(r/2,right) ~= 0
                values(r,c) = abs(tiles(pos(r/2,right),4));
            end
        elseif c == vsize(2) % last column
            if pos(r/2,left) ~= 0
                values(r,c) = abs(tiles(pos(r/2,left),2));
            end
        else % middle columns
            if pos(r/2,right) ~= 0 && pos(r/2,left) ~= 0
                values(r,c) = abs(tiles(pos(r/2,right),4)-...
                    tiles(pos(r/2,left),2));
            elseif pos(r/2,right) ~= 0 && pos(r/2,left) == 0
                values(r,c) = abs(tiles(pos(r/2,right),4));
            elseif pos(r/2,right) == 0 && pos(r/2,left) ~= 0
                values(r,c) = abs(tiles(pos(r/2,left),2));
            end
        end
    end
end
end

function h = drawResults(values,b)
%Define x coordinates of patches
X = repmat([  0:b(2)-1  -0.5:b(2)  ;
    0.5:b(2)       0:b(2)  ;
    1:b(2)     0.5:b(2)+1;
    0.5:b(2)       0:b(2)   ],1,b(1));
X = [X [  0:b(2)-1  ;
    0.5:b(2)    ;
    1:b(2)    ;
    0.5:b(2)]];

%Define y coordinates of patches
Y =    repmat([b(1)     b(1)-.5  b(1)     b(1)+.5]',1,b(2));
Y = [Y repmat([b(1)-.5  b(1)-1   b(1)-.5  b(1)   ]',1,b(2)+1)];
Y = repmat(Y,1,b(1));
start = 1;
for i = 1:b(1)
    Y(:,start:i*(b(2)*2+1)) = Y(:,start:i*(b(2)*2+1)) - (i-1);
    start = i*(b(2)*2+1) + 1;
end
Y = [Y repmat([0 -.5 0 .5]',1,b(2))];

%Define colors of patches
C = values';
C = C(:)';
C(isnan(C)) = [];

h.p = patch(X,Y,C);
set(h.p,'linestyle','none');
axis equal
axis off
set(gcf,'color','k')
colormap copper
if max(values(:)) < 1e-8
    set(gca,'climmode','manual','clim',[0 1]);
end

hold on
h.m = mesh(0:b(2),0:b(1),zeros(b(1)+1,b(2)+1));
set(h.m,'facealpha',0,'linewidth',2,'edgecolor','w')
hold off
end


text:
% Implement 1D CFAR using lagging cells on the given noise and target scenario.

% Close and delete all currently open figures
close all;

% Data_points
Ns = 1000;

% Generate random noise
s=abs(randn(Ns,1));

%Targets location. Assigning bin 100, 200, 300 and 700 as Targets with the amplitudes of 8, 9, 4, 11.
s([100 ,200, 350, 700])=[8 15 7 13];

%plot the output
plot(s);

% TODO: Apply CFAR to detect the targets by filtering the noise.

% 1. Define the following:
% 1a. Training Cells
T = 12;
% 1b. Guard Cells
G = 4;

% Offset : Adding room above noise threshold for desired SNR
offset=5;

% Vector to hold threshold values
threshold_cfar = [];

%Vector to hold final signal after thresholding
signal_cfar = [];

% 2. Slide window across the signal length
for i = 1:(Ns-(G+T+1))

    % 2. - 5. Determine the noise threshold by measuring it within the training cells

    noise_level =sum(s(i:i+T-1));

    % 6. Measuring the signal within the CUT

    threshold = (noise_level/T)*offset;
    threshold_cfar=[threshold_cfar,{threshold}];

    signal=s(i+T+G);

    % 8. Filter the signal above the threshold
    if (signal<threshold)
        signal=0;
    end
    signal_cfar = [signal_cfar, {signal}];
end

% plot the filtered signal
plot (cell2mat(signal_cfar),'g--');

% plot original sig, threshold and filtered signal within the same figure.
figure,plot(s);
hold on,plot(cell2mat(circshift(threshold_cfar,G)),'r--','LineWidth',2)
hold on, plot (cell2mat(circshift(signal_cfar,(T+G))),'g--','LineWidth',4);
legend('Signal','CFAR Threshold','detection')

author: williamhyin


text:
function HT = SoftThreshold(input,H);
% HT = input;
% for i = 1:numel(input)
%     if abs(input(i)) < H
%         HT(i) = 0;
%     end
%     if abs(input(i)) > H
%         HT(i) = sign(input(i))*(abs(input(i))-H);
%     end
% end
HT = sign(input).*max(abs(input)-H,0);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% DEMROBOTWIRELESSFGPLVM1 Wireless Robot data from University of Washington, without dynamics and without back constraints.

% FGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'robotWireless';
experimentNo = 1;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = fgplvmOptions('ftc');
latentDim = 2;
d = size(Y, 2);

model = fgplvmCreate(latentDim, d, Y, options);

% Optimise the model.
iters = 1000;
display = 1;

model = fgplvmOptimise(model, display, iters);

% Save the results.
modelWriteResult(model, dataSetName, experimentNo);

if exist('printDiagram') & printDiagram
  lvmPrintPlot(model, lbls, dataSetName, experimentNo);
end

% Load the results and display dynamically.
lvmResultsDynamic(model.type, dataSetName, experimentNo, 'vector')

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function v = dq2vel(dq)

% DQ2VEL     Transforms a dual quaternion point velocity representation
%            into a vector representation.
%
%        V = DQ2VEL(DQ) transforms the dual quaternion representation of
%        a point velocity (in the euclidean space) into a vector v, which
%        represents the point velocity coordinates.
%        DQ is either a vector of size 8 or an array of size 8*N (each
%        column represents a point velocity dual quaternion) where N is the
%        number of points. V is a vector of size 3 or an array of size 3*N
%        depending on the input format.
%
% See also VEL2DQ, DQ2POS, DQ2LINE, DQ2LINEVEL

sdq = size(dq);
if sdq == [1 8]
    dq = dq.';
    sdq = size(dq);
end

% wrong format
if sdq(1) ~= 8
    error('DualQuaternion:dquat2vel:wrongsize',...
        '%d rows in the DQ array. It should be 8.',sdq(1));
end

% extraction of the point velocity coordinates
v = dq(6:8,:);


text:
function sSubject = db_parse_subject( subjectsDir, subjectSubDir, sizeProgress )
% DB_PARSE_SUBJECT: Parse a subject directory.
%
% USAGE:  sSubject = db_parse_subject(subjectsDir, subjectSubDir, sizeProgress);
%
% INPUT:
%     - subjectsDir   : Top-level subjects directory (protocol/anatomy/)
%     - subjectSubDir : Study subdirectory
%     - sizeProgress : The full process increments the available progress bar by this amount
% OUTPUT:
%     - sSubject : array of Subject structures (one entry for each subdirectory with a brainstormsubject file),
%                  or [] if no brainstormsubject.mat file was found in directory.

% NOTES:
%    In a subject subdirectory there is :
%        - One and only one 'brainstormsubject*.mat' file,
%        - All the files associated with this subject :
%             - Surfaces (cortex, scalp, inner skull, outer skull, ...)
%             - Anatomical MRI
%    All the .MAT files in the directory are associated with the subject.
%    The files that are pointed by the 'Anatomy' and surfaces fields
%        ('Scalp','Cortex',...) in the brainstormsubject.mat file are the
%        current files for each category.
%
% WARNING: This function ignores the <DirDefaultSubject> directories.
%          These directories must be parsed independently.
%
% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2008-2021

%% ===== PARSE INPUTS =====
if ~file_exist(bst_fullfile(subjectsDir, subjectSubDir))
    error('Path in argument does not exist, or is not a directory');
end
if (nargin < 3) || isempty(sizeProgress)
    sizeProgress = [];
end

%% ===== LOOK FOR SUBJECT DEFINITION =====
% List 'brainstormsubject*.mat' files in the subjectSubDir directory
subjectFiles = dir(bst_fullfile(subjectsDir, subjectSubDir, 'brainstormsubject*.mat'));
% If there is no subject definition file : ignore the directory
if isempty(subjectFiles)
    subjMat = [];
% If there is more that one subject definition file : warning and ignore the directory
elseif (length(subjectFiles) > 1)
    subjMat = [];
    bst_error(sprintf('There is more than one brainstormsubject file in directory ''%s'' : ignoring directory.', bst_fullfile(subjectsDir, subjectSubDir)), ...
             'Parsing subject directory', 0);
% Else there is only one file in subjectFiles list
else
    subjectFile = subjectFiles(1);
    % Open file and get subject description
    try
        subjMat = load(bst_fullfile(subjectsDir, subjectSubDir, subjectFile.name));
    catch
        warning('Brainstorm:CannotOpenFile', 'Cannot open file ''%s'' : ignoring subject', bst_fullfile(subjectsDir, subjectSubDir, subjectFile.name));
        subjMat = [];
    end
end

%% ===== READ SUBJECT DESCRIPTION =====
% Copy fields of loaded subject matrix into the sSubject structure (output variable)
% IMPORTANT : All the filenames in the subject description matrix
%             must be relative to the SUBJECTS directory.
if ~isempty(subjMat)
    % Define structure templates
    sSubject = db_template('Subject');
    % Store subject's .MAT filename
    sSubject.FileName = bst_fullfile(subjectSubDir, subjectFile.name);
    sSubject.Name     = subjectSubDir;
    % Comments
    if isfield(subjMat, 'Comments')
        sSubject.Comments = subjMat.Comments;
    end
    % DateOfAcquisition
    if isfield(subjMat, 'DateOfAcquisition')
        sSubject.DateOfAcquisition = subjMat.DateOfAcquisition;
    end
    % UseDefaultAnat
    if isfield(subjMat, 'UseDefaultAnat')
        sSubject.UseDefaultAnat = subjMat.UseDefaultAnat;
    end
    % UseDefaultChannel
    if isfield(subjMat, 'UseDefaultChannel')
        sSubject.UseDefaultChannel = subjMat.UseDefaultChannel;
    end
else
    % Define structure templates
    sSubject = repmat(db_template('Subject'), 0);
end

%% ===== READ ALL FILES IN FOLDER =====
% List all '*.mat' files in the subjectSubDir directory
allFiles = dir(bst_fullfile(subjectsDir, subjectSubDir, '*'));
% Exclude files starting with a '.'
dirFiles = repmat(allFiles,0);
for iFile = 1:length(allFiles)
    if (allFiles(iFile).name(1) ~= '.')
        dirFiles(end+1,1) = allFiles(iFile);
    end
end
if isempty(dirFiles)
    return;
end
% Progress bar?
isProgressBar = bst_progress('isvisible') && ~isempty(sizeProgress);
startValue = bst_progress('get');
% Process all the files
for iFile = 1:length(dirFiles)
    % Increment progress bar
    if isProgressBar
        bst_progress('set',  round(startValue + iFile/length(dirFiles) * sizeProgress));
    end
    % Process sub-directories recursively (excluding the 'DirDefaultSubject' directories)
    if (dirFiles(iFile).isdir) && ~strcmpi(dirFiles(iFile).name, bst_get('DirDefaultSubject'))
        % Append subjects in the subdirectory to the sSubject array
        sSubject = [sSubject, db_parse_subject(subjectsDir, bst_fullfile(subjectSubDir, dirFiles(iFile).name))];
    % Files (only if there is a brainstormsubject file in current directory)
    elseif ~isempty(subjMat)
        % Reconstruct filename
        filenameRelative = bst_fullfile(subjectSubDir, dirFiles(iFile).name);
        filenameFull     = bst_fullfile(subjectsDir, subjectSubDir, dirFiles(iFile).name);
        % Skip "_openmeeg.mat" files (related with the BEM computation)
        if ~isempty(strfind(filenameRelative, '_openmeeg.mat'))
            continue;
        end
        % Determine filetype
        fileType = file_gettype(filenameFull);
        if isempty(fileType)
            continue;
        end
        % Process file
        switch(fileType)
            case 'subjectimage'
                % Try to load Anatomy info
                tempAnatomy = io_getAnatomyInfo(filenameRelative);
                % If an anatomy was identified : add it to the subject structure
                if ~isempty(tempAnatomy)
                    sSubject(1).Anatomy(end+1) = tempAnatomy;
                end
            case {'scalp', 'cortex', 'innerskull', 'outerskull', 'fibers', 'fem', 'tess'}
                tempSurface = io_getSurfaceInfo(filenameRelative);
                if ~isempty(tempSurface)
                    switch(fileType)
                        case 'scalp',       tempSurface.SurfaceType = 'Scalp';
                        case 'cortex',      tempSurface.SurfaceType = 'Cortex';
                        case 'innerskull',  tempSurface.SurfaceType = 'InnerSkull';
                        case 'outerskull',  tempSurface.SurfaceType = 'OuterSkull';
                        case 'fibers',      tempSurface.SurfaceType = 'Fibers';
                        case 'fem',         tempSurface.SurfaceType = 'FEM';
                        case 'tess',        tempSurface.SurfaceType = 'Other';
                    end
                    sSubject(1).Surface(end+1) = tempSurface;
                end
        end
    end
end
% Set progress bar to the end
if isProgressBar
    bst_progress('set', round(startValue + sizeProgress));
end

%% ===== DEFAULT ANATOMY/SURFACES =====
% Sort anatomy files: first T1 and MRI, then volume atlases
if ~isempty(sSubject) && (length(sSubject(1).Anatomy) >= 2)
    iNoAtlas = find(cellfun(@(c)isempty(strfind(c, '_volatlas')), {sSubject(1).Anatomy.FileName}));
    if (length(iNoAtlas) >= 2)
        iT1 = find(~cellfun(@(c)isempty([strfind(lower(c), '_t1'), strfind(lower(c), '_mri')]), {sSubject(1).Anatomy(iNoAtlas).FileName}));
        iNoAtlas = [iNoAtlas(iT1), setdiff(iNoAtlas, iNoAtlas(iT1))];
    end
    iSort = [iNoAtlas, setdiff(1:length(sSubject(1).Anatomy), iNoAtlas)];
    sSubject(1).Anatomy = sSubject(1).Anatomy(iSort);
end
% If there is a brainstormsubject.mat
if ~isempty(subjMat)
    % The brainstormsubject.mat can define what are the defaults files for the different
    % file categories : Anatomy, Scalp, Cortex, InnerSkull, OuterSkull, Fibers, FEM
    % Now we are going to read those brainstormsubject fields look for
    % the pointed files in the sSubject structure

    % ==== ANATOMY ====
    % By default : use the first anatomy in list (which is not a volume atlas)
    if ~isempty(sSubject(1).Anatomy)
        sSubject(1).iAnatomy = 1;
    else
        sSubject(1).iAnatomy = [];
    end
    % Anatomy field : filename of the current/default Anatomy for this subject (subjectimage.mat)
    % => If there are more than one MRI in subject directory,
    %    this one will be used by default for all the processes.
    if (isfield(subjMat, 'Anatomy') && ~isempty(subjMat.Anatomy))
        ind = find(file_compare({sSubject(1).Anatomy.FileName}, subjMat.Anatomy), 1);
        if ~isempty(ind)
            % Reorder anatomy entries
            sSubject(1).Anatomy = sSubject(1).Anatomy([ind, setdiff(1:length(sSubject(1).Anatomy), ind)]);
            sSubject(1).iAnatomy = 1;
        end
    end

    % ==== SURFACE ====
    % Sort surfaces by categories
    subjectSurfaces = db_surface_sort(sSubject(1).Surface);
    % Select one surface in each category
    for surfaceCatergory = {'Scalp', 'Cortex', 'InnerSkull', 'OuterSkull', 'Fibers', 'FEM'}
        % By default : use the last surface in list
        if ~isempty(subjectSurfaces.(surfaceCatergory{1}))
            sSubject(1).(['i' (surfaceCatergory{1})]) = subjectSurfaces.(['Index' surfaceCatergory{1}])(end);
        else
            sSubject(1).(['i' (surfaceCatergory{1})]) = [];
        end
        % If a default was defined in the brainstormsubject*.mat
        if (isfield(subjMat, surfaceCatergory{1}) && ~isempty(subjMat.(surfaceCatergory{1})) && ischar(subjMat.(surfaceCatergory{1})))
            ind = find(file_compare({sSubject(1).Surface.FileName}, subjMat.(surfaceCatergory{1})), 1);
            if (ind > 0)
                sSubject(1).(['i' (surfaceCatergory{1})]) = ind;
            end
        end
    end
end

%% ===================================================================================
%  === HELPER FUNCTIONS ==============================================================
%  ===================================================================================
    % Load an Anatomy from a relative filename
    % Return an Anatomy structure :
    %    |- FileName
    %    |- Comment
    % or an empty structure if anatomy is not or badly defined.
    function Anatomy = io_getAnatomyInfo(relativeFilename)
        Anatomy = repmat(db_template('Anatomy'), 0);
        % Check if the file exists, and load Comment field
        if file_exist(bst_fullfile(subjectsDir, relativeFilename))
            try
                % Load MRI 'Comment' field from SUBJECTIMAGE.MAT
                warning off MATLAB:load:variableNotFound;
                anatomyMat = load(bst_fullfile(subjectsDir, relativeFilename), 'Comment');
                warning on MATLAB:load:variableNotFound;
                % If 'Comment' not defined or empty
                if (length(fieldnames(anatomyMat)) == 0)
                    [fpath, basename] = bst_fileparts(relativeFilename);
                    Anatomy(1).Comment = basename;
                else
                    Anatomy(1).Comment = anatomyMat.Comment;
                end
                Anatomy(1).FileName = relativeFilename;
            catch
                % An error occured during the 'load' operation
                warning('Brainstorm:CannotOpenFile', 'Cannot open anatomy file ''%s''.', bst_fullfile(subjectsDir, relativeFilename));
            end
        else
            % File does not exist
            warning('Brainstorm:FileNotFound', 'Anatomy file ''%s'' was not found.', bst_fullfile(subjectsDir, relativeFilename));
        end
    end

    % Load a surface from a relative filename
    % surfaceType = {'Scalp', 'Cortex', 'InnerSkull', 'OuterSkull', 'Fibers', 'FEM'}
    % Return a Surface structure :
    %    |- FileName
    %    |- Comment
    % or an empty structure if surface is not or badly defined
    function Surface = io_getSurfaceInfo(relativeFilename)
        Surface = repmat(db_template('Surface'), 0);
        % Check if the file exists, and load Comment field
        if file_exist(bst_fullfile(subjectsDir, relativeFilename))
            try
                % Load Surface 'Comment' field from SUBJECTIMAGE.MAT
                surfaceMat = load(bst_fullfile(subjectsDir, relativeFilename), 'Comment');
                % If 'Comment' not defined or empty
                if (isempty(surfaceMat))
                    Surface(1).Comment = '';
                elseif iscell(surfaceMat.Comment)
                    Surface(1).Comment = surfaceMat.Comment{1};
                elseif ischar(surfaceMat.Comment)
                    Surface(1).Comment = surfaceMat.Comment;
                end
                Surface(1).FileName = relativeFilename;
            catch
                % An error occured during the 'load' operation
                warning('Brainstorm:CannotOpenFile', 'Cannot open surface file ''%s''. Adding ''.bak'' to the filename.', relativeFilename);
                fullfilename = bst_fullfile(subjectsDir, relativeFilename);
                file_move(fullfilename, [fullfilename, '.bak']);
            end
        else
            % File does not exist
            warning('Brainstorm:FileNotFound', 'Surface file ''%s'' was not found.', bst_fullfile(subjectsDir, relativeFilename));
        end
    end
end

author: brainstorm-tools


text:
function varargout = plotcoeffs(A, varargin)
%PLOTCOEFFS   PLOTCOEFFS for CHEBMATRIX objects.
%   PLOTCOEFFS(A) plots the Chebyshev coefficients CHEBMATRIX object A. If A
%   contains only CHEBFUN and DOUBLE objects, A is converted to a QUASIMATRIX,
%   and CHEBFUN/PLOTCOEFFS() is called. In this case PLOTCOEFFS(A, S) allows
%   various line types, plot symbols, and colors to be used, where S is a
%   character string. See CHEBFUN/PLOTCOEFFS() for further details.
%
%   If A contains inf x inf blocks, an error is thrown.
%
% See also CHEBFUN/PLOTCOEFFS.

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

% Deal with an empty input:
if ( isempty(A) )
    [varargout{1:nargout}] = plot([]);
    return
end

s = cellfun(@(b) min(size(b)), A.blocks);
isQuasi = all(isfinite(s(:)));

if ( ~isQuasi )
    % If A contains inf x inf blocks throw an error:
    error('CHEBFUN:CHEBMATRIX:plotcoeffs:notQuasi', ...
        'PLOTCOEFFS does not support CHEBMATRIX objects of size INFxINF.');
end

% If A contains only CHEBFUN or DOUBLE, convert it to a quasimatrix and call
% CHEBFUN/PLOTCOEFFS():
A.blocks = reshape(A.blocks, 1, numel(A.blocks));
A = quasimatrix(A.blocks);
[varargout{1:nargout}] = plotcoeffs(A, varargin{:});

end

author: chebfun


text:
function [Jul1,Jul2]=UTC2TDB(Jul1,Jul2,deltaT,clockLoc)
%%UTC2TDB Convert from universal coordinated time (UTC) to barycentric
%         dynamical time (TDB) to an accuracy of nanoseconds (if deltaT
%         is accurate) using the routines from the International
%         Astronomical Union's library that do not require external
%         ephemeris data.
%
%INPUTS: Jul1,Jul2 Two parts of a pseudo-Julian date given in UTC. The
%                  units of the date are days. The full date is the sum of
%                  both terms. The date is broken into two parts to provide
%                  more bits of precision. It does not matter how the date
%                  is split.
%           deltaT An optional parameter specifying the offset between UTC
%                  and UT1 in seconds. If this parameter is omitted or an
%                  empty matrix is passed, then the value of the function
%                  deltaTTUT1 with an appropriate offset will be used.
%         clockLoc An optional 3X1 vector specifying the location of the
%                  clock in WGS-84 ECEF Cartesian [x;y;z] coordinates with
%                  units of meters. Due to relativistic effects, clocks
%                  that are synchronized with respect to UTC are not
%                  synchronized with respect to TDB. If this parameter is
%                  omitted, then a clock at the center of the Earth is
%                  used and the precision declines to microseconds.
%
%OUTPUTS:Jul1,Jul2 Two parts of a Julian date given in TDB.
%
%This function converts the date to a terrestrial time and then calls
%TT2TDB.
%
%Many temporal coordinate systems standards are compared in [1].
%
%REFERENCES:
%[1] D. F. Crouse, "An Overview of Major Terrestrial, Celestial, and
%    Temporal Coordinate Systems for Target Tracking," Formal Report, Naval
%    Research Laboratory, no. NRL/FR/5344--16-10,279, 10 Aug. 2016, 173
%    pages.
%
%March 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3)
    deltaT=[];
end

if(nargin<4)
    clockLoc=[];
end

[Jul1,Jul2]=UTC2TT(Jul1,Jul2);
[Jul1,Jul2]=TT2TDB(Jul1,Jul2,deltaT,clockLoc);

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:
%SHOWPOSE Display a camera icon in 3D
%
%   h = cam.showpose(T)
%
%  Create a new camera at pose T, and return the graphics handle.
%
% The camera is depicted as a pyramid with the apex as the camera
% origin and the base plane normal to the optical axis.
% The sides are colored red, green, blue corresponding to the X, Y, Z axes
% respectively.

function h = showpose(c, T, sz)

    if nargin < 2
        T = eye(4,4);
    end
    if nargin < 3
        % get the overall scale factor from the existing graph
        sz = [get(gca, 'Xlim'); get(gca, 'Ylim'); get(gca, 'Zlim')];
        sz = max(sz(:,2)-sz(:,1));
        sz = sz / 20;
    end

    % define pyramid dimensions from the size parameter
    w = sz;
    l = sz*2;

    % define the vertices of the camera
    vertices = [
         0    0    0
         w/2  w/2  l
        -w/2  w/2  l
        -w/2 -w/2  l
         w/2 -w/2  l
        ];

    ud.vertices = vertices;

    % create the camera
    vertices = transformp(T, vertices')';

    % the first index for each face controls the face color
    faces = [
        1 2 5 NaN
        2 1 3 NaN
        3 4 1 NaN
        4 1 5 NaN
        5 2 3 4
        %2 3 4 5
        ];

    colors = [
        1 0 0       % R
        0 1 0       % G
        1 0 0       % R
        0 1 0       % G
        0 0 1       % B
        ];

    if isempty(c.camview)
        % create the camera view
        h = patch('Vertices', vertices, ...
            'Faces', faces, ...
            'FaceVertexCData', colors, ...
            'FaceColor','flat', ...
            'UserData', ud);

        set(h, 'FaceAlpha', 1.0);
        c.camview = h;
        c.camview_sz = sz;
    else
        set(c.camview, 'Vertices', vertices);
    end
end

author: petercorke

repo: machinevision-toolbox-matlab


text:
% SP_TO_VTK: Export to VTK format for plotting (store data in
%            binary base64 encoded format).
%
%  sp_to_vtk (u, space, geometry, npts, filename,
%             fieldname, [option], [precision])
%  sp_to_vtk (u, space, geometry, pts, filename,
%             fieldname, [option], [precision])
%
% INPUT:
%
%     u:          vector of dof weights
%     space:      object representing the space of discrete functions (see sp_scalar)
%     geometry:   geometry structure (see geo_load)
%     npts:       number of points along each parametric direction where to evaluate
%     pts:        cell array with the coordinates along each parametric direction of the points where to evaluate
%     filename:   name of the output file.
%     fieldname:  how to name the saved variable in the vtk file
%     option:     accepted options are 'value' (default), 'gradient',
%                 and for vectors also 'curl', 'divergence'
%     precision:  accepted options are 32 (default), or 64, other
%                 values are silently ignored
%
%
% OUTPUT:
%
%    a vtk structured mesh file named <filename> is produced, data are
%    base64-encoded and appended at the end of the file to save disk
%    space and write time.
%
% Copyright (C) 2009, 2010, 2014 Carlo de Falco
% Copyright (C) 2011, 2012 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function sp_to_vtk_raw (u, space, geometry, npts, filename, fieldname, varargin)

  if (nargin == 6)
    option = 'value';
  else
    option = varargin{1};
  end

  if (~(numel (varargin) >= 2) || ~(varargin{2} == 64))
    precision = 32;
  else
    precision = 64;
  end

  [eu, F] = sp_eval (u, space, geometry, npts, option);

  if (space.ncomp == 1 || ~strcmp(option, 'gradient'))
    msh_to_vtk_raw (F, eu, filename, fieldname, precision);
  else
    error ('sp_to_vtk: For vector fields, the gradient cannot be saved')
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function hash = xVOChash_init(strs)
% From the PASCAL VOC 2011 devkit

hsize=4999;
hash.key=cell(hsize,1);
hash.val=cell(hsize,1);

for i=1:numel(strs)
    s=strs{i};
    h=mod(str2double(s([4 6:end])),hsize)+1;
    j=numel(hash.key{h})+1;
    hash.key{h}{j}=strs{i};
    hash.val{h}(j)=i;
end

author: ShaoqingRen


function [ status ] = issues( k )
%ISSUES M2T Test cases related to issues
%
% Issue-related test cases for matlab2tikz
%
% See also: ACID, matlab2tikz_acidtest
    testfunction_handles = {
                             @scatter3Plot3
                           };

    numFunctions = length( testfunction_handles );

    if (k<=0)
        status = testfunction_handles;
        return;  % This is used for querying numFunctions.

    elseif (k<=numFunctions)
        status = testfunction_handles{k}();
        status.function = func2str(testfunction_handles{k});

    else
        error('issues:outOfBounds', ...
              'Out of bounds (number of testfunctions=%d)', numFunctions);
    end

end

% =========================================================================
function [stat] = scatter3Plot3()
    stat.description = 'Scatter3 plot with 2 colors';
    stat.issues = 292;

    hold on;
    x = sin(1:5);
    y = cos(3.4 *(1:5));
    z = x.*y;
    scatter3(x,y,z,150,...
             'MarkerEdgeColor','none','MarkerFaceColor','k');
    scatter3(-x,y,z,150,...
             'MarkerEdgeColor','none','MarkerFaceColor','b');
end

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


text:
%FFNC Feed-forward neural net classifier back-end
%
% 	[W,HIST] = FFNC (ALG,A,UNITS,ITER,W_INI,T,FID)
%
% INPUT
% 	ALG   Training algorithm: 'bpxnc' for back-propagation (default), 'lmnc'
%         for Levenberg-Marquardt
% 	A     Training dataset
% 	UNITS Array indicating number of units in each hidden layer (default: [5])
% 	ITER  Number of iterations to train (default: inf)
% 	W_INI Weight initialisation network mapping (default: [], meaning
%         initialisation by Matlab's neural network toolbox)
% 	T     Tuning set (default: [], meaning use A)
%   FID   File ID to write progress to (default [], see PRPROGRESS)
%
% OUTPUT
% 	W     Trained feed-forward neural network mapping
% 	HIST  Progress report (see below)
%
% DESCRIPTION
% This function should not be called directly, but through one of its
% front-ends, BPXNC or LMNC. Uses the Mathworks' Neural Network toolbox.
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, BPXNC, LMNC, NEURC, RNNC, RBNC, PRPROGRESS

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

% $Id: ffnc.m,v 1.10 2009/07/26 18:43:12 duin Exp $

function [w,hist] = ffnc(alg,a,units,max_iter,w_ini,tt,fid)

		% Settings for the various training algorithms.

  checktoolbox('nnet');

	if (strcmp(alg,'bpxnc'))
		mapname = 'BP-NeuralNet';
	elseif (strcmp(alg,'lmnc'))
		mapname = 'LM-NeuralNet';
	else
		error('illegal training algorithm specified');
	end;

	% Check arguments

	if (nargin < 7), fid = []; end;
	if (nargin < 6) | (isempty(tt))
		prwarning(2,'no tuning set supplied, using training set for tuning (risk of overfit)');
		if (nargin < 2), t = []; else, t = a; end;
    tt = []; % preserve input for regopt calls
  else
    t = tt;
	end
	if (nargin < 5) | (isempty(w_ini))
		prwarning(3,'no initialisation supplied, using Nguyen-Widrow random initialisation');
		w_ini = [];
	end
	if (nargin < 4) | (isempty(max_iter))
		prwarning(3,'no maximum number of iterations supplied, assuming infinite');
		max_iter = inf;
	end
	if (nargin < 3) | (isempty(units))
		prwarning(2,'no network architecture specified, assuming one hidden layer of 5 units');
		units = 5;
	end
	if (nargin < 2) | (isempty(a))
		w = prmapping(alg,{units,max_iter,w_ini,t,fid});
		w = setname(w,mapname);
		hist = [];
		return
	end

	a = setprior(a,getprior(a));
	t = setprior(a,getprior(t));

  if isnan(units) % optimize complexity parameter: number of neurons
		defs = {5,[],[],[],[]};
		parmin_max = [1,30;0,0;0,0;0,0;0,0];
		[w,hist] = regoptc(a,alg,{units,max_iter,w_ini,tt,fid},defs,[1],parmin_max,testc([],'soft'),0);
		return
  end

	% Training target values.

	prwarning (4, 'using training targets 0.9/0.1');
	target_high	= 0.9;
	target_low	= 0.1;

	% Check whether the dataset is valid.

	islabtype(a,'crisp');
	isvaldfile(a,1,2); 							% At least 1 object per class, 2 classes
	a = testdatasize(a);
	t = testdatasize(t);
	iscomdset(a,t);   							% Check whether training and tuning set match
	%a = setprior(a,getprior(a));
	%t = setprior(a,getprior(t));

	[m,k,c] = getsize(a);
	lablist = getlablist(a);

  nodipimage('set');

	% Standard training parameters.

	disp_freq   = inf;
	err_goal 		= 0.02/m;						% Mean-squared error goal, stop if reached
	trnsf_fn	  = 'logsig';					% Transfer function
  perf_fn     = 'mse';            % Performance function

	% Settings for the different training algorithms.

	tp.show   = disp_freq;
	tp.time   = inf;
	tp.goal   = err_goal;

  %DXD apparently this is now required:
  tp.showWindow = 0;
  tp.showCommandLine = 0;
	if (strcmp(alg,'bpxnc'))
		trnalg 					= 'traingdx';
		lrnalg 					= 'learngdm';
		burnin       		= 500;			% Never stop training before this many iters
		tp.epochs    		= min(50,max_iter); 	% Iteration unit
		tp.lr        		= 0.01;			% BP, initial value for adaptive learning rate
		tp.lr_inc    		= 1.05;			% BP, multiplier for increasing learning rate
		tp.lr_dec    		= 0.7;			% BP, multiplier for decreasing learning rate
		tp.mc        		= 0.95;			% BP, momentum
		tp.max_perf_inc = 1.04;			% BP, error ratio
		tp.min_grad  		= 1e-6;			% BP, minimum performance gradient
		tp.max_fail  		= 5;				% BP, maximum validation failures
		speed           = 10000;    % waitbar speed
	elseif (strcmp(alg,'lmnc'))
		trnalg 					= 'trainlm';
		lrnalg 					= 'learngdm';
		burnin       		= 50;				% Never stop training before this many iters
		tp.epochs 	 		= min(1,max_iter); 		% Iteration unit
		%tp.mem_reduc 		= 1;				% Trade-off between memory & speed
		tp.max_fail  		= 1; 				% LM, maximum validation failures
		tp.min_grad  		= 1e-6;			% LM, minimum gradient, stop if reached
		tp.mu  	    		= 0.001;		% LM, initial value for adaptive learning rate
		tp.mu_inc    		= 10;				% LM, multiplier for increasing learning rate
		tp.mu_dec    		= 0.1;			% LM, multiplier for decreasing learning rate
		tp.mu_max    		= 1e10;			% LM, maximum learning rate
		speed           = 100;      % waitbar speed
	end;

	% Scale each feature to the range [0,1].
	prwarning(3,'scaling such that training set features have range [0,1]');
	ws = scalem(a,'domain'); a_scaled = a*ws; t_scaled = t*ws;

	% Set number of network outputs: 1 for 2 classes, c for c > 2 classes.
	if (c == 2), cout = 1; else cout = c; end

	% Create target matrix: row C contains a high value at position C,
	% the correct class, and a low one for the incorrect ones (place coding).
	if (cout > 1)
    target = target_low * ones(c,c) + (target_high - target_low) * eye(c);
	else
		target = [target_high; target_low];
	end

	% Create the target arrays for both datasets.
	target_a = target(getnlab(a),:)';
	target_t = target(getnlab(t),:)';

	% Create the network layout: K inputs, N(:) hidden units, COUT outputs.
	numlayers = length(units)+1; numunits = [k,units(:)',cout];
	transfer_fn = cellstr(char(ones(numlayers,1)*trnsf_fn));

	% Create network and set training parameters. The network is initialised
	% by the Nguyen-Widrow rule by default.

	warning('OFF','NNET:Obsolete')
	net = newff(ones(numunits(1),1)*[0 1],numunits(2:end),...
							transfer_fn,trnalg,lrnalg,perf_fn);
    net.trainParam = tp;

	% If an initial network is specified, use its weights and biases.

	if (~isempty(w_ini))
		% Use given initialisation.
		[data,lab,type_w] = get(w_ini,'data','labels','mapping_file');
		if (strcmp(type_w,'sequential'))
			a_scaled = a*data{1}; t_scaled = t*data{1}; ws = data{1};
			[data,lab,type_w] = get(data{2},'data','labels','mapping_file');
		end

		% Check whether the mapping's dimensions are the same as the network's.
		[kw,cw] = size(w_ini); net_ini = data{1};
		if (~strcmp(type_w,'neurc')) | (kw ~= k) | (cw ~= c) | ...
				(net.numInputs  ~= net_ini.numInputs) | ...
				(net.numLayers  ~= net_ini.numLayers) | ...
				(net.numOutputs ~= net_ini.numOutputs) | ...
				any(net.biasConnect ~= net_ini.biasConnect) | ...
				any(net.inputConnect ~= net_ini.inputConnect) | ...
				any(net.outputConnect ~= net_ini.outputConnect)
			error('incorrect size initialisation network supplied')
		end
		% Check whether the initialisation network was trained on the same data.
		[dummy1,nlab,dummy2] = renumlab(lablist,lab);
		if (max(nlab) > c)
			error('initialisation network should be trained on same classes')
		end
		net.IW = net_ini.IW; net.LW = net_ini.LW; net.b = net_ini.b;
	end

	% Initialize loop

	opt_err = inf; opt_iter = inf; opt_net = net;
	iter = 0; this_iter = 1; hist = [];

	% Loop while
	% - training has not gone on for longer than 50 iterations or 2 times the
	%   number of iterations for which the error was minimal, and
	% - the number of iterations does not exceed the maximum
	% - the actual training function still performed some iterations

	prprogress(fid,'%s: neural net, %i units: \n',alg,units);
	prprogress(fid,'%i %5.3f %5.3f\n',0,1,1);

	s = sprintf('%s: neural net, %i units',alg,units);
	prwaitbar(100,s);
  while ((iter <= 2*opt_iter) | (iter < burnin)) & ...
				(iter < max_iter) & (this_iter > 0) & (opt_err > 0)
		prwaitbar(100,100-100*exp(-iter/speed));
		% Call TRAIN, from Matlab's NN toolbox.

		prwarning(4,'[%d] calling NNETs train', iter);
		%net.trainParam.mu = min(net.trainParam.mu_max*0.9999,net.trainParam.mu);
		%net.trainParam.mem_reduc = 1;
    %net.efficiency.memoryReduction = 1;
    if verLessThan('nnet','7.0')
      net.trainParam.mem_reduc = 1;
    else
      net.efficiency.memoryReduction = 1;
    end
		net.trainParam.showWindow = 0;
		net.trainParam.showCommandLine = 0;

		[net,tr] = train(net,+a_scaled',target_a);
		this_iter = length(tr.epoch)-1; iter = iter + this_iter;

		% Copy current learning rate as the one to start with for the next time.
		if (strcmp(alg,'bpxnc'))
			net.trainParam.lr = tr.lr(end);
		else
			net.trainParam.mu = tr.mu(end);
		end;

		% Map train and tuning set.
  	w = prmapping('neurc','trained',{net},lablist,k,c);
  	w = setname(w,mapname);

		% Calculate mean squared errors (MSE).
  	out_a = a_scaled*w; out_t = t_scaled*w;
		mse_a = mean(mean(((out_a(:,1:cout))-target_a').^2,1));
		mse_t = mean(mean(((out_t(:,1:cout))-target_t').^2,1));

		% Calculate classification errors.
		e_a = testc(a_scaled,w); e_t = testc(t_scaled,w);

		% If this error is minimal, store the iteration number and weights.
  	if (e_t < opt_err)
  		opt_err = e_t; opt_iter = iter; opt_net = net;
  	end
		w1 = cell2mat(net.IW); w1 = w1(:);
		%w2 = cell2mat(net.LW'); bugfix, doesnot work for multilayer networks
		%w2 = w2(:);
		netLW = net.LW(:);
		w2 = [];
		for j=1: length(netLW)
			ww = netLW{j};
			w2 = [w2; ww(:)];
		end
  	hist = [hist; iter e_a e_t mse_a mse_t ...
						mean([w1; w2].^2)];
		prprogress(fid,'%i %5.3f %5.3f\n',iter,e_t,opt_err);

  end
	prwaitbar(0);

	% Create mapping.

  w = ws*prmapping('neurc','trained',{opt_net},lablist,k,c);
  w = setname(w,mapname);
	w = setcost(w,a);

  nodipimage('reset');

return

%NODIPIMAGE Remove or reset DIPIMAGE path
%
%   NODIPIMAGE(FLAG)
%
% INPUT
%  FLAG  'set' or 'reset'. Removes DIPIMAGE from the path ('set')or adds
%        it again if it was removed ('reset')
%
% DESCRIPTION
% DIPIMAGE and the NNET toolbox both contain a routine MSE. If MSE is found
% in the DIPIMAGE path this path is removed ('set') or re-established.

function nodipimage(flag)

persistent DIPPATH

if strcmp(flag,'set')

  DIPPATH = fileparts(which('mse'));
  [ff,pp] = fileparts(DIPPATH);
  if strcmp(pp,'dipimage')
    rmpath(DIPPATH);
  else
    DIPPATH = [];
  end

elseif strcmp(flag,'reset')

  if ~isempty(DIPPATH)
    addpath(DIPPATH);
  end

else

  error('Unknown option')

end

author: marianux

repo: ecg-kit

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


text:
  function fatrix2_check_dim1(x, dim)
%|function fatrix2_check_dim1(x, dim)
%|
%| check that array x has dimensions "dim"
%| account for the "trailing 1" of size() for 1d column vectors

if nargin < 2, ir_usage, end

siz = size(x);
if numel(siz) == numel(dim) + 1 && siz(end) == 1
	siz = siz(1:end-1);
end
if ~isequal(siz, dim)
	fail('siz "%s" vs dim "%s" mismatch', num2str(size(x)), num2str(dim))
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
classdef MW11 < PROBLEM
% <multi> <real> <large/none> <constrained>
% Constrained benchmark MOP proposed by Ma and Wang

%------------------------------- Reference --------------------------------
% Z. Ma and Y. Wang, Evolutionary constrained multiobjective optimization:
% Test suite construction and performance comparisons. IEEE Transactions on
% Evolutionary Computation, 2019, 23(6): 972-986.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 15; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values and constraint violations
        function Population = Evaluation(obj,varargin)
            X = varargin{1};
            X = max(min(X,repmat(obj.upper,size(X,1),1)),repmat(obj.lower,size(X,1),1));
            g = 1 + sum(2*(X(:,obj.M:end) + (X(:,obj.M-1:end-1) - 0.5).^2 - 1).^2,2);
            PopObj(:,1) = g.*X(:,1)*sqrt(1.9999);
            PopObj(:,2) = g.*sqrt(2 - (PopObj(:,1)./g).^2);
            PopCon(:,1) = -(3 - PopObj(:,1).^2 - PopObj(:,2)).*(3 - 2*PopObj(:,1).^2 - PopObj(:,2));
            PopCon(:,2) = (3 - 0.625*PopObj(:,1).^2 - PopObj(:,2)).*(3 - 7*PopObj(:,1).^2 - PopObj(:,2));
            PopCon(:,3) = -(1.62 - 0.18*PopObj(:,1).^2 - PopObj(:,2)).*(1.125 - 0.125*PopObj(:,1).^2 - PopObj(:,2));
            PopCon(:,4) = (2.07 - 0.23*PopObj(:,1).^2 - PopObj(:,2)).*(0.63 - 0.07*PopObj(:,1).^2 - PopObj(:,2));
            Population  = SOLUTION(X,PopObj,PopCon,varargin{2:end});
            obj.FE      = obj.FE + length(Population);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1)  = (0:1/(N-1):1)';
            R(:,2)  = 1 - R(:,1);
            R       = R./repmat(sqrt(sum(R.^2,2)/2),1,2);
            c1      = (3 - R(:,1).^2 - R(:,2)).*(3 - 2*R(:,1).^2 - R(:,2));
            c2      = (3 - 0.625*R(:,1).^2 - R(:,2)).*(3 - 7*R(:,1).^2 - R(:,2));
            c3      = (1.62 - 0.18*R(:,1).^2 - R(:,2)).*(1.125 - 0.125*R(:,1).^2 - R(:,2));
            c4      = (2.07 - 0.23*R(:,1).^2 - R(:,2)).*(0.63 - 0.07*R(:,1).^2 - R(:,2));
            invalid = c1<0 | c2>0 | c3<0 | c4>0;
            while any(invalid)
                R(invalid,:) = R(invalid,:).*1.001;
                R(any(R>2.2,2),:) = [];
                c1      = (3 - R(:,1).^2 - R(:,2)).*(3 - 2*R(:,1).^2 - R(:,2));
                c2      = (3 - 0.625*R(:,1).^2 - R(:,2)).*(3 - 7*R(:,1).^2 - R(:,2));
                c3      = (1.62 - 0.18*R(:,1).^2 - R(:,2)).*(1.125 - 0.125*R(:,1).^2 - R(:,2));
                c4      = (2.07 - 0.23*R(:,1).^2 - R(:,2)).*(0.63 - 0.07*R(:,1).^2 - R(:,2));
                invalid = c1<0 | c2>0 | c3<0 | c4>0;
            end
            R = [R;1,1];
            R = R(NDSort(R,1)==1,:);
        end
        %% Generate the feasible region
        function R = GetPF(obj)
            [x,y] = meshgrid(linspace(0,2.1,400));
            z     = nan(size(x));
            fes1  = -(3-x.^2-y).*(3-2*x.^2-y) <= 0;
            fes2  = (3-0.625*x.^2-y).*(3-7*x.^2-y) <= 0;
            fes3  = -(1.62-0.18*x.^2-y).*(1.125-0.125*x.^2-y) <= 0;
            fes4  = (2.07-0.23*x.^2-y).*(0.63-0.07*x.^2-y) <= 0;
            z(fes1 & fes2 & fes3 & fes4 & x.^2+y.^2>=2) = 0;
            R = {x,y,z};
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
classdef InverterFactory < handle

    properties (Access = private)
        tensor
        inverter
    end

    methods (Access = public)

        function inverter = create(obj,tensor)
            obj.init(tensor);
            obj.createInverter()
            inverter = obj.getInverter();
        end
    end

    methods (Access = private)

        function init(obj,tensor)
            obj.tensor = tensor;
        end

        function createInverter(obj)
            if obj.isVoigt()
                obj.inverter = VoigtTensorInverter(obj.tensor);
            elseif obj.isFourthOrderTensor() && ~obj.isSymmetric()
                obj.inverter = TensorInverter(obj.tensor);
            elseif obj.isFourthOrderTensor() && obj.isSymmetric()
                obj.inverter = SymmetricTensorInverter(obj.tensor);
            else
             error('Not admitted object to be Inverted')
            end
        end

        function inv = getInverter(obj)
            inv = obj.inverter;
        end

        function isVoigt = isVoigt(obj)
            isVoigt   = strcmp(obj.tensor.getRepresentation(),'voigt');
        end

        function itIs = isFourthOrderTensor(obj)
            itIs   = strcmp(obj.tensor.getOrder(),'fourth');
        end

        function isSymmetric = isSymmetric(obj)
            symmetrizer = FourthOrderSymmetrizer();
            isSymmetric = symmetrizer.isSymmetric(obj.tensor.getValue());
        end

    end

end