text:
function uudecode = i4_to_uudecode ( i )

%*****************************************************************************80
%
%% I4_TO_UUDECODE returns the I-th character in the UUDECODE encoding.
%
%  Example:
%
%    I  UUDECODE
%
%    1  '`'
%    2  '!'
%    3  '"'
%   ..
%   64  '_'
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 March 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the index of the character.
%    1 <= I <= 64.
%
%    Output, character UUDECODE, the requested character.
%
  s = '`!"#$%&''()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_';

  if ( 1 <= i && i <= 64 )
    uudecode = s(i);
  else
    uudecode = ' ';
  end

  return
end

repo: jburkardt-m


text:
function [ Gnew] = gsp_graph_sparsify(G,epsilon)
%GSP_GRAPH_SPARSIFY sparsify a graph using Spielman-Srivastava algorithm
%   Usage: Gnew = gsp_graph_sparsify(G,epsilon);
%
%   Input parameters
%       G       : Graph structure or Laplacian matrix
%       epsilon : Sparsification parameter
%
%   Ouput parameters:
%       Gnew    : New sparsified graph or new laplacian
%
%   This function sparsifies a graph using Spielman-Srivastava algorithm.
%   Note that *epsilon* should be between $1/\sqrt(N)$ and $1$.
%
%   Example:::
%
%         epsilon = 0.5;
%         param.distribute = 1;
%         nnparam.k = 20;
%         param.nnparam = nnparam;
%         G = gsp_sensor(256,param);
%         G2 = gsp_graph_sparsify(G,epsilon);
%         figure(1)
%         gsp_plot_graph(G);
%         title('Original graph')
%         figure(2)
%         gsp_plot_graph(G2);
%         title('Sparsified graph')
%
%   References: spielman2011graph rudelson1999random rudelson2007sampling

% Author: David Shuman, Nathanael Perraudin
% Date  : 22 June 2014
% Testing: test_sparsify

% Test the input parameters
if isstruct(G)
    if isfield(G,'lap_type')
        if ~strcmp(G.lap_type,'combinatorial')
            error('Not implemented yet');
        end
    end
    L = G.L;
else
    L = G;
end

N=size(L,1);

if N<3
    error('GSP_GRAPH_SPARSIFY: Cannot sparsify a graph with less than 3 nodes');
end

% Epsilon should be between 1/sqrt(N) and 1, with a larger epsilon leading to a sparser graph
if ( (epsilon <= 1/sqrt(N)) || (epsilon >1) )
    error('GSP_GRAPH_SPARSIFY: Epsilon out of required range');
end

% Compute resistance distances
resistance_distances = gsp_resistance_distances(L);

% Get the weight matrix, and check if the original graph is connected
if isstruct(G)
    W = G.W;
else
    W = diag(diag(L)) - L;
end
W(W<1e-10) = 0;
W = sparse(W);
original_connected=gsp_check_connectivity_undirected(W);
if (original_connected==0)
    warning('Original graph not connected before sparsification');
end

% Initialize the probability distribution that will be used to select edges in the sparsified graph
[start_nodes,end_nodes,weights]=find(tril(W));
weights=max(0,weights);
Re=max(0,resistance_distances(sub2ind(size(resistance_distances),start_nodes,end_nodes)));
Pe=weights.*Re;
Pe=Pe/sum(Pe);

max_tries=10; % maximum number of tries to get a connected graph

for i=1:max_tries

    % Set Q
    C0=1/30; % Rudelson, 1996 Random Vectors in the Isotropic Position (too hard to figure out actual C0)
    C= 4*C0; % Rudelson and Vershynin, 2007, Thm. 3.1
    q=round(9*C^2*N*log(N)/(epsilon^2));

    % Choose random edges in the new graph according the probability
    % distribution initialized above
    results=gendist(Pe',q,1);
    spin_counts=hist(results,1:length(Pe'));
    per_spin_weights=weights./(q*Pe);

    % Tally the new weights and form the new graph
    new_weights=spin_counts'.*per_spin_weights;

    sparserW=sparse(start_nodes,end_nodes,new_weights,N,N);
    sparserW=sparserW+sparserW';
    sparserL=diag(sum(sparserW))-sparserW;
    sparserA=sign(sparserW);

    % Check if new graph is connected. If not, reduce epsilon and try again
    new_graph_connected=gsp_check_connectivity_undirected(sparserA);
    if new_graph_connected
        break;
    elseif i==max_tries
        warning('Despite attempts to reduce epsilon, sparsified graph is disconnected');
    else
        epsilon=epsilon-(epsilon-1/sqrt(N))/2;
    end
end

if isstruct(G)
    if ~G.directed
        sparserW = (sparserW + sparserW')/2;
        sparserL = (sparserL + sparserL')/2;
    end
    Gnew=G;
    Gnew.W = sparserW;
    Gnew.L = sparserL;
    Gnew.A = sparserA;
    Gnew.d = diag(sparserL);
    Gnew.Ne = nnz(Gnew.W)/2;
else
    Gnew = sparse(sparserL);
end

end

repo: gspbox


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function call_price_fft = CallPricingFFT(model,n,S,K,T,r,d,varargin)

lnS = log(S);
lnK = log(K);

%optAlpha = optimalAlpha(model,lnS,lnK,T,r,d,varargin{:});
optAlpha = .75;

DiscountFactor = exp(-r*T);

%-------------------------
%--- FFT Option Pricing --
%-------------------------
% from: Option Valuation Using the Fast Fourier Transform,
%       Peter Carr, March 1999, pp 10-11
%-------------------------

% predefined parameters
FFT_N = 2^n;                               % must be a power of two (2^14)
FFT_eta = 0.05;                             % spacing of psi integrand

% effective upper limit for integration (18)
% uplim = FFT_N * FFT_eta;

FFT_lambda = (2 * pi) / (FFT_N * FFT_eta);  %spacing for log strike output (23)
FFT_b = (FFT_N * FFT_lambda) / 2;           % (20)

uvec = 1:FFT_N;
%log strike levels ranging from lnS-b to lnS+b
ku = - FFT_b + FFT_lambda * (uvec - 1);     %(19)

jvec = 1:FFT_N;
vj = (jvec-1) * FFT_eta;

% optimal alpha illustration (payoff independent)
% alpharange = -3:0.1:8;
% resultrangef = zeros(1,length(alpharange));
% resultrangef1 = zeros(1,length(alpharange));
% eps = 0.000001;
% for n = 1:length(alpharange)
%     resultrangef(n)= (-alpharange(n) * log(K) + log(psialpha(model,alpharange(n),lnS,T,r,d,varargin{:})));
%     resultrangef1(n)= (-(alpharange(n) + eps) * log(K) ...
%                         + log(psialpha(model,alpharange(n)+eps,lnS,T,r,d,varargin{:})));
%     resultrangef1(n)= (resultrangef1(n) - resultrangef(n)) / eps;
% end
% plot(alpharange,resultrangef); hold on; plot(alpharange,resultrangef1, 'g'); hold off;

%applying FFT
tmp = DiscountFactor * psi(model,vj,optAlpha,lnS,T,r,d,varargin{:}) .* exp(1i * vj * (FFT_b)) * FFT_eta;
tmp = (tmp / 3) .* (3 + (-1).^jvec - ((jvec - 1) == 0) );   %applying simpson's rule
cpvec = real(exp(-optAlpha .* ku) .* fft(tmp) / pi);        %call price vector resulting in equation 24

indexOfStrike = floor((lnK + FFT_b)/FFT_lambda + 1);
iset = max(indexOfStrike)+1:-1:min(indexOfStrike)-1;
xp = ku(iset);
yp = cpvec(iset);
call_price_fft = real(interp1(xp,yp,lnK));

end

%analytical formula for zhi in equation ( 6 ) of Madan's paper
function ret = psi(model,v,alpha,varargin)
  ret = exp(feval(@CharacteristicFunctionLib, model, v - (alpha + 1) * 1i,varargin{:})) ./ (alpha.^2 + alpha - v.^2 + 1i * (2 * alpha + 1) .* v);
end

% function ret = psialpha(model,alpha,varargin)
%   ret = exp(feval(@CharacteristicFunctionLib, model, - (alpha + 1) * 1i,varargin{:}))./ (alpha.^2 + alpha);
% end

author: Sable

repo: mcbench-benchmarks


text:
function gen_mask_patch_cat_idx_for_super_res(patchified_mask)
    persistent mask_idx;
    if(isempty(mask_idx))
        global mem;
        [v, idx] = max(patchified_mask);
        idx_all = unique(idx);
        mask_idx = {};
        for m = 1:length(idx_all)
            mask_idx{m} = find(idx ==idx_all(m));
        end
        mem.mask_idx = mask_idx;
    end
end

author: jimmy-ren


text:
function [theta, J_history] = gradientDescentMulti(X, y, theta, alpha, num_iters)
%GRADIENTDESCENTMULTI Performs gradient descent to learn theta
%   theta = GRADIENTDESCENTMULTI(x, y, theta, alpha, num_iters) updates theta by
%   taking num_iters gradient steps with learning rate alpha

% Initialize some useful values
m = length(y); % number of training examples
J_history = zeros(num_iters, 1);

for iter = 1:num_iters

    % ====================== YOUR CODE HERE ======================
    % Instructions: Perform a single gradient step on the parameter vector
    %               theta.
    %
    % Hint: While debugging, it can be useful to print out the values
    %       of the cost function (computeCostMulti) and gradient here.
    %

    predictions =  X * theta;
    updates = X' * (predictions - y);
    theta = theta - alpha * (1/m) * updates;

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

    % Save the cost J in every iteration
    J_history(iter) = computeCostMulti(X, y, theta);

end

end

author: vkosuri

repo: CourseraMachineLearning


text:
function varargout = gradientSpect(f, dn, dim, deriv_order)
%GRADIENTSPECT Calculate the gradient using a Fourier spectral method.
%
% DESCRIPTION:
%       gradientSpect calculates the gradient of an n-dimensional input
%       matrix using the Fourier collocation spectral method. The gradient
%       for singleton dimensions is returned as 0.
%
%       Examples:
%       1. 1D Sinusoid:
%
%           % compute gradient of a 2 period sinusoid
%           x = pi/20:pi/20:4*pi;
%           y = sin(x);
%           dydx = gradientSpect(y, pi/20);
%
%           % plot gradient and error compared to analytical solution
%           subplot(2, 1, 1), plot(x, cos(x), 'k-', x, dydx, 'bx');
%           axis tight;
%           title('dy/dx');
%           subplot(2, 1, 2), plot(x, cos(x) - dydx, 'k-');
%           axis tight;
%           title('Relative Error');
%
%       2. Modified 2D Mathworks gradient example (x and y reversed):
%
%           [x, y] = meshgrid(-2:.2:2, -2:.2:2);
%           z = x .* exp(-x.^2 - y.^2);
%           [px,py] = gradientSpect(z, [.2, .2]);
%           contour(z), hold on, quiver(py, px), hold off;
%
% USAGE:
%       fx = gradientSpect(f, dx)
%       fx = gradientSpect(f, dx, [], deriv_order)
%       fn = gradientSpect(f, dn, dim)
%       fn = gradientSpect(f, dn, dim, deriv_order)
%       [fx, fy] = gradientSpect(f, dn)
%       [fx, fy] = gradientSpect(f, dn, [], deriv_order)
%       [fx, fy, fz, ...] = gradientSpect(f, dn)
%       [fx, fy, fz, ...] = gradientSpect(f, dn, [], deriv_order)
%
% INPUTS:
%       f           - matrix or vector to find the gradient of
%       dn          - array of values for the grid point spacing in each
%                     dimension. If a value for dim is given, dn is the
%                     spacing in dimension dim.
%
% OPTIONAL INPUTS:
%       dim         - optional input to specify a single dimension over
%                     which to compute the gradient for n-dimension input
%                     functions
%       deriv_order - order of the derivative to compute, e.g., use 1 to
%                     compute df/dx, 2 to compute df^2/dx^2, etc.
%                     (default = 1)
%
%
% OUTPUTS:
%       fx, fy, ... - gradient in the each dimension, where x corresponds
%                     to dim = 1, y corresponds to dim = 2 etc
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 25th January 2012
%       last update - 27th August 2012
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
%
% See also gradient, gradientFD

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

% check for user input for deriv_order
if nargin < 4  || isempty(deriv_order)
    deriv_order = 1;
elseif ~((deriv_order > 0) && (deriv_order == round(deriv_order)))
    error('input for deriv_order must be an integer > 0');
end

% check for user input for dim
if (nargin < 3) || isempty(dim)
    dim = 0;
end

% get size of the input function
sz = size(f);

% check if input is 1D or user defined input dimension is given
if dim || (nargout == 1 && (max(sz) == prod(sz)))

    % check if a single dn value is given, if not, extract the required
    % value
    if numel(dn) ~= 1
        dn = dn(dim);
    end

    % get the grid size along the specified dimension, or the longest
    % dimension if 1D
    if nargout == 1 && (max(sz) == prod(sz))
        [Nx, dim] = max(sz);
    else
        Nx = sz(dim);
    end

    % get the wavenumbers
    kx = getWavenumbers(Nx, dn, dim);

    % calculate derivative and assign output
    varargout{1} = ifft( bsxfun(@times, (1i*kx).^deriv_order, fft(f, [], dim)) , [], dim, 'symmetric');

else

    % check for the correct number of dn values
    if length(dn) ~= length(sz)
        error([num2str(length(sz)) ' values for dn must be specified for a ' num2str(length(sz)) '-dimensional input matrix']);
    end

    % check for the correct number of outputs
    if nargout ~= length(sz)
        error(['Incorrect number of output arguments for ' num2str(length(sz)) '-dimensional input matrix']);
    end

    % calculate the gradient for each non-singleton dimension
    for dim = 1:ndims(f)
        if sz(dim) > 1

            % get the wavenumbers
            kx = getWavenumbers(sz(dim), dn(dim), dim);

            % calculate derivative and assign output
            varargout{dim} = ifft( bsxfun(@times, (1i*kx).^deriv_order, fft(f, [], dim)) , [], dim, 'symmetric');

        else
            % assign the derivate for singleton dimensions to be 0
            varargout{dim} = 0;
        end
    end
end

function kx = getWavenumbers(Nx, dx, dim)
% subfuction to compute the wavenumber vector

% calculate the wavenumbers
if rem(Nx, 2) == 0
    % grid dimension has an even number of points
    nx = ((-Nx/2:Nx/2-1)/Nx).';
else
    % grid dimension has an odd number of points
    nx = ((-(Nx-1)/2:(Nx-1)/2)/Nx).';
end
kx = ifftshift((2*pi/dx).*nx);

% permute to be in the correction direction for use with bsxfun
if dim == 1
    kx = reshape(kx, Nx, 1);
else
    kx = reshape(kx, [ones(1, dim - 1), Nx]);
end

repo: Aero-matlab


text:
function [ inlierPointsLeft, inlierPointsRight ] = extractBinnedFeatures( viLeftImage, viRightImage )

inlierPointsLeft = [];
inlierPointsRight = [];
       % Binning

        uBin = 1:floor(size(viLeftImage,2)/6):size(viLeftImage,2);
        vBin = 1:floor(size(viLeftImage,1)/2):size(viLeftImage,1);
        uBinSize = diff([uBin,size(viLeftImage,2)]);
        vBinSize = diff([vBin,size(viLeftImage,1)]);
        [UBIN, VBIN] = meshgrid(uBin, vBin);
        [UBINSIZE, VBINSIZE] = meshgrid(uBinSize, vBinSize);
        UBIN = UBIN(:);
        VBIN = VBIN(:);
        UBINSIZE = UBINSIZE(:);
        VBINSIZE = VBINSIZE(:);

        for b = 1:size(UBIN,1)
            roiVec = [UBIN(b), VBIN(b), UBINSIZE(b), VBINSIZE(b)];

            %Detect strongest corners
            leftPoints = detectSURFFeatures(viLeftImage,'ROI',roiVec);
            rightPoints = detectSURFFeatures(viRightImage,'ROI',roiVec);

            %leftPoints = leftPoints.selectStrongest(50);
            %rightPoints = rightPoints.selectStrongest(50);

            %Extract features and stereo match
           [featuresLeft, validLeftPoints] = extractFeatures(viLeftImage, leftPoints);
           [featuresRight, validRightPoints] = extractFeatures(viRightImage, rightPoints);

            indexPairs = matchFeatures(featuresLeft, featuresRight);
            matchedPointsLeft = validLeftPoints(indexPairs(:, 1), :);
            matchedPointsRight = validRightPoints(indexPairs(:, 2), :);

            inliers = abs((matchedPointsLeft.Location(:, 2) - matchedPointsRight.Location(:, 2))) <= 1 & abs((matchedPointsLeft.Location(:, 1) - matchedPointsRight.Location(:, 1))) > 3;

            inlierPointsLeft = [inlierPointsLeft; matchedPointsLeft(inliers).Location];
            inlierPointsRight = [inlierPointsRight; matchedPointsRight(inliers).Location];
        end

end

repo: MSCKF


text:
function [ X ] = LineFaceIntersection( faceX, faceN, lineX, lineD )
%LINEFACEINTERSECTION Intersection of a plane and a line
%   faceX: a point on the plane
%   faceN: the normal direction of the plane
%   lineX: a point on the line
%   lineD: direction of the line
%

% A = dot(faceN,lineD,2)/lineD(1);
% B = -dot(faceX,faceN,2) + (dot(lineX(2:3),faceN(2:3),2)) - lineX(1)/lineD(1)*(dot(lineD(2:3),faceN(2:3),2));
A = sum(faceN.*lineD)/lineD(1);
B = -sum(faceX.*faceN) + (sum(lineX(2:3).*faceN(2:3),2)) - lineX(1)/lineD(1)*(sum(lineD(2:3).*faceN(2:3),2));

x = -B/A;

y = (x-lineX(1))/lineD(1)*lineD(2) + lineX(2);
z = (x-lineX(1))/lineD(1)*lineD(3) + lineX(3);
X = [x y z];
end

author: zouchuhang

repo: LayoutNet


text:
function [transprob, termprob] = remove_hhmm_end_state(A)
% REMOVE_END_STATE Infer transition and termination probabilities from automaton with an end state
% [transprob, termprob] = remove_end_state(A)
%
% A(i,k,j) = Pr( i->j | Qps=k), where i in 1:Q, j in 1:(Q+1), and Q+1 is the end state
% This implements the equation in footnote 3 of my NIPS 01 paper,
% transprob(i,k,j) = \tilde{A}_k(i,j)
% termprob(k,j) = \tau_k(j)
%
% For the top level, the k index is missing.

Q = size(A,1);
toplevel = (ndims(A)==2);
if toplevel
  Qk = 1;
  A = reshape(A, [Q 1 Q+1]);
else
  Qk = size(A, 2);
end

transprob = A(:, :, 1:Q);
term = A(:,:,Q+1)'; % term(k,j) = P(Qj -> end | k)
termprob = term;
%termprob = zeros(Qk, Q, 2);
%termprob(:,:,2) = term;
%termprob(:,:,1) = 1-term;

for k=1:Qk
  for i=1:Q
    for j=1:Q
      denom = (1-termprob(k,i));
      denom = denom + (denom==0)*eps;
      transprob(i,k,j) = transprob(i,k,j) / denom;
    end
  end
end

if toplevel
  termprob = squeeze(termprob);
  transprob = squeeze(transprob);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function cerr = c_err(pk, pc, t, ts, tt, pp)
    %c_err calculate c error
    qsum=pk*((1/(tt+pc)^pp)-(1/(ts+pc)^pp));
    psum=sum(1./(t+pc));
    cerr=qsum+pp*psum;
end

author: CelsoReyes


text:
classdef TestGetFontScaleFromHeight
    %TestGetFontScaleFromHeight

    methods (Static)
        function test_simple
            x = cv.getFontScaleFromHeight(20);
            validateattributes(x, {'numeric'}, {'scalar', 'real'});
        end

        function test_options
            x = cv.getFontScaleFromHeight(20, ...
                'FontFace','HersheyComplex', 'FontStyle','Italic', 'Thickness',2);
            validateattributes(x, {'numeric'}, {'scalar', 'real'});
        end

        function test_error_argnum
            try
                cv.getFontScaleFromHeight();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


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

%*****************************************************************************80
%
%% HYPER_1F1_VALUES returns some values of the hypergeometric function 1F1.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      fx = Hypergeometric1F1 [ a, b, x ]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 March 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Cambridge University Press, 1999,
%    ISBN: 0-521-64314-7,
%    LC: QA76.95.W65.
%
%    Daniel Zwillinger, editor,
%    CRC Standard Mathematical Tables and Formulae,
%    30th Edition,
%    CRC Press, 1996,
%    ISBN: 0-8493-2479-3,
%    LC: QA47.M315.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0
%    before the first call.  On each call, the routine increments N_DATA by 1,
%    and returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real A, B, X, the parameters.
%
%    Output, real FX, the value of the function.
%
  n_max = 24;

  a_vec = [ ...
    -2.500, ...
    -0.500, ...
     0.500, ...
     2.500, ...
    -2.500, ...
    -0.500, ...
     0.500, ...
     2.500, ...
    -2.500, ...
    -0.500, ...
     0.500, ...
     2.500, ...
     0.825, ...
     1.100, ...
     1.650, ...
     3.300, ...
     0.825, ...
     1.100, ...
     1.650, ...
     3.300, ...
     0.825, ...
     1.100, ...
     1.650,...
     3.300 ];
  b_vec = [ ...
     3.3, ...
     1.1, ...
     1.1, ...
     3.3, ...
     3.3, ...
     1.1, ...
     1.1, ...
     3.3, ...
     3.3, ...
     1.1, ...
     1.1, ...
     3.3, ...
     6.7, ...
     6.7, ...
     6.7, ...
     6.7, ...
     6.7, ...
     6.7, ...
     6.7, ...
     6.7, ...
     6.7, ...
     6.7, ...
     6.7, ...
     6.7 ];
  fx_vec = [ ...
     0.81879926689265186854, ...
     0.88283984828032972070, ...
     1.1245023764952626690, ...
     1.2101049301639599598, ...
     0.12723045536781567174, ...
     0.12326016871544045107, ...
     2.3297954665128293051, ...
     3.3890020264468009733, ...
    -0.18819510282516768874, ...
    -1.0764203806547022727, ...
     5.7521824680907968433, ...
     9.9998567403304086593, ...
     1.0317208964319891384, ...
     1.0424867029249952040, ...
     1.0643112000949092012, ...
     1.1321844369742336326, ...
     1.2328402688568452181, ...
     1.3200654482027340732, ...
     1.5104811522310825217, ...
     2.2307520785940524365, ...
     1.5197286298183137741, ...
     1.7364938170250847619, ...
     2.2492330307668135926, ...
     4.6377737119178965298 ];
  x_vec = [ ...
     0.25, ...
     0.25, ...
     0.25, ...
     0.25, ...
     1.55, ...
     1.55, ...
     1.55, ...
     1.55, ...
     2.85, ...
     2.85, ...
     2.85, ...
     2.85, ...
     0.25, ...
     0.25, ...
     0.25, ...
     0.25, ...
     1.55, ...
     1.55, ...
     1.55, ...
     1.55, ...
     2.85, ...
     2.85, ...
     2.85, ...
     2.85 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

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

  return
end

author: johannesgerer


text:
function printTensor(tensor,order)
%-fanDTasia ToolBox------------------------------------------------------------------
% This Matlab script is part of the fanDTasia ToolBox: a Matlab library for Diffusion
% Weighted MRI (DW-MRI) Processing, Diffusion Tensor (DTI) Estimation, High-order
% Diffusion Tensor Analysis, Tensor ODF estimation, Visualization and more.
%
% A Matlab Tutorial on DW-MRI can be found in:
% http://www.cise.ufl.edu/~abarmpou/lab/fanDTasia/tutorial.php
%
%-CITATION---------------------------------------------------------------------------
% If you use this software please cite the following work:
% A. Barmpoutis and B.C. Vemuri, "A Unified Framework for Estimating Diffusion Tensors
% of any order with Symmetric Positive-Definite Constraints",
% In the Proceedings of ISBI, 2010
%
%-DESCRIPTION------------------------------------------------------------------------
% This function prints in the command line the tensor coefficients of a given tensor
% (in 3 variables) of a specific order.
%
%-USE--------------------------------------------------------------------------------
% printTensor(tensor,order);
%
% tensor: is a vector that contains all the unique coefficients of a symmetric tensor
% order: is the order of the tensor
%
%-DISCLAIMER-------------------------------------------------------------------------
% You can use this source code for non commercial research and educational purposes
% only without licensing fees and is provided without guarantee or warrantee expressed
% or implied. You cannot repost this file without prior written permission from the
% authors. If you use this software please cite the following work:
% A. Barmpoutis and B.C. Vemuri, "A Unified Framework for Estimating Diffusion Tensors
% of any order with Symmetric Positive-Definite Constraints", In Proc. of ISBI, 2010.
%
%-AUTHOR-----------------------------------------------------------------------------
% Angelos Barmpoutis, PhD
% Computer and Information Science and Engineering Department
% University of Florida, Gainesville, FL 32611, USA
% abarmpou at cise dot ufl dot edu
%------------------------------------------------------------------------------------

    c=1;
    for i=0:order
		for j=0:order-i
			fprintf(1,'D%d%d%d = %.8f\n',i,j,order-i-j,tensor(c));
			c=c+1;
        end
    end


text:
function [modifiedMap, xmlMap, model2] = setColorToMapRxns (mapFileName, model, v1, rxnsFluxTask, rxnLineWidth, rxnLineColor, molFillColor)
% USAGE:
% [modifiedMap, xmlMap, model2] =
%   setColorToMapRxns (mapFileName, model, v1, rxnsFluxTask, rxnLineWidth, rxnLineColor, molFillColor)
%
% INPUTS:
% mapFileName           relative path to XML map (e.g. 'yeastMap1.0.xml')
% model                 model structure
% v1                    vector of fluxes
% rxnsFluxTask          1 = Flux
%                       2 = Flux essentiality; map reactions causing any
%                           effect on growth (i.e. obj function)
%                       3 = Flux essentiality; map reactions causing no
%                           growth(i.e. obj function) when deleted, i.e.
%                           essential reactions
%                       4 = Flux essentiality; map reactions with not any
%                           effect on growth (i.e. obj function) when
%                           deleted
%                       5 = Omics
% ADDITIONAL INPUTS:
% rxnLineWidth          reaction line width; default = 5
% rxnLineColor          reaction line color; default = Pathways; write
%                       'Subsystems' to color five subsystems only (central,
%                       AA, lipids, cofactors, nuclotides and energy).
%                       'Subsystems' option is the deault option at task 5
% molFillColor          molecule fill color; default = true (i.e. white)
%                       molecule fill color; true (i.e. like CellDesigner default)
% OUTPUTS:
% modifiedMap           map with MATLAB structure, and modified reactions' line color and width
% xmlMap                xml structure of the original map obtained from mapFileName
% model2                model with modified subsystems

if ~(exist('transformXML2Map')==2)
    error('COBRA Toolbox is required, go to https://opencobra.github.io/cobratoolbox/')
end

if nargin < 5 || rxnLineWidth == 0
    rxnLineWidth = 0;
    rxnLineColor = '';
    molFillColor = true;
elseif nargin < 6
    rxnLineColor = '';
    molFillColor = true;
elseif nargin < 7
    molFillColor = true;
end
if rxnsFluxTask == 5
    rxnLineColor = 'Subsystems';
end

[xmlMap, mapMap] = transformXML2Map(mapFileName);
if molFillColor == true
    modifiedMap = unifyMetabolicMapCD(mapMap); % molColor whrite and rxnColor = ligthgray
elseif molFillColor == false
    modifiedMap = mapMap; % like CellDesigner default
end

%%%%%%%% Gettin reactions line width with normalized fluxes %%%%%%%%%%%%%
model2 = model; % to not overwrite model
flxs(:, 1) = v1; % to not overwrite v1
absFlxs = abs(flxs);
rxnLineWidthInMap = absFlxs / max(absFlxs); %normalize fluxes line width
if rxnsFluxTask == 1 && rxnLineWidth == 0
    rxnLineWidthInMap(rxnLineWidthInMap > 0.1 & rxnLineWidthInMap <=  1) = 50;
    rxnLineWidthInMap(rxnLineWidthInMap > 0.01 & rxnLineWidthInMap <= 0.1) = 25;
    rxnLineWidthInMap(rxnLineWidthInMap > 0.001 & rxnLineWidthInMap <= 0.01) = 16;
    rxnLineWidthInMap(rxnLineWidthInMap > 0.0001 & rxnLineWidthInMap <= 0.001) = 10;
    rxnLineWidthInMap(rxnLineWidthInMap > 0.00001 & rxnLineWidthInMap <= 0.0001) = 6;
elseif rxnsFluxTask == 2 && rxnLineWidth == 0
    rxnLineWidthInMap(rxnLineWidthInMap > 0.90 & rxnLineWidthInMap <=  0.95) =  6;
    rxnLineWidthInMap(rxnLineWidthInMap > 0.80 & rxnLineWidthInMap <= 0.90) =  8;
    rxnLineWidthInMap(rxnLineWidthInMap > 0.50 & rxnLineWidthInMap <= 0.80) = 11;
    rxnLineWidthInMap(rxnLineWidthInMap > 0.25 & rxnLineWidthInMap <= 0.50) = 15;
    rxnLineWidthInMap(rxnLineWidthInMap >= 0.00 & rxnLineWidthInMap <= 0.25) = 20;
elseif rxnsFluxTask == 1 && rxnLineWidth ~= 0.0
    rxnLineWidthInMap(rxnLineWidthInMap == 1.0) = rxnLineWidth;
elseif rxnsFluxTask == 2 && rxnLineWidth ~= 0.0
    rxnLineWidthInMap(rxnLineWidthInMap < 0.95) = rxnLineWidth;
elseif rxnsFluxTask == 3 && rxnLineWidth == 0.0
    rxnLineWidthInMap(rxnLineWidthInMap == 0.0) = 5;
elseif rxnsFluxTask == 3 && rxnLineWidth ~= 0.0
    rxnLineWidthInMap(rxnLineWidthInMap == 0.0) = rxnLineWidth;
elseif rxnsFluxTask == 4 && rxnLineWidth == 0.0
    rxnLineWidthInMap(rxnLineWidthInMap > 0.95 & rxnLineWidthInMap <= 1.0) = 5.0;
elseif rxnsFluxTask == 4 && rxnLineWidth ~= 0.0
    rxnLineWidthInMap(rxnLineWidthInMap > 0.95 & rxnLineWidthInMap <= 1.0) = rxnLineWidth;
elseif rxnsFluxTask == 5 && rxnLineWidth == 0.0
    cutValue = 2.0; % transcriptomics minimum cut value of log change
    rxnLineWidthInMap = absFlxs;
    rxnLineWidthInMap(rxnLineWidthInMap > 5.0) = 30;
    rxnLineWidthInMap(rxnLineWidthInMap > 4.0 & rxnLineWidthInMap <= 5.0) = 25;
    rxnLineWidthInMap(rxnLineWidthInMap > 3.0 & rxnLineWidthInMap <= 4.0) = 20;
    rxnLineWidthInMap(rxnLineWidthInMap > 2.0 & rxnLineWidthInMap <= 3.0) = 15;
    rxnLineWidthInMap(rxnLineWidthInMap > 1.0 & rxnLineWidthInMap <= 2.0) = 10;
    rxnLineWidthInMap(rxnLineWidthInMap >= 0.5 & rxnLineWidthInMap <= 1.0) =  5;
end
for i = 1:length(modifiedMap.rxnName)
    a = find(ismember(model2.rxns, modifiedMap.rxnName{i}));
    if isempty(a)
        modifiedMap.rxnWidth{i} = 1;
    else
        modifiedMap.rxnWidth{i} = rxnLineWidthInMap(a);
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%% Cofactors not included in the analyses %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% Compartments included in the analyses %%%%%%%%%%%%%%%%%%%%%%%%%

% propossed cofactors list for PART A
cofactorsListA ={'h'; 'h2o'; 'nad'; 'nadh'; 'nadp'; 'nadph';...
    'fad'; 'fadh'; 'ppi'; 'adp'; 'gtp'; 'gdp'; 'gmp'; ...
    'utp'; 'udp'; 'ump'; 'ctp'; 'cdp'; 'cmp'};

% propossed cofactors list for PART B
cofactorsListB ={'h'; 'h2o'; 'co2'; 'coa'; 'accoa'; 'nad'; 'nadh'; 'nadp'; 'nadph';...
    'fad'; 'fadh'; 'pi'; 'ppi'; 'atp'; 'adp'; 'amp'; 'gtp'; 'gdp'; 'gmp';...
    'utp'; 'udp'; 'ump'; 'ctp'; 'cdp'; 'cmp'; 'gthox'; 'gthrd_c'; 'amet_c';...
    'ahcys_c'; 'glu__L'; 'gln__L'; 'akg'; 'nh4'};

% compartments list
compartmentsList = {'_e'; '_c'; '_m'; '_n'; '_x'; '_v'; '_r'; '_g'};

% include compartment to cofactor lists A and B
cofactorsListWcompA = cell(length(cofactorsListA), 1);
c = 1; % general counter
for i = 1:length(cofactorsListA)
    for j = 1:length(compartmentsList)
        cofactorsListWcompA{c, 1} = horzcat(cofactorsListA{i},compartmentsList{j});
        c = c + 1;
    end
end
cofactorsListWcompB = cell(length(cofactorsListB), 1);
c = 1;
for i = 1:length(cofactorsListB)
    for j = 1:length(compartmentsList)
        cofactorsListWcompB{c, 1}= horzcat(cofactorsListB{i},compartmentsList{j});
        c = c + 1;
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%% PART A %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% getting transport reactions
[nrxns,~] = size(model2.rxns);
transRxns = cell(nrxns, 1);
for j = 1:size(model2.rxns)
    if true(contains(model2.subSystems {j}, 'S_Transport'))
        transRxns{j, 1} = model2.rxns{j};
    end
end
transRxns = transRxns(~cellfun('isempty',transRxns));

% getting subStystems without transport
SubSystems = (unique(getModelSubSystems(model2)));
[nSubSystems, ~] = size(SubSystems);
subSystemsNoTrans = cell(nSubSystems, 1);
for i = 1:nSubSystems
    if ~contains(SubSystems{i,1},'S_Transport')&& ...
            ~contains(SubSystems{i, 1},'Extracellular exchange') &&...
            ~contains(SubSystems{i, 1},'S_tRNA_charging')
        subSystemsNoTrans{i, 1} = SubSystems{i, 1};
    else
        subSystemsNoTrans{i, 1} = {};
    end
end
subSystemsNoTrans = subSystemsNoTrans(~cellfun('isempty',subSystemsNoTrans));

% defining subSystems of trasnport reactions to be colored, these
% transport reactions complete reaction conections in subSystems which
% require conexions among diferents compartments
[nTransRxns,~] = size(transRxns);
[nSubSystemsNoTrans, ~] = size (subSystemsNoTrans);
totalRxnsToIncludeA = cell(nTransRxns,1);
for i = 1:nSubSystemsNoTrans
    subSystemRxns = findRxnsFromSubSystem(model2, subSystemsNoTrans(i));
    metsSubSystemRxns = findMetsFromRxns(model2, subSystemRxns);
    [nMetsSubSystemRxns, ~] = size (metsSubSystemRxns);
    % getting unique compartments
    metsCompartment = cell(nMetsSubSystemRxns,1);
    for j = 1:nMetsSubSystemRxns
        str = metsSubSystemRxns{j};
        [~, s] = size(str);
        metsCompartment{j,1} = str(s-1:s);
    end
    metsCompartment = (unique(metsCompartment));
    [nUniqueMetsComp, ~] = size(metsCompartment);

    if nUniqueMetsComp > 1
        % getting duplicated mets like 'pyr_c' and 'pyr_m' which are 'pyr'
        % and pyr after erasing the compartment signature
        metsFromSubSystemRxnsNoComp = cell(nMetsSubSystemRxns,1);
        for j = 1:nMetsSubSystemRxns
            str = metsSubSystemRxns{j};
            [z, s] = size(str);
            metsFromSubSystemRxnsNoComp{j, 1} = str(1:s-2); % erasing comp signature
        end
        [~, o, ~] = unique(metsFromSubSystemRxnsNoComp,'first');
        indexToDupes = find(not(ismember(1:numel(metsFromSubSystemRxnsNoComp),o)));
        [~, nIndexToDupes] = size(indexToDupes);
        duplicatedMets = cell(nIndexToDupes,1);
        for j = 1:nIndexToDupes
            duplicatedMets(j, 1) = metsFromSubSystemRxnsNoComp (indexToDupes(1, j));
        end

        % getting reactions from duplicated mets
        dulicatedRxnMets = cell(nIndexToDupes * nUniqueMetsComp,1);
        c = 1; % duplicatedRxnMets counter
        for j = 1:nIndexToDupes
            str = duplicatedMets{j, 1};
            if ~ismember(str, cofactorsListA)
                for k = 1:nUniqueMetsComp
                    str1 = metsCompartment{k};
                    if ~ismember(str1, compartmentsList(2)) % compartment '_c'
                        dulicatedRxnMets{c, 1} = {horzcat(str, '_c'), horzcat(str, str1)};
                        c = c + 1;
                    end
                end
            end
        end
        dulicatedRxnMets = dulicatedRxnMets(~cellfun('isempty',dulicatedRxnMets));
        [ndulicatedRxnMets, ~] = size(dulicatedRxnMets);
        % from duplicated reactions, get reactions to subscribe a reactions
        % subsystem from non transport subsystems
        rxnsToInclude = cell(ndulicatedRxnMets, 1);
        for j = 1:ndulicatedRxnMets
            metsToFindInTrasnportRxns = dulicatedRxnMets{j, 1};
            for k = 1:nTransRxns
                rxnFormula = printRxnFormula(model2, transRxns{k, 1}, false);
                if true(contains(string(rxnFormula), metsToFindInTrasnportRxns (1, 1)))
                    if true(contains(rxnFormula, metsToFindInTrasnportRxns (1, 2)))
                        rxnsToInclude{j, 1} = transRxns{k, 1};
                    end
                end
            end
        end
        rxnsToInclude = rxnsToInclude(~cellfun('isempty',rxnsToInclude));
        rxnsToInclude = unique(rxnsToInclude);
        [s, ~] = size(rxnsToInclude);
        totalRxnsToIncludeA = [totalRxnsToIncludeA; rxnsToInclude];

        for j = 1:s
            rxnID = findRxnIDs(model2, rxnsToInclude{j});
            model2.subSystems(rxnID) = {subSystemsNoTrans(i)};
        end

    end
end
totalRxnsToIncludeA = totalRxnsToIncludeA(~cellfun('isempty',totalRxnsToIncludeA));
totalRxnsToIncludeA = unique(totalRxnsToIncludeA);

%%%%%%%%%%%%%%%%%%%%%%%%%%%% PART B %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Get subsystems category. This is based on the number of conexions that a
% subsystem includes in the genome scale metabolic model2
subsLength ={};
for i = 1:length(subSystemsNoTrans)
    rxnsSubs = findRxnsFromSubSystem(model2,subSystemsNoTrans(i));
    metsSubsRxns = findMetsFromRxns(model2, rxnsSubs);
    metsSubsRxnsNoCof = {};
    for j = 1:length(metsSubsRxns)
        if ~ismember(metsSubsRxns(j),cofactorsListWcompB)
            metsSubsRxnsNoCof{j,1} = char(metsSubsRxns(j));
        end
    end
    metsSubsRxnsNoCof = metsSubsRxnsNoCof(~cellfun('isempty',metsSubsRxnsNoCof));

    metsRxns = findRxnsFromMets(model2, metsSubsRxnsNoCof);
    metsRxnsSubs = {};
    for j = 1:length(metsRxns)
        rxnID = findRxnIDs(model2,metsRxns(j));
        metsRxnsSubs{j,1} = string(model2.subSystems(rxnID));
    end
    metsRxnsSubs = unique(string(metsRxnsSubs));
    subsLength{i,1} = length(metsRxnsSubs);
    subsLength{i,2} = subSystemsNoTrans{i};
end

% Get  remaining transport reacctions, namely, reactions not included in
% PART A

[nTransRxns, ~] = size(transRxns);
remainingTransRxns = cell(nTransRxns, 1);
for i = 1:nTransRxns
    if ~ismember(transRxns(i),totalRxnsToIncludeA)
        remainingTransRxns{i,1} = transRxns{i, 1};
    end
end
remainingTransRxns = remainingTransRxns(~cellfun('isempty',remainingTransRxns));

% Assign remaining transport reactions
nonIncludSubs = {'S_tRNA_charging', 'S_Other',...
    'Biomass and maintenance functions', };
totalRxnsToIncludeAB = totalRxnsToIncludeA;
remTransRxnsWnewSubs = {};
remTransRxnsSubs = {};
rxnsToInclude = cell(length(remainingTransRxns), 1);
subSToInclude = cell(length(remainingTransRxns), 1);
for i = 1:length(remainingTransRxns)
    metsRemTransRxns = findMetsFromRxns(model2, remainingTransRxns(i));
    metsRxns = {};
    for j = 1:length(metsRemTransRxns)
        if true(contains(metsRemTransRxns{j},'_c') &&...
                ~ismember(metsRemTransRxns(j), cofactorsListWcompB))
            metsRxns = findRxnsFromMets(model2,metsRemTransRxns(j));
        end
    end
    metsRxns = metsRxns(~cellfun('isempty',metsRxns));

    rxnSubs = {};
    for j = 1:length(metsRxns)
        rxnID = findRxnIDs(model2, metsRxns(j));
        rxnSubs{j,1} = model2.subSystems{rxnID};
    end
    rxnSubs = rxnSubs(~cellfun('isempty',rxnSubs));

    r = [];
    rankMax = 0;
    subSToUseStr = cell(1,1);
    for j = 1:length(rxnSubs)
        if true(~ismember(string(metsRxns(j)),totalRxnsToIncludeA))
            for k = 1:length(subsLength)
                if true(and(ismember(rxnSubs{j}, subsLength{k,2}),...
                        ~ismember(rxnSubs{j}, nonIncludSubs)))
                    r(k, 1) = subsLength{k, 1};
                    subSToUseStr{1, 1} = string(subsLength(k,2));
                    if k > 1
                        if true(rankMax < r(k, 1))
                            rankMax = r(k, 1);
                            subSToUseStr{1, 1} = string(subsLength(k,2));
                        end
                    end
                end
            end
        end
    end
    if ~isempty(subSToUseStr(1,1))
        rxnsToInclude{i} = remainingTransRxns{i};
        subSToInclude{i} = char(subSToUseStr{1, 1});
        remTransRxnsWnewSubs{i,1} = remainingTransRxns(i);
        remTransRxnsSubs{i,1} = subSToUseStr(1, 1);
        totalRxnsToIncludeAB = [totalRxnsToIncludeAB; rxnsToInclude];
    end
end
totalRxnsToIncludeAB = totalRxnsToIncludeAB(~cellfun('isempty',totalRxnsToIncludeAB));
totalRxnsToIncludeAB = unique(totalRxnsToIncludeAB);

for j = 1:length(rxnsToInclude)
    rxnID = findRxnIDs(model2, rxnsToInclude{j});
    model2.subSystems(rxnID) = {subSToInclude(j)};
end

%%%%%%%%%%%%%%%%%%%%%%%% ASIGN SUBSYSTEMS' COLOR %%%%%%%%%%%%%%%%%%%%%%%%%%%

% Selection of reactions to include in the map
if rxnsFluxTask == 1
    for i = 1:length(model2.rxns)
        if absFlxs(i)==0 % zero flux reactions are not included
            model2.subSystems(i) = {{'NonIncluded'}};  % included at the end
        end
    end
elseif rxnsFluxTask == 2
    for i = 1:length(model2.rxns)
        if absFlxs(i) > 0.95 * max(absFlxs)% reactions causing growth (i.e. obj function) defect when deleting
            model2.subSystems(i) = {{'NonIncluded'}};   % included at the end
        end
    end
elseif rxnsFluxTask == 3
    for i = 1:length(model2.rxns)
        if absFlxs(i) > 0.0 % zero flux reactions only, i.e. essential reations
            model2.subSystems(i) = {{'NonIncluded'}};   % included at the end
        end
    end
elseif rxnsFluxTask == 4
    for i = 1:length(model2.rxns)
        if absFlxs(i) < 0.95 * max(absFlxs)% reactions fluxes causing no growth (i.e. obj function) defect when deleting
            model2.subSystems(i) = {{'NonIncluded'}};   % included at the end
        end
    end
elseif rxnsFluxTask == 5
    for i = 1:length(model2.rxns)
        if absFlxs(i) < cutValue % zero flux reactions are not included
            model2.subSystems(i) = {{'NonIncluded'}};  % included at the end
        end
    end
end

%See colors letter-code at:
%la.mathworks.com/matlabcentral/mlc-downloads/downloads/submissions/24497/versions/2/screenshot.PNG
%[xmlMap, mapMap] = transformXML2Map('YeastMapUnified1.0b.xml');
centralPathsList = {'S_GlycolysisGluconeogenesis', 'S_Pentose_Phosphate_Pathway',...
    'S_Pyruvate_Metabolism', 'S_Citric_Acid_Cycle', 'S_Taurine_Metabolism',...
    'S_Starch_and_Sucrose_Metabolism', 'S_Xylose_Metabolism',...
    'S_Fructose_and_Mannose_Metabolism', 'S_Galactose_metabolism',...
    'S_Arabinose_Metabolism', 'S_Alternate_Carbon_Metabolism', 'S_Methane_Metabolism'};
aminoacPathsList = {'S_Glutamate_metabolism', 'S_Glutamine_Metabolism',...
    'S_Alanine_and_Aspartate_Metabolism', 'S_Histidine_Metabolism',...
    'S_Threonine_and_Lysine_Metabolism', 'S_Glycine_and_Serine_Metabolism',...
    'S_Tyrosine__Tryptophan__and_Phenylalanine_Metabolism',...
    'S_Arginine_and_Proline_Metabolism', 'S_Cysteine_Metabolism',...
    'S_Methionine_Metabolism', 'S_Valine__Leucine__and_Isoleucine_Metabolism',...
    'S_Other_Amino_Acid_Metabolism', 'S_Asparagine_metabolism', 'S_Nitrogen_Metabolism'};
cofactorsPathsList = {'S_NAD_Biosynthesis', 'S_Pantothenate_and_CoA_Biosynthesis',...
    'S_Folate_Metabolism', 'S_Riboflavin_Metabolism', 'S_Quinone_Biosynthesis',...
    'S_Porphyrin_and_Chlorophyll_Metabolism'};
lipidsPathsList = {'S_Fatty_Acid__Biosynthesis', 'S_Fatty_Acid_Degradation',...
    'S_Fatty_Acid_Metabolism', 'S_Glycerolipid_Metabolism', 'S_Glycoprotein_Metabolism',...
    'S_Phospholipid_Biosynthesis', 'S_Phospholipid_Metabolism',...
    'S_Sphingolipid_Metabolism', 'S_Sterol_Metabolism'};
nucleotidsPathsList = {'S_Purine_and_Pyrimidine_Biosynthesis',...
    'S_Nucleotide_Salvage_Pathway', 'S_Thiamine_Metabolism', 'S_Pyridoxine_Metabolism'};
enrgyPathsList = {'Biomass and maintenance functions', 'S_Anaplerotic_reactions',...
    'S_Oxidative_Phosphorylation'};
colors = createColorsMap;

if isempty(rxnLineColor)||rxnLineColor == "Pathways"
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_GlycolysisGluconeogenesis', 'DARKBLUE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Pentose_Phosphate_Pathway', 'STEELBLUE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Pyruvate_Metabolism', 'DODGERBLUE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Citric_Acid_Cycle', 'MEDIUMBLUE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Taurine_Metabolism', 'MEDIUMBLUE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Starch_and_Sucrose_Metabolism', 'BLUE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Xylose_Metabolism', 'STEELBLUE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Fructose_and_Mannose_Metabolism', 'DODGERBLUE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Galactose_metabolism', 'MEDIUMBLUE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Arabinose_Metabolism', 'BLUE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Alternate_Carbon_Metabolism', 'DARKBLUE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Methane_Metabolism', 'DARKBLUE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Taurine_Metabolism', 'DARKBLUE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Anaplerotic_reactions', 'CRIMSON');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Complex_Alcohol_Metabolism', 'MAROON');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Oxidative_Phosphorylation', 'RED');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Purine_and_Pyrimidine_Biosynthesis', 'PINK');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Nucleotide_Salvage_Pathway', 'HOTPINK');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Thiamine_Metabolism', 'MAGENTA');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Pyridoxine_Metabolism', 'PINK');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_NAD_Biosynthesis', 'PLUM');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Pantothenate_and_CoA_Biosynthesis', 'VIOLET');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Folate_Metabolism', 'DARKORCHID');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Riboflavin_Metabolism', 'PURPLE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Quinone_Biosynthesis', 'MEDIUMORCHID');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Porphyrin_and_Chlorophyll_Metabolism', 'PLUM');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Glutamate_metabolism', 'ORANGE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Glutamine_Metabolism', 'DARKORANGE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Alanine_and_Aspartate_Metabolism', 'CORAL');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Histidine_Metabolism', 'GOLD');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Threonine_and_Lysine_Metabolism', 'ORANGE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Glycine_and_Serine_Metabolism', 'DARKORANGE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Tyrosine__Tryptophan__and_Phenylalanine_Metabolism', 'CORAL');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Arginine_and_Proline_Metabolism', 'GOLD');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Cysteine_Metabolism', 'ORANGE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Methionine_Metabolism', 'DARKORANGE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Valine__Leucine__and_Isoleucine_Metabolism', 'CORAL');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Other_Amino_Acid_Metabolism', 'GOLD');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Asparagine_metabolism', 'ORANGE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Nitrogen_Metabolism', 'DARKORANGE');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Fatty_Acid__Biosynthesis', 'FORESTGREEN');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Fatty_Acid_Degradation', 'MEDIUMSEAGREEN');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Fatty_Acid_Metabolism', 'DARKGREEN');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Glycerolipid_Metabolism', 'SEAGREEN');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Glycoprotein_Metabolism', 'FORESTGREEN');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Phospholipid_Biosynthesis', 'MEDIUMSEAGREEN');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Phospholipid_Metabolism', 'DARKGREEN');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Sphingolipid_Metabolism', 'SEAGREEN');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Sterol_Metabolism', 'FORESTGREEN');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Other', 'BLACK');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'Biomass and maintenance functions', 'RED');
elseif rxnLineColor == "Subsystems"
    modifiedMap = colorSubsystem(modifiedMap, model2, centralPathsList, 'BLUE'); %mapMap to fluxMap
    modifiedMap = colorSubsystem(modifiedMap, model2, aminoacPathsList, 'ORANGE');
    modifiedMap = colorSubsystem(modifiedMap, model2, cofactorsPathsList, 'ORCHID');
    modifiedMap = colorSubsystem(modifiedMap, model2, lipidsPathsList, 'LIME');
    modifiedMap = colorSubsystem(modifiedMap, model2, nucleotidsPathsList, 'MAGENTA');
    modifiedMap = colorSubsystem(modifiedMap, model2, enrgyPathsList, 'RED');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Complex_Alcohol_Metabolism', 'GOLDENROD');
    modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Other', 'GOLDENROD');
elseif ~(rxnLineColor == "Subsystems" || rxnLineColor == "Payhways" || isempty(rxnLineColor))
    for i = 1:length(modifiedMap.rxnName)
        a = find(ismember(model2.rxns, modifiedMap.rxnName{i}));
        if ~isempty(a) & string(model2.subSystems(a)) ~= "NonIncluded"
            modifiedMap.rxnColor(i, 1) = {colors(rxnLineColor)};
        end
    end
end
if rxnLineColor == "Subsystems" || rxnLineColor == "Payhways" || isempty(rxnLineColor)
    h2o = {'H2Ot', 'H2Otm',	'H2Otn', 'H2Otv', 'H2Otp', 'H2Oter'};
    for i = 1:length(modifiedMap.rxnName)
        a = find(ismember(model2.rxns, modifiedMap.rxnName{i}));
        for j = 1:length(h2o)
            if ~isempty(a) & string(modifiedMap.rxnName{i}) == string(h2o(j))...
                    & string(model2.subSystems(a)) ~= "NonIncluded"
                modifiedMap.rxnColor(i, 1) = {colors('CYAN')};
            end
        end
    end
    o2 = {'O2t', 'O2tm', 'O2ter', 'O2tp'};
    for i = 1:length(modifiedMap.rxnName)
        a = find(ismember(model2.rxns, modifiedMap.rxnName{i}));
        for j = 1:length(o2)
            if ~isempty(a) & string(modifiedMap.rxnName{i}) == string(o2(j))...
                    & string(model2.subSystems(a)) ~= "NonIncluded"
                modifiedMap.rxnColor(i, 1) = {colors('CYAN')};
            end
        end
    end
    Pi = {'PIt2r', 'PIt2m', 'PIt5m', 'PIt2n', 'PIt2v', 'PIt2p', 'ATPtm_H', 'ATPtp_H', ...
        'ATP2tp_H', 'ATPS3g', 'ATPS3v', 'UTPtm', 'CTPtm', 'GTPt2m', 'ATPS'};
    for i = 1:length(modifiedMap.rxnName)
        a = find(ismember(model2.rxns, modifiedMap.rxnName{i}));
        for j = 1:length(Pi)
            if ~isempty(a) & string(modifiedMap.rxnName{i}) == string(Pi(j))...
                    & string(model2.subSystems(a)) ~= "NonIncluded"
                modifiedMap.rxnColor(i, 1) = {colors('RED')};
            end
        end
    end
    co2 = {'CO2t', 'CO2tg', 'CO2tm', 'CO2tn', 'HCO3tn', 'CO2tv', 'CO2tp'};
    for i = 1:length(modifiedMap.rxnName)
        a = find(ismember(model2.rxns, modifiedMap.rxnName{i}));
        for j = 1:length(co2)
            if ~isempty(a) & string(modifiedMap.rxnName{i}) == string(co2(j))...
                    & string(model2.subSystems(a)) ~= "NonIncluded"
                modifiedMap.rxnColor(i, 1) = {colors('SLATEGRAY')};
            end
        end
    end
    aa = {'MMETt2', 'GLUt2r', 'GLNt2r', 'GLNt7', 'GLUt7', 'LCYSTintv', 'NH4t', 'TCHOLAabcv'};
    for i = 1:length(modifiedMap.rxnName)
        a = find(ismember(model2.rxns, modifiedMap.rxnName{i}));
        for j = 1:length(aa)
            if ~isempty(a) & string(modifiedMap.rxnName{i}) == string(aa(j))...
                    & string(model2.subSystems(a)) ~= "NonIncluded"
                modifiedMap.rxnColor(i, 1) = {colors('ORANGE')};
            end
        end
    end
    cenMet = {'COAtim', 'NAt3_1', 'MALS'};
    for i = 1:length(modifiedMap.rxnName)
        a = find(ismember(model2.rxns, modifiedMap.rxnName{i}));
        for j = 1:length(cenMet)
            if ~isempty(a) & string(modifiedMap.rxnName{i}) == string(cenMet(j))...
                    & string(model2.subSystems(a)) ~= "NonIncluded"
                modifiedMap.rxnColor(i, 1) = {colors('BLUE')};
            end
        end
    end
    lip = {'G3PIt'};
    for i = 1:length(modifiedMap.rxnName)
        a = find(ismember(model2.rxns, modifiedMap.rxnName{i}));
        for j = 1:length(lip)
            if ~isempty(a) & string(modifiedMap.rxnName{i}) == string(lip(j))...
                    & string(model2.subSystems(a)) ~= "NonIncluded"
                modifiedMap.rxnColor(i, 1) = {colors('LIME')};
            end
        end
    end
end

if rxnsFluxTask == 5 % change color to dark-color in negative values
    for i = 1:length(modifiedMap.rxnName)
        valueskeys = values(colors)';
        valueskeys (:,2) = keys(colors)';
        a = find(ismember(valueskeys, modifiedMap.rxnColor{i, 1}));
        b = find(ismember(model2.rxns, modifiedMap.rxnName{i}));
        if ~isempty(b) & v1(b) < 0 & string(model2.subSystems(b)) ~= "NonIncluded"
            if string(valueskeys (a,2)) == 'LIME'
                modifiedMap.rxnColor(i, 1) = {colors('GREEN')};
            else
                rxnColor =strcat('DARK', char(valueskeys (a(1,1),2)));
                if string(rxnColor) == 'DARKAQUA'
                    rxnColor = 'DARKCYAN';
                elseif string(rxnColor) == 'DARKFUCHSIA'
                    rxnColor = 'DARKMAGENTA';
                end
                modifiedMap.rxnColor(i, 1) = {colors(rxnColor)};
            end
        end
    end
end

end

author: SysBioChalmers

repo: RAVEN


text:
function value = p14_gu ( u )

%*****************************************************************************80
%
%% P14_GU computes the auxilliary function G(U).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real U, the argument of the function.
%
%    Output, real VALUE, the value of the function.
%
  beta = 0.1;
  gamma = 0.1;

  value = u / ( beta * ( u + gamma ) );

  return
end


text:
function [ap_results05,ap_results07,ap_resultsCOCO,mean_recall_per_IoU] = print_detection_evaluation_metrics(...
    all_bbox_gt, abbox_dets, classes)
%
% This file is part of the code that implements the following paper:
% Title      : "LocNet: Improving Localization Accuracy for Object Detection"
% Authors    : Spyros Gidaris, Nikos Komodakis
% Institution: Universite Paris Est, Ecole des Ponts ParisTech
% ArXiv link : http://arxiv.org/abs/1511.07763
% code       : https://github.com/gidariss/LocNet
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2016 Spyros Gidaris
%
% Title     : "LocNet: Improving Localization Accuracy for Object Detection"
% ArXiv link: http://arxiv.org/abs/1511.07763
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

% evaluate the average precision for IoU >=0.5
ap_results05 = evaluate_average_precision_pascal(all_bbox_gt, abbox_dets, classes );
fprintf('AP with IoU >= 0.5\n');
printAPResults(classes, ap_results05);
% evaluate the average precision for IoU >=0.7
ap_results07 = evaluate_average_precision_pascal(all_bbox_gt, abbox_dets, classes, 'minoverlap', 0.7);
fprintf('AP with IoU >= 0.7\n');
printAPResults(classes, ap_results07);
% evaluate the coco style average precision
ap_resultsCOCO = evaluate_average_precision_pascal(all_bbox_gt, abbox_dets, classes, 'coco_style', true);
fprintf('AP with coco style\n');
printAPResults(classes, ap_resultsCOCO);

num_classes = length(classes);
all_bbox_gt_no_difficults = remove_difficult_ground_truth_bboxes(all_bbox_gt);
all_bbox_gt_per_class = get_per_class_ground_truth_bboxes(all_bbox_gt_no_difficults, num_classes);

ave_recall_per_class = zeros(num_classes, 1);
recall_per_class     = cell(1,num_classes);

% compute the recall per IoU threshold and the average recall of the set of
% detections for each class separately
for j = 1:num_classes
    [ave_recall_per_class(j), recall_per_class{j}, thresholds] = ...
        compute_ave_recall_of_bbox(abbox_dets{j}, all_bbox_gt_per_class{j});
    recall_per_class{j} = recall_per_class{j}(:);
    recall_per_class{j} = recall_per_class{j}(thresholds>=0.5);
end
thresholds = thresholds(thresholds>=0.5);
fprintf('Average Recall of detections :\n');
printAPResults(classes, ave_recall_per_class);
fprintf('Recall per IoU of detections (averaged over the classes):\n');
mean_recall_per_IoU = mean(cell2mat(recall_per_class),2); % average over the classes
thresholds_string = strsplit(num2str(thresholds));
assert(length(mean_recall_per_IoU) == length(thresholds_string));
printAPResults(thresholds_string, mean_recall_per_IoU);
end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
function F=tstatis(X,y,P,N)

if nargin<4;N=3;end
if nargin<3;P=1000;end

[M Nvar]=size(X);
t0=tstatis1(X,y);
k1=find(y==1);
k2=find(y~=1);
T=zeros(P,Nvar);
TMC=zeros(N,Nvar);
Q=floor(M*0.9);

%++++ Reliability index
i=1;
while i<=N
  index=randperm(M);   % MCS
  cal=index(1:Q);
  tmc=tvalue(X(cal,:),y(cal));
  TMC(i,:)=tmc;
  fprintf('The %dth MCS finished.\n',i);
  i=i+1;
end
RI=abs(mean(TMC)./std(TMC));

%+++ Permutation test for p value computation.
i=1;
while i<=P
  index=randperm(M);   % permutation
  yr=y(index);
  if sum(abs(yr-y))==0;continue;end
  tr=tvalue(X,yr);
  T(i,:)=tr;
  fprintf('The %dth permutation finished.\n',i);
  i=i+1;
end

p=zeros(1,Nvar);
for i=1:Nvar
   p(i)=length(find(abs(T(:,i))>abs(t0(i))))/P;
end

%+++ Output
F.t0=abs(t0)/max(abs(t0));
F.TMC=TMC;
F.RI=RI;
F.tr=T;
F.p=p;

author: viggin

repo: domain-adaptation-toolbox


text:
classdef NaiveClustererX < ClustererX
% NAIVECLUSTERERX Class
%
% Summary of NaiveClustererX:
% This is a class implementation of a naive clusterer that generates clusters
% of tracks sharing common measurements.
%
% NaiveClustererX Properties:
%   + NumMeasDims - The number of observation dimensions.
%
% NaiveClustererX Methods:
%   + NaiveClustererX  - Constructor method
%   + cluster - Perform clustering and generate a list of clusters
%
% (+) denotes public properties/methods
%
% See also SystematicResamplerX
    properties
    end

    properties (SetAccess = immutable)
        NumMeasDims
    end

    methods
        function this = NaiveClustererX(varargin)
        % ELLIPSOIDALGATERX Constructor method
        %
        % Usage
        % -----
        % * nv = NaiveClustererX() returns a NaiveClustererX object
        %
        % See also NaiveClustererX/cluster

        end

        function [ClusterList,UnassocTrackInds] = cluster(this,ValidationMatrix)
        % CLUSTER Perform naive clustering to generate clusters of tracks
        % sharing common measurements.
        %
        % Parameters
        % ----------
        % ValidationMatrix: matrix
        %   A (Nt x Nm) validation matrix, Nt being the number of tracks and
        %   Nm being the number of measurements, where each element (t,m) is
        %   a binary variable (0 or 1), representing whether measurement m
        %   fell in the gate of target t.
        %
        % Returns
        % -------
        % ClusterList: cell vector
        %   A (1 x Nc) cell vector, where each cell represents one of Nc Cluster
        %   objects, with the following fields:
        %       - MeasIndList: A list of the indices of all measurements
        %                     contained within the cluster.
        %       - TrackIndList: A list of indices of all tracks belonging
        %                       to the cluster.
        %
        % UnassocTrackInds: column vector
        %  A (1 x Nu) column vector, where each element contains the index
        %  (as ordered in ValidationMatrix) of any tracks that have not
        %  been associated to any measurements. As such, 0<= Nu <= Nt.
        %
        % Usage
        % -----
        % * [ClusterList,UnassocTracks] = cluster(this,ValidationMatrix)
        %   returns a list of clusters ClusterList and a list of unassociated
        %   track indices UnassocTracks (corresponding to the row indices of
        %   ValidationMatrix).
        %   ClusterList is a list of Cluster objects, where each cluster
        %   object has two properties:
        %       - MeasIndList: A list of the indices of all measurements
        %                     contained within the cluster.
        %       - TrackIndList: A list of indices of all tracks belonging
        %                       to the cluster.
        %
        % See also NaiveClustererX/NaiveClustererX

            % Initiate parameters
            NumTracks = size(ValidationMatrix,1); % Number of measurements

            % Form clusters of tracks sharing measurements
            UnassocTrackInds = [];
            ClusterList = [];
            ClusterObj.MeasIndList = [];
            ClusterObj.TrackIndList = [];

            % Iterate over all tracks
            for trackInd=1:NumTracks
                % Extract valid measurement indices
                validMeasInds = find(ValidationMatrix(trackInd,:));

                % If there exist valid measurements
                if (~isempty(validMeasInds))

                    % Check if matched measurements are members of any clusters
                    NumClusters = numel(ClusterList);
                    matchedClusterIndFlags = zeros(1, NumClusters);
                    for ClusterInd=1:NumClusters
                        if (sum(ismember(validMeasInds, ClusterList(ClusterInd).MeasIndList)))
                            matchedClusterIndFlags(ClusterInd) = 1; % Store matched cluster ids
                        end
                    end

                    NumMatchedClusters = sum(matchedClusterIndFlags);
                    matchedClusterInds = find(matchedClusterIndFlags);

                    % If only matched with a single cluster, join.
                    switch(NumMatchedClusters)
                        case(1)
                            ClusterList(matchedClusterInds).TrackIndList = union(ClusterList(matchedClusterInds).TrackIndList, trackInd);
                            ClusterList(matchedClusterInds).MeasIndList = union(ClusterList(matchedClusterInds).MeasIndList, validMeasInds);
                        case(0)
                            ClusterList(end+1).TrackIndList = trackInd;
                            ClusterList(end).MeasIndList = validMeasInds;
                            %ClusterList(end+1) = ClusterObj;
                        otherwise
                            % Start from last cluster, joining each one with the previous
                            %   and removing the former.
                            for matchedClusterInd = NumMatchedClusters-1:-1:1
                                ClusterList(matchedClusterInds(matchedClusterInd)).TrackIndList = ...
                                    union(ClusterList(matchedClusterInds(matchedClusterInd)).TrackIndList, ...
                                        ClusterList(matchedClusterInds(matchedClusterInd+1)).TrackIndList);
                                ClusterList(matchedClusterInds(matchedClusterInd)).MeasIndList = ...
                                    union(ClusterList(matchedClusterInds(matchedClusterInd)).MeasIndList, ...
                                        ClusterList(matchedClusterInds(matchedClusterInd+1)).MeasIndList);
                                ClusterList(matchedClusterInds(matchedClusterInd+1)) = [];
                            end
                            % Finally, join with associated track.
                            ClusterList(matchedClusterInds(matchedClusterInd)).TrackIndList = ...
                                union(ClusterList(matchedClusterInds(matchedClusterInd)).TrackIndList, trackInd);
                            ClusterList(matchedClusterInds(matchedClusterInd)).MeasIndList = ...
                                union(ClusterList(matchedClusterInds(matchedClusterInd)).MeasIndList, validMeasInds);
                    end
                else
                    ClusterList(end+1).TrackIndList = trackInd;
                    ClusterList(end).MeasIndList = [];
                    %ClusterList(end+1) = ClusterObj;
                    UnassocTrackInds = [UnassocTrackInds trackInd];
                end
                this.ClusterList = ClusterList;
                this.UnassocTrackInds = UnassocTrackInds;
            end
        end
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function [V, policy, iter, cpu_time] = mdp_policy_iteration_modified(P, R, discount, epsilon, max_iter)

% mdp_policy_iteration_modified    Resolution of discounted MDP
%                                  with modified policy iteration algorithm
% Arguments -------------------------------------------------------------
% Let S = number of states, A = number of actions
%   P(SxSxA) = transition matrix
%              P could be an array with 3 dimensions or
%              a cell array (1xA), each cell containing a matrix (SxS) possibly sparse
%   R(SxSxA) or (SxA) = reward matrix
%              R could be an array with 3 dimensions (SxSxA) or
%              a cell array (1xA), each cell containing a sparse matrix (SxS) or
%              a 2D array(SxA) possibly sparse
%   discount = discount rate in ]0, 1]
%              beware to check conditions of convergence for discount = 1.
%   epsilon  = epsilon-optimal policy search, upper than 0,
%              optional (default : 0.01)
%   max_iter = maximum number of iteration to be done in the inner loop,
%              upper than 0, optional (default: 10)
% Evaluation -------------------------------------------------------------
%   V(S)     = value function
%   policy(S)= epsilon-optimal policy
%   iter     = number of main iterations
%   cpu_time = used CPU time
%--------------------------------------------------------------------------
% In verbose mode, at each iteration, displays the variation of V

% MDPtoolbox: Markov Decision Processes Toolbox
% Copyright (C) 2009  INRA
% Redistribution and use in source and binary forms, with or without modification,
% are permitted provided that the following conditions are met:
%    * Redistributions of source code must retain the above copyright notice,
%      this list of conditions and the following disclaimer.
%    * Redistributions in binary form must reproduce the above copyright notice,
%      this list of conditions and the following disclaimer in the documentation
%      and/or other materials provided with the distribution.
%    * Neither the name of the <ORGANIZATION> nor the names of its contributors
%      may be used to endorse or promote products derived from this software
%      without specific prior written permission.
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
% ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
% IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
% INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
% BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
% OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
% OF THE POSSIBILITY OF SUCH DAMAGE.

cpu_time = cputime;

global mdp_VERBOSE;

% check of arguments
if discount <= 0 || discount > 1
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: Discount rate must be in ]0; 1]')
    disp('--------------------------------------------------------')
elseif nargin > 4 && epsilon <= 0
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: epsilon must be upper than 0')
    disp('--------------------------------------------------------')
elseif nargin > 5 && max_iter <= 0
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: The maximum number of iteration must be upper than 0')
    disp('--------------------------------------------------------')
else

    if discount == 1
        disp('-------------------------------------------------------')
        disp('MDP Toolbox WARNING: check conditions of convergence.')
        disp('With no discount, convergence is not always assumed.')
        disp('--------------------------------------------------------')
    end;

    if iscell(P); S = size(P{1},1); else S = size(P,1); end;

    PR = mdp_computePR(P,R);

    % initialization of optional arguments
    if nargin < 5; max_iter = 10; end;
    if nargin < 4; epsilon = 0.01; end;

    % computation of threshold of variation for V for an epsilon-optimal policy
    if discount ~= 1
        thresh = epsilon * (1-discount)/discount;
    else
        thresh = epsilon;
    end;

    if discount == 1
        V = zeros(S,1);
    else
        V = 1/(1-discount)*min(min(PR))*ones(S,1);
    end;

    if mdp_VERBOSE; disp('  Iteration  V_variation'); end;

    iter = 0;
    is_done = false;
    while ~is_done

        iter = iter + 1;

        [Vnext, policy] = mdp_bellman_operator(P,PR,discount,V);
        %[Ppolicy, PRpolicy] = mdp_computePpolicyPRpolicy(P, PR, policy);

        variation = mdp_span(Vnext - V);
        if mdp_VERBOSE;
             disp(['      ' num2str(iter,'%5i') '         ' num2str(variation)]);
        end;

        V=Vnext;
        if variation < thresh
            is_done = true;
        else
	    is_verbose = false;
            if mdp_VERBOSE; mdp_VERBOSE = 0; is_verbose = true; end;
            V = mdp_eval_policy_iterative(P, PR, discount, policy, V, epsilon, max_iter);
            if is_verbose; mdp_VERBOSE = 1; end;
        end;
    end;
end;

cpu_time = cputime - cpu_time;


text:
function model = Initialization_Discrete(Seqs, step, order)

switch nargin
    case 1
        Num = zeros(length(Seqs), 1);
        Tmax = zeros(length(Seqs), 1);
        for i = 1:length(Seqs)
            Num(i) = length(Seqs(i).Time);
            Tmax(i) = Seqs(i).Time(end) + eps;
        end
        model.h = mean(Tmax./Num);
        model.k = floor(min(0.7*Tmax/model.h));
    case 2
        model.h = step;
        Tmax = zeros(length(Seqs), 1);
        for i = 1:length(Seqs)
            Tmax(i) = Seqs(i).Time(end) + eps;
        end
        model.k = floor(min(0.7*Tmax/model.h));
    case 3
        model.h = step;
        model.k = order;
end

D = zeros(length(Seqs),1);
for i = 1:length(Seqs)
    D(i) = max(Seqs(i).Mark);
end
model.D = max(D);
model.Tmax = [];

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
% MAIN - Pendulum
%
% Demonstrates simple swing-up for a single pendulum with a torque motor.
%

clc; clear;
addpath ../../

% Physical parameters of the pendulum
p.k = 1;  % Normalized gravity constant
p.c = 0.1;  % Normalized damping constant

% User-defined dynamics and objective functions
problem.func.dynamics = @(t,x,u)( dynamics(x,u,p) );
problem.func.pathObj = @(t,x,u)( u.^2 );

% Problem bounds
problem.bounds.initialTime.low = 0;
problem.bounds.initialTime.upp = 0;
problem.bounds.finalTime.low = 0.5;
problem.bounds.finalTime.upp = 2.5;

problem.bounds.state.low = [-2*pi; -inf];
problem.bounds.state.upp = [2*pi; inf];
problem.bounds.initialState.low = [0;0];
problem.bounds.initialState.upp = [0;0];
problem.bounds.finalState.low = [pi;0];
problem.bounds.finalState.upp = [pi;0];

problem.bounds.control.low = -5; %-inf;
problem.bounds.control.upp = 5; %inf;

% Guess at the initial trajectory
problem.guess.time = [0,1];
problem.guess.state = [0, pi; pi, pi];
problem.guess.control = [0, 0];

% Select a solver:
problem.options.method = 'trapezoid';
problem.options.defaultAccuracy = 'medium';

% Solve the problem
soln = optimTraj(problem);
t = soln.grid.time;
q = soln.grid.state(1,:);
dq = soln.grid.state(2,:);
u = soln.grid.control;

% Plot the solution:
figure(1); clf;

subplot(3,1,1)
plot(t,q)
ylabel('q')
title('Single Pendulum Swing-Up');

subplot(3,1,2)
plot(t,dq)
ylabel('dq')

subplot(3,1,3)
plot(t,u)
ylabel('u')

author: MatthewPeterKelly

repo: OptimTraj


text:
function hRD = plugin_rd_tx
% Reference design definition

% Call the common reference design definition function
hRD = AnalogDevices.adrv9361z7035.common.plugin_rd('ccbox_lvds','Tx');
AnalogDevices.adrv9361z7035.common.add_io(hRD, 'Tx', 'ccbox_lvds');

author: analogdevicesinc


text:
function linpack_d_test18 ( )

%*****************************************************************************80
%
%% TEST18 tests DPOFA and DPOSL.
%
%  Discussion:
%
%    DPOFA factors a positive definite symmetric matrix,
%    and DPOSL can solve a factored linear system.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 20;
  lda = n;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST18\n' );
  fprintf ( 1, '  For a positive definite symmetric matrix,\n' );
  fprintf ( 1, '  DPOFA computes the LU factors.\n' );
  fprintf ( 1, '  DPOSL solves a factored linear system.\n' );
  fprintf ( 1, '  The matrix size is N = %d\n', n );
%
%  Set the matrix A.
%
  a(1:n,1:n) = 0.0;

  for i = 1 : n
    a(i,i) = 2.0;
    if ( 1 < i )
      a(i,i-1) = -1.0;
    end
    if ( i < n )
      a(i,i+1) = -1.0;
    end
  end
%
%  Set the right hand side.
%
  x = zeros ( n, 1 );

  for i = 1 : n
    x(i) = i;
  end

  b(1:n) = a(1:n,1:n) * x(1:n);
%
%  Factor the matrix.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Factor the matrix.\n' );

  [ a, info ] = dpofa ( a, lda, n );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Error, DPOFA returns INFO = %d\n', info );
    return
  end
%
%  Solve the linear system.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solve the linear system.\n' );

  b = dposl ( a, lda, n, b );
%
%  Print the result.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The first and last five entries of the solution:\n' );
  fprintf ( 1, '  (Should be 1,2,3,4,5,...,n-1,n.)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    if ( i <= 5 | n-5 < i )
      fprintf ( 1, '  %6d  %14f\n', i, b(i) );
    end
    if ( i == 5 )
      fprintf ( 1, '  ......  ..............\n' );
    end
  end

  return
end


text:
function [SuperFiber, fg] = dtiComputeSuperFiberRepresentation(fg, clusterlabels, numNodes, M)
%
% [SuperFiber, fg] = dtiComputeSuperFiberRepresentation(fg, [clusterlabels], [numNodes])
%
% Superfibers are a single curve (fiber) that summarizes all of the fibers
% in a single labeled cluster.  For example, one might compute the super
% fiber associated with the arcuate.
%
% Create SuperFiber representation of all the fibers assigned to one
% cluster. If numNodes is not provided, it is assumed that they all are
% resampled to same number of Nodes as in fg.fibers{1}.
%
%  INPUTS:
%    fg            - fiber group structure
%    clusterLabels - Possible labels of the fiber subgroups
%    numNodes      - Number of samples (nodes) along the fiber
%    M             - 'mean' or 'median' for calculating tract core
%
%  OUTPUTS:
%   The parameters in the returned representation are:
%   SuperFiber.fibermeans: means (for every coordinate, every node)
%   SuperFiber.fibervarcovs: low diag var/cov matrix for coordinates
%                          (every node)): [var(x), cov(x, y),
%                          cov(x, z), var(y), cov(y, z), var(z)]
%   SuperFiber.n is number of fibers included in this cluster.
%   fg - returned containes flipped and resampled fibers that produced the
%        superfiber.
%
% Note: while computing superFiber representation, some fibers (or even
% majority of fibers) may have to have been  flipped to group first and
% last points together. The 'fg' returned containes flipped and resampled
% fibers that produced the superfiber.
%
% Algorithm:
%
% All fibers in the group (e.g., optic tract) are resampled to represent
% equal numbers of nodes. We determine the mean and 3D covariance of the
% voxel positions representing that node. We define the core fiber as the
% set of mean positions of the sample nodes.
%
% In subsequent routines, we determine the weight contributed by each voxel
% to the diffusivity measure using the covariance matrix. Specifically, we
% calculate the Mahalanobis distance of each voxel from the mean. If a
% voxel position is X, the mean position is X0, and the covariance matrix
% of the 3D positions is C, then the Mahalanobis distance, d, is sqrt(delta
% Cinv delta), where X-X0 is delta.
%
% The assgined weight is the inverse of the distance, d. This procedure
% assigns more weight to fibers close to the core of the bundle, and the
% meaning of close is derived from the distribution of the points at each
% node as captured by the covariance matrix.
%
%  WEB RESOURCES:
%   mrvBrowseSVN('dtiComputeSuperFiberRepresentation');
%   http://white.stanford.edu/newlm/index.php/Diffusion_properties_along_tr
%   ajectory
%
%  EXAMPLE USAGE:
%   % Create superfiber
%       fg = mtrImportFibers('fiberGroupName.pdb');
%       numNodes = 30;
%       [SuperFiber, sfg] = dtiComputeSuperFiberRepresentation(fg, [], numNodes);
%       sFg = dtiNewFiberGroup([fg.name '_superFiber'],[0 0 155],[],[],SuperFiber.fibers);
%   % Save the superfiber
%       mtrExportFibers(sFg,sFg.name);
%
%  HISTORY:
%   2008       ER wrote it
%   08/11/2009 ER added checks that all start points and all end points are
%              grouped respectively.
%   09.13.2011 LMP Added comments, example usage, web resources and
%              input/output descriptions

if(notDefined('numNodes'))
    numNodes=length(fg.fibers{1});
    if length(unique(cellfun(@length, fg.fibers)))>1
        % The numNodes should be the same across all fiber groups, or
        % The user should provide the numNodes for each of the fiber groups
        error('Either provide numNodes to resample the fg to, or number of nodes in the existing fg should be same across the fibers');
    end
end
if notDefined('M')
    M='mean';
elseif ~strcmp(M,'mean') && ~strcmp(M,'median')
    error('M must be either mean or median')
end
% If cluster labels is not sent in, then we assign cluster labels all the
% value of one according to the size of fg.fibers.  FG is a fiber group,
% and it consists of a set of cells that define the nodes of fibers.  This
% assignment says that all of the fibers in the fiber group are in the same
% cluster.  That cluster is labeled cluster 1.
if(notDefined('clusterlabels')), clusterlabels = ones(size(fg.fibers));end

%Check that the all the starting points and all the end points are grouped.
%This routine also resamples the fibers to numNodes samples.
fg = dtiReorientFibers(fg, numNodes);

% The super name
SuperFiber.name=[fg.name '_SuperFiber'];

%The number of fibers in the fiber group
nfibers = size(fg.fibers, 1);
% This should be the same numNodes we sent in above
numNodes= size(fg.fibers{1}, 2);

% curves is the 3D representation of each node for each fiber.
curves  = zeros(3, numNodes, nfibers);
for ii=1:nfibers
    curves(:, :, ii) = fg.fibers{ii};
end

% For every cluster, there is a superfiber representation
% The clusters are labeled as numbers, so max(clusterlabels) is the number
% of clusters.
for clust = 1:max(clusterlabels)

    % For this cluster, find the number of fibers.
    SuperFiber.n(clust, 1) = size(find(clusterlabels==clust), 1);

    % For every node in the numNode list, compute the mean and covariance
    % cloud of the nodes in 3-space.
    for node = 1:numNodes

        %Find the mean position of the 3D coordinates at this node, for
        %this cluster.
        if strcmp(M,'mean')
            SuperFiber.fibers{clust, 1}(:, node) = mean(curves(:, node, clusterlabels==clust), 3);
        elseif strcmp(M,'median')
            SuperFiber.fibers{clust, 1}(:, node) = median(curves(:, node, clusterlabels==clust), 3);
        end
        % Find the lower diagonal matrix of the cloud of points at this
        % node.  Say, for cluster 1 we find the 3D values at this node.  We
        % permute the coordinates to because (help, explain why and what it
        % is matched to), we compute the covariance matrix of these 3D
        % points. The function tril zeros the upper diagonal entries.
        % Later, ER stores the lower triangular part of the symmetric)
        % matrix.
        %
        % This is a very complicated representation that is probably
        % unpacked later in other functions.  It would be nice to specify
        % in the header what functions need to be aware of this
        % representation.
        LowDiagonalMatrixThisNode = ...
            tril(cov(permute(curves(:, node, clusterlabels==clust), [3 1 2])));

        if size(find(clusterlabels==clust), 1)==1
            %ONLY ONE FIBER IN THE GROUP
            LowDiagonalMatrixThisNode = zeros(3);
        end

        %These are low diagonal elements of the 3x3 covariance matrix.
        SuperFiber.fibervarcovs{clust, 1}(:, node) = ...
            LowDiagonalMatrixThisNode([1:3 5:6 9]');

    end %node

end %clust

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [c,relres,iter]=franaiter(F,f,varargin)
%FRANAITER  Iterative analysis
%   Usage:  c=franaiter(F,f);
%           [c,relres,iter]=franaiter(F,f,...);
%
%   Input parameters:
%         F       : Frame.
%         f       : Signal.
%         Ls      : Length of signal.
%   Output parameters:
%         c       : Array of coefficients.
%         relres  : Vector of residuals.
%         iter    : Number of iterations done.
%
%   `c=franaiter(F,f)` computes the frame coefficients *c* of the signal *f*
%   using an iterative method such that perfect reconstruction can be
%   obtained using |frsyn|. `franaiter` always works, even when |frana|
%   cannot generate perfect reconstruction coefficients.
%
%   `[c,relres,iter]=franaiter(...)` additionally returns the relative
%   residuals in a vector *relres* and the number of iteration steps *iter*.
%
%   **Note:** If it is possible to explicitly calculate the canonical dual
%   frame then this is usually a much faster method than invoking
%   `franaiter`.
%
%   `franaiter` takes the following parameters at the end of the line of
%   input arguments:
%
%     'tol',t      Stop if relative residual error is less than the
%                  specified tolerance. Default is 1e-9 (1e-5 for single precision)
%
%     'maxit',n    Do at most n iterations.
%
%     'pg'        Solve the problem using the Conjugate Gradient
%                  algorithm. This is the default.
%
%     'pcg'        Solve the problem using the Preconditioned Conjugate Gradient
%                  algorithm.
%
%     'print'      Display the progress.
%
%     'quiet'      Don't print anything, this is the default.
%
%   Examples
%   --------
%
%   The following example shows how to rectruct a signal without ever
%   using the dual frame:::
%
%      f=greasy;
%      F=frame('dgtreal','gauss',40,60);
%      [c,relres,iter]=franaiter(F,f,'tol',1e-14);
%      r=frsyn(F,c);
%      norm(f-r)/norm(f)
%      semilogy(relres);
%      title('Conversion rate of the CG algorithm');
%      xlabel('No. of iterations');
%      ylabel('Relative residual');
%
%   See also: frame, frana, frsyn, frsyniter

% AUTHORS: Peter L. Søndergaard

complainif_notenoughargs(nargin,2,'FRANAITER');
complainif_notvalidframeobj(F,'FRANAITER');

tolchooser.double=1e-9;
tolchooser.single=1e-5;

definput.keyvals.Ls=[];
definput.keyvals.tol=tolchooser.(class(f));
definput.keyvals.maxit=100;
definput.flags.alg={'cg','pcg'};
definput.keyvals.printstep=10;
definput.flags.print={'quiet','print'};

[flags,kv,Ls]=ltfatarghelper({'Ls'},definput,varargin);

%% ----- step 1 : Verify f and determine its length -------
% Change f to correct shape.
[f,~,Ls,W,dim,permutedsize,order]=assert_sigreshape_pre(f,[],[],upper(mfilename));

F=frameaccel(F,Ls);
L=F.L;

%% -- run the iteration

A=@(x) F.frsyn(F.frana(x));

% An explicit postpad is needed for the pcg algorithm to not fail
f=postpad(f,L);

if flags.do_pcg
    d=framediag(F,L);
    M=spdiags(d,0,L,L);

    [fout,flag,~,iter,relres]=pcg(A,f,kv.tol,kv.maxit,M);
else

    [fout,flag,~,iter,relres]=pcg(A,f,kv.tol,kv.maxit);
end;

c=F.frana(fout);

if nargout>1
    relres=relres/norm(fout(:));
end;

%% --- cleanup -----

permutedsize=[size(c,1),permutedsize(2:end)];

c=assert_sigreshape_post(c,dim,permutedsize,order);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function info = xif_read_header(filename)
% function for reading header of XIF HDllab/ATL ultrasound ( .xif )
% volume file
%
% info = xif_read_header(filename);
%
% examples:
% 1,  info=xif_read_header()
% 2,  info=xif_read_header('volume.xif');

if(exist('filename','var')==0)
    [filename, pathname] = uigetfile('*.xif', 'Read xif-file');
    filename = [pathname filename];
end

%Open file.
fid = fopen(filename, 'r');

%Abort if error opening file (does not exist, etc.)
if (fid == -1)
    error('Error opening file: %s, aborting...', filename);
end

% Store file name of volume
info.Filename=filename;

%Run through the XML header, reading in scanlines, samples etc.
s = [];
while (size(strfind(s, '#SOH')) == 0)
    s = fgets(fid);
    if (strfind(s, 'echoNumLines'))
       % <dataItem
       %   name="echoNumLines"
       %   type="INT"
       %   value="105"/>
        s = fgets(fid);
        s = fgets(fid);
        info.Scanlines = str2num(cell2mat(regexp(s, '(\d+)', 'match')));
    end
    if (strfind(s, 'echoNumDisplaySamples'))
        s = fgets(fid);
        s = fgets(fid);
        info.Samples = str2num(cell2mat(regexp(s, '(\d+)', 'match')));
    end
    if (strfind(s, 'numFrames'))
        s = fgets(fid);
        s = fgets(fid);
        info.Frames = str2num(cell2mat(regexp(s, '(\d+)', 'match')));
    end
    if (strfind(s, 'frameRate'))
        s = fgets(fid);
        s = fgets(fid);
        info.Framerate = str2num(cell2mat(regexp(s, '(\d+)', 'match')));
    end
    if (strfind(s, 'echoScan_linearWidth'))
        s = fgets(fid);
        s = fgets(fid);
        %Different regexp, as we need to match a decimal number.
        info.Width = str2num(cell2mat(regexp(s, '(\d+\.\d+)', 'match')));
    end
    if (strfind(s, 'echoDepth_twodDepth'))
        s = fgets(fid);
        s = fgets(fid);
        %Same as above.
        info.Height = str2num(cell2mat(regexp(s, '(\d+\.\d+)', 'match')));
    end
end
info.Dimensions=[info.Samples, info.Scanlines, info.Frames];
info.PixelDimensions=[0 0 0];
fclose(fid);

author: tomdoel

repo: pulmonarytoolkit


text:
%% Implementation of Progressive Switching Median Filter
%% Base Paper : Zhou Wang and David Zhang, "Progressive Switching Median
%% Filter for the Removal of Impulse Noise from Highly Corrupted Images",
%% IEEE Trans. on Cir. and Sys., vol. 46, no. 1, Jan. 1999.
%% Function Y = PSMF(x)
%% input    x = Image is corrupted by Salt & Pepper Noise
%%
%%  Example: Y = PSMF(x);
%%      Posted date   : 16 - 10 - 2008
%%      Modified date :
%%
%% Developed By : K.Kannan (kannan.keizer@gmail.com)
%%                & Jeny Rajan (jenyrajan@gmail.com)
%%                  Medical Imaging Research Group (MIRG), NeST,
%%                  Trivandrum.
%% Progressive Switching Median Filter
function Y = PSMF(x)
x = double(x);
WF = 3; ND = 3;T = 40;a = 65;b = -50;
M = medfilt2(x,[3 3]);
N = abs(x - M);
N(N>T)=0;
N = N ~= 0;
N = double(N);
R = sum(N(:))/(size(x,1) * size(x,2));
if R <= 0.25
    WD = 3;
else
    WD = 5;
end
TD = a + (b * R);
z = IMPDET(x,ND,WD,TD);
Y = NF(x,z,WF);

%% Impulse Detection
function F1 = IMPDET(x,ND,WD,TD)
X = x;
M = medfilt2(X,[WD WD]);
D = abs(X - M);
F = zeros(size(x));
F(D>=TD)=1;
F1 = F;
X(F1==F)=X(F1==F);
X(F1~=F)=M(F1~=F);
for i = 1:ND-1
    M = medfilt2(X,[WD WD]);
    F1(abs(X - M)<TD)=F(abs(X - M)<TD);
    F1((X - M)>=TD)=1;
    X(F1==F)=X(F1==F);
    X(F1~=F)=M(F1~=F);
    F = F1;
end
return;

%% Noise Filtering
function Y = NF(x,f,WF)
g = f;
Y = x;
Y1 = Y;
g1 = g;
s = sum(g(:));
while s ~= 0
    M = medfilt2(Y,[WF WF]);
    Y1(g==1)=M(g==1);
    g1(Y~=Y1)=0;
    Y = Y1;
    g = g1;
    s1 = sum(g(:));
    if s1 ~= s
        s = s1;
    else
        s = 0;
    end
end
return;

author: Sable


text:
%DETECTCHARUCODIAMOND  Detect ChArUco Diamond markers
%
%     [diamondCorners, diamondIds] = cv.detectCharucoDiamond(img, markerCorners, markerIds, squareMarkerLengthRate)
%     [...] = cv.detectCharucoDiamond(..., 'OptionName',optionValue, ...)
%
% ## Input
% * __img__ input image necessary for corner subpixel (8-bit grayscale or
%   color).
% * __markerCorners__ cell array of detected marker corners from
%   cv.detectMarkers function `{{[x,y],..}, ..}`.
% * __markerIds__ vector of marker ids in `markerCorners` (0-based).
% * __squareMarkerLengthRate__ rate between square and marker length:
%   `squareMarkerLengthRate = squareLength/markerLength`. The real units are
%   not necessary.
%
% ## Output
% * __diamondCorners__ output list of detected diamond corners (4 corners per
%   diamond). The order is the same than in marker corners: top left, top
%   right, bottom right and bottom left. Similar format than the corners
%   returned by cv.detectMarkers (e.g `{{[x y],..}, ..}`).
% * __diamondIds__ ids of the diamonds in `diamondCorners` (0-based). The id
%   of each diamond is in fact a vector of four integers, so each diamond has
%   4 ids, which are the ids of the aruco markers composing the diamond.
%
% ## Options
% * __CameraMatrix__ Optional 3x3 camera calibration matrix
%   `A = [fx 0 cx; 0 fy cy; 0 0 1]`.
% * __DistCoeffs__ Optional vector of camera distortion coefficients
%   `[k1,k2,p1,p2,k3,k4,k5,k6,s1,s2,s3,s4]` of 4, 5, 8 or 12 elements.
%
% This function detects Diamond markers from the previous detected ArUco
% markers. The diamonds are returned in the `diamondCorners` and `diamondIds`
% parameters. If camera calibration parameters are provided, the diamond
% search is based on reprojection. If not, diamond search is based on
% homography. Homography is faster than reprojection but can slightly reduce
% the detection rate.
%
% See also: cv.detectMarkers, cv.interpolateCornersCharuco,
%  cv.drawDetectedDiamonds
%


text:
function [h] = m_legend(varargin)
% M_LEGEND Add a legend to a map. This function does not have the
%         full functionality of LEGEND but avoids some of the problems
%         encountered using that function in M_MAP. Feel free to
%         add code!
%
%    [legend_handle] = M_LEGEND(HANDLES,String1,String2,String3,...)
%
%    puts a legend on the plot containing the handles in the vector HANDLES
%    using the specified strings as labels for the corresponding handles.
%    By holding the mouse button down, the legend can be moved around the
%    plot.

% Original Author:: Deirdre Byrne, dbyrne@umeoce.maine.edu 00/06/23

% This software is provided "as is" without warranty of any kind.

global MAP_PROJECTION MAP_VAR_LIST

if isempty(MAP_PROJECTION)
  disp('No Map Projection initialized - call M_PROJ first!');
  return;
end

if nargin < 2
  help m_legend
  return
end

handles = varargin{1};
handles = handles(:);
varargin = varargin(:);
s = size(varargin,1);
if s-1 ~= size(handles,1)
  disp('must have same number of legends as handles')
  return
end

varargin = varargin(2:s);

% get axis location
ax=gca;
units = get(ax,'units');
p = get(ax,'pos');
c = get(ax,'color');
if strcmp(c,'none')
  c = get(gcf,'color');
end
if strcmp(c,'none')
  c = [0.702 0.702 0.702];
end

% calculate initial size for legend
xlen = (p(3))/5;
ylen = (p(4))/7;

% calculate position for legend
% default is to put it in bottom right corner
axpos(1) = p(1)+p(3) - xlen;
axpos(2) = p(2);
axpos(3:4) = [xlen ylen];

ButtonDownFcn = 'moveaxis';
Tag = 'legend';

h=axes('units',units,'position',axpos, ...
    'visible','on', ...
    'color',c, ...
    'xcolor','k', ...
    'ycolor','k', ...
    'tag',Tag, ...
    'Buttondownfcn',ButtonDownFcn, ...
    'xlim',[0 1], ...
    'xtick',[], ...
    'ylim',[0 1], ...
    'ytick',[], ...
    'nextplot','add', ...
    'box','on');

lh = length(handles);

% set width based on legend strings
axw = 4;
for i = 1:lh
  l = length(varargin{i});
  axw = max(axw,l);
end
axw = max(2*axw,axw+10);

sf = 1;
set(h,'units','characters');
axpos = get(h,'pos');
axh = (sf*lh+1);
set(h,'position',[(axpos(1)+axpos(3) - axw) axpos(2) axw axh])

for i = 1:lh
  xp = [0.05 0.2 0.35];
  iy = (i*sf-1)*0.8/lh + 0.25;
  yp = [iy iy iy];
  if strcmp(get(handles(i),'type'),'line')
    mark = get(handles(i),'marker');
    msize = get(handles(i),'markersize');
    linest = get(handles(i),'linestyle');
    linew = get(handles(i),'linewidth');
    clr = get(handles(i),'color');
    if ~strcmp(linest,'none')
      l1(i) = plot(xp,yp,'marker','none', ...
	  'color',clr,'linewidth',linew);
    end
    l2(i) = plot(xp(2),yp(2),'marker',mark, ...
	'linewidth',linew, 'markersize', msize, ...
	'color',clr);
    t(i) = text(0.5, iy, varargin{i});
  end
end

if nargout == 0
  clear h
end


text:
function h = noise3dPlot(noise, plotType, CX, CY, CZ, CZVar, varargin)

% NOISE3DPLOT Draw a 3D or contour plot for the relevant noise model.
% FORMAT
% DESC draws a 3D or contour plot for the relevant noise model.
% ARG noise : the noise structure for which the plot is required.
% ARG plotType : string containing the name of the plotting function (for example mesh, contour).
% ARG X : the input X data in the form of a 'mesh' matrix.
% ARG Y : the input Y data in the form of a 'mesh' matrix.
% ARG mu : the input mean in the form of a 'mesh' matrix.
% ARG varSigma : the input variance in the form of a 'mesh' matrix.
% ARG P1, P2, P3 ... : optional additional arguments for the given plot type.
% RETURN h : the gradients of the diagonal with respect to each element
% of X. The returned matrix has the same dimensions as X.
%
% SEEALSO : noiseParamInit, noise3dPlot,
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005

% NOISE

functionName = [noise.type 'Noise3dPlot'];
if exist(functionName) == 2
  fhandle = str2func(functionName);
  h = fhandle(noise, plotType, CX, CY, CZ, CZVar, varargin{:});
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function x = p30_start ( n )

%*****************************************************************************80
%
%% P30_START returns a starting point for optimization for problem 30.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables X.
%
%    Output, real X(N), a starting point for the optimization.
%
  x = [ -1.0, 1.0 ]';

  return
end

repo: jburkardt-m


text:
function plotregionData(varargin)
%PLOTREGIONDATA   SINGFUN does not support PLOTREGIONDATA.

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

error('CHEBFUN:SINGFUN:plotregionData:notSupported', ['SINGFUN does not' ...
    'support PLOTREGIONDATA']);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function W = randInitializeWeights(L_in, L_out)
%RANDINITIALIZEWEIGHTS Randomly initialize the weights of a layer with L_in
%incoming connections and L_out outgoing connections
%   W = RANDINITIALIZEWEIGHTS(L_in, L_out) randomly initializes the weights
%   of a layer with L_in incoming connections and L_out outgoing
%   connections.
%
%   Note that W should be set to a matrix of size(L_out, 1 + L_in) as
%   the column row of W handles the "bias" terms
%

% You need to return the following variables correctly
W = zeros(L_out, 1 + L_in);

% ====================== YOUR CODE HERE ======================
% Instructions: Initialize W randomly so that we break the symmetry while
%               training the neural network.
%
% Note: The first row of W corresponds to the parameters for the bias units
%

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

end

author: vugsus

repo: coursera-machine-learning


text:
function p = lbfgsbprod(H,g)
%lbfgshprod  computes products with the L-BFGS matrix B.
%
%   p = lbfgsbprod(H,g)  returns  p = B*g.
%
%   The product is computed using the factoriation
%   [(9.15),p.231] described in Nocedal and Wright, 1999.

%   See also lbfgsadd, lbfgsdel, lbfgsupdate, lbfgsinit.

% $Id$

% ----------------------------------------------------------------------
% Explicit matrix formulation
% ----------------------------------------------------------------------

%jMax = H.jMax;
%jNew = H.jNew;
%jOld = H.jOld;

% |---|---|---|---|---|
%       ^   ^
%       |   |
%      old new
%   4   5   1   2   3    1=newest, 2=2nd newest,... 5=oldest

%if false
%   Sk = H.S(:,[jNew-1:-1:1,jMax:-1:jNew]);
%   Yk = H.Y(:,[jNew-1:-1:1,jMax:-1:jNew]);
%   Lk = (Sk' * Yk) .* (repmat((1:jMax)',1,jMax) > repmat((1:jMax),jMax,1)) ;
%   Dk = diag(sum(Sk.*Yk));

%   deltak = H.delta;

%   M = [deltak*Sk'*Sk, Lk; Lk', -Dk];

%   p = [deltak*Sk'; Yk'] * g;
%   p = M \ p;
%   p = [deltak*Sk, Yk] * p;
%   p = deltak * g - p;
%end

% ----------------------------------------------------------------------
% Formulation directly based on arrays in permuted form :-)
% ----------------------------------------------------------------------

if ~isempty(H.ML)

   % This code works much faster for larger problems
   v1 = H.delta * (g' * H.S)';
   v2 = (g' * H.Y)';
   p  = [v1(H.valid); v2(H.valid)];

%  v1 = H.delta * (g' * H.S(:,H.valid))';
%  v2 = (g' * H.Y(:,H.valid))';
%  q = [v1; v2];

   p = H.MU \ (H.ML \ p); % H.M \ p (optionally: use linsolve)

%  v1 = H.S(:,H.valid) * p(1:H.rank) * H.delta;
%  v2 = H.Y(:,H.valid) * p(H.rank+1:2*H.rank);
%  q = v1 + v2;

   % This code works much faster for larger problems
   pe = zeros(size(H.S,2),1); pe(H.valid) = p(1:H.rank);
   v1 = H.S * (pe * H.delta);
   pe = zeros(size(H.Y,2),1); pe(H.valid) = p(H.rank+1:2*H.rank);
   v2 = H.Y * pe;
   p  = v1 + v2;
else
   p = 0;
end

p = H.delta * g - p;

author: mpf


text:
%SETLABTYPE Reset label type of dataset
%
%    A = SETLABTYPE(A,TYPE,LABELS)
%
% The label type of the dataset A is converted to TYPE ('crisp','soft' or
% 'targets'). A conversion of the dataset fields 'nlab', 'lablist' and
% 'targets' is made where necessary. If given, LABELS replaces the labels
% or targets of A.
%
% EXAMPLE
% a = prdataset(rand(10,5)); % create dataset of 10 objects and 5 features
% a = setlabtype(a,'soft',rand(10,1)); % give it random soft labels

author: marianux

repo: ecg-kit

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