text:
function theVec=vech(H,offDiagScalFact,omitDiagonal)
%%VECH The vector half operator. Given a square 2D matrix H, this returns
%      the entries that are on or below the main diagonal of H, listed
%      column-by-column, unless omitDiagonal=true, in which case the main
%      diagonal is omitted too.
%
%INPUTS: H An nXn matrix. It can be real or complex.
% offDiagScalFact An optional scalar factor by which the off-diagonal
%          elements are scaled. If this parameter is omitted or an empty
%          matrix is passed, then the default value of 1 is used.
% omitDiagonal An optional parameter indicating whether the elements on the
%          main diagonal should be omitted. The default if omitted or an
%          empty matrix is passed is false.
%
%OUTPUTS: theVec A (n*(n+1)/2)X1 matrix holding the lower-triangular part
%                of the matrix, listed column-by-column. If omitDiagonal is
%                true, then this is a (n*(n+1)/2-n)X1 vector, as the
%                elements on the main diagonal were omitted.
%
%The opposite of this function is vech2Mat.
%
%December 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(omitDiagonal))
    omitDiagonal=false;
end

if(nargin<2||isempty(offDiagScalFact))
    offDiagScalFact=1;
end

n=size(H,1);

if(omitDiagonal)
    numEls=n*(n+1)/2-n;
    theVec=zeros(numEls,1);
    if(numEls==0)
        return;
    end
    curStart=1;
    for curCol=1:n
        %Off-diagonal elements
        span=(curStart):(curStart+n-curCol-1);
        theVec(span)=H((curCol+1):end,curCol)*offDiagScalFact;

        curStart=curStart+n-curCol;
    end
else
    theVec=zeros(n*(n+1)/2,1);

    curStart=1;
    for curCol=1:n
        %The diagonal element
        theVec(curStart)=H(curCol,curCol);

        %Off-diagonal elements
        span=(curStart+1):(curStart+n-curCol);
        theVec(span)=H((curCol+1):end,curCol)*offDiagScalFact;

        curStart=curStart+n-curCol+1;
    end
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%   AUTHORSHIP
%   Primary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause
%
function [x, y]=clockwise(xy)
x=xy(:,1);
y=xy(:,2);
cx = mean(x);
cy = mean(y);
a = atan2(y - cy, x - cx);
[~, order] = sort(a);
x = x(order);
y = y(order);
%assert(all(ismember([x y],xy,'rows')));
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function Mh = hmxFusion(Mh)
%+========================================================================+
%|                                                                        |
%|         OPENHMX - LIBRARY FOR H-MATRIX COMPRESSION AND ALGEBRA         |
%|           openHmx is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : hmxFusion.m                                   |
%|    #    |   VERSION    : 0.40                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.03.2018                                    |
%| ( === ) |   SYNOPSIS   : Fusion and recompress for full, sparse and    |
%|  `---'  |                low-rank leaves                               |
%+========================================================================+

%%% H-Matrix
if (Mh.typ == 0)
    % Check and convert leaf data
    typ = zeros(1,4);
    for i = 1:4
        % Children data
        typ(i) = Mh.chd{i}.typ;
        dat    = Mh.chd{i}.dat;
        tol    = Mh.chd{i}.tol;
        dim    = size(Mh.chd{i});

        % Full or sparse matrix
        if (typ(i) == 2)
            % Non zeros terms
            n = nnz(dat);

            % Sparse
            if issparse(dat)
                % Compress
                if (n == 0)
                    A    = zeros(dim(1),0);
                    B    = zeros(0,dim(2));
                    flag = 1;

                % Full
                elseif (n > 1/4*prod(dim))
                    dat  = full(dat);
                    flag = 2;

                % Sparse
                else
                    flag = 0;
                end

            % Full
            else
                % Compress
                if (n == 0)
                    A    = zeros(dim(1),0);
                    B    = zeros(0,dim(2));
                    flag = 1;

                % Sparse
                elseif (n <= 1/4*prod(dim))
                    dat  = sparse(double(dat));
                    flag = 2;

                % Compress
                else
                    [A,B,flag] = hmxSVD(dat,tol);
                end
            end

            % Update
            if (flag == 1)
                Mh.chd{i}.dat = {A,B};
                Mh.chd{i}.typ = 1;
                typ(i)        = 1;

            elseif (flag == 2)
                Mh.chd{i}.dat = dat;
            end
        end
    end

    % Low-rank fusion (QRSVD)
    if (sum(typ==1) == 4)
        % Rank for each leaf
        rk = zeros(1,4);
        nk = 0;
        for i = 1:4
            rk(i) = size(Mh.chd{i}.dat{1},2);
            nk    = nk + sum(size(Mh.chd{i}))*rk(i);
        end

        % Low rank matrix
        A = zeros(size(Mh,1),sum(rk),class(Mh.pos{1}));
        B = zeros(sum(rk),size(Mh,2),class(Mh.pos{2}));
        j = 0;
        for i = 1:4
            A(Mh.row{i},j+(1:rk(i))) = Mh.chd{i}.dat{1};
            B(j+(1:rk(i)),Mh.col{i}) = Mh.chd{i}.dat{2};
            j = j + rk(i);
        end

        % Recompression
        [A,B] = hmxQRSVD(A,B,Mh.tol);

        % Update
        if (sum(size(Mh))*size(A,2) <= nk)
            Mh.typ = 1;
            Mh.row = cell(1,4);
            Mh.col = cell(1,4);
            Mh.chd = cell(1,4);
            Mh.dat = {A,B};
        end
    end

%%% Unavalaible case
else
    error('hmxFusion.m : unavailable case')
end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function version = om_getgccversion
% checks for gcc compiler version (works if superior to gcc3)
% Copyright (C) 2010-2017, OpenMEEG developers

tmpdir = pwd;
cd /tmp
[junk,tname] = fileparts(tempname);
txtfile  = [tname '.txt'];
dos(['gcc -v >& ' txtfile]);
efid = fopen(txtfile);

tmp = ''; cnt = 1;
vec = [];
while ~isnumeric(tmp)
    tmp = fgetl(efid);
    vec{cnt} = tmp;
    cnt = cnt+1;
end
fclose(efid);
delete(txtfile);
cd(tmpdir);
tmp = deblank(vec{cnt-2});
num = findstr('gcc version ',tmp);
version = str2num(tmp(num+11:num+12));
if isempty(version) && findstr('clang',vec{2}) > 0
    version = 5;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function vw = makeROIfromSelectedVoxels(vw,name,select,color, comments)
%
% vw = makeROIfromSelectedVoxels(vw,[name],[select],[color], [comments])
%
% Makes an ROI consisting of all coords according to cothresh. phWindow,
% and mapWindow.
%
% name:     name (string) for the ROI
% select:   if non-zero, chooses the new ROI as the selectedROI (default=1)
% color:    sets color for drawing the ROI (default 'b').
% comments: test string.

% jw, 1/16/2010

mrGlobals;

if notDefined('vw'),        vw      = getCurView;   end
if notDefined('name'),      name    = [];           end
if notDefined('select'),    select  = [];           end
if notDefined('color'),     color   = [];           end
if notDefined('comments'),  comments= [];           end

% Create an ROI with all coordinates
coords = allCoords(vw);
vw = newROI(vw,name,select,color,coords, comments);

% Restrict the coordinates according to cothresh. phWindow, and mapWindow.
vw = restrictROIfromMenu(vw);

return

function coords = allCoords(vw)

viewType = viewGet(vw, 'viewType');

switch lower(viewType)
    case {'gray', 'volume'}
        coords = viewGet(vw, 'coords');
    case 'inplane'
        sz =  viewGet(vw,'anatsize');
        [a b c] = ind2sub(sz, 1:prod(sz));
        coords = [a; b; c];
    case 'flat'
        leftcoords  = viewGet(vw, 'coords', 'left');
        rightcoords = viewGet(vw, 'coords', 'right');
        coords      = [leftcoords rightcoords];
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [model] = createPeriplasmaticSpace(model,microbeID,infoFile)
% Part of the DEMETER pipeline. This function creates a periplasmatic space
% for refined reconstructions if it is appropriate for the organism. The
% periplasmatic space is created by by retrieving all extracellular
% metabolites and adding a third compartment.
%
% USAGE
%       [model] = createPeriplasmaticSpace(model,microbeID,infoFile)
%
% INPUT
% model             COBRA model structure
% microbeID:        ID of the reconstructed microbe that serves as the
%                   reconstruction name and to identify it in input tables
% infoFile:         Table with taxonomic and gram staining information on
%                   microbes to reconstruct
%
% OUTPUT
% model             COBRA model structure
%
% AUTHOR:
%       - Almut Heinken, 03/2020

% get the information on taxonomy and gram staining to find out if a periplasmatic space
% should be added
phylCol=find(strcmp(infoFile(1,:),'Phylum'));
if ~isempty(find(strcmp(infoFile(:,1),microbeID)))
    phylum=infoFile{find(strcmp(infoFile(:,1),microbeID)),phylCol};
    genCol=find(strcmp(infoFile(1,:),'Genus'));
    genus=infoFile{find(strcmp(infoFile(:,1),microbeID)),genCol};

    [~,gramCol,~]=intersect(infoFile(1,:),{'Gram Staining','Gram_Staining'});
    gramStatus=infoFile(find(strcmp(infoFile(:,1),microbeID)),gramCol);
    if length(gramStatus)>1
        gramStatus = gramStatus{1};
    end
        if (strcmp(gramStatus,'Gram-') || strcmp(phylum,'Deinococcus-Thermus')) && ~any(strcmp(phylum,{'Euryarchaeota','Crenarchaeota','Thaumarchaeota','Tenericutes','Chloroflexi','Chloroflexota'})) && ~any(strcmp(genus,{'Acidaminobacter','Gracilibacter','Intestinibacillus','Mageeibacillus','Massiliomicrobiota','Flavonifractor','Pseudoflavonifractor'}))

            % get all extracellular metabolites
            exMets=model.mets(find(contains(model.mets,'[e]')));

            % Add periplasmatic metabolites, metNames and metFormulas
            metNames = {};
            metFormulas = {};
            metCharges = [];

            for i=1:length(exMets)
                pMets{i} = strrep(exMets{i},'[e]','[p]');
                metNames{i} = model.metNames{find(strcmp(model.mets,exMets{i}))};
                metFormulas{i} = model.metFormulas{find(strcmp(model.mets,exMets{i}))};
                metCharges(i) = model.metCharges(find(strcmp(model.mets,exMets{i})));
            end

            modelNew=model;
            % Add all new periplasmatic metabolites
            pMets=setdiff(pMets,model.mets);
            for i=1:length(pMets)
                modelNew = addMetabolite(modelNew, pMets{i}, 'metName', metNames{i}, 'metFormula', metFormulas{i}, 'Charge', metCharges(i));
            end

            % Get all transport reactions associated with the exchange metabolites
            rxnsToAdd = {};
            rxnNames = {};
            rxnFormulas = {};
            rxnLB = [];
            rxnUB = [];
            rxnSubsystem = {};

            % convert extracellular transport reactions to transport reactions from
            % periplasmatic space to cytosol
            for i=1:length(exMets)
                transpRxns = findRxnsFromMets(model,exMets{i});
                % remove exchange reactions
                transpRxns(find(strncmp(transpRxns,'EX_',3)))=[];
                for j=1:length(transpRxns)
                    rxnsToAdd{end+1} = [transpRxns{j} 'pp'];
                    rxnNames{end+1} = [model.rxnNames{find(strcmp(model.rxns,transpRxns{j}))} ', periplasmatic'];
                    rxnLB(end+1) = model.lb(find(strcmp(model.rxns,transpRxns{j})));
                    rxnUB(end+1) = model.ub(find(strcmp(model.rxns,transpRxns{j})));
                    form = printRxnFormula(model,transpRxns{j});
                    rxnFormulas{end+1} = strrep(form{1},'[e]','[p]');
                    rxnSubsystem{end+1} = 'Transport, periplasmatic';
                end
                modelNew = removeRxns(modelNew,transpRxns, 'metFlag', false);
            end

            % add transport reactions to transport reactions from
            % extracellular to periplasmatic space
            for i=1:length(exMets)
                rxnsToAdd{end+1} = [upper(strrep(exMets{i},'[e]','')) 'tex'];
                rxnNames{end+1} = [model.metNames{find(strcmp(model.mets,exMets{i}))} ' diffusion extracellular to periplasm'];
                rxnLB(end+1) = -1000;
                rxnUB(end+1) = 1000;
                rxnFormulas{end+1} = [exMets{i} ' <=> ' strrep(exMets{i},'[e]','[p]')];
                rxnSubsystem{end+1} = 'Transport, extracellular';
            end

            % add all new reactions
            for i=1:length(rxnsToAdd)
                modelNew = addReaction(modelNew,rxnsToAdd{i},'reactionName',rxnNames{i},...
                    'reactionFormula',rxnFormulas{i},'lowerBound',rxnLB(i),'upperBound',rxnUB(i),'subSystem',rxnSubsystem{i});
            end

            model = convertOldStyleModel(modelNew);
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function K = curvature_2D(X,Y,degree)
% CURVATURE_2D estimate curvature of a curve defined by points
%
%   K = CURVATURE_2D(X, Y, DEGREE)
%
%   DEGREE is the degree of the polynom used for the interpolation
%

%
%   Created by Alexandre Gramfort on 2008-12-15.
%   Copyright (c) 2007 Alexandre Gramfort. All rights reserved.
%

% $Id: curvature_2D.m 2 2009-06-16 19:24:10Z gramfort $
% $LastChangedBy: gramfort $
% $LastChangedDate: 2009-06-16 15:24:10 -0400 (Mar, 16 jui 2009) $
% $Revision: 2 $

if nargin<3
    % default values
    degree = 5; % degree of polynum used for interpolation
end

% Find a parametrization of the curve
t = cumsum(sqrt(diff(X).^2 + diff(Y).^2));
t = [0,t];
t = t / t(end); % normalize between 0 and 1

% compute coefficients of interpolation functions
x0 = polyfit(t, X, degree);
y0 = polyfit(t, Y, degree);

% compute coefficients of first and second derivatives. In the case of a
% polynom, it is possible to compute coefficient of derivative by
% multiplying with a matrix.
derive = diag(degree:-1:0);
xp = circshift(x0*derive, [0 1]);
yp = circshift(y0*derive, [0 1]);
xs = circshift(xp*derive, [0 1]);
ys = circshift(yp*derive, [0 1]);

% compute values of first and second derivatives for needed points
xprime = polyval(xp, t);
yprime = polyval(yp, t);
xsec = polyval(xs, t);
ysec = polyval(ys, t);

% compute value of curvature
K = (xprime.*ysec - xsec.*yprime)./ ...
    power(xprime.*xprime + yprime.*yprime, 3/2);


text:
% Count the number of self-loops in the graph
%
% INPUT: adjacency matrix, nxn
% OUTPUT: integer, number of self-loops
%
% Note: in the adjacency matrix representation loops appear as non-zeros on the diagonal
% GB: last updated, Sep 20 2012

function sl=selfLoops(adj)

sl=sum(diag(adj));

author: aeolianine

repo: octave-networks-toolbox


text:
function grid = getGrid(~, N, dom)
%GETGRID   Returns a grid correspoding to a SPINOP3 object.
%
% See also SPINOP3.

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

xx = trigpts(N, dom(1:2));
yy = trigpts(N, dom(3:4));
zz = trigpts(N, dom(5:6));
[xx, yy, zz] = meshgrid(xx, yy, zz);
grid = {xx; yy; zz};

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ f, i ] = fibonacci_floor ( n )

%*****************************************************************************80
%
%% FIBONACCI_FLOOR returns the largest Fibonacci number less than or equal to N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the positive integer whose Fibonacci "floor" is desired.
%
%    Output, integer F, the largest Fibonacci number less than or equal to N.
%
%    Output, integer I, the index of the F.
%
  if ( n <= 0 )

    i = 0;
    f = 0;

  else

    i = floor ( log ( 0.5 * ( 2 * n + 1 ) * sqrt ( 5.0 ) ) ...
      / log ( 0.5 * ( 1.0 + sqrt ( 5.0 ) ) ) );

    f = fibonacci_direct ( i );

    if ( n < f )
      i = i - 1;
      f = fibonacci_direct ( i );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function img = imMergeLabels(img, lbls, varargin)
%IMMERGELABELS Merge regions in a labeled image
%
%   Deprecated: replaced by the "imMergeRegions" function.
%
%   Usage:
%   LBL2 = imMergeLabels(IMG, LABELS);
%   IMG is a label image, LABELS are the labels of the regions to be
%   merged. IMG can be 2D or 3D.
%   The indices of regions are set to the same value, and the boundary
%   between the regions is also merged to the new region, preserving
%   topology of adjacent regions.
%   The function uses morphologial operation, processing time can be
%   consuming.
%
%   LBL2 = imMergeLabels(IMG, LABELS, SE);
%   Specify the structuring element used for morphological detection of the
%   boundary.
%
%   Example:
%   lbl = [1 1 0 3 3;1 1 0 3 3;1 0 2 0 3;0 2 2 2 0];
%   lbl2 = imMergeLabels(lbl, [1 2])
%   lbl2 =
%        1     1     0     3     3
%        1     1     0     3     3
%        1     1     1     0     3
%        1     1     1     1     0
%
%   See Also
%     imMergeRegions
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2007-08-07
% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

% simply call the new function
img = imMergeRegions(img, lbls, varargin{:});

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [node,elem] = cuboidmesh(cuboid,h)
%% CUBEHEXMESH uniform mesh of cuboid
%
% [node,elem] = cuboidmesh([x0,x1,y0,y1,z0,z1],h) generates a uniform mesh of the
% cuboid [x0,x1]*[y0,y1]*[z0,z1] with mesh size h.
%
% Example
%
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

x0 = cuboid(1); x1 = cuboid(2);
y0 = cuboid(3); y1 = cuboid(4);
z0 = cuboid(5); z1 = cuboid(6);

[x,y,z] = meshgrid(x0:h:x1,y0:h:y1,z0:h:z1);

node = [x(:),y(:),z(:)];

%% Generate elements
nx = size(x,1) - 1; % number of cells in x-direction
ny = size(y,2) - 1; % number of cells in y-direction
nz = size(z,2) - 1; % number of cells in z-direction

elem = zeros(nx*ny*nz,8);
cellidx = 1:nx*ny*nz;

[i, j, k] = ind2sub([nx ny nz],cellidx);  % index of cells in subscript form

s =[ nx+1, ny+1, nz+1];

elem(cellidx,1) = sub2ind3(s,i,j,k);
elem(cellidx,2) = sub2ind3(s,i,j+1,k);
elem(cellidx,3) = sub2ind3(s,i+1,j+1,k);
elem(cellidx,4) = sub2ind3(s,i+1,j,k);
elem(cellidx,5) = sub2ind3(s,i,j,k+1);
elem(cellidx,6) = sub2ind3(s,i,j+1,k+1);
elem(cellidx,7) = sub2ind3(s,i+1,j+1,k+1);
elem(cellidx,8) = sub2ind3(s,i+1,j,k+1);

    function idx = sub2ind3(siz,i,j,k)
        nr = siz(1); nc = siz(2); nv = siz(3);
        if (max(j)>nc) || (max(i)>nr) || (max(k)>nv)
            error(message('MATLAB:mysub2ind:IndexOutOfRange'));
        end
        idx = (k-1)*nr*nc + (j-1)*nr+i;
    end

end

repo: ifem


text:
function y = r8vec_sqstb ( n, x )

%*****************************************************************************80
%
%% R8VEC_SQSTB computes a "slow" quarter sine transform backward of an R8VEC.
%
%  Discussion:
%
%    This routine is provided for illustration and testing.  It is inefficient
%    relative to optimized routines that use fast Fourier techniques.
%
%    For 0 <= I <= N-1,
%
%      Y(I) = -2 Sum ( 1 <= J <= N-1 ) X(J) * sin ( PI * J * (I+1/2) / N )
%             - X(N) * cos ( pi * I )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    William Briggs, Van Emden Henson,
%    The Discrete Fourier Transform,
%    SIAM, 1995,
%    QA403.5 B75
%
%  Parameters:
%
%    Input, integer N, the number of data values.
%
%    Input, real X(N), the data sequence.
%
%    Output, real Y(N), the transformed data.
%
  y(1:n) = 0.0;

  for i = 1 : n

    for j = 1 : n - 1
      theta = 0.5 * pi * j * ( 2 * i - 1 ) / n;
      y(i) = y(i) - 2.0 * x(j) * sin ( theta  );
    end

    theta = pi * ( i - 1 );
    y(i) = y(i) - x(n) * cos ( theta );

  end

  return
end


text:
%ISCOMDSET Test whether datasets are compatible
%
% 	N = ISCOMDSET(A,B,CLAS);
%
% INPUT
%   A    Input argument, to be tested on dataset
%   B    Input argument, to be tested on compatibility with A
%   CLAS 1/0, test on equal classes (1) or don't test (0)
%        (optional; default 1)
%
% OUTPUT
%   N    1/0 if A and B are / are not compatible datasets
%
% DESCRIPTION
% The function ISCOMDSET tests whether A and B are compatible
% datasets, i.e. have the same features and the same classes.
%
% SEE ALSO (PRTools Guide)
% ISDATASET, ISMAPPING, ISDATAIM, ISFEATIM, ISVALDFILE, ISVALDSET

function n = iscomdset(a,b,clas)

	if nargin < 3, clas = 1; end

	if isempty(b) % return of second dataset empty (i.e. not supplied)
		return
	end

	if nargout == 0
		isdataset(a);
		isdataset(b);
	else
		n = isdataset(a) & isdataset(b);
	end

	featlaba = setstr(getfeatlab(a));
	featlabb = setstr(getfeatlab(b));

	nf = strcmp(featlaba,featlabb);

	if ~nf
		if nargout == 0
			error([newline 'Datasets for training and testing/tuning should' newline ...
					'have the same features in the same order.'])
		else
			n = 0;
		end
	end

	if clas
		lablista = setstr(getlablist(a));
		lablistb = setstr(getlablist(b));

		no = strcmp(lablista,lablistb);

		if ~no
			if nargout == 0
				error([newline 'Datasets for training and testing/tuning should' newline ...
					'have the same class labels in the same order.'])
			else
				n = 0;
			end
		end
	end

return

author: marianux

repo: ecg-kit

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


text:
function [beta_gibbs,F_gibbs,gamma_gibbs,L_gibbs,phi_gibbs,sigma_gibbs,lambda_t_gibbs,sigma_t_gibbs,sbar,favar,It,Bu]=...
    favar_stvol2gibbs(Xbart,yt,beta0,omega0,alpha0,delta0,gamma0,zeta0,f0,upsilon0,betahat,sigmahat,I_o,omega,T,n,q,It,Bu,pick,pickf,favar,data_endo,lags)

%% preliminaries
% initialise variables
nfactorvar=favar.nfactorvar;
numpc=favar.numpc;
favarX=favar.X(:,favar.plotX_index);
onestep=favar.onestep;
Sigma=bear.nspd(favar.Sigma);
Ll=favar.L;
favar_X=favar.X;
% load priors
L0=favar.L0*eye(n);
a0=favar.a0;
b0=favar.b0;
% sigmahat=(1/T)*(EPS'*EPS);

% preallocation
Ll_gibbs=zeros(size(Ll(:),1),It-Bu);
R2_gibbs=zeros(size(favarX,2),It-Bu);

if onestep==0 %static factors in this case
    FY=data_endo;
    pbstring='two-step'; %string for the progress bar
    % elseif onestep==1
    %     pbstring='one-step'; %string for the progress bar
end

% preliminary elements for the algorithm
% compute alphabar
alphabar=T+alpha0;

% initiate the Gibbs sampler
% initiate the counting of iterations
count=1;
pickcount=1;
% initiate the record matrices and cells
beta_gibbs=[];
F_gibbs=[];
L_gibbs=[];
phi_gibbs=[];
sigma_gibbs=[];
lambda_t_gibbs={};
sigma_t_gibbs={};

% step 1: determine initial values for the algorithm

% initial value for beta
beta=betahat;
% initial value for f_2,...,f_n
% obtain the triangular factorisation of sigmahat
[Fhat,Lambdahat]=bear.triangf(sigmahat);
% obtain the inverse of Fhat
[invFhat]=bear.invltod(Fhat,n);
% create the cell storing the different vectors of invF
Finv=cell(n,1);
% store the vectors
for ii=2:n
    Finv{ii,1}=invFhat(ii,1:ii-1);
end
% initial values for L_1,...,L_n
L=zeros(T,n);
% initial values for gamma_1,...,gamma_n
gamma=0.85*ones(1,n);
% initial values for G_1,...,G_n
G=cell(n,1);
for ii=1:n
    G{ii,1}=speye(T)-sparse(diag(gamma(1,ii)*ones(T-1,1),-1));
end
% initial values for phi_1,...,phi_n
phi=ones(1,n);

% step 2: determine the sbar values and Lambda
sbar=diag(Lambdahat);
Lambda=sparse(diag(sbar));

% step 3: recover the series of initial values for lambda_1,...,lambda_T and sigma_1,...,sigma_T
lambda_t=repmat(diag(sbar),1,1,T);
sigma_t=repmat(sigmahat,1,1,T);

% create a progress bar
hbar = bear.parfor_progressbar(It,['Progress of the Gibbs sampler (',pbstring,').']);

%% run the Gibbs sampler
while count<=It

    % step 4: draw beta from its conditional posterior
    % first compute the summations required for omegabar and betabar
    summ1=zeros(q,q);
    summ2=zeros(q,1);
    % run the summation
    for jj=1:T
        prodt=Xbart{jj,1}'/sigma_t(:,:,jj);
        summ1=summ1+prodt*Xbart{jj,1};
        summ2=summ2+prodt*yt(:,:,jj);
    end
    % then obtain the inverse of omega0
    invomega0=diag(1./diag(omega0));
    % obtain the inverse of omegabar
    invomegabar=summ1+invomega0;
    % recover omegabar
    C=chol(bear.nspd(invomegabar),'Lower')';
    invC=C\speye(q);
    omegabar=invC*invC';
    % recover betabar
    betabar=omegabar*(summ2+invomega0*beta0);
    % finally, draw beta from its posterior
    beta=betabar+chol(bear.nspd(omegabar),'lower')*randn(q,1);

    % step 5: draw the series f_2,...,f_n from their conditional posteriors
    % recover first the residuals
    for jj=1:T
        epst(:,:,jj)=yt(:,:,jj)-Xbart{jj,1}*beta;
    end
    % then draw the vectors in turn
    for jj=2:n
        % first compute the summations required for upsilonbar and fbar
        summ1=zeros(jj-1,jj-1);
        summ2=zeros(jj-1,1);
        % run the summation
        for kk=1:T
            prodt=epst(1:jj-1,1,kk)*exp(-L(kk,jj));
            summ1=summ1+prodt*epst(1:jj-1,1,kk)';
            summ2=summ2+prodt*epst(jj,1,kk)';
        end
        summ1=(1/sbar(jj,1))*summ1;
        summ2=(-1/sbar(jj,1))*summ2;
        % then obtain the inverse of upsilon0
        invupsilon0=diag(1./diag(upsilon0{jj,1}));
        % obtain upsilonbar
        invupsilonbar=summ1+invupsilon0;
        C=chol(bear.nspd(invupsilonbar));
        invC=C\speye(jj-1);
        upsilonbar=full(invC*invC');
        % recover fbar
        fbar=upsilonbar*(summ2+invupsilon0*f0{jj,1});
        % finally draw f_i^(-1)
        Finv{jj,1}=fbar+chol(bear.nspd(upsilonbar),'lower')*randn(jj-1,1);
    end
    % recover the inverse of F
    invF=eye(n);
    for jj=2:n
        invF(jj,1:jj-1)=Finv{jj,1};
    end
    % eventually recover F
    F=bear.invltod(invF,n);
    % then update sigma
    sigma=F*Lambda*F';

    % step 6: draw the series gamma_1,...,gamma_n from their conditional posteriors
    % draw the parameters in turn
    for jj=1:n
        % estimate zetabar
        zetabar=1/((1/phi(1,jj))*L(1:T-1,jj)'*L(1:T-1,jj)+1/zeta0);
        % estimate zetabar
        gammabar=zetabar*((1/phi(1,jj))*L(2:T,jj)'*L(1:T-1,jj)+gamma0/zeta0);
        % draw the value gamma_i
        gamma(1,jj)=gammabar+zetabar^0.5*randn;
        % obtain G_i
        G{jj,1}=speye(T)-sparse(diag(gamma(1,jj)*ones(T-1,1),-1));
    end

    % step 7: draw the series phi_1,...,phi_n from their conditional posteriors
    % draw the parameters in turn
    for jj=1:n
        % estimate deltabar
        deltabar=L(:,jj)'*G{jj,1}'*I_o*G{jj,1}*L(:,jj)+delta0;
        % draw the value phi_i
        phi(1,jj)=bear.igrandn(alphabar/2,deltabar/2);
    end

    % step 8: draw the series lambda_i,t from their conditional posteriors, i=1,...,n and t=1,...,T
    % consider variables in turn
    for jj=1:n
        % consider periods in turn
        for kk=1:T
            % a candidate value will be drawn from N(lambdabar,phibar)
            % the definitions of lambdabar and phibar varies with the period, thus define them first
            % if the period is the first period
            if kk==1
                lambdabar=(gamma(1,jj)*L(2,jj))/(1/omega+gamma(1,jj)^2);
                phibar=phi(1,jj)/(1/omega+gamma(1,jj)^2);
                % if the period is the final period
            elseif kk==T
                lambdabar=gamma(1,jj)*L(T-1,jj);
                phibar=phi(1,jj);
                % if the period is any period in-between
            else
                lambdabar=(gamma(1,jj)/(1+gamma(1,jj)^2))*(L(kk-1,jj)+L(kk+1,jj));
                phibar=phi(1,jj)/(1+gamma(1,jj)^2);
            end
            % now draw the candidate
            cand=lambdabar+phibar^0.5*randn;
            % compute the acceptance probability
            prob=bear.mhprob2(jj,cand,L(kk,jj),sbar(jj,1),epst(:,1,kk),Finv{jj,1});
            % draw a uniform random number
            draw=rand;
            % keep the candidate if the draw value is lower than the prob
            if draw<=prob
                L(kk,jj)=cand;
                % if not, just keep the former value
            end
        end
    end
    % then recover the series of matrices lambda_t and sigma_t
    for jj=1:T
        lambda_t(:,:,jj)=diag(sbar).*diag(exp(L(jj,:)));
        sigma_t(:,:,jj)=F*lambda_t(:,:,jj)*F';
    end

    %% Sample Sigma and Ll (static)
    [Sigma,Ll]=bear.favar_SigmaL(Sigma,Ll,nfactorvar,numpc,onestep,n,favar_X,FY,a0,b0,T,lags,L0);

    %% record phase
    % if the burn-in sample phase is not yet over
    if count<=Bu
        % simply add 1 to the iteration count
        count=count+1;
        % on the other hand, if the burn-in sample phase is over
    elseif count>Bu
        % adding one iteration to the count will depend on wether post-burn selection applies
        % if there is no post burn selection
        if pick==0
            % record the results
            beta_gibbs(:,count-Bu)=beta;
            F_gibbs(:,:,count-Bu)=F;
            gamma_gibbs(count-Bu,:)=gamma;
            L_gibbs(:,:,count-Bu)=L;
            phi_gibbs(count-Bu,:)=phi;
            sigma_gibbs(:,count-Bu)=sigma(:);
            % save the factors and loadings (keep the notation in the code consistent, although - except L - they don't change)
            Ll_gibbs(:,count-Bu)=Ll(:);

            % compute R2 (Coefficient of Determination) for plotX variables (keep the notation in the code consistent, although they don't change)
            R2=bear.favar_R2(favarX,FY);
            R2_gibbs(:,count-Bu)=R2(:);
            for jj=1:T
                lambda_t_gibbs{jj,1}(:,:,count-Bu)=lambda_t(:,:,jj);
                sigma_t_gibbs{jj,1}(:,:,count-Bu)=sigma_t(:,:,jj);
            end
            % then add one to the count
            count=count+1;
            % if there is post burn selection, only one draw over 'fpick' draws will be retained
        elseif pick==1
            % if the iteration does not correspond to fpick, don't record the results, don't increase the regular count, but do increase pickcount by 1, and do record the acceptance rate of the Metropolis-Hastings step
            if pickcount~=pickf
                pickcount=pickcount+1;
                % on the other hand, if the iteration does correspond to fpick
            elseif pickcount==pickf
                % do record the results
                beta_gibbs(:,count-Bu)=beta;
                F_gibbs(:,:,count-Bu)=F;
                gamma_gibbs(count-Bu,:)=gamma;
                L_gibbs(:,:,count-Bu)=L;
                phi_gibbs(count-Bu,:)=phi;
                sigma_gibbs(:,count-Bu)=sigma(:);
                % save the factors and loadings (keep the notation in the code consistent, although - except L - they don't change)
                Ll_gibbs(:,count-Bu)=Ll(:);

                % compute R2 (Coefficient of Determination) for plotX variables (keep the notation in the code consistent, although they don't change)
                R2=bear.favar_R2(favarX,FY);
                R2_gibbs(:,count-Bu)=R2(:);
                for jj=1:T
                    lambda_t_gibbs{jj,1}(:,:,count-Bu)=lambda_t(:,:,jj);
                    sigma_t_gibbs{jj,1}(:,:,count-Bu)=sigma_t(:,:,jj);
                end
                % then increase the regular count by 1 and re-initialise pickcount
                count=count+1;
                pickcount=1;
            end
        end
    end

    % update progress by one iteration
    hbar.iterate(1);

end

% in case we have thinning of the draws,
thin=abs(round(favar.thin)); % should be a positive integer
if thin~=1
    beta_gibbs=beta_gibbs(:,thin:thin:end);
    F_gibbs=F_gibbs(:,thin:thin:end);
    L_gibbs=L_gibbs(:,thin:thin:end);
    phi_gibbs=phi_gibbs(:,thin:thin:end);
    sigma_gibbs=sigma_gibbs(:,thin:thin:end);
    Ll_gibbs=Ll_gibbs(:,thin:thin:end);
    R2_gibbs=R2_gibbs(:,thin:thin:end);
    for jj=1:T
        lambda_t_gibbs=lambda_t_gibbs(:,:,thin:thin:end);
        sigma_t_gibbs=sigma_t_gibbs(:,:,thin:thin:end);
    end
    It=(1/thin)*It;
    Bu=(1/thin)*Bu;
end

% save in favar structure
favar.L_gibbs=Ll_gibbs;
favar.R2_gibbs=R2_gibbs;

close(hbar);   %close progress bar

author: european-central-bank

repo: BEAR-toolbox


text:
% EEGLAB Toolbox to process EEG data
% Version - see eeg_getversion()
%
% eeglab() - Matlab graphic user interface environment for
%   electrophysiological data analysis incorporating the ICA/EEG toolbox
%   (Makeig et al.) developed at CNL / The Salk Institute, 1997-2001.
%   Released 11/2002- as EEGLAB (Delorme, Makeig, et al.) at the Swartz Center
%   for Computational Neuroscience, Institute for Neural Computation,
%   University of California San Diego (http://sccn.ucsd.edu/).
%   User feedback welcome: email eeglab@sccn.ucsd.edu
%
% Authors: Arnaud Delorme, Scott Makeig, with substantial contributions
%   from Colin Humphries, Sigurd Enghoff, Tzyy-Ping Jung, plus contributions
%   from Tony Bell, Te-Won Lee, Luca Finelli and other contributors.
%
% Description:
%   EEGLAB is Matlab software for processing continuous or event-related
%   EEG or other physiological data. It is designed for use by both novice and
%   expert Matlab users. In normal use, the EEGLAB graphic interface calls
%   graphic functions via pop-up function windows. The EEGLAB history mechanism
%   can save the resulting Matlab calls to disk for later incorporation into
%   Matlab scripts.  A single data structure ('EEG') containing all dataset
%   parameters may be accessed and modified directly from the Matlab commandline.
%   EEGLAB now recognizes "plugins," sets of EEGLAB functions linked to the EEGLAB
%   main menu through an "eegplugin_[name].m" function (Ex. >> help eeplugin_besa.m).
%

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [r0,u0,p0] = Euler_IC1d(x,input)
% Load the IC of a classical 2D Riemann Problem configuration.
% By Manuel Diaz 2012.10.24.
% In the notation we take advantage of the matlab array notation as follows
%
% prop = [prop_left , prop_right]
%
% Notation:
% u   = Velocity in x direction
% p   = Pressure
% rho = Density
% r   = Fugacity
% E   = Enerty
% t   = temperature
%
% Based on:
% http://wonka.physics.ncsu.edu/pub/VH-1/testpage/ and
% http://sitemaker.umich.edu/anand/files/riemann_shock-tube.pdf
% See also my routine CFD/Riemann.m to compute exact solutions.
%
%% Initial Physical Properties per case:
switch input
    case{1} % Configuration 1, Sod's Problem
        fprintf('Case 1: Sods problem \n');
        p   = [1    0.1  ];
        u   = [0.75 0    ];
        rho = [1    0.125];

    case{2} % Configuration 2, Left Expansion and right strong shock
        fprintf('Case 2: Left Expansion and right strong shock \n');
        p   = [1000 0.1  ];
        u   = [0    0    ];
        rho = [3    2    ];

    case{3} % Configuration 3, Right Expansion and left strong shock
        fprintf('Case 3: Right Expansion and left strong shock \n');
        p   = [7    10   ];
        u   = [0    0    ];
        rho = [1    1    ];

    case{4} % Configuration 4, Double Shock
        fprintf('Case 4: Double Shock \n');
        p   = [450  45   ];
        u   = [20   -6   ];
        rho = [6    6    ];

    case{5} % Configuration 5, Double Expansion
        fprintf('Case 5: Double Expansion \n');
        p   = [40   40   ];
        u   = [-2   2    ];
        rho = [1    2.5  ];

    case{6} % Configuration 6, Cavitation
        fprintf('Case 6: Cavitation \n');
        p   = [0.4  0.4  ];
        u   = [-20  20   ];
        rho = [1    1    ];

    otherwise
        error('Available cases: 1, 2, 3, 4, 5 and 6');

end
% Print
fprintf('\n');
fprintf('density (L): %1.3f\n',rho(1));
fprintf('velocity(L): %1.3f\n',u(1));
fprintf('Presure (L): %1.3f\n',p(1));
fprintf('\n');
fprintf('density (R): %1.3f\n',rho(2));
fprintf('velocity(R): %1.3f\n',u(2));
fprintf('Presure (R): %1.3f\n',p(2));
fprintf('\n');

%% Load Selected case Initial condition:
% number of points required
 nx = length(x);

% Parameters of regions dimensions
x_middle = ceil(nx/2);
l_1 = 1:x_middle; l_2 = x_middle+1:nx;

% Pre-Allocate variables
r0 = zeros(1,nx);
u0 = zeros(1,nx);
p0 = zeros(1,nx);

% Initial Condition for our 2D domain
% Fugacity
r0(l_1) = rho(1); % region 1
r0(l_2) = rho(2); % region 2
% Velovity in x
u0(l_1) = u(1); % region 1
u0(l_2) = u(2); % region 2
% temperature
p0(l_1) = p(1); % region 1
p0(l_2) = p(2); % region 2

repo: Aero-matlab


text:
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http ://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

function A = buildSecondDerivZMatrix(grid_size)

% d/dz for every entry in the cube except the first and last slice.
m = grid_size(1) * grid_size(2) * (grid_size(3) - 2);
n = grid_size(1) * grid_size(2) * grid_size(3);
e = ones(m, 1);
interior = spdiags([e, -2 * e, e], ...
    [0, grid_size(1) * grid_size(2), 2 * grid_size(1) * grid_size(2)], ...
    m, n);

boundary_z1 = makeBoundaryZ1(n);
boundary_zend = makeBoundaryZEnd(n);

% The matrix for a full cube.
cube = [boundary_z1; interior; boundary_zend];

% Repeat the cube for the last 2 dimensions.
A = sparse(0, 0);
for v = 1:grid_size(5)
    for u = 1:grid_size(4)
        A = blkdiag(A, cube);
    end
end

% Boundary conditions for the first slice.
% n is the number of columns for the full matrix (the number of variables
% in the grid).
function A = makeBoundaryZ1(n)
    mm = grid_size(1) * grid_size(2);
    nn = grid_size(1) * grid_size(2) * 2;
    e = ones(mm, 1);
    B = spdiags([-e, e], [0, grid_size(1) * grid_size(2)], mm, nn);
    % Concat zero columns to the right so we can vertical concat with the full
    % matrix later.
    A = [B, sparse(mm, n - nn)];
end

% Boundary conditions for the last slice.
function A = makeBoundaryZEnd(n)
    mm = grid_size(1) * grid_size(2);
    nn = grid_size(1) * grid_size(2) * 2;
    e = ones(mm, 1);
    B = spdiags([e, -e], [0, grid_size(1) * grid_size(2)], mm, nn);
    % Concat zero columns to the right so we can vertical concat with the full
    % matrix later.
    A = [sparse(mm, n - nn), B];
end

end

author: mahmoudnafifi


text:
function J = imPad( I, pad, type )
% Pad an image along its four boundaries.
%
% Similar to Matlab's padarray, with the following differences:
%  (1) limited to padding along height and width
%  (2) input format allows for separate padding along each dimension
%  (3) padding values may be negative, in which case performs *cropping*
%  (4) optimized (speedup can be significant, esp. for small arrays)
%
% The amount of padding along each of the four boundaries (referred to as
% T/B/L/R) is determined by the parameter "pad" as follows:
%  if(numel(pad)==1): T=B=L=R=pad
%  if(numel(pad)==2): T=B=pad(1), L=R=pad(2)
%  if(numel(pad)==4): T=pad(1), B=pad(2), L=pad(3), R=pad(4)
%
% USAGE
%  J = imPad( I, pad, type )
%
% INPUTS
%  I      - [hxwxk] input image (single, double or uint8 array)
%  pad    - pad or crop amount: 1, 2, or 4 element vector (see above)
%  type   - pad value or 'replicate', 'symmetric', 'circular'
%
% OUTPUTS
%  J      - [T+h+B x L+w+R x k] padded image
%
% EXAMPLE
%  I=imread('peppers.png'); pad=[10 20]; type=50;
%  tic, J1=imPad(I,pad,type); toc
%  tic, J2=padarray(I,pad,type,'both'); toc
%  figure(1); im(J1); isequal(J1,J2)
%
% See also padarray
%
% Piotr's Computer Vision Matlab Toolbox      Version 3.00
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

J = imPadMex( I, pad, type );

%%% OLD Matlab code - slower (although still faster than padarray)
% [h,w,~]=size(I); p=pad; k=length(p);
% if(k==1), p=[p p p p]; elseif(k==2), p=[p(1) p(1) p(2) p(2)]; end
% if( length(type)==1 )
%   J=imPad(I,p,'replicate'); v=feval(class(I),type);
%   J(1:p(1),:,:)=v; J(end-p(2)+1:end,:,:)=v;
%   J(:,1:p(3),:)=v; J(:,end-p(4)+1:end,:)=v;
%   return;
% elseif( strcmp(type,'replicate') )
%   rs = [uint32(ones(1,p(1))) 1:h h*ones(1,p(2))];
%   cs = [uint32(ones(1,p(3))) 1:w w*ones(1,p(4))];
% elseif( strcmp(type,'symmetric') )
%   rs = uint32([1:h h:-1:1]); rs=rs(mod(-p(1):h+p(2)-1,2*h)+1);
%   cs = uint32([1:w w:-1:1]); cs=cs(mod(-p(3):w+p(4)-1,2*w)+1);
% elseif( strcmp(type,'circular') )
%   rs = uint32(1:h); rs=rs(mod(-p(1):h+p(2)-1,h)+1);
%   cs = uint32(1:w); cs=cs(mod(-p(3):w+p(4)-1,w)+1);
% end
% J = I(rs,cs,:);

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
% DEMVOWELSISOMAP Model the vowels data with a 2-D FGPLVM using RBF kernel.

% FGPLVM

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

dataSetName = 'vowels';
experimentNo = 1;
ind
% load data
[Y, lbls] = lvmLoadData(dataSetName);
ind = randperm(size(Y,1));
model.X = isomapEmbed(Y, 2);

save('demVowelsIsomap', 'model');

figure,
hold on

ax = gca;
lvmTwoDPlot(model.X, lbls, getSymbols(size(lbls, 2)));
xLim = [min(model.X(:, 1)) max(model.X(:, 1))]*1.1;
yLim = [min(model.X(:, 2)) max(model.X(:, 2))]*1.1;
set(ax, 'xLim', xLim);
set(ax, 'yLim', yLim);

set(ax, 'fontname', 'arial');
set(ax, 'fontsize', 20);

errors = fgplvmNearestNeighbour(model, lbls);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%% Data import
%Copyright (c) 2011, The MathWorks, Inc.

% Import tumor-weight profile for drugs A & B
ds = dataset('xlsfile', 'Data.xls') ;

% Convert drug type column to nominal
ds.Drug = nominal(ds.Drug)               ;

%% Curve fitting

% Initial estimates for [L0, L1, k1, k2_A and k2_B]
p0 = [0.11 , .149 , 1.836, 0.0154, 0.00732] ;

% Estimate parameters
fh   = @(p , t ) objFcn(p , t , ds.Drug) ;% Function handle
tic
pFit = lsqcurvefit(fh , p0 , ds.Time , ds.TumorWeight ) ;
toc
%% Display results & plot

disp([' Drug-independent parameters: '          , ...
        num2str(pFit(1)), ' (L0), '             , ...
        num2str(pFit(2)), ' (L1) and '          , ...
        num2str(pFit(3)), ' (k1) '              ]) ;

disp([' Drug-dependent parameters: ' , ...
        num2str(pFit(4)), ' (k2 - Drug A) and '              , ...
        num2str(pFit(5)), ' (k2 - Drug B)'   ]) ;

% Plot observations
figure; hold on  ;
idx_A = ds.Drug == 'A' ;
plot(ds.Time(idx_A) , ds.TumorWeight(idx_A) , 'ro')  ;
plot(ds.Time(~idx_A), ds.TumorWeight(~idx_A), 'bo')  ;

% Plot predictions
yPred = fh(pFit, ds.Time) ;
plot(ds.Time(idx_A)  , yPred(idx_A) , 'r:')  ;
plot(ds.Time(~idx_A) , yPred(~idx_A), 'b:')  ;

xlabel('Time (days)')
ylabel('Tumor Weight (milligram)')
title('Tumor growth profile')
legend({'Drug A', 'Drug B'})

author: Sable

repo: mcbench-benchmarks


text:
function f = not(f)
%~   CHEBFUN logical NOT.
%   NOT(F) returns a CHEBFUN which evaluates to zero at all points where F is
%   nonzero and one otherwise.

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

% Add breaks at the roots (since these will take the value 1 in the output).
f = addBreaksAtRoots(f);

% Tolerance:
vs = vscale(f);
vs(vs < eps) = 1;
tol = eps*vs;

% Loop over the FUNs:
for k = 1:numel(f.funs)
    f.funs{k} = not(f.funs{k});
end

% pointValues:
f.pointValues = abs(f.pointValues) < tol;

end


text:
function [y,a]=plane_project2(x,E)
    [N,M]=size(x);
    p=size(E,2);
    a=zeros(p,M);
    ct=E(:,1);
    Ep=E(:,2:p)-ct*ones(1,p-1);
    a(2:p,:)=Ep\(x-ct*ones(1,M));
    a(1,:)=ones(1,M)-sum(a(2:p,:),1);
    y=E*a;
end

author: zhouyuanzxcv

repo: Hyperspectral


text:
	function f=def_funct(t,t1,a1,t2,a2,a3)
			% **************************************************
			% Functie MATLAB destinata generarii
			% unei functii de forma
			%
			%        a1, pentru t<t1
			% f(t)=  a2, pentru t>t2
			%        a3, pentru t>=t1 si t<=t2
			%
			% Variabile de intrare:
			% - t intervalul total de definitie al functiei
			% - t1, t2 limitele de definitie
			% - a1, a2, a3 valorile functiei
			%   in diferitele intervale
			% Variabila de iesire: vectorul f
			% **************************************************

			% Determinarea indicilor pentru care t<t1
	k1=find(t<t1);
			% Determinarea indicilor pentru care t>t2
	k3=find(t>t2);
			% Generarea tuturor indicilor
	k2=1:length(t);
			% Generarea indicilor din mijloc
			% prin extragerea celor pentru care t<t1 si t>t2
	k2([k1 k3])=[];
			% Incarcarea vectorului rezultat
			% (utilizand tricul lui Tony):
			% - prima parte cu valorile a1
	f(k1)=a1(1,ones(1,length(k1)));
			% - partea din mijloc cu valorile a3
	f(k2)=a3(1,ones(1,length(k2)));
			% - ultima parte cu valorile a2
	f(k3)=a2(1,ones(1,length(k3)));


text:
% download.m

close all;
clear;
clc;

%% dataset
fprintf('###### Dataset downlowd ######\n\n');
% local folder
local_folder_name = 'dataset';

if ~exist('dataset', 'dir')
    mkdir('dataset');
end

% remote url, path, and filelist
site_url = 'http://www.kasailab.com/';
site_path = 'public/github/ClassifierToolbox/dataset/';
filename_array = {'ORL_Face_img.mat', 'ORL_Face_img_cov.mat', 'AR_Face_img_27x20.mat', ...
                    'AR_Face_img_60x43.mat','Brodatz_texture_img_small_set.mat', ...
                    'USPS.mat', 'MNIST.mat', ...
                    'test_cov.mat'};
dataset_num = length(filename_array);

% download files
for i = 1 : dataset_num
    file_name = filename_array{i};
    filename_full = sprintf('%s%s%s', site_url, site_path, file_name);
    filename_full_local = sprintf('%s/%s', local_folder_name, file_name);

    fprintf('# Dataset downlaoding [%d/%d]\n', i, dataset_num);

    if ~exist(filename_full_local, 'file')
        fprintf('  * file:\t "%s"\n', file_name);
        fprintf('  * from:\t "%s%s"\n', site_url, site_path);
        fprintf('  * into:\t "%s"\n', filename_full_local);
        fprintf('  * ........ ');
        websave(filename_full_local, filename_full);
        fprintf('done\n\n');
    else
        fprintf('  * ........ %s already exists. Skip downloading.\n\n', file_name);
    end
end

fprintf('# %d dataset files have been successfully downloaded.\n\n\n', dataset_num);

%% libraries
fprintf('###### Library downlowd ######\n\n');

if ~exist('lib', 'dir')
    mkdir('lib');
end

cd lib/;

lib_num = 4;
lib_cell = cell(1,lib_num);

lib_struct.name = 'NMFLibrary-master';
lib_struct.url = 'https://github.com/hiroyuki-kasai/NMFLibrary/archive/master.zip';
lib_cell{1} = lib_struct;

lib_struct.name = 'SparseGDLibrary-master';
lib_struct.url = 'https://github.com/hiroyuki-kasai/SparseGDLibrary/archive/master.zip';
lib_cell{2} = lib_struct;

lib_struct.name = 'DSK-master';
lib_struct.url = 'https://github.com/seuzjj/DSK/archive/master.zip';
lib_cell{3} = lib_struct;

lib_struct.name = 'KMeans_SPD_Matrices';
lib_struct.url = 'https://jp.mathworks.com/matlabcentral/mlc-downloads/downloads/submissions/46343/versions/2/download/zip';
lib_cell{4} = lib_struct;

% download files
for i = 1 : lib_num
    file_name = lib_cell{i}.name;
    file_url =  lib_cell{i}.url;

    fprintf('# Library downlaoding [%d/%d]\n', i, lib_num);

    if ~exist(file_name, 'dir')
        fprintf('  * file:\t "%s"\n', file_name);
        fprintf('  * from:\t "%s"\n', file_url);
        fprintf('  * into:\t "%s"\n', 'download.zip');
        fprintf('  * ........ ');
        websave('download.zip', file_url);
        fprintf('done.\n');

        fprintf('# unzip download.zip ... ');
        unzip download.zip;
        fprintf('done.\n');

        fprintf('# delete download.zip ... ');
        delete('download.zip');
        fprintf('done.\n\n');
    else
        fprintf('  * ........ %s already exists. Skip downloading.\n\n', file_name);
    end
end

cd ..;

fprintf('# %d library files have been successfully downloaded and installed.\n', lib_num);

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function XYdata = add2XYdata(XYdata,newx,newy)
%This function takes points that have been simulated and adds them to the
%existing data structure
%XYdata is a structure with n elements of .x and .y. n is the number of
% design points tested. .x is a 1xd where d is the number of dimensions. .y
% is 1xr where r is the number of replicates for the design point.

n= size(newx,1);

%If the structure does not exist, create it
if (isempty(XYdata))
    XYdata.x = [];
    XYdata.y = [];
end

%Convert structure to a matrix of points
x = reshape([XYdata.x],length(XYdata(1).x),[])';

%Loop through all of the points to be added
for i=1:n
    %Does the points exist already in the array?
    %TODO this can be searched all at once.
    I = find(ismember(x,newx(i,:),'rows')==1);
    %If the point exists, add the response
    if (I>0)
        XYdata(I).y=[XYdata(I).y,newy(i)];
    %else if the point does not exist, add it
    else
        %Add the points to a temp structure
        temp.x = newx(i,:);
        temp.y = newy(i);
        %Add the point in case there is a future point with the same
        % location.
        x=[x;newx(i,:)];

        %If the first array is empty, add point to the first location
        if (isempty(XYdata(1).x))
            XYdata(1) = temp;
        else
            %else, build the array
            XYdata = [XYdata temp];
        end
    end
end


text:
%% Runme!
% Coded by Manuel Diaz, f99543083
% NTU, IAM, 2012.04.01.

% Exercise 1.a
HWproblem1adIso

% Exercise 1.b
HWproblem1bdIso

% Exercise 2
HWproblem2dIso

% Exercise 3
HWproblem3dIso

% Exercise 4
HWproblem4dIso


text:
% Fast convolution, replacement for both conv2 and convn.
%
% See conv2 or convn for more information on convolution in general.
%
% This works as a replacement for both conv2 and convn.  Basically,
% performs convolution in either the frequency or spatial domain, depending
% on which it thinks will be faster (see below). In general, if A is much
% bigger then B then spatial convolution will be faster, but if B is of
% similar size to A and both are fairly big (such as in the case of
% correlation), convolution as multiplication in the frequency domain will
% tend to be faster.
%
% The shape flag can take on 1 additional value which is 'smooth'.  This
% flag is intended for use with smoothing kernels.  The returned matrix C
% is the same size as A with boundary effects handled in a special manner.
% That is instead of A being zero padded before being convolved with B;
% near the boundaries a cropped version of the matrix B is used, and the
% results is scaled by the fraction of the weight found in  the cropped
% version of B.  In this case each dimension of B must be odd, and all
% elements of B must be positive.  There are other restrictions on when
% this flag can be used, and in general it is only useful for smoothing
% kernels.  For 2D filtering it does not have much overhead, for 3D it has
% more and for higher dimensions much much more.
%
% For optimal performance some timing constants must be set to choose
% between doing convolution in the spatial and frequency domains, for more
% info see timeConv below.
%
% USAGE
%  C = convnFast( A, B, [shape] )
%
% INPUTS
%  A       - d dimensional input matrix
%  B       - d dimensional matrix to convolve with A
%  shape   - ['full'] 'valid', 'full', 'same', or 'smooth'
%
% OUTPUTS
%  C       - result of convolution
%
% EXAMPLE
%
% See also CONV2, CONVN

% Piotr's Image&Video Toolbox      Version 2.0
% Written and maintained by Piotr Dollar    pdollar-at-cs.ucsd.edu
% Please email me if you find bugs, or have suggestions or questions!

function C = convnFast( A, B, shape )

if( nargin<3 || isempty(shape)); shape='full'; end
if(isempty(strmatch(shape, char({'same', 'valid', 'full', 'smooth'}))))
  error( 'convnFast: unknown shape flag' ); end

shapeorig = shape;
smoothFlag = (strcmp(shape,'smooth'));
if( smoothFlag ); shape = 'same'; end;

% get dimensions of A and B
ndA = ndims(A);  ndB = ndims(B); nd = max(ndA,ndB);
sizA = size(A); sizB = size(B);
if (ndA>ndB); sizB = [sizB ones(1,ndA-ndB)]; end
if (ndA<ndB); sizA = [sizA ones(1,ndB-ndA)]; end

% ERROR CHECK if smoothflag
if( smoothFlag )
  if( ~all( mod(sizB,2)==1 ) )
    error('If flag==''smooth'' then must have odd sized mask');
  end;
  if( ~all( B>0 ) )
    error('If flag==''smooth'' then mask must have >0 values.');
  end;
  if( any( (sizB-1)/2>sizA ) )
    error('B is more then twice as big as A, cannot use flag==''smooth''');
  end;
end

% OPTIMIZATION for 3D conv when B is actually 2D - calls (spatial) conv2
% repeatedly on 2D slices of A.  Note that may need to rearange A and B
% first and use recursion. The benefits carry over to convnBound
% (which is faster for 2D arrays).
if( ndA==3 && ndB==3 && (sizB(1)==1 || sizB(2)==1) )
  if (sizB(1)==1)
    A = permute( A, [2 3 1]);  B = permute( B, [2 3 1]);
    C = convnFast( A, B, shapeorig );
    C = permute( C, [3 1 2] );
  elseif (sizB(2)==1)
    A = permute( A, [3 1 2]);  B = permute( B, [3 1 2]);
    C = convnFast( A, B, shapeorig );
    C = permute( C, [2 3 1] );
  end
  return;
elseif( ndA==3 && ndB==2 )
  C1 = conv2( A(:,:,1), B, shape );
  C = zeros( [size(C1), sizA(3)] ); C(:,:,1) = C1;
  for i=2:sizA(3); C(:,:,i) = conv2( A(:,:,i), B, shape ); end
  if (smoothFlag)
    for i=1:sizA(3)
      C(:,:,i) = convnBound(A(:,:,i),B,C(:,:,i),sizA(1:2),sizB(1:2));
    end
  end
  return;
end

% get predicted time of convolution in frequency and spatial domain
% constants taken from timeConv
sizfft = 2.^ceil(real(log2(sizA+sizB-1))); psizfft=prod(sizfft);
frequenPt = 3 * 1e-7 * psizfft * log(psizfft);
if (nd==2)
  spatialPt = 5e-9 * sizA(1) * sizA(2) * sizB(1) * sizB(2);
else
  spatialPt = 5e-8 * prod(sizA) * prod(sizB);
end

% perform convolution
if ( spatialPt < frequenPt )
  if (nd==2)
    C = conv2( A, B, shape );
  else
    C = convn( A, B, shape );
  end
else
  C = convnFreq( A, B, sizA, sizB, shape );
end;

% now correct boundary effects (if shape=='smooth')
if( ~smoothFlag ); return; end;
C = convnBound( A, B, C, sizA, sizB );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate boundary values for C in spatial domain
function C = convnBound( A, B, C, sizA, sizB )
nd = length(sizA);
radii = (sizB-1)/2;

% flip B appropriately (conv flips B)
for d=1:nd; B = flipdim(B,d); end

% get location that need to be updated
% this is the LEAST efficient part (fixing is annoying though)
inds = {':'}; inds = inds(:,ones(1,nd));
Dind = zeros( sizA );
for d=1:nd
  inds1 = inds; inds1{ d } = 1:radii(d);
  inds2 = inds; inds2{ d } = sizA(d)-radii(d)+1:sizA(d);
  Dind(inds1{:}) = 1;  Dind(inds2{:}) = 1;
end
Dind = find( Dind );
Dndx = ind2sub2( sizA, Dind );
nlocs = length(Dind);

% get cuboid dimensions for all the boundary regions
sizeArep = repmat( sizA, [nlocs,1] );
radiiRep = repmat( radii, [nlocs,1] );
Astarts = max(1,Dndx-radiiRep);
Aends = min( sizeArep, Dndx+radiiRep);
Bstarts = Astarts + (1-Dndx+radiiRep);
Bends = Bstarts + (Aends-Astarts);

% now update these locations
vs = zeros( 1, nlocs );
if( nd==2 )
  for i=1:nlocs % accelerated for 2D arrays [fast]
    Apart = A( Astarts(i,1):Aends(i,1), Astarts(i,2):Aends(i,2) );
    Bpart = B( Bstarts(i,1):Bends(i,1), Bstarts(i,2):Bends(i,2) );
    v = (Apart.*Bpart); vs(i) = sum(v(:)) ./ sum(Bpart(:));
  end
elseif( nd==3 ) % accelerated for 3D arrays
  for i=1:nlocs
    Apart = A( Astarts(i,1):Aends(i,1), Astarts(i,2):Aends(i,2), ...
      Astarts(i,3):Aends(i,3) );
    Bpart = B( Bstarts(i,1):Bends(i,1), Bstarts(i,2):Bends(i,2), ...
      Bstarts(i,3):Bends(i,3) );
    za = sum(sum(sum(Apart.*Bpart))); zb=sum(sum(sum(Bpart)));
    vs(1,i) = za./zb;
  end
else % general case [slow]
  extract=cell(1,nd);
  for i=1:nlocs
    for d=1:nd; extract{d} = Astarts(i,d):Aends(i,d); end
    Apart = A( extract{:} );
    for d=1:nd; extract{d} = Bstarts(i,d):Bends(i,d); end
    Bpart = B( extract{:} );
    v = (Apart.*Bpart); vs(i) = sum(v(:)) ./ sum(Bpart(:));
  end
end
C( Dind ) = vs * sum(B(:));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Convolution as multiplication in the frequency domain
function C = convnFreq( A, B, sizA, sizB, shape )
siz = sizA + sizB - 1;

% calculate correlation in frequency domain
Fa = fftn(A,siz);
Fb = fftn(B,siz);
C = ifftn(Fa .* Fb);

% make sure output is real if inputs were both real
if(isreal(A) && isreal(B)); C = real(C); end

% crop to size
if(strcmp(shape,'valid'))
  C = arrayToDims( C, max(0,sizA-sizB+1 ) );
elseif(strcmp(shape,'same'))
  C = arrayToDims( C, sizA );
elseif(~strcmp(shape,'full'))
  error('unknown shape');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function used to calculate constants for prediction of convolution in the
% frequency and spatial domains.  Method taken from normxcorr2.m
% May need to reset K's if placing this on a new machine, however, their
% ratio should be about the same..
function K = timeConv() %#ok<DEFNU>
mintime = 4;

switch 3
  case 1  % conv2  [[empirically K = 5e-9]]
    % convolution time = K*prod(size(a))*prod(size(b))
    siza = 30;  sizb = 200;
    a = ones(siza);  b = ones(sizb);
    t1 = cputime;  t2 = t1; k = 0;
    while (t2-t1)<mintime;
      disc = conv2(a,b); k = k + 1; t2 = cputime; %#ok<NASGU>
    end
    K = (t2-t1)/k/siza^2/sizb^2;

  case 2  % convn  [[empirically K = 5e-8]]
    % convolution time = K*prod(size(a))*prod(size(b))
    siza = [10 10 10];  sizb = [30 30 10];
    a = ones(siza);  b = ones(sizb);
    t1 = cputime;  t2 = t1;  k = 0;
    while (t2-t1)<mintime;
      disc = convn(a,b); k = k + 1; t2 = cputime; %#ok<NASGU>
    end
    K = (t2-t1)/k/prod(siza)/prod(sizb);

  case 3 % fft (one dimensional) [[empirically K = 1e-7]]
    % fft time = K * n log(n)  [if n is power of 2]
    % Works fastest for powers of 2.  (so always zero pad until have
    % size of power of 2?).  2 dimensional fft has to apply single
    % dimensional fft to each column, and then signle dimensional fft
    % to each resulting row.  time = K * (mn)log(mn).  Likewise for
    % highter dimensions.  convnFreq requires 3 such ffts.
    n = 2^nextpow2(2^15);
    vec = complex(rand(n,1),rand(n,1));
    t1 = cputime;  t2 = t1;  k = 0;
    while (t2-t1) < mintime;
      disc = fft(vec); k = k + 1; t2 = cputime; %#ok<NASGU>
    end
    K = (t2-t1) / k / n / log(n);
end

author: kristinbranson

repo: JAABA


text:
function [result] = ft_test_run(varargin)

% FT_TEST_RUN

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

narginchk(1, inf);
command = varargin{1};
assert(isequal(command, 'run') || isequal(command, 'inventorize'));
varargin = varargin(2:end);

optbeg = find(ismember(varargin, {'dependency', 'dccnpath', 'maxmem', 'maxwalltime', 'upload', 'sort', 'returnerror'}));
if ~isempty(optbeg)
  optarg   = varargin(optbeg:end);
  varargin = varargin(1:optbeg-1);
else
  optarg = {};
end

% varargin contains the file (or files) to test
% optarg contains the command-specific options

% get the optional input arguments
dependency  = ft_getopt(optarg, 'dependency', {});
hasdccnpath = ft_getopt(optarg, 'dccnpath');  % default is handled below
maxmem      = ft_getopt(optarg, 'maxmem', inf);
maxwalltime = ft_getopt(optarg, 'maxwalltime', inf);
upload      = ft_getopt(optarg, 'upload', 'yes'); % this will be set to 'no' in case FieldTrip version is not clean
sortarg     = ft_getopt(optarg, 'sort', 'alphabetical');
returnerror = ft_getopt(optarg, 'returnerror', 'no');

if ischar(dependency)
  % this should be a cell-array
  dependency = {dependency};
end

if isempty(hasdccnpath)
  % true when central storage is available, false otherwise
  hasdccnpath = exist(dccnpath('/home/common/matlab/fieldtrip/data/'), 'dir');
end

if ischar(maxwalltime)
  % it is probably formatted as HH:MM:SS, convert to seconds
  maxwalltime = str2walltime(maxwalltime);
end

if ischar(maxmem)
  % it is probably formatted as XXmb, or XXgb, convert to bytes
  maxmem = str2mem(maxmem);
end

% get the version and the path
[revision, ftpath] = ft_version;

%% determine the list of functions to test
if ~isempty(varargin) && exist(varargin{1}, 'file')
  functionlist = varargin;
else
  addpath(fullfile(ftpath, 'test'));
  d = dir(fullfile(ftpath, 'test', 'test_*.m'));
  functionlist = {d.name}';
  for i=1:numel(functionlist)
    functionlist{i} = functionlist{i}(1:end-2); % remove the extension
  end
end

%% determine the list of files to test
filelist = cell(size(functionlist));
for i=1:numel(functionlist)
  filelist{i} = which(functionlist{i});
end

fprintf('considering %d test scripts\n', numel(filelist));

%% make a subselection based on the filters
dep  = true(size(filelist));
file = false(size(filelist)); % by default ignore file reads
mem  = zeros(size(filelist));
tim  = zeros(size(filelist));

for i=1:numel(filelist)
  fid = fopen(filelist{i}, 'rt');
  str = fread(fid, [1 inf], 'char=>char');
  fclose(fid);
  line = tokenize(str, 10);

  if ~isempty(dependency)
    dep(i) = false;
  else
    dep(i) = true;
  end

  for k=1:numel(line)
    for j=1:numel(dependency)
      % search for the dependencies in each of the test functions
      [s, e] = regexp(line{k}, sprintf('%% DEPENDENCY.*%s.*', dependency{j}), 'once', 'start', 'end');
      if ~isempty(s)
        dep(i) = true;
      end
    end

    if ~istrue(hasdccnpath)
      % search for the occurence of the DCCNPATH function in each of the test functions
      [s, e] = regexp(line{k}, 'dccnpath', 'once', 'start', 'end');
      if ~isempty(s)
        file(i) = true;
      end
    end

    [s, e] = regexp(line{k}, '% WALLTIME.*', 'once', 'start', 'end');
    if ~isempty(s)
      s = s + length('% WALLTIME'); % strip this part
      tim(i) = str2walltime(line{k}(s:e));
    end

    [s, e] = regexp(line{k}, '% MEM.*', 'once', 'start', 'end');
    if ~isempty(s)
      s = s + length('% MEM'); % strip this part
      mem(i) = str2mem(line{k}(s:e));
    end
  end % for each line

end % for each function/file

fprintf('%3d scripts are excluded due to the dependencies\n',                  sum(~dep));
fprintf('%3d scripts are excluded due to loading files from the DCCN path\n',  sum(file));
fprintf('%3d scripts are excluded due to the requirements for memory\n',       sum(mem>maxmem));
fprintf('%3d scripts are excluded due to the requirements for walltime \n',    sum(tim>maxwalltime));

% make selection of test scripts, remove scripts that exceed walltime or memory
sel                  = true(size(filelist));
sel(tim>maxwalltime) = false;
sel(mem>maxmem)      = false;
sel(dep==false)      = false;
sel(file==true)      = false;

% make the subselection of functions to test
functionlist = functionlist(sel);
tim = tim(sel);
mem = mem(sel);

switch (sortarg)
  case {'alphabet' 'alphabetic' 'alphabetical'}
    [functionlist, indx] = sort(functionlist);
    tim = tim(indx);
    mem = mem(indx);
  case {'mem' 'memory'}
    [mem, indx]  = sort(mem);
    tim          = tim(indx);
    functionlist = functionlist(indx);
  case {'walltime' 'time'}
    [tim, indx]  = sort(tim);
    mem          = mem(indx);
    functionlist = functionlist(indx);
  case 'random'
    indx = randperm(numel(functionlist));
    functionlist = functionlist(indx);
    tim          = tim(indx);
    mem          = mem(indx);
  otherwise
    ft_error('incorrect specification for ''sort''');
end

% this will hold all results
result = [];

if strcmp(command, 'inventorize')
  % do not execute anything, only return the function names
  for i=1:numel(functionlist)
    result(i).functionname = functionlist{i};
    result(i).mem          = mem(i);
    result(i).walltime     = tim(i);
  end
  fprintf('selecting %d test scripts\n', numel(functionlist));
  return
else
  fprintf('executing %d test scripts\n', numel(functionlist));
end

% uploading results to the dashboard from a work-in-progress version is not supported
if istrue(upload) && ~istrue(ft_version('clean'))
  warning('FieldTrip version is not clean, not uploading results to the dashboard')
  upload = 'no';
end

%% run over all tests
for i=1:numel(functionlist)

  close all
  fprintf('================================================================================\n');;
  fprintf('=== evaluating %s\n', functionlist{i});
  try
    stopwatch = tic;
    eval(functionlist{i});
    passed = true;
    runtime = round(toc(stopwatch));
    fprintf('=== %s PASSED in %d seconds\n', functionlist{i}, runtime);
  catch me
    passed = false;
    runtime = round(toc(stopwatch));
    % show the error with the stack trace
    fprintf('Error using %s (line %d)\n', me.stack(1).name, me.stack(1).line);
    disp(me.message)
    for j=2:numel(me.stack)
      fprintf('Error in %s (line %d)\n', me.stack(j).name, me.stack(j).line);
    end
    fprintf('=== %s FAILED in %d seconds\n', functionlist{i}, runtime);
  end
  close all

  result(i).matlabversion    = version('-release');
  result(i).fieldtripversion = revision;
  result(i).branch           = ft_version('branch');
  result(i).arch             = computer('arch');
  result(i).hostname         = gethostname;
  result(i).user             = getusername;
  result(i).passed           = passed;
  result(i).runtime          = runtime;
  result(i).functionname     = functionlist{i};
  result(i).mem              = mem(i);
  result(i).walltime         = tim(i);

  if istrue(upload)
    try
      % the weboptions function is available in 2014b onward, but behaves inconsistently
      options = weboptions('MediaType','application/json');
    catch
      options = [];
    end
    url = 'https://dashboard.fieldtriptoolbox.org/api/';
    webwrite(url, result(i), options);
  end

  if strcmp(returnerror, 'immediate') && ~passed
    % report only on the tests completed sofar and give an error
    printstruct_as_table(result);
    rethrow(me);
  end

end

if strcmp(returnerror, 'final') && any(~[result.passed])
  % report on the complete batch and give an error
  printstruct_as_table(result);
  failed = find(~[result.passed]);
  error('%d of the test scripts failed', numel(failed));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function walltime = str2walltime(str)
str = lower(strtrim(str));
walltime = str2double(str);
if isnan(walltime)
  str = strtrim(str);
  hms = sscanf(str, '%d:%d:%d'); % hours, minutes, seconds
  walltime = 60*60*hms(1) + 60*hms(2) + hms(3);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mem = str2mem(str)
str = lower(strtrim(str));
mem = str2double(str);
if isnan(mem)
  if ~isempty(regexp(str, '[0-9]*kb', 'once'))
    mem = str2double(str(1:end-2)) * 2^10;
  elseif ~isempty(regexp(str, '[0-9]*mb', 'once'))
    mem = str2double(str(1:end-2)) * 2^20;
  elseif ~isempty(regexp(str, '[0-9]*gb', 'once'))
    mem = str2double(str(1:end-2)) * 2^30;
  elseif ~isempty(regexp(str, '[0-9]*tb', 'once'))
    mem = str2double(str(1:end-2)) * 2^40;
  else
    mem = str2double(str);
  end
end

author: fieldtrip

repo: fieldtrip


% global parameters (modified based on the original MDP code)
function opt = globals()

opt.root = pwd;

% path for MOT benchmark
mot_paths = {'data/', ''};
for i = 1:numel(mot_paths)
    if exist(mot_paths{i}, 'dir')
        opt.mot = mot_paths{i};
        break;
    end
end

opt.mot2d = 'MOT16';
opt.results_dir = 'results/';
opt.mot2d_train_seqs = {'MOT16-02', 'MOT16-04', 'MOT16-05', 'MOT16-09', 'MOT16-10', 'MOT16-11', 'MOT16-13'};
opt.mot2d_train_nums = [600, 1050, 837, 525, 654, 900, 750];
opt.mot2d_test_seqs = {'MOT16-01', 'MOT16-03', 'MOT16-06', 'MOT16-07', 'MOT16-08', 'MOT16-12', 'MOT16-14'};
opt.mot2d_test_nums = [450, 1500, 1194, 500, 625, 900, 750];

addpath(fullfile(opt.root, 'devkit', 'utils'));
addpath([opt.root '/3rd_party/libsvm-3.20/matlab']);
addpath([opt.root '/3rd_party/Hungarian']);

% tracking parameters
opt.aspect_ratio_thre = 0.6;    % aspect ratio threshold in target association
opt.distance_thre = 2;        % distance threshold in target association, multiple of the width of target
opt.det_overlap_thre = 0.5;        % overlap with detection
opt.tracking_score_thre = 0.18;
opt.termination_thre = 2;
opt.association_score_thre = 0.6;
opt.exit_thre = 0.3;
opt.initialization_thre = 0.25;
opt.overlap_sup = 0.7;

% scene context
opt.backgrounds = cell(1,7);
opt.backgrounds{1} = [240, 570, 540, 760; 1, 1, 1920, 400]; % [x1, y1, x2, y2]
opt.backgrounds{2} = [1, 1, 155, 460; 1610, 670, 1710, 770; 700, 445, 760, 520; 1270, 210, 1360, 300; 1760, 900, 1830, 1040; 1455, 770, 1515, 890];
opt.backgrounds{7} = [1, 1, 1920, 450];
opt.h_max = cell(1,7);
opt.h_max{2} = 315;
opt.h_min = cell(1,7);
opt.fps = [30, 30, 14, 30, 30, 30, 25];

% tracking state: START -> ACTIVATED -> TRACKED -> LOST -> STOP
opt.STATE_START = 4;
opt.STATE_ACTIVATED = 2;
opt.STATE_TRACKED = 1;
opt.STATE_LOST = 3;
opt.STATE_STOP = 5;


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

%*****************************************************************************80
%
%% CHI_SQUARE_CDF_VALUES returns some values of the Chi-Square CDF.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      Needs["Statistics`ContinuousDistributions`"]
%      dist = ChiSquareDistribution [ df ]
%      CDF [ dist, x ]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer A, the parameter of the function.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 21;

  a_vec = [ ...
     1,  2,  1,  2, ...
     1,  2,  3,  4, ...
     1,  2,  3,  4, ...
     5,  3,  3,  3, ...
     3,  3, 10, 10, ...
    10 ];

  fx_vec = [ ...
     0.7965567455405796E-01, ...
     0.4987520807317687E-02, ...
     0.1124629160182849E+00, ...
     0.9950166250831946E-02, ...
     0.4729107431344619E+00, ...
     0.1812692469220181E+00, ...
     0.5975750516063926E-01, ...
     0.1752309630642177E-01, ...
     0.6826894921370859E+00, ...
     0.3934693402873666E+00, ...
     0.1987480430987992E+00, ...
     0.9020401043104986E-01, ...
     0.3743422675270363E-01, ...
     0.4275932955291202E+00, ...
     0.6083748237289110E+00, ...
     0.7385358700508894E+00, ...
     0.8282028557032669E+00, ...
     0.8883897749052874E+00, ...
     0.1721156299558408E-03, ...
     0.3659846827343712E-02, ...
     0.1857593622214067E-01 ];

  x_vec = [ ...
     0.01E+00, ...
     0.01E+00, ...
     0.02E+00, ...
     0.02E+00, ...
     0.40E+00, ...
     0.40E+00, ...
     0.40E+00, ...
     0.40E+00, ...
     1.00E+00, ...
     1.00E+00, ...
     1.00E+00, ...
     1.00E+00, ...
     1.00E+00, ...
     2.00E+00, ...
     3.00E+00, ...
     4.00E+00, ...
     5.00E+00, ...
     6.00E+00, ...
     1.00E+00, ...
     2.00E+00, ...
     3.00E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

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

  return
end


text:
function x = VBA_spm_invFcdf(F,v,w)
% Inverse Cumulative Distribution (CDF) of F (Fisher-Snedecor) distribution
% FORMAT x = spm_invFpdf(F,df)
% FORMAT x = spm_invFpdf(F,v,w)
%
% F  - CDF (lower tail p-value)
% df - Degrees of freedom, concatenated along last dimension
%      Eg. Scalar (or column vector) v & w. Then df=[v,w];
% v  - Shape parameter 1 /   numerator degrees of freedom (v>0)
% w  - Shape parameter 2 / denominator degrees of freedom (w>0)
% x  - F-variate   (F has range [0,Inf) )
%__________________________________________________________________________
%
% spm_Fcdf implements the inverse Cumulative Distribution Function
% for the F-distribution.
%
% Definition:
%--------------------------------------------------------------------------
% The CDF F(x) of the F distribution with degrees of freedom v & w,
% defined for positive integer degrees of freedom v & w, is the
% probability that a realisation of an F random variable X has value
% less than x F(x)=Pr{X<x} for X~F(v,w). The F-distribution is defined
% for v>0 & w>0, and for x in [0,Inf) (See Evans et al., Ch16).
%
% Variate relationships: (Evans et al., Ch16 & 37)
%--------------------------------------------------------------------------
% The square of a Student's t variate with w degrees of freedom is
% distributed as an F-distribution with [1,w] degrees of freedom.
%
% For X an F-variate with v,w degrees of freedom, w/(w+v*X^2) has
% distribution related to a Beta random variable with shape parameters
% w/2 & v/2, as described below.
%
% Algorithm:
%--------------------------------------------------------------------------
% Using the routine spm_invBcdf for the Beta distribution, with
% appropriate parameters:  The CDF of the F-distribution with v,w
% degrees of freedom is related to the incomplete beta function by:
%       Pr(X<x) = 1 - betainc(w/(w+v*x^2),w/2,v/2)
% See Abramowitz & Stegun, 26.6.2; Press et al., Sec6.4 for
% definitions of the incomplete beta function. The relationship is
% easily verified by substituting for w/(w+v*x^2) in the integral of the
% incomplete beta function.
%
%
% References:
%--------------------------------------------------------------------------
% Evans M, Hastings N, Peacock B (1993)
%       "Statistical Distributions"
%        2nd Ed. Wiley, New York
%
% Abramowitz M, Stegun IA, (1964)
%       "Handbook of Mathematical Functions"
%        US Government Printing Office
%
% Press WH, Teukolsky SA, Vetterling AT, Flannery BP (1992)
%       "Numerical Recipes in C"
%        Cambridge
%
%__________________________________________________________________________
% Copyright (C) 1993-2011 Wellcome Trust Centre for Neuroimaging

% Andrew Holmes
% $Id: spm_invFcdf.m 4182 2011-02-01 12:29:09Z guillaume $

%-Format arguments, note & check sizes
%--------------------------------------------------------------------------
if nargin<2, error('Insufficient arguments'), end

%-Unpack degrees of freedom v & w from single df parameter (v)
if nargin<3
    vs = size(v);
    if prod(vs)==2
        %-DF is a 2-vector
        w = v(2); v = v(1);
    elseif vs(end)==2
        %-DF has last dimension 2 - unpack v & w
        nv = prod(vs);
        w  = reshape(v(nv/2+1:nv),vs(1:end-1));
        v  = reshape(v(1:nv/2)   ,vs(1:end-1));
    else
        error('Can''t unpack both df components from single argument')
    end
end

%-Check argument sizes
ad = [ndims(F);ndims(v);ndims(w)];
rd = max(ad);
as = [[size(F),ones(1,rd-ad(1))];...
      [size(v),ones(1,rd-ad(2))];...
      [size(w),ones(1,rd-ad(3))]];
rs = max(as);
xa = prod(as,2)>1;
if sum(xa)>1 && any(any(diff(as(xa,:)),1))
    error('non-scalar args must match in size');
end

%-Computation
%--------------------------------------------------------------------------
%-Initialise result to zeros
x = zeros(rs);

%-Only defined for F in [0,1] & strictly positive v & w.
% Return NaN if undefined.
md = ( F>=0  &  F<=1  &  v>0  &  w>0 );
if any(~md(:))
    x(~md) = NaN;
    warning('Returning NaN for out of range arguments');
end

%-Special cases: x=0 when F=0, x=Inf when F=1
x(md & F==1) = Inf;

%-Compute where defined & not special case
Q  = find( md  &  F>0  &  F<1 );
if isempty(Q), return, end
if xa(1), QF=Q; else QF=1; end
if xa(2), Qv=Q; else Qv=1; end
if xa(3), Qw=Q; else Qw=1; end

%-Compute
bQ   = VBA_spm_invBcdf(1-F(QF),w(Qw)/2,v(Qv)/2);
x(Q) = (w(Qw)./bQ -w(Qw))./v(Qv);

author: MBB-team

repo: VBA-toolbox


text:
function runBinEnergy(energy, fieldSize, numBin, extraBin, startBin, endBin, IncreaseFluence, OnlyHorn, UseFlatFilter, FFA, FFB, FFDistance, ParticleType, DoseNormFlag, nhist, r, c, s)
% 'runBinEnergy'
% JC 1/29/06
% How to use
%> runBinEnergy 18 20 4 100000000
%
% Written by JC. Dec 2005
% Run DPM for the different energy bins, for optimize the photon spectrum later.
%"input" is the .MAT file have DPMIN, CTscan, and indV
%DPMIN - DPM input structure
%CTscan - DPM scan structure
%indV = sub2ind(getUniformizedSize(planC), r, c, s); see DPMInfluence.m
%line 72
%DPMIN_CTscan_indV DMPIN, CTScan, indV
%"number" is how many bins we're going to divide the whole energy range
% Typical value for "number" is 20
%nhis - number of histories
%
% JC Oct 18, 2006  use target to flattening filter distance as 12.5cm instead of 9.5cm);
% JC Feb 26, 2007
    % Add  DoseNormFlag
    % DPMIN.NumParticles = nhist*fieldSize*fieldSize;
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

setappdata(0, 'usenativesystemdialogs', false)
currentDir = cd;

[FileName,path] = uigetfile('*.mat','Select MAT file containing DPMIN, CTscan, & indV, p');

if path == 0
    errordlg('File Should exist');
    error('File Should exist');
end

cd(path);
load (FileName, 'DPMIN', 'CTscan', 'indV')
cd(currentDir);

if (class(energy)== 'char')
    energy = str2num(energy);numBin = str2num(numBin);extraBin = str2num(extraBin);nhist = str2num(nhist);
    r = str2num(r); c = str2num(c); s = str2num(s); fieldSize = str2num(fieldSize);
end

[FileName,path] = uigetfile('*.*','Select dpm.dll to execute');

if path == 0
    errordlg('File Should exist');
    error('File Should exist');
end

% need ./local to run dpm. if no. crush.
% so creat it
mkdir local

rand('state',sum(100*clock));
%for i = 0: numBin + extraBin - 1,

% Add input ParticleType as 0 for photon or -1 for electron
DPMIN.ParticleType = str2num(ParticleType);
DPMIN.SourceEnergy = energy;
% DPMIN.NumParticles = nhist;
DPMIN.NumParticles = nhist*fieldSize*fieldSize;

source = DPMIN.SourceXYZ;
% scale 40cm/100cm back from the iso center.
% 0.5*0.4 = (half field size on the plus/minus sides; 0.5;
%           (distance from source-to-jaw: 40cm / isodistance 100cm = 0.4)
%  == 0.2
DPMIN.ClmtCorner1 = [source(1)-0.2*fieldSize    source(2)+40    source(3)+0.2*fieldSize];
DPMIN.ClmtCorner2 = [source(1)+0.2*fieldSize    source(2)+40    source(3)+0.2*fieldSize];
DPMIN.ClmtCorner3 = [source(1)+0.2*fieldSize    source(2)+40    source(3)-0.2*fieldSize];

% DPMIN.FlatFilterA = -0.21300000000000;
% DPMIN.FlatFilterB= -0.10000000000000;
% DPMIN.FlatFilterDist = 12.50000000000000;
DPMIN.FlatFilterA = str2num(FFA);
DPMIN.FlatFilterB= str2num(FFB);
DPMIN.FlatFilterDist = str2num(FFDistance);

DPMIN.IsoDistance = 100;
DPMIN.IncreaseFluence = str2num(IncreaseFluence);
DPMIN.OnlyHorn = str2num(OnlyHorn);
DPMIN.UseFlatFilter = str2num(UseFlatFilter)
DPMIN.DoseNormFlag = str2num(DoseNormFlag);

if(DPMIN.ParticleType == 0)

    for i = str2num(startBin): str2num(endBin)
        DPMIN.BinEnergy = [i*energy/numBin (i+1)*energy/numBin]
        DPMIN.UsePhotSpectrum = 0;
        DPMIN.OutputError = 0;
        DPMIN.RandSeeds = [rand*1000 rand*10000];

        if (DPMIN.ParticleType == 0)
            DPMIN.Prefix = 'pre4phot'
        elseif (DPMIN.ParticleType == -1)
            DPMIN.Prefix = 'pre4elec'
        else
            error('DPMIN.ParticleType has to be 0 or -1')
        end

        cd(path);
        tic; doseV = dpm(DPMIN, CTscan); toc
        cd(currentDir);

        dose3D = zeros(r,c,s);
        %   Change the convention of the filename, using the upper bound of the
        %   bin as the file name, instead of the lower bound.
        %    filename = ['doseV_', num2str(i*energy/numBin),'MV.mat'];
        dose3D(indV) = doseV;
        filename = ['dose3D_', num2str((i+1)*energy/numBin),'MV.mat'];
        if (DPMIN.OnlyHorn == 1)
            filename =  ['dose3D_Horn_', num2str((i+1)*energy/numBin),'MV.mat'];
		end

		% Jan 11, 2008. After the extension, UseFlatFilter can be 0, 1, or 2.
        if (DPMIN.UseFlatFilter ~= 0)
            filename =  ['dose3D_FF_', num2str((i+1)*energy/numBin),'MV.mat'];
        end

        save (filename, 'dose3D');

    end

elseif (DPMIN.ParticleType == -1)
    DPMIN.Prefix = 'pre4elec';
    DPMIN.UseFlatFilter = 0

    cd(path);
    tic; doseV = dpm(DPMIN, CTscan); toc
    cd(currentDir);

    dose3D = zeros(r,c,s);
    dose3D(indV) = doseV;
    filename = ['dose3D_elec.mat'];
    save(filename,'dose3D');

else
    error('DPMIN.ParticleType has to be 0 or -1')
end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% LFReadESLF.m - Read ESLF-encoded light fields
%
% Usage:
%   LF = LFReadESLF( FName )
%   LF = LFReadESLF( FName, [LensletSize_pix], [LoadAlpha], [<imread options>] )
%
% This function reads the ESLF files commonly used to store light fields, e.g. those at
% lightfields.stanford.edu.
%
% Inputs:
%   FName: path to input file
%
% Optional Inputs:
%   LensletSize_pix: number of pixels per lenslet, default 14
%
%   LoadAlpha: default true; if true, aattempts to load an alpha channel from the ESLF and use as
%             a weight channel for the LF; useful for indicating empty subaperture images; use in
%             conjunction with png or other formats with alpha channels
%
%  <imread options>: additional options are passed to imread; useful examples: 'BackgroundColor'
%
% Outputs:
%   LF: a 4D light field with index order [t,s,v,u,c], where s,t are horizontal and vertical
%   subaperture index; u,v are horizontal and vertical pixel index; and c is colour channel.
%
% User guide: LFToolbox.pdf
% See also: LFWriteESLF, LFReadGantryArray, LFUtilDecodeLytroFolder

% Copyright (c) 2013-2020 Donald G. Dansereau

function LF = LFReadESLF( FName, LensletSize_pix, LoadAlpha, varargin )

LensletSize_pix = LFDefaultVal('LensletSize_pix', [14,14]);
LoadAlpha = LFDefaultVal('LoadAlpha', true);
NChans = 3;

if( LoadAlpha )
	[Img,ColorMap, Alpha] = imread( FName, varargin{:} ); % note: requesting Alpha sets bg to black
	if( ~isempty(Alpha) )
		Img(:,:,NChans+1) = Alpha;
		NChans = NChans + 1;
	end
else
	[Img,ColorMap] = imread( FName, varargin{:} );
end

if( ~isempty(ColorMap) )  % indexed colour image
	Img = ind2rgb(Img, ColorMap); % todo[optimization] may wish to convert to uint8
end

ImgSize = size(Img(:,:,1));

LFSize(3:4) = ceil(ImgSize./LensletSize_pix);
PadSize = LFSize(3:4) .* LensletSize_pix;
PadAmt = PadSize-ImgSize;

LFSize(1:2) = PadSize ./ LFSize(3:4);
LFSize(5) = NChans;

ImgPad = padarray(Img, PadAmt, 0,'post');

LF = reshape(ImgPad, LFSize([1,3,2,4,5]));
LF = permute(LF, [1,3,2,4,5]);

repo: LFToolbox


text:
classdef TestDrawKeypoints
    %TestDrawKeypoints

    properties (Constant)
        img = fullfile(mexopencv.root(),'test','fruits.jpg');
    end

    methods (Static)
        function test_1
            im = imread(TestDrawKeypoints.img);
            detector = cv.FeatureDetector('ORB');
            kpts = detector.detect(im);
            out = cv.drawKeypoints(im, kpts);
            validateattributes(out, {class(im)}, {'size',size(im)});
        end

        function test_2
            im = imread(TestDrawKeypoints.img);
            detector = cv.FeatureDetector('ORB');
            kpts = detector.detect(im);
            out = cv.drawKeypoints(im, kpts, ...
                'OutImage',im, 'Color',[255 0 0], 'DrawRichKeypoints',true);
            validateattributes(out, {class(im)}, {'size',size(im)});
        end

        function test_error_argnum
            try
                cv.drawKeypoints();
                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 [x_max y_max A]=crit_interp_g(y,x);
%fits a gaussian to three points: y=[y(x(1)) y(x(2)) y(x(3))]. Returns the
%position (x_max) and value (y_max) of the interpolated critical point
%(peak or trough). Things go sideways (complex) if y is mixed sign and
%rounding errors can cause spurious complex results if y is negative.
%If x is omitted, it is assumed to be [-1 0 1].
%y=A(1)*exp(-A(2)*(x-A(3)).^2)

%    Copyright Travis Wiens 2009 travis.mlfx@nutaksas.com

if nargin<2
    x=[-1 0 1];
end

lny=log(y);

%lny=1/denom*(a*x^2+b*x+c)
a =(x(3) * (lny(2) - lny(1)) + x(2) * (lny(1) - lny(3)) + x(1) * (lny(3) - lny(2)));
b =(x(3)*x(3) * (lny(1) - lny(2)) + x(2)*x(2) * (lny(3) - lny(1)) + x(1)*x(1) * (lny(2) - lny(3)));
c =(x(2) * x(3) * (x(2) - x(3)) * lny(1) + x(3) * x(1) * (x(3) - x(1)) * lny(2) + x(1) * x(2) * (x(1) - x(2)) * lny(3));

%y=A*exp(-B*(x-C)^2);
C=-b/(2*a);
x_max=C;

if nargout>1
    denom = (x(1) - x(2)) * (x(1) - x(3)) * (x(2) - x(3));
    A=exp(c/denom-b*b/(4*a*denom));
    y_max=A;
    if nargout>2
        B=-a/denom;
        A=[A B C];
    end
end


text:
% Test of subdivision of a base polyhedron.

mesh_types = {'tetra','oct','ico'};
nmesh = length(mesh_types);
nsub = 3;

clf;
for i=1:nmesh
    mesh_type = mesh_types{i};
    [vertex,face] = compute_base_mesh(mesh_type);
    for s=0:nsub
        subplot(nmesh,nsub+1,s+1+(nsub+1)*(i-1));
        plot_mesh(vertex,face);
        lighting flat;
        if s~=nsub
			options.spherical = 1;
            options.relaxation = 3;
            [vertex,face] = perform_mesh_subdivision(vertex,face,1,options);
        end
    end
end

author: gpeyre

repo: matlab-toolboxes


text:
function [ vX, mX ] = SolveLsLInfProx( mA, vB, paramLambda, numIterations )
% ----------------------------------------------------------------------------------------------- %
%[ vX, mX ] = SolveLsL1Prox( mA, vB, lambdaFctr, numIterations )
% Solve L Infinity Regularized Least Squares Using Proximal Gradient (PGM) Method.
% Input:
%   - mA                -   Input Matirx.
%                           The model matrix.
%                           Structure: Matrix (m X n).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - vB                -   input Vector.
%                           The model known data.
%                           Structure: Vector (m X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - paramLambda       -   Parameter Lambda.
%                           The L Infinity Regularization parameter.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: (0, inf).
%   - numIterations     -   Number of Iterations.
%                           Number of iterations of the algorithm.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range {1, 2, ...}.
% Output:
%   - vX                -   Output Vector.
%                           Structure: Vector (n X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References
%   1.  Wikipedia PGM - https://en.wikipedia.org/wiki/Proximal_gradient_method.
% Remarks:
%   1.  Using vanilla PGM.
% Known Issues:
%   1.  A
% TODO:
%   1.  B
% Release Notes:
%   -   1.0.000     23/08/2017
%       *   First realease version.
% ----------------------------------------------------------------------------------------------- %

mAA = mA.' * mA;
vAb = mA.' * vB;
vX  = pinv(mA) * vB; %<! Dealing with "Fat Matrix"

stopThr = 1e-6;

stepSize = 1 / (2 * (norm(mA, 2) ^ 2));
% stepSize = 1 / sum(mA(:) .^ 2); %<! Faster to calculate, conservative (Hence slower)

mX = zeros([size(vX, 1), numIterations]);
mX(:, 1) = vX;

for ii = 2:numIterations

    vG = (mAA * vX) - vAb;
    vV = vX - (stepSize * vG);
    vX = vV - (stepSize * paramLambda * ProjectL1Ball(vV / (paramLambda * stepSize), 1, stopThr));

    mX(:, ii) = vX;

end

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function [number_of_objectives, number_of_decision_variables, min_range_of_decesion_variable, max_range_of_decesion_variable] = objective_description_function()

%% function [number_of_objectives, number_of_decision_variables, min_range_of_decesion_variable, max_range_of_decesion_variable] = objective_description_function()
% This function is used to completely describe the objective functions and
% the range for the decision variable space etc. The user is prompted for
% inputing the number of objectives, numebr of decision variables, the
% maximum and minimum range for each decision variable and finally the
% function waits for the user to modify the evaluate_objective function to
% suit their need.

%  Copyright (c) 2009, Aravind Seshadri
%  All rights reserved.

% g = sprintf('Input the number of objective: ');
% Obtain the number of objective function
number_of_objectives = input('Input the number of objective: '); % modified by zzb
if number_of_objectives < 2
    error('This is a multi-objective optimization function hence the minimum number of objectives is two');
end
% g = sprintf('\nInput the number of decision variables: ');
% Obtain the number of decision variables
number_of_decision_variables = input('Input the number of decision variables: '); % modified by zzb
min_range_of_decesion_variable = input('Input the array of minimum value for decision variable: '); % modified by zzb
while(length(min_range_of_decesion_variable) ~= number_of_decision_variables)
    min_range_of_decesion_variable = input('The size is wrong. \nInput the array of minimum value for decision variable: ');
end
max_range_of_decesion_variable = input('Input the array of maximum value for decision variable: '); % modified by zzb
while(length(max_range_of_decesion_variable) ~= number_of_decision_variables)
    max_range_of_decesion_variable = input('The size is wrong. \nInput the array of maximum value for decision variable: ');
end
clc;
% for i = 1 : number_of_decision_variables
%     clc
%     g = sprintf('\nInput the minimum value for decision variable %d : ', i);
%     % Obtain the minimum possible value for each decision variable
%     min_range_of_decesion_variable(i) = input(g);
%     g = sprintf('\nInput the maximum value for decision variable %d : ', i);
%     % Obtain the maximum possible value for each decision variable
%     max_range_of_decesion_variable(i) = input(g);
%     clc
% end
open('D:\ProgramFiles\MATLAB\toolbox\genetic\NSGA-II\evaluate_objective.m');  % modified by zzb
g = sprintf('Now edit the function named "evaluate_objective" appropriately to match your needs.\nMake sure that the number of objective functions and decision variables match your numerical input. \nMake each objective function as a corresponding array element. \nAfter editing do not forget to save. \nPress "c" and enter to continue... ');
% Prompt the user to edit the evaluate_objective function and wait until
% 'c' is pressed.
x = input(g, 's');
if isempty(x)
    x = 'x';
end
while x ~= 'c'
    clc;
    x = input(g, 's');
    if isempty(x)
        x = 'x';
    end
end

author: sfvsfv

repo: Mathematical-modeling


text:
%  Figure 10.28      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
% fig10_28.m is a script to generate Fig. 10.28,
% the step response of the collocated
% design for the satellite with PD compensation
% response at theta-2
clf;
m=[1, 0.1]; k0=[0, 0.091] ; d0=[0, 0.0036]; k1=[0, 0.4];
[f,g,h,j]=twomass(m,k0,d0); [f1,g,h,j] = twomass(m,k1,d0);
h1=[0, 0, 1, 0];

nc1=0.25*[2, 1];
dc1=[1/40, 1];

[ac,bc,cc,dc]=tf2ss(nc1, dc1);
[aol,bol,col,dol]= series(ac, bc,cc,dc,f,g,h1,j);
[acl]=aol-bol*col;
ccl2=[h,0*cc]
[aol1,bol1,col1,dol1] = series(ac,bc,cc,dc,f1,g,h1,j);
acl1=aol1-bol1*col1;
t=0:0.5:100;
sys=ss(acl,bol,ccl2,dol);
[y,t]=step(sys,t);
plot(t,y);
xlabel('Time (sec)');
ylabel('Amplitude');
grid on;
title('Fig. 10.28: Response at \theta_2 of the collocated design')