text:
function ret = checkSolver(instr,err)
%CHECKSOLVER  See if a solver is installed on your PC
%
%   checkSolver() prints a list of OPTI interfaced solvers, their version
%   numbers and and whether they are available on your PC.
%
%   checkSolver('matrix') prints a solver vs problem type matrix for
%   determining which solvers and setup for which problems.
%
%   checkSolver('config') prints a solver vs config matrix for
%   determining which options are configured for which solvers.
%
%   checkSolver(type) prints a list of solvers for a given problem, 'type',
%   together with a table of constraint and derivative information. See
%   below for heading information. Examples are 'lp', 'milp', etc.
%
%   ok = checkSolver(solver) returns true if the passed solver is available
%   on your PC. Examples are 'clp', 'ipopt', etc.
%
%   checkSolver(solver) displays information about the solver (if
%   available) and configuration specifics (no output args).
%
%   solver = checkSolver(bestSolver) returns the best available solver as a
%   string for a given problem type. Examples are 'best_lp', 'best_milp',
%   etc.
%
%   Problem Type table headings are as follows:
%
%   BD - Bounds
%   LI - Linear Inequality (* indicates a conversion will take place)
%   LE - Linear Equality   (* indicates a conversion will take place)
%   QI - Quadratic Inequality
%   NI - Nonlinear Inequality
%   NE - Nonlinear Equality
%   SP - Supports Sparse Matrices
%   D1 - Requires First Derivatives (Gradient) (* indicates optional)
%   D2 - Requires Second Derivatives (Hessian) (* indicates optional)
%   GL - Global Solver

%   Copyright (C) 2011-2013 Jonathan Currie (I2C2)

global SLE LP MILP BILP QP QCQP MIQP MIQCQP SDP MISDP SNLE SCNLE NLS UNO NLP MINLP DESC PTYPES

%Ordered lists of 'best' solvers
SLE = {'matlab','MUMPS'};
LP = {'CPLEX','clp','scip','qsopt','ooqp','glpk','lp_solve','MATLAB','lipsol','csdp','dsdp','SeDuMi'};
MILP = {'CPLEX','cbc','scip','glpk','lp_solve',};
BILP = {'CPLEX','cbc','scip','glpk','lp_solve','MATLAB'};
QP = {'CPLEX','clp','ooqp','scip','ipopt','MATLAB'};
QCQP = {'CPLEX','scip','ipopt'};
MIQP = {'CPLEX','scip','bonmin'};
MIQCQP = {'CPLEX','scip','bonmin'};
SDP = {'csdp','dsdp','SeDuMi'};
MISDP = {};
SNLE = {'nl2sol','hybrj','mkltrnls','lmder','levmar','MATLAB','IPOPT','nomad','nlopt'};
SCNLE = {'nl2sol','mkltrnls','levmar','IPOPT','nomad','pswarm','nlopt'};
NLS = {'nl2sol','mkltrnls','lmder','levmar','MATLAB','IPOPT','nomad','pswarm','nlopt'};
UNO = {'IPOPT','matlab','m1qn3','NLOPT','nomad','scip','GMATLAB'};
NLP = {'IPOPT','MATLAB','FILTERSD','NLOPT','LBFGSB','nomad','PSwarm','scip','GMATLAB'};
MINLP = {'BONMIN','nomad','scip','GMATLAB'};
%Description List
DESC = containers.Map({'cplex','glpk','lp_solve','clp','matlab','mumps','ooqp','ipopt','nlopt','qsopt','cbc','bonmin','lbfgsb','levmar','lmder','nl2sol','hybrj',...
                       'mkltrnls','pswarm','gmatlab','nomad','scip','m1qn3','filtersd','dsdp','csdp','sedumi','lipsol'},...
                      {'IBM ILOG CPLEX','GNU Linear Programming Kit','LP_Solve','COIN-OR Linear Programming',...
                       'Optimization Toolbox: ','a MUltifrontal Massively Parallel sparse direct Solver','Object-Orientated Quadratic Programming',...
                       'Interior-Point Optimization','Nonlinear Optimization','QSopt Library','COIN-OR Branch and Cut',...
                       'Basic Open-source Nonlinear Mixed INteger programming','Limited Memory Broyden-Fletcher-Goldfarb-Shanno Bounded Optimization',...
                       'Levenberg-Marquardt Nonlinear Least Squares in C/C++','MINPACK Levenberg-Marquardt','Adaptive Nonlinear Least Squares','MINPACK Powell Hybrid',...
                       'Intel MKL Trust Region Nonlinear Least Squares','Pattern and Particle Swarm Global Optimizer','Global Optimization Toolbox: ',...
                       'Nonlinear Optimization with the MADS algorithm','Solving Constraint Integer Programs','Large-Scale Unconstrained L-BFGS Minimization',...
                       'FilterSD: Filter & Trust Region Nonlinear Optimization','Dual Scaling Interior Point SDP Solver','Predictor-Corrector Primal-Dual SDP Solver',...
                       'Self-Dual-Minimization','Linear programming Interior-Point SOLvers'});
%Problem Type List
PTYPES = {'sle','lp','milp','bilp','qp','qcqp','miqp','miqcqp','sdp','misdp','snle','scnle','nls','uno','nlp','minlp'};

%Default Input Args
if(nargin < 2), err = 1; end %default to create an error
if(nargin < 1), instr = []; end
%Default Internal Args
ret = [];

if(isempty(instr))
    printSolvers();     %print a list of solvers and availability
else
    %Check if we have a known solver
    allS = getAllSolvers();
    for i = 1:length(allS)
        if(strcmpi(instr,allS{i}))
            ret = check(instr,err);
            if(nargout == 0 && ret), printSolverInfo(instr); end
            return;
        end
    end

    %If we got here, check for a known problem type
    for i = 1:length(PTYPES)
        if(strcmpi(instr,PTYPES{i}))
            ret = findAllForProb(instr,nargout);
            return;
        end
    end

    %Still not found, check standard options and alternative names for solvers
    switch(lower(instr))
        %Options
        case 'matrix'
            printSolverMatrix;
        case 'config'
            printConfigMatrix;
        case 'ptypes'
            ret = PTYPES;
        case 'ver'
            ret = genSolverList();
        %Alternative Solver Names
        case 'auto'
            ret = 1; %assume we can find a solver
        case 'optim'
            ret = check('matlab',err);
            if(nargout == 0 && ret), printSolverInfo('matlab'); end
        case {'mldivide','fminsearch'}
            ret = 1; %base matlab
        case 'l-bfgs-b'
            ret = check('lbfgsb',err);
            if(nargout == 0 && ret), printSolverInfo('lbfgsb'); end
        case 'lpsolve'
            ret = check('lp_solve',err);
            if(nargout == 0 && ret), printSolverInfo('lp_solve'); end
        case {'lmdif','lm_der','lm_dif'}
            ret = check('lmder',err);
            if(nargout == 0 && ret), printSolverInfo('lmder'); end
        case 'nl2sno'
            ret = check('nl2sol',err);
            if(nargout == 0 && ret), printSolverInfo('nl2sol'); end
        case 'hybrd'
            ret = check('hybrj',err);
            if(nargout == 0 && ret), printSolverInfo('hybrd'); end
        case 'mkltr'
            ret = check('mkltrnls',err);
            if(nargout == 0 && ret), printSolverInfo('mkltrnls'); end
        case 'baron'
            ret = check('baron',err);

        %Check for all, best
        otherwise
            if(strfind(instr,'best'))
                ret = findBestSolver(instr);
            elseif(strfind(instr,'all'))
                ret = findAllForProb(instr,nargout);
            else
                error('Unknown Solver or Option: %s',instr); %not found
            end
    end
end

function ok = check(name,err)
%Check a particular solver is installed
switch(name)
    case 'MATLAB'
        c = which('fmincon'); %using ver is quite slow
        if(isempty(c))
            if(err)
                error('The MATLAB Optimization Toolbox must be installed');
            else
                ok = 0;
            end
        else
            ok = 1;
        end
        return
    case {'GMATLAB','gmatlab'}
        c = which('patternsearch'); %using ver is quite slow
        if(isempty(c))
            if(err)
                error('The MATLAB Global Optimization Toolbox must be installed');
            else
                ok = 0;
            end
        else
            ok = 1;
        end
        return
    case {'matlab'} %always installed (mldivide)
        ok = 1;
        return
    case {'CPLEX','cplex'} %m file only
        c = which(name);
        %Problems with OPTI + R2013a x64 (really can't find the issue currently [19/5/13])
        if(~isempty(c))
            if(~isempty(strfind(computer,'64')))
                v = ver('matlab');
                if(isfield(v,'Release') && strcmpi(v.Release,'(R2013a)'))
                    c = []; %force not available (due to crashes in cplexlink)
                end
%                 cp = Cplex; %still problems crashing! arrh
%                 solVer = cp.getVersion;
%                 clear cp;
%                 if(strcmp(solVer,'12.5.0.0'))
%
%                 end
            end
        end
    case {'MOSEK','mosek'} %find mex file
        c = which(['mosekopt.' mexext]);
    case {'sedumi','SEDUMI','SeDuMi'}
        %Check matlab file
        c = which('sedumi');
        %Check mex file exists as well
        if(isempty(which(['eigK.' mexext])))
            c = [];
        %Check for v1.31 or higher
        elseif(~isempty(c))
            [~,major,minor] = getVerFromVerFile('sedumi');
            if(major == 1 && minor < 31)
                c = [];
            elseif(major < 1)
                c = [];
            end
        end
    case {'BARON','baron'}
        c = which('baron.m');
    case {'LIPSOL','lipsol'}
        c = which(['ls_blkslv.' mexext]);
    otherwise
        c = which([name '.' mexext]); %bug fix 22/2/12
end

if(isempty(c))
    if(err)
        throwAsCaller(MException('opti:checkSolver','OPTI cannot find a valid version of %s. Ensure it is installed on your system.',name));
    else
        ok = 0;
    end
else
    ok = 1;
end

function str = findBestSolver(prob)
%Find the best available solver for a given problem type

global SLE LP MILP BILP QP QCQP MIQP MIQCQP SDP MISDP SNLE SCNLE NLS UNO NLP MINLP PTYPES %#ok<NUSED>

s = regexp(prob,'_','split');
if(length(s) < 2)
    error('When requesting the best solver it should be of the form best_xxx, e.g. best_lp');
end

%Remove 'D' from dynamic problems
if(s{2}(1) == 'D')
    s{2} = s{2}(2:end);
end

fnd = 0;
for i = 1:length(PTYPES)
    if(strcmpi(s{2},PTYPES{i}))
        try
            list = eval(upper(s{2}));
            fnd = 1;
        catch %#ok<CTCH>
            fnd = 0; %global var doesn't exist?
        end
        break;
    end
end
%If not found, we don't have a solver for this problem
if(~fnd)
    error('Could not find a problem type match for type %s',upper(s{2}));
%Otherwise return the best solver available
else
    str = [];
    for i = 1:length(list)
        if(check(list{i},0))
            str = list{i};
            break;
        end
    end
    if(isempty(str))
        error('You do not have a solver which can solve a %s explicitly',upper(s{2}));
    end
end

function list = findAllForProb(prob,nout)
%Find all available solvers for a given problem type

global SLE LP MILP BILP QP QCQP MIQP MIQCQP SDP MISDP SNLE SCNLE NLS UNO NLP MINLP PTYPES %#ok<NUSED>

%Backwards compatibility check
if(strfind(prob,'all_'))
    s = regexp(prob,'_','split');
    prob = s{2};
end

fnd = 0;
for i = 1:length(PTYPES)
    if(strcmpi(prob,PTYPES{i}))
        try
            plist = eval(upper(prob));
            fnd = 1;
        catch %#ok<CTCH>
            fnd = 0; %global var doesn't exist?
        end
        break;
    end
end
%If found keep availale solvers
if(fnd)
    j = 1; list = [];
    for i = 1:length(plist)
        if(check(plist{i},0))
            list{j} = plist{i}; %#ok<AGROW>
            j = j + 1;
        end
    end
%Check if we want all...
else
    if(strcmpi(prob,'all'))
        list = getAllSolvers();
        return;
    else
        error('Unknown problem type: %s\nNote the API has changed in v1.34 to just the problem type, e.g. checkSolver(''lp'')',prob);
    end
end

%If not returning anything, just print
if(~nout)
    printSelSolvers(list,upper(prob))
    list = [];
end

function list = getAllSolvers()
%Create a list of all solvers OPTI is interfaced to (not neccesarily found however).
global SLE LP MILP BILP QP QCQP MIQP MIQCQP SDP MISDP SNLE SCNLE NLS UNO NLP MINLP

fulllist = [SLE LP MILP BILP QP QCQP MIQP MIQCQP SDP MISDP SNLE SCNLE NLS UNO NLP MINLP];
%Remove duplicates
fulllist = unique(lower(fulllist));
%Sort & Return
list = sort(upper(fulllist));

function printSolvers()
%Print each solver and if it exists

fprintf('\n------------------------------------------------\n')
fprintf('OPTI AVAILABLE SOLVERS:\n\n');
list = genSolverList();
for i = 1:length(list)
    fprintf('%s\n',list{i});
end

fprintf('------------------------------------------------\n')

function list = genSolverList()
%Generate a display compatible list of all solvers, their availability and version number
solv = getAllSolvers();
list = cell(length(solv),1);
for i = 1:length(solv)
    [avail,solVer] = getSolverVer(solv{i});
    list{i} = sprintf('%-10s       %-13s  %s',[solv{i} ':'],avail,solVer);
end

function printSelSolvers(solv,type)
%Print selected solvers

fprintf('\n------------------------------------------------\n')
fprintf(['OPTI ' upper(type) ' SOLVERS:\n']);
if(isempty(solv))
    fprintf('None Available\n');
else
    [~,hd] = getSolverInfo(solv{1},type);
    len = length(strtrim(hd));
    fprintf(' %s\n',hd); j = 1;
    for i = 1:length(solv)
        if(check(solv{i},0)) %only print if the solver is available
            fprintf(['(%2d) %-10s  %-'  num2str(len+4) 's %s\n'],j,upper(solv{i}),getSolverInfo(solv{i},type),getDesc(solv{i},type));
            j = j + 1;
        end
    end
end
fprintf('------------------------------------------------\n')

function printSolverMatrix()
%Print solver vs problem matrix

global SLE LP BILP MILP QP QCQP MIQP MIQCQP SDP MISDP SNLE SCNLE NLS UNO NLP MINLP

fprintf('\n------------------------------------------------\n')
fprintf('OPTI SOLVER vs PROBLEM MATRIX:\n\n');
fprintf(' %-8s  %-5s %-4s %-5s %-6s %-4s %-5s %-5s %-4s %-4s %-4s %-5s %-6s %-5s %-5s %-4s %-4s\n','Solver','SLE','LP','BILP','MILP','QP','QCQP','MIQP','MIQCQP','SDP','MISDP','SNLE','SCNLE','NLS','UNO','NLP','MINLP');

solv = getAllSolvers();
list = {SLE LP BILP MILP QP QCQP MIQP MIQCQP SDP MISDP SNLE SCNLE NLS UNO NLP MINLP};

for i = 1:length(solv)
    fprintf('%-8s |',solv{i});
    sok = checkSolver(solv{i},0);
    msok = check('MATLAB',0);
    for j = 1:length(list)
        if(any(strcmpi(solv{i},list{j})))
            if(strcmpi(solv{i},'matlab')) %matlab has optimization toolbox + base packages
                if(msok || any(strcmp({'matlab'},list{j})))
                    fprintf('  x  |');
                else
                    fprintf('  o  |');
                end
            else
                if(sok)
                    fprintf('  x  |');
                else
                    fprintf('  o  |');
                end
            end
        else
            fprintf('     |');
        end
    end
    fprintf('\n');
end

fprintf('------------------------------------------------\n')

function printConfigMatrix()
%Print solver vs config matrix

%[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf]
fprintf('\n------------------------------------------------\n')
fprintf('OPTI SOLVER vs CONFIG MATRIX:\n\n');
fprintf('   %-11s %-7s %-7s %-7s %-7s %-7s %-6s %-8s %-7s %-7s %-6s %-7s\n','Solver','miter','meval','mnode','mtime','tolr','tola','tolint','disp','sOpts','ctrlc','iterfnc');

solv = getAllSolvers();

for i = 1:length(solv)
    info = optiSolverInfo(solv{i},'');
    %Separate out variables we use
    miter   = getMarker(info.opt.miter);
    meval   = getMarker(info.opt.meval);
    mnode   = getMarker(info.opt.mnode);
    mtime   = getMarker(info.opt.mtime);
    tolr    = getMarker(info.opt.tolr);
    tola    = getMarker(info.opt.tola);
    tolint  = getMarker(info.opt.tolint);
    disp    = getMarker(info.opt.disp);
    sOpts   = getMarker(info.opt.opts);
    ctrlc   = getMarker(info.opt.ctrlc);
    iterf   = getMarker(info.opt.iterf);

    %Print line
    fprintf('%-11s  |   %1s   |   %1s   |   %1s   |   %1s   |   %1s   |   %1s   |   %1s   |   %1s   |   %1s   |   %1s   |   %1s   |\n',solv{i},miter,meval,mnode,mtime,tolr,tola,tolint,disp,sOpts,ctrlc,iterf);
end

fprintf('------------------------------------------------\n')

function str = getDesc(solv,type)
%Return solver description
global DESC

str = DESC(lower(solv));
if(strcmp(solv,'MATLAB'))
    str = [str getMatlabDesc(type)];
elseif(strcmp(solv,'GMATLAB'))
    str = [str getGMatlabDesc(type)];
elseif(strcmp(solv,'matlab'))
    str = getMatlabDesc(type);
end

function str = getMatlabDesc(type)
%Return MATLAB solver

switch(lower(type))
    case 'sle'
        str = 'mldivide';
    case 'lp'
        str = 'linprog';
    case 'qp'
        str = 'quadprog';
    case 'bilp'
        str = 'bintprog';
    case 'nlp'
        str = 'fmincon';
    case 'snle'
        str = 'fsolve';
    case 'nls'
        str = 'lsqnonlin';
    case 'uno'
        str = 'fminsearch';
    otherwise
        str = '';
end

function str = getGMatlabDesc(type)
%Return Global MATLAB Solver

switch(lower(type))
    case {'uno','nlp'}
        str = 'patternsearch';
    case 'minlp'
        str = 'ga';
    otherwise
        str = '';
end

function [avail,solVer] = getSolverVer(solver)
%Return currently installed solver version

ok = check(solver,0);
if(ok)
    avail = 'Available';
    switch(lower(solver))
        %Interfaced solvers
        case 'cplex'
            c = Cplex;
            solVer = ['v' c.getVersion];
            clear c;
        case 'matlab'
            c = ver('optim');
            if(isfield(c,'Version'))
                solVer = ['v' c.Version];
            else
                solVer = '';
            end
            clear c;
        case 'gmatlab'
            c = ver('globaloptim');
            if(isfield(c,'Version'))
                solVer = ['v' c.Version];
            else
                solVer = '';
            end
            clear c;
        case 'mosek'
            solVer = ''; %unknown for the moment
        case 'sedumi'
            solVer = getVerFromVerFile('sedumi');
        case 'lipsol'
            solVer = getVerFromVerFile('lipsol');
        case 'opti'
            c = optiver;
            solVer = ['v' num2str(c)];
        %Supplied Solvers
        otherwise
            try
                v = feval(lower(solver));
                if(~strcmp(v,''))
                    solVer = ['v' v];
                else
                    solVer = '';
                end
            catch %#ok<CTCH>
                solVer = '';
            end
    end
else
    avail = 'Not Available';
    solVer = '';
end

function [str,hd] = getSolverInfo(solver,type)
%Generate a string of solver information for a particular solver

%Get Solver Information
info = optiSolverInfo(solver,type);
%Separate out variables we use
bnd     = getMarker(info.con.bnd);
lineq   = getMarker(info.con.lineq);
leq     = getMarker(info.con.leq);
qineq   = getMarker(info.con.qineq);
sdp     = getMarker(info.con.sdp);
nlineq  = getMarker(info.con.nlineq);
nleq    = getMarker(info.con.nleq);
der1    = getMarker(info.der1);
der2    = getMarker(info.der2);
glb     = getMarker(info.global);
sp      = getMarker(info.sparse);

lhs = '%-17s';
switch(lower(type))
    case {'lp','milp','qp','miqp'}
        hd = sprintf([lhs 'BD LI LE SP'],'');
        str = sprintf(' %-2s %-2s %-2s %-2s',bnd,lineq,leq,sp);
    case 'bilp'
        hd = sprintf([lhs 'LI LE SP'],'');
        str = sprintf(' %-2s %-2s %-2s',lineq,leq,sp);
    case {'qcqp','miqcqp'}
        hd = sprintf([lhs 'BD LI LE QI SP'],'');
        str = sprintf(' %-2s %-2s %-2s %-2s',bnd,lineq,leq,qineq,sp);
    case {'sdp','misdp'}
        hd = sprintf([lhs 'BD LI LE QI SD SP'],'');
        str = sprintf(' %-2s %-2s %-2s %-2s %-2s',bnd,lineq,leq,qineq,sdp,sp);
    case 'nls'
        hd = sprintf([lhs 'BD LI LE D1'],'');
        str = sprintf(' %-2s %-2s %-2s %-2s',bnd,lineq,leq,der1);
    case {'nlp','minlp'}
        hd = sprintf([lhs 'BD LI LE QI NI NE SP D1 D2 GL'],'');
        str = sprintf(' %-2s %-2s %-2s %-2s %-2s %-2s %-2s',bnd,lineq,leq,qineq,nlineq,nleq,sp,der1,der2,glb);
    otherwise %no constraints
        hd = '';
        str = '';
end

function str = getMarker(value)
switch(value)
    case 1
        str = 'x';
    case 0
        str = '';
    case -1
        str = '*';
end

%Print Solver Information
function printSolverInfo(solver)
fprintf('\n-----------------------------------------------------------\n');
%Print heading
[~,solVer] = getSolverVer(solver);
fprintf('OPTI SOLVER: %s [%s]\n',upper(solver),solVer);
fprintf('-----------------------------------------------------------\n');
%Get solver info
info = optiSolverInfo(solver);
%Print Special Info
if(info.der1 || info.der2 || info.sparse || info.global || info.multialg)
   fprintf('Solver Features:\n');
   if(info.parallel == 1), fprintf('\t- Multi-Threaded (Parallel) Solver\n'); end
   if(info.multialg == 1), fprintf('\t- Contains Multiple Algorithms\n'); end
   if(info.whitebox == 1), fprintf('\t- White-Box NL Solver Which Exploits Problem Structure\n'); end
   if(info.global == 1)
       fprintf('\t- Solves Global Optimization Problems\n');
   elseif(info.global == -1)
       fprintf('\t- Some Algorithms Solve Global Optimization Problems\n');
   end
   if(info.con.int), fprintf('\t- Solves Mixed Integer Problems\n'); end
   if(info.sparse), fprintf('\t- Solves Sparse Problems\n'); end
   if(info.der1 == 1)
       fprintf('\t- Requires First Derivatives\n');
   elseif(info.der1 == -1)
       if(info.multialg==1)
            fprintf('\t- Some Algorithms Require First Derivatives\n');
       else
            fprintf('\t- Can Accept First Derivatives (Or Solver Will Approximate)\n');
       end
   end
   if(info.der2 == 1)
       fprintf('\t- Requires Second Derivatives\n');
   elseif(info.der2 == -1)
       fprintf('\t- Can Accept Second Derivatives (Or Solver Will Approximate)\n');
   end
   fprintf('-----------------------------------------------------------\n');
end
%Print Constraint Info
if(info.con.con)
    fprintf('Supported Constraints:\n');
    printCon('Decision Variable Bounds',info.con.bnd);
    printCon('Linear Inequalities',info.con.lineq);
    printCon('Linear Equalities',info.con.leq);
    printCon('Quadratic Inequalities',info.con.qineq);
    printCon('Quadratic Equalities',info.con.qeq);
    printCon('Semidefinite (Linear Matrix Inequalities)',info.con.sdp);
    printCon('Nonlinear Inequalities',info.con.nlineq);
    printCon('Nonlinear Equalities',info.con.nleq);
    printCon('Integer Variables',info.con.int);
    fprintf('-----------------------------------------------------------\n');
end
%See if we have any unique info
u=struct2cell(info.unq); uu = 0;
for i = 1:length(u)
    if(ischar(u{i}) || any(u{i} ~= 0))
        uu = 1;
        break;
    end
end
%Print Unique Info
if(uu)
    fprintf('Special Requirements / Options:\n');
    if(info.unq.x0 == 1)
        fprintf('\t- Requires an initial solution guess via x0\n');
    elseif(info.unq.x0 == -1)
        fprintf('\t- Can accept an initial solution guess via x0\n');
    end
    switch(info.unq.qpform)
        case 'sym'
            fprintf('\t- QP H Matrix must be symmetric\n');
        case 'tril'
            fprintf('\t- QP H Matrix must be symmetric lower triangular\n');
        case 'triu'
            fprintf('\t- QP H Matrix must be symmetric upper triangular\n');
    end
    switch(info.unq.hform)
        case 'sym'
            fprintf('\t- Nonlinear Hessian (if supplied) must be symmetric\n');
        case 'tril'
            fprintf('\t- Nonlinear Hessian (if supplied) must be symmetric lower triangular\n');
        case 'triu'
            fprintf('\t- Nonlinear Hessian (if supplied) must be symmetric upper triangular\n');
    end
    if(info.unq.jstr), fprintf('\t- Requires the structure of the Jacobian\n'); end
    if(info.unq.hstr), fprintf('\t- Requires the structure of the Hessian (if supplied)\n'); end
    fprintf('-----------------------------------------------------------\n');
end

%Print info from MEX file
switch(lower(solver))
    case {'cplex','mosek','matlab','gmatlab','sedumi','lipsol'}
        %Nothing
    case {'filtersd'}
        fprintf('MEX FILE DETAILS:');
        eval([lower(solver) 'sp']);
    otherwise
        fprintf('MEX FILE DETAILS:');
        eval(lower(solver));
end

function printCon(con,val)
if(val == 1)
    fprintf('\t- %s\n',con);
elseif(val == -1)
    fprintf('\t- %s (but may be converted internally)\n',con);
end

function [solVer,major,minor] = getVerFromVerFile(solverP)
solVer = ''; major = 0; minor = 0;
%Find sedumi.m and Version .txt
p = which(solverP); p2 = which('Version.txt','-all');
if(~isempty(p) && ~isempty(p2))
    %Find final file sep in solver path
    ind = strfind(p,filesep);
    %For each version path, see if we can find a matching solver path
    for i = 1:length(p2)
        ind2 = strfind(p2{i},filesep);
        %compare paths up to find file sep
        if(strcmp(p(1:ind(end)),p2{i}(1:ind2(end))))
            %open and reading version
            fid = fopen(p2{i});
            if(fid > 0)
                v = fscanf(fid,'%d.%d');
                fclose(fid);
                major = v(1); minor = v(2);
                solVer = sprintf('v%d.%d',v(1),v(2));
            else
                fclose(fid);
            end
            break;
        end
    end
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [varargout] = read_nex5(filename, varargin)

% READ_NEX5 reads header or data from a Nex Technologies *.nex5 file,
% which is a file containing action-potential (spike) timestamps and waveforms
% (spike channels), event timestamps (event channels), and continuous
% variable data (continuous A/D channels).
%
% LFP and spike waveform data that is returned by this function is
% expressed in microVolt.
%
% Use as
%   [hdr] = read_nex5(filename)
%   [dat] = read_nex5(filename, ...)
%   [dat1, dat2, dat3, hdr] = read_nex5(filename, ...)
%
% Optional arguments should be specified in key-value pairs and can be
%   header      structure with header information
%   feedback    0 or 1
%   tsonly      0 or 1, read only the timestamps and not the waveforms
%   channel     number, or list of numbers (that will result in multiple outputs)
%   begsample   number (for continuous only)
%   endsample   number (for continuous only)
%
% See also READ_NEX5_HEADER
%
% Copyright (C) 2020 Robert Oostenveld, Alex Kirillov
%
% 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$

% parse the optional input arguments
hdr       = ft_getopt(varargin, 'header');
channel   = ft_getopt(varargin, 'channel');
feedback  = ft_getopt(varargin, 'feedback', false);
tsonly    = ft_getopt(varargin, 'tsonly', false);
begsample = ft_getopt(varargin, 'begsample', 1);
endsample = ft_getopt(varargin, 'endsample', inf);

% start with empty return values and empty data
varargout = {};

if isempty(hdr)
  if feedback, fprintf('reading header from %s\n', filename); end
  hdr = read_nex5_header(filename);
  if hdr.FileHeader.NumVars<1
    ft_error('no channels present in file');
  end
end

% use Matlab for automatic byte-ordering
fid = fopen_or_error(filename, 'r', 'ieee-le');

for i=1:length(channel)
  chan = channel(i);
  vh = hdr.VarHeader(chan);
  clear buf
  status = fseek(fid, vh.DataOffset, 'bof');
  if status < 0;  ft_error('error with fseek');  end
  switch vh.Type
    case 0
      % Neurons, only timestamps
      buf.ts = Nex5ReadTimestamps(fid, vh);

    case 1
      % Events, only timestamps
      buf.ts = Nex5ReadTimestamps(fid, vh);

    case 2
      % Interval variables
      buf.begs = Nex5ReadTimestamps(fid, vh);
      buf.ends = Nex5ReadTimestamps(fid, vh);

    case 3
      % Waveform variables
      buf.ts = Nex5ReadTimestamps(fid, vh);
      if ~tsonly
        if vh.ContDataType == 0
          buf.dat = fread(fid, [vh.NumberOfDataPoints vh.Count], 'int16');
        else
          buf.dat = fread(fid, [vh.NumberOfDataPoints vh.Count], 'float32');
        end
        % convert the AD values to miliVolt, subsequently convert from miliVolt to microVolt
        buf.dat = (buf.dat * vh.ADtoUnitsCoefficient + vh.UnitsOffset) * 1000;
      end

    case 4
      % Population vector
      ft_error('population vectors are not supported');

    case 5
      % Continuously recorded variables
      buf.ts = Nex5ReadTimestamps(fid, vh);
      if vh.ContIndexOfFirstPointInFragmentDataType == 0
        buf.indx = fread(fid, [1 vh.Count], 'uint32');
      else
        buf.indx = fread(fid, [1 vh.Count], 'uint64');
      end
      if vh.Count>1 && (begsample~=1 || endsample~=inf)
        ft_error('reading selected samples from multiple AD segments is not supported');
      end
      if ~tsonly
        numsamples = min(endsample - begsample + 1, vh.NumberOfDataPoints);
        status = fseek(fid, (begsample-1)*2, 'cof');
        if status < 0;  ft_error('error with fseek');  end
        if vh.ContDataType == 0
          buf.dat  = fread(fid, [1 numsamples], 'int16');
        else
          buf.dat  = fread(fid, [1 numsamples], 'float32');
        end
        % convert the AD values to miliVolt, subsequently convert from milliVolt to microVolt
        buf.dat = (buf.dat * vh.ADtoUnitsCoefficient + vh.UnitsOffset) * 1000;
      end

    case 6
      % Markers
      buf.ts = Nex5ReadTimestamps(fid, vh);
      for j=1:vh.NumberOfMarkerFields
        buf.MarkerNames{j,1} = fread(fid, [1 64], 'uint8=>char');
        if vh.MarkerDataType == 0
          for k=1:vh.Count
            buf.MarkerValues{j,k} = fread(fid, [1 vh.MarkerLength], 'uint8=>char');
          end
        else
          buf.MarkerValues{j} = fread(fid, vh.Count, 'uint32');
        end
      end

    otherwise
      ft_error('incorrect channel type');
  end % switch channel type

  % return the data of this channel
  varargout{i} = buf;
end % for channel

% always return the header as last
varargout{end+1} = hdr;

fclose(fid);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction to read nex5 timestamps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ts = Nex5ReadTimestamps(fid, varHeader)
if varHeader.TimestampDataType == 0
  ts = fread(fid, [1 varHeader.Count], 'int32');
else
  ts = fread(fid, [1 varHeader.Count], 'int64');
end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function plot_events(wevent, filepattern, figureOutDirectory)
numEvents = numel(wevent);
for eventNumber=1:numEvents
    fprintf('- plotting event %d of %d\n',eventNumber,numEvents );
    plot_panels(wevent{eventNumber});
    ah=get(gcf,'Children');
    title(ah(7),sprintf('Event %d',eventNumber));
    outfile = sprintf('%s/%s%03d.png',figureOutDirectory, filepattern, eventNumber);
    feval('print', '-dpng', outfile);
    close
end

author: geoscience-community-codes

repo: GISMO


text:
function out = CO_Embed2_Basic(y,tau)
% CO_Embed2_Basic Point density statistics in a 2-d embedding space.
%
% Computes a set of point-density statistics in a plot of y_i against y_{i-tau}.
%
% INPUTS:
% y, the input time series.
%
% tau, the time lag (can be set to 'tau' to set the time lag the first zero
%                       crossing of the autocorrelation function).
%
% Outputs include the number of points near the diagonal, and similarly, the
% number of points that are close to certain geometric shapes in the y_{i-tau},
% y_{tau} plot, including parabolas, rings, and circles.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

if nargin < 2
    tau = 1;
end

doPlot = false; % plot outputs to a figure

if strcmp(tau,'tau')
	% Make tau the first zero crossing of the autocorrelation function
    tau = CO_FirstCrossing(y,'ac',0,'discrete');
end

xt = y(1:end-tau); % part of the time series
xtp = y(1+tau:end); % time-lagged time series
N = length(y) - tau; % Length of each time series subsegment

% Points in a thick bottom-left -- top-right diagonal
out.updiag01 = sum(abs(xtp-xt) < 0.1)/N;
out.updiag05 = sum(abs(xtp-xt) < 0.5)/N;

% Points in a thick bottom-right -- top-left diagonal
out.downdiag01 = sum(abs(xtp+xt) < 0.1)/N;
out.downdiag05 = sum(abs(xtp+xt) < 0.5)/N;

% Ratio of these
out.ratdiag01 = out.updiag01/out.downdiag01;
out.ratdiag05 = out.updiag05/out.downdiag05;

% In a thick parabola concave up
out.parabup01 = sum(abs(xtp-xt.^2) < 0.1)/N;
out.parabup05 = sum(abs(xtp-xt.^2) < 0.5)/N;

% In a thick parabola concave down
out.parabdown01 = sum(abs(xtp+xt.^2) < 0.1)/N;
out.parabdown05 = sum(abs(xtp+xt.^2) < 0.5)/N;

% In a thick parabola concave up, shifted up 1
out.parabup01_1 = sum(abs(xtp-(xt.^2+1)) < 0.1)/N;
out.parabup05_1 = sum(abs(xtp-(xt.^2+1)) < 0.5)/N;

% In a thick parabola concave down, shifted up 1
out.parabdown01_1 = sum(abs(xtp+(xt.^2-1)) < 0.1)/N;
out.parabdown05_1 = sum(abs(xtp+(xt.^2-1)) < 0.5)/N;

% In a thick parabola concave up, shifted down 1
out.parabup01_n1 = sum(abs(xtp-(xt.^2-1)) < 0.1)/N;
out.parabup05_n1 = sum(abs(xtp-(xt.^2-1)) < 0.5)/N;

% In a thick parabola concave down, shifted down 1
out.parabdown01_n1 = sum(abs(xtp+(xt.^2+1)) < 0.1)/N;
out.parabdown05_n1 = sum(abs(xtp+(xt.^2+1)) < 0.5)/N;

% RINGS (points within a radius range)
out.ring1_01 = sum(abs(xtp.^2+xt.^2-1) < 0.1)/N;
out.ring1_02 = sum(abs(xtp.^2+xt.^2-1) < 0.2)/N;
out.ring1_05 = sum(abs(xtp.^2+xt.^2-1) < 0.5)/N;

% CIRCLES (points inside a given circular boundary)
out.incircle_01 = sum(xtp.^2+xt.^2 < 0.1)/N;
out.incircle_02 = sum(xtp.^2+xt.^2 < 0.2)/N;
out.incircle_05 = sum(xtp.^2+xt.^2 < 0.5)/N;
out.incircle_1 = sum(xtp.^2+xt.^2 < 1)/N;
out.incircle_2 = sum(xtp.^2+xt.^2 < 2)/N;
out.incircle_3 = sum(xtp.^2+xt.^2 < 3)/N;
out.medianincircle = median([out.incircle_01, out.incircle_02, out.incircle_05 ...
                            out.incircle_1, out.incircle_2, out.incircle_3]);
out.stdincircle = std([out.incircle_01, out.incircle_02, out.incircle_05 ...
                        out.incircle_1, out.incircle_2, out.incircle_3]);

if doPlot
    figure('color','w'); box('on');
    plot(xt,xtp,'.k');
    hold on
    r = (xtp.^2+xt.^2 < 0.2);
    plot(xt(r),xtp(r),'.g')
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [params, names] = diagKernExtractParam(kern)

% DIAGKERNEXTRACTPARAM Extract parameters from the DIAG kernel structure.
% FORMAT
% DESC extracts parameters from the diagonal noise covariance function
% kernel structure into a vector of parameters for optimisation.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel structure, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
%
% FORMAT
% DESC extracts parameters and parameter names from the diagonal noise covariance function
% kernel structure.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel structure, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
% RETURN names : cell array of strings containing names for each
% parameter.
%
% SEEALSO diagKernParamInit, diagKernExpandParam, kernExtractParam, scg, conjgrad
%
% COPYRIGHT : Neil D. Lawrence, 2011
%
% KERN

params = kern.variance;
if nargout > 1
  names = {'variance'};
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [t,amps,data,aux] = read_intan_data(filename)

% [t,amps,data,aux] = read_intan_data
%
% Opens file selection GUI to select and then read data from an Intan
% amplifier data file (*.int).
%
% t = time vector (in seconds)
% amps = vector listing active amplifier channels
% data = matrix of electrode-referred amplifier signals (in microvolts)
% aux = matrix of six auxiliary TTL input signals
%
% Example usage:
%  >> [t,amps,data,aux] = read_intan_data;
%  >> plot(t,data(:,1));
%
% Version 1.1, June 26, 2010
% (c) 2010, Intan Technologies, LLC
% For more information, see http://www.intantech.com
% For updates and latest version, see http://www.intantech.com/software.html
%
% 06-22-10 Added GUI file selection and optimized: Craig Patten, Plexon, Inc.

% use MATLAB predefined gui uigetfile to select the file(s) to analyze
if ~exist('filename','var')
    [file, path, filterindex] = uigetfile('*.int','Select a .int file','MultiSelect', 'off');
    filename = [path,file];
end

fid = fopen(filename, 'r');

% Read first three header bytes encoding file version
for i=1:3
    header(i) = fread(fid, 1, 'uint8');
end

if (header(1) ~= 128)
    error('Improper data file format.');
end

if (header(2) ~= 1 || header(3) ~= 1)
    warning('Data file version may not be compatible with this m-file.');
end

% Now see which amplifier channels are saved in this file.
for i=1:64
    amp_on(i) = fread(fid, 1, 'uint8');
end

num_amps = sum(amp_on);

% Create a list of amplifier channels in this file.
amps = zeros(1,num_amps);
index = 1;
for i=1:64
    if (amp_on(i) == 1)
        amps(index) = i;
        index = index + 1;
    end
end

% Now search for the end of the file to find out the length of the data.
% t_count = 0;
% while (~feof(fid))
%    fread(fid, 1+4*num_amps, 'uint8');
%    t_count = t_count + 1;
% end
% t_count = t_count - 1;
% t_max = t_count/25000;

%-----------------------------------
% replace above code with a more efficient method CDP 06-24-10
s = dir(filename);
filesize = s.bytes;
t_count = (filesize - 67)/(num_amps*4 + 1);
t_max = t_count/25000;
%-----------------------------------

% print channel (singular) when there is only one channel! CDP 06-24-10
if num_amps == 1;
    fprintf(1, '\nData file contains %0.2f seconds of data from %d amplifier channel.\n', t_max, num_amps);
    fprintf(1, 'Channel: ');
else
    fprintf(1, '\nData file contains %0.2f seconds of data from %d amplifier channels.\n', t_max, num_amps);
    fprintf(1, 'Channels: ');
end

for i=1:num_amps
    fprintf(1, '%d ', amps(i));
end
fprintf(1, '\n\n');

% Pre-allocate large data matrices.
aux = zeros(t_count,6,'uint8');
t = (0:1:(t_count-1))/25000;
t = t';

%--------------------------------------
% Replace code code below with much faster code CDP 06-24-10
% Go back to the beginning of the file...
frewind(fid);

% ...skip the header this time...
fread(fid, 3+64, 'uint8');

% allocate space to read the entire file
data2 = zeros((filesize-67),1,'uint8');
% read the entire file
data2 = fread(fid,(filesize-67),'uint8=>uint8');

% extract the digital data
aux_data = data2((num_amps*4)+1:num_amps*4+1:filesize-67);

% extract individual bits
aux = [bitget(aux_data,6),bitget(aux_data,5),bitget(aux_data,4),bitget(aux_data,3),bitget(aux_data,2),bitget(aux_data,1)];
clear aux_data;

% delete the digital data
data2((num_amps*4)+1:num_amps*4+1:filesize-67) = [];

% convert the remaining data from bytes to single
data2 = typecast(data2,'single');

data = zeros(t_count,num_amps);
% de-mux the channels
for ind = 1:num_amps
    data(:,ind) = data2(ind:num_amps:length(data2));
end
%---------------------------------------

% % Go back to the beginning of the file...
% frewind(fid);
%
% % ...skip the header this time...
% fread(fid, 3+64, 'uint8');
%
% % ...and read all the data.
% fprintf(1, 'Reading data...  (This may take a while.)\n\n');
% for i=1:t_count
%     for j=1:num_amps
%         data(i,j) = double(fread(fid, 1, 'float32'));
%     end
%
%     aux_byte = fread(fid, 1, 'uint8');
%
%     % Decode auxiliary TTL inputs
%     if aux_byte >= 32
%         aux(i,6) = 1;
%         aux_byte = aux_byte - 32;
%     end
%     if aux_byte >= 16
%         aux(i,5) = 1;
%         aux_byte = aux_byte - 16;
%     end
%     if aux_byte >= 8
%         aux(i,4) = 1;
%         aux_byte = aux_byte - 8;
%     end
%     if aux_byte >= 4
%         aux(i,3) = 1;
%         aux_byte = aux_byte - 4;
%     end
%     if aux_byte >= 2
%         aux(i,2) = 1;
%         aux_byte = aux_byte - 2;
%     end
%     if aux_byte >= 1
%         aux(i,1) = 1;
%         aux_byte = aux_byte - 1;
%     end
% end

% Close file, and we're done.
fclose(fid);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function V=covnw(data,nlag,demean)
% Long-run covariance estimation using Newey-West (Bartlett) weights
%
% USAGE:
%   [V] = covnw(DATA)
%   [V] = covnw(DATA,NLAG,DEMEAN)
%
% INPUTS:
%   DATA   - T by K vector of dependent data
%   NLAG   - Non-negative integer containing the lag length to use.  If empty or not included,
%              NLAG=min(floor(1.2*T^(1/3)),T) is used
%   DEMEAN - Logical true or false (0 or 1) indicating whether the mean should be subtracted when
%              computing the covariance
%
% OUTPUTS:
%   V      - A K by K covariance matrix estimated using Newey-West (Bartlett) weights
%
% COMMENTS:
%
% EXAMPLES:
%   Simulate an AR(1)
%       y = armaxfilter_simulate(1000,0,1,.9);
%   Newey-West covariance with automatic BW selection
%       lrcov = covnw(y)
%   Newey-West covariance with 10 lags
%       lrcov = covnw(y, 10)
%   Newey-West covariance with 10 lags and no demeaning
%       lrcov = covnw(y, 10, 0)
%
% See also COVVAR

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 5/1/2007

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T=size(data,1);
if nargin==1
    nlag=min(floor(1.2*T^(1/3)),T);
    demean=true;
elseif nargin==2
    demean=true;
end
if isempty(nlag)
    nlag=min(floor(1.2*T^(1/3)),T);
end
if isempty(demean)
    demean=true;
end
if ~ismember(demean,[0 1])
    error('DEMEAN must be either logical true or false.')
end
if floor(nlag)~=nlag || nlag<0
    error('NLAG must be a non-negative integer.')
end
if ndims(data)>2
    error('DATA must be a T by K matrix of data.')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if demean
    data=data-repmat(mean(data),T,1);
end

% NW weights
w=(nlag+1-(0:nlag))./(nlag+1);
% Start the covariance
V=data'*data/T;
for i=1:nlag
    Gammai=(data((i+1):T,:)'*data(1:T-i,:))/T;
    GplusGprime=Gammai+Gammai';
    V=V+w(i+1)*GplusGprime;
end

author: bashtage

repo: mfe-toolbox


text:
% Example 8.7: Floorplan generation test script
% Section 8.8.1/2, Boyd & Vandenberghe "Convex Optimization"
% Joelle Skaf 12/04/05
%
% Rectangles aligned with the axes need to be place in the smallest
% possible bounding box. No overlap is allowed. Each rectangle to be placed
% can be reconfigured, within some limits.
% In the current problem, 60 rectangles are to be place. We are given 2
% acyclic graphs H and V (for horizontal and vertical) that specify the
% relative positioning constraints of those rectangles.
% We are also given minimal areas for the rectangles and aspect ratio
% constraints

% input data
load data_floorplan_60;
rho = 1;
Amin = 100*ones(1,n);

[W, H, w, h, x, y] = floorplan(adj_H, adj_V, rho, Amin,ones(60,1)*0.5,ones(60,1)*2);
fill([0; W; W; 0],[0;0;H;H],[1 1 1]);           % bounding box
hold on
for i=1:n
    fill([x(i); x(i)+w(i); x(i)+w(i); x(i)],[y(i);y(i);y(i)+h(i);y(i)+h(i)],0.90*[1 1 1]);
    hold on;
    text(x(i)+w(i)/2, y(i)+h(i)/2,int2str(i));
end
axis([0 W 0 H]);
axis equal; axis off;

author: yu-jiang

repo: radpbook

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


text:
function demo(no_compile)

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2011-2012 Ross Girshick
% Copyright (C) 2008, 2009, 2010 Pedro Felzenszwalb, Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

startup;

if ~exist('no_compile', 'var')
  fprintf('compiling the code...');
  compile;
  fprintf('done.\n\n');
end

load('VOC2007/car_final');
model.vis = @() visualizemodel(model, ...
                  1:2:length(model.rules{model.start}));
test('000034.jpg', model, 1);

load('INRIA/inriaperson_final');
model.vis = @() visualizemodel(model, ...
                  1:2:length(model.rules{model.start}));
test('000061.jpg', model, 2);

load('VOC2007/person_grammar_final');
model.class = 'person grammar';
model.vis = @() visualize_person_grammar_model(model, 6);
test('000061.jpg', model, 2);

load('VOC2007/bicycle_final');
model.vis = @() visualizemodel(model, ...
                  1:2:length(model.rules{model.start}));
test('000084.jpg', model, 1);

function test(imname, model, num_dets)
cls = model.class;
fprintf('///// Running demo for %s /////\n\n', cls);

% load and display image
im = imread(imname);
clf;
image(im);
axis equal;
axis on;
title('input image');
disp('input image');
disp('press any key to continue'); pause;
disp('continuing...');

% load and display model
model.vis();
disp([cls ' model visualization']);
disp('press any key to continue'); pause;
disp('continuing...');

% detect objects
tic;
[ds, bs] = imgdetect(im, model, -1);
toc;
top = nms(ds, 0.5);
top = top(1:min(length(top), num_dets));
ds = ds(top, :);
bs = bs(top, :);
clf;
if model.type == model_types.Grammar
  bs = [ds(:,1:4) bs];
end
showboxes(im, reduceboxes(model, bs));
title('detections');
disp('detections');
disp('press any key to continue'); pause;
disp('continuing...');

if model.type == model_types.MixStar
  % get bounding boxes
  bbox = bboxpred_get(model.bboxpred, ds, reduceboxes(model, bs));
  bbox = clipboxes(im, bbox);
  top = nms(bbox, 0.5);
  clf;
  showboxes(im, bbox(top,:));
  title('predicted bounding boxes');
  disp('bounding boxes');
  disp('press any key to continue'); pause;
end

fprintf('\n');

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
% Copyright 2018 Marc René Schädler
%
% This file is part of the mobile hearing aid prototype project
% The the mobile hearing aid prototype project 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.
%
% The mobile hearing aid prototype project 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 the mobile hearing aid prototype project. If not, see http://www.gnu.org/licenses/.

function [signal, fs] = gensweep(type, variable, frequency)
  fs = 48000; % Hz
  reference_level = 93.979; % dB SPL
  sweep_duration = 0.250; % s
  sweep_width = [0.95 1.05];
  flank_duration = 0.010; % s
  sweep_samples = round(fs.*sweep_duration);
  flank_samples = round(fs.*flank_duration);

  % Generate stimulus
  if type > 0
    signal = sinesweepphase(sweep_samples, sweep_width.*frequency./fs, rand(1).*2.*pi);
    signal = normalize(signal, variable - reference_level);
    signal = flank(signal, flank_samples);
  else
    signal = zeros(sweep_samples,1);
  end
end

author: m-r-s

repo: hearingaid-prototype


text:
function inoutsig = gaindb(inoutsig,gn,varargin)
%GAINDB  Increase/decrease level of signal
%   Usage:  outsig = gaindb(insig,gn);
%
%   `gaindb(insig,gn)` increases the energy level of the signal by *gn*
%   dB.
%
%   If *gn* is a scalar, the whole input signal is scaled.
%
%   If *gn* is a vector, each column is scaled by the entries in
%   *gn*. The length of *gn* must match the number of columns.
%
%   `gaindb(insig,gn,dim)` scales the signal along dimension *dim*.
%
%   See also: rms

%   AUTHOR: Peter L. Søndergaard, 2009

% ------ Checking of input parameters ---------

if nargin<2
  error('%s: Too few input parameters.',upper(mfilename));
end;

if ~isnumeric(inoutsig)
  error('%s: insig must be numeric.',upper(mfilename));
end;

if ~isnumeric(gn)
  error('%s: gn must be numeric.',upper(mfilename));
end;

definput.keyvals.dim=[];
[flags,kv]=ltfatarghelper({'dim'},definput,varargin);

% ------ Computation --------------------------

if isscalar(gn)
  inoutsig = inoutsig*10^(gn/20);
else
  if isvector(gn)
    M=length(gn);

    [inoutsig,L,Ls,W,dim,permutedsize,order]=...
        assert_sigreshape_pre(inoutsig,[],kv.dim,upper(mfilename));

    if M~=W
      error('%s: Length of gn and signal size must match.',upper(mfilename));
    end;

    for ii=1:W
      inoutsig(:,ii)=inoutsig(:,ii)*10^(gn(ii)/20);
    end;

    inoutsig=assert_sigreshape_post(inoutsig,kv.dim,permutedsize,order);

  else
    if ~isnumeric(gn)
      error('%s: gn must be a scalar or vector.',upper(mfilename));
    end;
  end;
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


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

function test_searchlight_singlethread()
    opt=struct();
    opt.progress=false;
    helper_test_searchlight(opt);

function test_searchlight_matlab_multithread()
    has_function=@(x)~isempty(which(x));
    has_parallel_toolbox=all(cellfun(has_function, {'gcp','parpool'}));

    if ~has_parallel_toolbox
        cosmo_notify_test_skipped('Matlab parallel toolbox not available');
        return;
    end

    warning_state=cosmo_warning();
    warning_resetter=onCleanup(@()cosmo_warning(warning_state));
    cosmo_warning('off');

    opt=struct();
    opt.progress=false;
    opt.nproc=2;
    helper_test_searchlight(opt);

function test_searchlight_octave_multithread()
    has_function=@(x)~isempty(which(x));
    has_parallel_toolbox=all(cellfun(has_function, {'parcellfun',...
                                                    'pararrayfun'}));

    if ~has_parallel_toolbox
        cosmo_notify_test_skipped('Octave parallel toolbox not available');
        return;
    end

    warning_state=cosmo_warning();
    warning_resetter=onCleanup(@()cosmo_warning(warning_state));
    cosmo_warning('off');

    opt=struct();
    opt.progress=false;
    opt.nproc=2;
    helper_test_searchlight(opt);

function helper_test_searchlight(opt)

    ds=cosmo_synthetic_dataset('size','normal');
    m=any(abs(ds.samples)>3,1);
    ds=cosmo_slice(ds,~m,2);
    ds=cosmo_dim_prune(ds);

    measure=@(x,a) cosmo_structjoin('samples',size(x.samples,2));
    nh=cosmo_spherical_neighborhood(ds,'radius',2,'progress',0);

    m=cosmo_searchlight(ds,nh,measure,opt);

    assertEqual(m.samples,[8 12 10 9 12 10 16 13 12 17 14 15 ...
                            13 11 15 14 10 9 14 11 5 7 6 7]);
    assertEqual(m.fa.i,ds.fa.i);
    assertEqual(m.fa.j,ds.fa.j);
    assertEqual(m.fa.k,ds.fa.k);
    assertEqual(m.a,ds.a);

    nh2=cosmo_spherical_neighborhood(ds,'count',17,'progress',0);
    m=cosmo_searchlight(ds,nh2,measure,opt);
    assertEqual(m.samples,[17 17 17 17 17 17 17 17 17 17 18 16 ...
                                17 17 16 15 17 17 17 17 17 17 17 17]);

    measure=@cosmo_correlation_measure;

    nh3=cosmo_spherical_neighborhood(ds,'radius',2,...
                                cosmo_structjoin('progress',0));
    m=cosmo_searchlight(ds, nh3, measure,...
                            'center_ids',[4 21],opt);

    assertVectorsAlmostEqual(m.samples, [0.9742,-.0273]...
                                        ,'relative',.001);
    assertEqual(m.fa.i,[1 1]);
    assertEqual(m.fa.j,[2 1]);
    assertEqual(m.fa.k,[1 5]);

    sa=struct();
    sa.time=(1:6)';
    sdim=struct();
    sdim.values={10:15};
    sdim.labels={'time'};

    nh4=cosmo_spherical_neighborhood(ds,'radius',0,'progress',false);
    measure2=@(x,opt)cosmo_structjoin('samples',mean(x.samples,2),...
                                       'sa',sa,...
                                       'a',cosmo_structjoin('sdim',sdim));
    m2=cosmo_searchlight(ds,nh4,measure2,opt);
    assertEqual(m2.sa,sa);
    assertEqual(m2.a.sdim,sdim);
    assertElementsAlmostEqual(m2.samples,ds.samples);

function test_searchlight_partial_classification
    ds=cosmo_synthetic_dataset('nchunks',6);

    partitions=struct();
    train_msk=mod(ds.sa.chunks,2)==0;
    partitions.train_indices={find(train_msk)};
    partitions.test_indices={find(~train_msk)};

    measure=@cosmo_crossvalidation_measure;
    opt=struct();
    opt.classifier=@cosmo_classify_lda;
    opt.partitions=partitions;
    opt.progress=false;
    opt.output='winner_predictions';

    nh=cosmo_spherical_neighborhood(ds,'radius',1,'progress',false);
    res=cosmo_searchlight(ds,nh,measure,opt);

    assertEqual(res.samples([3 4 7 8 11 12],:),NaN(6,6))
    s=res.samples([1 2 5 6 9 10],:);
    assertTrue(all(s(:)==1 | s(:)==2));

function test_searchlight_exceptions
    aet=@(varargin)assertExceptionThrown(@()...
                            cosmo_searchlight(varargin{:},...
                            'progress',0),'');
    ds=cosmo_synthetic_dataset();
    nh=cosmo_spherical_neighborhood(ds,'radius',1,'progress',false);
    measure=@(x,opt)cosmo_structjoin('samples',mean(x.samples,2));

    aet(struct,nh,measure);
    aet(ds,ds,measure);
    aet(ds,measure,nh);

    measure_bad=@(x,opt)cosmo_structjoin('samples',mean(x.samples,1));
    aet(ds,nh,measure_bad);

function test_searchlight_progress()
    if cosmo_skip_test_if_no_external('!evalc')
        return;
    end

    ds=cosmo_synthetic_dataset();
    nh=cosmo_spherical_neighborhood(ds,'count',2,'progress',false);
    measure=@(x,opt)cosmo_structjoin('samples',mean(x.samples,2));
    f=@()cosmo_searchlight(ds,nh,measure);
    res=evalc('f();');
    assert(~isempty(strfind(res,'[####################]')));

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
classdef SMD8 < PROBLEM
% <multi> <real> <bilevel>
% Bilevel optimization problems proposed by Sinha, Malo, and Deb
% maxFElower --- 500 --- Maximum number of lower level function evaluations for each solution

%------------------------------- Reference --------------------------------
% A. Sinha, P. Malo, K. Deb, Test problem construction for single-objective
% bilevel optimization, Evolutionary Computation, 2014, 22(3): 439-477.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties(SetAccess = private)
        maxFElower; % Maximum number of lower level function evaluations for each solution
        DU;         % Number of decision variables of the upper level
        DL;         % Number of decision variables of the lower level
        C;       	% Number of upper constraints
        p;          % The length of xu1
        q;          % The length of xl1
        r;          % The length of xu2 and xl2
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.maxFElower = obj.ParameterSet(500);
            obj.M = 2;
            obj.C = 0;
            if isempty(obj.D); obj.D = 5; end
            obj.DU = floor(obj.D/2);
            obj.DL = obj.D - obj.DU;
            obj.r  = floor(obj.DU/2);
            obj.p  = obj.DU - obj.r;
            obj.q  = obj.DL - obj.r;
            obj.lower    = -5*ones(1,obj.D);
            obj.upper    = 10*ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate upper level and lower level objective values
        function PopObj = CalObj(obj,PopDec)
            xu1   = PopDec(:,1:obj.p);
            xu2   = PopDec(:,obj.p+1:obj.p+obj.r);
            xl1   = PopDec(:,obj.p+obj.r+1:obj.p+obj.r+obj.q);
            xl2   = PopDec(:,obj.p+obj.r+obj.q+1:end);
            term2 = sum((xl1(:,2:end)-xl1(:,1:end-1).^2).^2+(xl1(:,1:end-1)-1).^2,2);
            % Upper level function value
            PopObj(:,1) = 20 + exp(1) - 20*exp(-0.2*sqrt(1/obj.p*sum(xu1.^2,2))) - exp(1/obj.p*sum(cos(2*pi*xu1),2)) - term2 + sum(xu2.^2,2) - sum((xu2-xl2.^3).^2,2);
            % Lower level function value
            PopObj(:,2) = sum(abs(xu1),2) + term2 + sum((xu2-xl2.^3).^2,2);
        end
        %% Calculate lower level objective values
        function llPopulation = EvaluationLower(obj,varargin)
            PopDec            = obj.CalDec(varargin{1});
            PopObj            = obj.CalObj(PopDec);
            PopObj(:,1)       = nan;
            PopCon            = obj.CalCon(PopDec);
            PopCon(:,1:obj.C) = nan;
            llPopulation      = SOLUTION(PopDec,PopObj,PopCon,varargin{2:end});
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function fe2d_d_fast_test ( )

%*****************************************************************************80
%
%% FE2D_D_FAST_TEST tests the FE2D_D_FAST code.
%
%  Discussion:
%
%    This function sets all parameter values and initial condition information
%    necessary to execute the "fast" version of the fe2d_d algorithm.
%
%  Licensing:
%
%    Copyright (C) 2014 Marcus R. Garvie.
%    See 'mycopyright.txt' for details.
%
%  Modified:
%
%    28 April 2014
%
%  Author:
%
%    Marcus R. Garvie.
%
%  Reference:
%
%    Marcus R Garvie, John Burkardt, Jeff Morgan,
%    Simple Finite Element Methods for Approximating Predator-Prey Dynamics
%    in Two Dimensions using MATLAB,
%    Submitted to Bulletin of Mathematical Biology, 2014.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FE2D_D_FAST_TEST:\n' );
  fprintf ( 1, '  Test the FE2D_D_FAST function\n' );
  fprintf ( 1, '  which applies Dirichlet boundary conditions as it\n' );
  fprintf ( 1, '  approximates a solution to a predator-prey system.\n' );
%
%  Set the parameters.
%
  alpha = 0.4;
  beta = 2.0;
  gamma = 0.6;
  delta = 1.0;
%
%  Use T=150.0 for normal run.
%  Use T=0.50 for a "quick" run that might take 15 minutes of computing.
%
% T = 150.0;
  T = 0.50;
  delt = 1.0 / 384.0;

  t = tic;
  fe2d_d_fast ( alpha, beta, gamma, delta, T, delt, @u0f, @v0f, @guf, @gvf );
  t = toc ( t );

  fprintf ( 1, '  Execution took %10.2g minutes \n', t / 60.0 );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FE2D_D_FAST_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

function value = u0f ( x, y )

%*****************************************************************************80
%
%% U0F evaluates the initial condition for U.
%
%  Licensing:
%
%    Copyright (C) 2014 Marcus R. Garvie.
%    See 'mycopyright.txt' for details.
%
%  Modified:
%
%    26 April 2014
%
%  Author:
%
%    Marcus R. Garvie.
%
%  Parameters:
%
%    Input, real X, Y, a location in the region.
%
%    Output, real VALUE, the initial condition for U at (X,Y).
%
  value = 6.0 / 35.0 - 2.0E-07 * ( x - 0.1 * y - 225.0 ) * ( x - 0.1 * y - 675.0 );

  return
end

function value = v0f ( x, y )

%*****************************************************************************80
%
%% V0F evaluates the initial condition for V.
%
%  Licensing:
%
%    Copyright (C) 2014 Marcus R. Garvie.
%    See 'mycopyright.txt' for details.
%
%  Modified:
%
%    26 April 2014
%
%  Author:
%
%    Marcus R. Garvie.
%
%  Parameters:
%
%    Input, real X, Y, a location in the region.
%
%    Output, real VALUE, the initial condition for V at (X,Y).
%
  value = 116.0 / 245.0 - 3.0E-05 * ( x - 450.0 ) - 1.2E-04 * ( y - 150.0 );

  return
end

function value = guf ( x, y, t )

%*****************************************************************************80
%
%% GUF evaluates the Dirichlet boundary condition for U.
%
%  Licensing:
%
%    Copyright (C) 2014 Marcus R. Garvie.
%    See 'mycopyright.txt' for details.
%
%  Modified:
%
%    28 April 2014
%
%  Author:
%
%    Marcus R. Garvie.
%
%  Parameters:
%
%    Input, real X, Y, a location on the boundary.
%
%    Input, real T, the time.
%
%    Output, real VALUE, the prescribed value for U at (X,Y,T).
%
  value = 0.0;

  return
end
function value = gvf ( x, y, t )

%*****************************************************************************80
%
%% GVF evaluates the Dirichlet boundary condition for V.
%
%  Licensing:
%
%    Copyright (C) 2014 Marcus R. Garvie.
%    See 'mycopyright.txt' for details.
%
%  Modified:
%
%    28 April 2014
%
%  Author:
%
%    Marcus R. Garvie.
%
%  Parameters:
%
%    Input, real X, Y, a location on the boundary.
%
%    Input, real T, the time.
%
%    Output, real VALUE, the prescribed value for V at (X,Y,T).
%
  value = 0.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [fileList] = getFileList(inputDir, filterExp,varargin)
	% Gathers a list of files based on an input regular expression.
	% Biafra Ahanonu
	% started: 2013.10.08 [11:02:31]
	% inputs
		% inputDir - directory to gather files from and regexp filter for files
		% filterExp - regexp used to find files
	% outputs
		% file list, full path

	[fileList] = ciapkg.io.getFileList(inputDir, filterExp,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function xcell = cellK(x,K)
% xcell = cellK(x,K)
%
% CELLK  Stores SeDuMi cone K-vector in cell-array format.
%
% On output xcell.f and xcell.l are the free and >=0 components,
% xcell.q{k}, xcell.r{k} and xcell.s{k} contain the Lorentz,
% Rotated Lorentz, and PSD-components, resp.
% xcell.s{k} is a K.s(k) x K.s(k) matrix.
%
% See also eigK, eyeK, sedumi.

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
%

% Remark: the cell-array scheme is inspired by SDPT3.
% Unlike in SDPT3 though, the cells are grouped in
% fields like xcell.q and xcell.s.
if nargin < 2
    error('cellK needs 2 input arguments')
end
if ~isstruct(K)
    error('K should be a structure')
end
i = 1;
if isfield(K,'f')
    xcell.f = x(i:i+K.f-1);
    i = i + K.f;
end
if isfield(K,'l')
    xcell.l = x(i:i+K.l-1);
    i = i + K.l;
end
% ------------------------------------------------------------
% Lorentz: this only works OUTSIDE sedumi, not for internal storage
% ------------------------------------------------------------
if isfield(K,'q')
    for k = 1: length(K.q)
        xcell.q{k} = x(i:i+K.q(k)-1);
        i = i + K.q(k);
    end
end
if isfield(K,'r')
    for k = 1: length(K.r)
        xcell.r{k} = x(i:i+K.r(k)-1);
        i = i + K.r(k);
    end
end
% ------------------------------------------------------------
% PSD-blocks:
% ------------------------------------------------------------
if isfield(K,'s')
    if ~isfield(K,'rsdpN')
        for k = 1: length(K.s)
            xcell.s{k} = mat(x(i:i+K.s(k)^2-1));
            i = i + K.s(k)^2;
        end
    else
        % ------------------------------------------------------------
        % Only applicable for internal storage: complex-as-real-storage
        % ------------------------------------------------------------
        for k = 1: K.rsdpN
            xcell.s{k} = mat(x(i:i+K.s(k)^2-1));
            i = i + K.s(k)^2;
        end
        j = sqrt(-1);
        for k = K.rsdpN+1:length(K.s)
            xcell.s{k} = mat(x(i:i+K.s(k)^2-1));
            i = i + K.s(k)^2;
            xcell.s{k} = xcell.s{k} + j * mat(x(i:i+K.s(k)^2-1));
            i = i + K.s(k)^2;
        end
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function param = copulaparam(type,tau)
%COPULAPARAM Copula parameter, given Kendall's rank correlation.
%   RHO = COPULAPARAM('Gaussian',TAU) returns the linear correlation
%   parameter RHO corresponding to a Gaussian copula having Kendall's rank
%   correlation TAU.  If TAU is a scalar correlation coefficient, RHO is a
%   scalar correlation coefficient corresponding to a bivariate copula.  If
%   TAU is a P-by-P correlation matrix, RHO is a P-by-P correlation matrix
%   corresponding to a P-variate copula.
%
%   RHO = COPULAPARAM('t',TAU) returns the linear correlation parameter RHO
%   corresponding to a t copula having Kendall's rank correlation TAU.  If
%   TAU is a scalar correlation coefficient, RHO is a scalar correlation
%   coefficient corresponding to a bivariate copula.  If TAU is a P-by-P
%   correlation matrix, RHO is a P-by-P correlation matrix corresponding to
%   a P-variate copula.
%
%   ALPHA = COPULAPARAM(TYPE,TAU) returns the copula parameter ALPHA
%   corresponding to a bivariate Archimedean copula having Kendall's rank
%   correlation TAU.  TYPE is one of 'Clayton', 'Frank', or 'Gumbel'.
%
%   Example:
%      % Determine the linear correlation coefficient corresponding to a
%      % bivariate Gaussian copula having a rank correlation of -0.5
%      tau = -0.5
%      rho = copulaparam('gaussian',tau)
%
%      % Generate dependent beta random values using that copula
%      u = copularnd('gaussian',rho,100)
%      b = betainv(u,2,2)
%
%      % Verify that those pairs have a sample rank correlation approximately
%      % equal to tau
%      tau_sample = kendall(b)

%   Written by Peter Perkins, The MathWorks, Inc.
%   Revision: 1.0  Date: 2003/09/05
%   This function is not supported by The MathWorks, Inc.
%
%   Requires MATLAB R13.

if nargin < 2
    error('Requires two input arguments.');
end

switch lower(type)
case {'gaussian' 't'}
    if ((numel(tau) == 1) && (tau < -1 | 1 < tau)) || ((numel(tau) ~= 1) && ~iscor(tau))
        error('TAU must be a correlation coefficient between -1 and 1, or a positive semidefinite correlation matrix.');
    end
    param = sin(tau.*pi./2);

case {'clayton' 'frank' 'gumbel'}
    if (numel(tau) ~= 1) || (tau < -1 | 1 < tau)
        error('TAU must be a correlation coefficient between -1 and 1.');
    end
    switch lower(type)
    case 'clayton'
        if tau < 0
            error('TAU must be nonnegative for the Clayton copula.');
        end
        param = 2*tau ./ (1-tau);
    case 'frank'
        if tau == 0
            param = 0;
        elseif abs(tau) < 1
            % There's no closed form for alpha in terms of tau, so alpha has to be
            % determined numerically.
            warn = warning('off','MATLAB:fzero:UndeterminedSyntax');
            param = fzero(@frankRootFun,sign(tau),[],tau);
            warning(warn);
        else
            param = sign(tau).*Inf;
        end
    case 'gumbel'
        if tau < 0
            error('TAU must be nonnegative for the Gumbel copula.');
        end
        param = 1 ./ (1-tau);
    end

otherwise
    error('Unrecognized copula type: ''%s''',type);
end

function err = frankRootFun(alpha,targetTau)
if abs(alpha) < realmin
    tau = 0;
else
    tau = 1 + 4 .* (debye1(alpha)-1) ./ alpha;
end
err = tau - targetTau;

author: Sable

repo: mcbench-benchmarks


text:
function [execResult, errMessage] = mysql_dbexecute(dbconn, sqlCommand)
% MYSQL_DBEXECUTE
%
% Retrieves data from a database connection, using either the Matlab
% database toolbox, or a java JDBC4 connection (depending on the input
% connection class).

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

switch class(dbconn)
case {'database','database.jdbc.connection'}
    % Database connection opened using database toolbox:
    curs = exec(dbconn,sqlCommand);
    errMessage = curs.Message;

    % Close cursor object:
    close(curs);

    if isempty(errMessage)
        % Success!:
        execResult = 0;
    else
        execResult = [];
    end

case 'com.mysql.jdbc.JDBC4Connection'
    % Database connection opened using java commands:
    try
        dbstmt = dbconn.createStatement();
        execResult = dbstmt.execute(sqlCommand);
        errMessage = [];
    catch emsg
        errMessage = emsg.message;
        execResult = [];
    end

otherwise
    error('Unknown database connection class %s',class(dbconn));
end

end


text:
function b = isposdef(a)
% ISPOSDEF   Test for positive definite matrix.
%    ISPOSDEF(A) returns 1 if A is positive definite, 0 otherwise.
%    Using chol is much more efficient than computing eigenvectors.

% Written by Tom Minka

[R,p] = chol(a);
b = (p == 0);

author: andrewssobral

repo: lrslibrary


text:
classdef FixedThre < ThreAlgs.ThreAlg
    properties
        thre;
    end
    methods
        function obj = FixedThre(thre)
            obj.algName = 'Segmentation using a fixed threshold';
            obj.thre = thre;
        end
        function CM = segment(obj, DI)
            if size(DI, 3) > 1
                cdMap2d = Utilities.mergeAvg(DI);
            end
            CM = (cdMap2d > obj.thre);
        end
    end
end

author: Bobholamovic

repo: ChangeDetectionToolbox


text:
function [ alphas, betas, triX, mask, xmin, ymin, npix ] = InitialisePieceWiseAffine( triangulation, sourcePoints )
%INITIALISEPIECEWICEAFFINE Summary of this function goes here
%   Detailed explanation goes here

    triangulation = triangulation + 1;
    numPoints = size(sourcePoints, 1);

    numTris = size(triangulation, 1);

    alphas = zeros(size(triangulation, 1), 3);
    betas = zeros(size(triangulation, 1), 3);

    xs = sourcePoints(:,1);
    ys = sourcePoints(:,2);

	for i = 1:numTris

        j = triangulation(i, 1);
        k = triangulation(i, 2);
        l = triangulation(i, 3);

        c1 = ys(l) - ys(j);
        c2 = xs(l) - xs(j);
        c4 = ys(k) - ys(j);
        c3 = xs(k) - xs(j);

        c5 = c3*c1 - c2*c4;

        alphas(i, 1) = (ys(j) * c2 - xs(j) * c1) / c5;
        alphas(i, 2) = c1/c5;
        alphas(i, 3) = -c2/c5;

        betas(i, 1) = (xs(j) * c4 - ys(j) * c3)/c5;
        betas(i, 2) = -c4/c5;
        betas(i, 3) = c3/c5;
    end

    xmin = min(xs);
    ymin = min(ys);

    xmax = max(xs);
    ymax = max(ys);

	w = int32(xmax - xmin + 1);
    h = int32(ymax - ymin + 1);

    mask = zeros(h, w);
    triX = zeros(h, w);

    shape = [xs, ys];

    for i = 1:h
        for j = 1:w

            currTri = findTriangle(double([double(j)-1+xmin, double(i)-1+ymin])', triangulation, shape);
            if(currTri ~= -1)
                triX(i, j) =  currTri - 1;
                mask(i, j) = 1;
            else
                triX(i, j) = -1;
            end
        end
    end
    npix = sum(sum(mask));
end

function [tri] = findTriangle(point, tris, controlPoints)

    numTris = size(tris, 1);
    tri = -1;

    for i=1:numTris

        if(PointInTriangle(point, controlPoints(tris(i,1),:)', controlPoints(tris(i,2),:)', controlPoints(tris(i,3),:)'))
           tri = i;
           break;
        end

    end

end

function inTriangle = PointInTriangle(point, v1, v2, v3)

    inTriangle = SameSide(point, v1,v2,v3) && SameSide(point, v2,v1,v3) && SameSide(point, v3,v1,v2);

end

function sameSide = SameSide(toTest,v1,v2,v3)

    x0 = toTest(1);
    y0 = toTest(2);

    x1 = v1(1);
    x2 = v2(1);
    x3 = v3(1);

    y1 = v1(2);
    y2 = v2(2);
    y3 = v3(2);

    x = (x3-x2)*(y0-y2) - (x0-x2)*(y3-y2);
    y = (x3-x2)*(y1-y2) - (x1-x2)*(y3-y2);
    if(x*y >= 0)
        sameSide = 1;
    else
        sameSide = 0;
    end
%     cross1 = cross( v3 - v2, toTest - v2);
%     cross2 = cross( v3 - v2, v1 - v2);
%
%     sameSide = (cross1 * cross2') >= 0;
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function fpc = eci2fpc1(gast, reci, veci)

% convert inertial state vector to flight path coordinates

% input

%  gast = greenwich apparent sidereal time (radians)
%  reci = inertial position vector (kilometers)
%  veci = inertial velocity vector (kilometers/second)

% output

%  fpc(1) = east longitude (radians)
%  fpc(2) = geocentric declination (radians)
%  fpc(3) = flight path angle (radians)
%  fpc(4) = azimuth (radians)
%  fpc(5) = position magnitude (kilometers)
%  fpc(6) = velocity magnitude (kilometers/second)

% global

%  omega = inertial rotation rate (radians/second)

% Orbital Mechanics with MATLAB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global omega

% compute geocentric radius

rmag = norm(reci);

% compute earth-relative position vector

c(1, 1) = cos(gast);
c(1, 2) = sin(gast);
c(1, 3) = 0.0d0;

c(2, 1) = -sin(gast);
c(2, 2) = cos(gast);
c(2, 3) = 0.0d0;

c(3, 1) = 0.0d0;
c(3, 2) = 0.0d0;
c(3, 3) = 1.0d0;

recf = c * reci';

% add earth rotation effect

vtmp(1) = veci(1) + reci(2) * omega;

vtmp(2) = veci(2) - reci(1) * omega;

vtmp(3) = veci(3);

% compute relative velocity vector and magnitude

vecf = c * vtmp';

vrel = norm(vecf);

% compute east longitude and geocentric declination

xlong = atan3(recf(2), recf(1));

decl = asin(recf(3) / rmag);

% compute flight path angle and azimuth

c(1, 1) = -sin(decl) * cos(xlong);
c(1, 2) = -sin(decl) * sin(xlong);
c(1, 3) = cos(decl);

c(2, 1) = -sin(xlong);
c(2, 2) = cos(xlong);
c(2, 3) = 0.0d0;

c(3, 1) = -cos(decl) * cos(xlong);
c(3, 2) = -cos(decl) * sin(xlong);
c(3, 3) = -sin(decl);

vr = c * vecf;

fpa = asin(-vr(3) / vrel);

azimuth = atan3(vr(2), vr(1));

% flight path coordinates

fpc(1) = xlong;

fpc(2) = decl;

fpc(3) = fpa;

fpc(4) = azimuth;

fpc(5) = rmag;

fpc(6) = vrel;


text:
hydro = struct();

hydro = readNEMOH(hydro,'../Ellipsoid/');
% hydro = readWAMIT(hydro,'../../WAMIT/Ellipsoid/ellipsoid.out',[]);
% hydro = combineBEM(hydro); % Compare to WAMIT
hydro = radiationIRF(hydro,10,[],[],[],[]);
hydro = radiationIRFSS(hydro,[],[]);
hydro = excitationIRF(hydro,15,[],[],[],[]);
writeBEMIOH5(hydro)
plotBEMIO(hydro)

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function outsig=rampsignal(insig,varargin)
%RAMPSIGNAL  Ramp signal
%   Usage: outsig=rampsignal(insig,L);
%
%   `rampsignal(insig,L)` applies a ramp function of length *L* to the
%   beginning and the end of the input signal. The default ramp is a
%   sinusoide starting from zero and ending at one (also known as a cosine
%   squared ramp).
%
%   If *L* is scalar, the starting and ending ramps will be of the same
%   length. If *L* is a vector of length 2, the first entry will be used
%   for the rising ramp, and the second for the falling.
%
%   If the input is a matrix or an N-D array, the ramp will be applied
%   along the first non-singleton dimension.
%
%   `rampsignal(insig)` will use a ramp length of half the signal.
%
%   `rampsignal(insig,L,wintype)` will use another window for ramping. This
%   may be any of the window types from |firwin|. Please see the help on
%   |firwin| for more information. The default is to use a piece of the
%   Hann window.
%
%   `rampsignal` accepts the following optional parameters:
%
%     'dim',d   Apply the ramp along dimension d. The default value of []
%               means to use the first non-singleton dimension.
%
%   See also: rampdown, rampsignal, firwin

definput.import={'firwin'};
definput.keyvals.dim=[];
definput.keyvals.L=[];
[flags,kv]=ltfatarghelper({'L','dim'},definput,varargin);

[insig,L,Ls,W,dim,permutedsize,order]=assert_sigreshape_pre(insig,[],kv.dim,'RAMPSIGNAL');
% Note: Meaning of L has changed, it is now the length of the signal.

switch numel(kv.L)
 case 0
  L1=L/2;
  L2=L/2;
 case 1
  L1=kv.L;
  L2=kv.L;
 case 2
  L1=kv.L(1);
  L2=kv.L(2);
 otherwise
  error('%s: The length must a scalar or vector.',upper(mfilename));
end;

if rem(L1,1)~=0 || rem(L2,1)~=0
  error('The length of the ramp must be an integer.');
end;

if L<L1+L2
  error(['%s: The length of the input signal must be greater than the length of the ramps ' ...
         'combined.'],upper(mfilename));
end;

r1=rampup(L1,flags.wintype);
r2=rampdown(L2,flags.wintype);

ramp=[r1;ones(L-L1-L2,1);r2];

% Apply the ramp
for ii=1:W
  insig(:,ii)=insig(:,ii).*ramp;
end;

outsig=assert_sigreshape_post(insig,dim,permutedsize,order);


text:
clear all
uri = 'ip:192.168.86.35';

%% Turn on DDS
tx = adi.AD9361.Tx('uri',uri);
tx.DataSource = 'DDS';
tx.SamplingRate = 20e6;
tx.RFBandwidth = 20e6;
toneFreq = 4e6;
tx.DDSFrequencies = repmat(toneFreq,2,4);
tx.AttenuationChannel0 = -10;
tx();
pause(1);
%% Set up fastlock profiles
freqs = 2.4e9+(1:8).*1e6;
for f = 0:length(freqs)-1
    % Update LO
    tx.CenterFrequency = freqs(f+1);
    pause(1);
    % Save profile
    tx.setAttributeLongLong('altvoltage1','fastlock_store',int64(f),true,8);
    fprintf('Saving profile at LO %d\n',freqs(f+1));
end
% Set pin control mode
tx.setDebugAttributeBool('adi,tx-fastlock-pincontrol-enable',1);
tx.setAttributeLongLong('altvoltage1','fastlock_recall',int64(0),true,8);
%% Configure Hop mode
h = FrequencyHopper;
h.uri = uri;
h.DwellSamples = 20000;
h.HoppingEnable = true;
h();
%% Capture
rx = adi.AD9361.Rx('uri',uri);
rx.SamplingRate = 20e6;
rx.CenterFrequency = 2.4e9+0.5e6;
rx.RFBandwidth = 20e6;
rx.SamplesPerFrame = 1e6;
rx.channelCount = 2;
rx.kernelBuffersCount = 1;
for k=1:10
    valid = false;
    while ~valid
        [out, valid] = rx();
    end
end
rx.release();
%% Plot
spectrogram(double(out),128,80,100,20e6,'yaxis','centered')
view(-77,72)
shading interp
colorbar off

author: analogdevicesinc


text:
function [ resultg, resultk, error_est ] = sum_sub_gk ( func, a, b, nsub, ...
  ng, wg, nk, xk, wk )

%*****************************************************************************80
%
%% SUM_SUB_GK carries out a composite Gauss-Kronrod rule.
%
%  Discussion:
%
%    The integral:
%
%      integral ( a <= x <= b ) f(x) dx
%
%    The quadrature rule:
%
%      H = ( B - A ) / NSUB
%      XMID(J) = A + 0.5 * H * ( 2 * J - 1 )
%
%      Sum ( 1 <= J <= NSUB )
%        Sum ( 1 <= I <= NK )
%          WK(I) * F ( XMID(J) + 0.5 * H * XK(I) )
%
%    The Gauss-Legendre weights should be computed by LEGCOM or LEGSET.
%    The Kronrod abscissas and weights should be computed by KRONSET.
%
%    The orders of the Gauss-Legendre and Kronrod rules must satisfy
%    NK = 2 * NG + 1.
%
%    The Kronrod rule uses the abscissas of the Gauss-Legendre rule,
%    plus more points, resulting in an efficient and higher order estimate.
%
%    The difference between the Gauss-Legendre and Kronrod estimates
%    is taken as an estimate of the error in the approximation to the
%    integral.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, external FUNC, the name of the function which
%    evaluates the integrand.  The function must have the form
%      function func ( x ).
%
%    Input, real A, B, the lower and upper limits of integration.
%
%    Input, integer NSUB, the number of equal subintervals into
%    which the finite interval (A,B) is to be subdivided for
%    higher accuracy.  NSUB must be at least 1.
%
%    Input, integer NG, the order of the Gauss-Legendre rule.
%    NG must be at least 1.
%
%    Input, real WG(NG), the weights of the
%    Gauss-Legendre rule.
%
%    Input, integer NK, the order of the Kronrod rule.
%    NK must be at least 1.
%
%    Input, real XK(NK), the abscissas of the
%    Kronrod rule.
%
%    Input, real WK(NK), the weights of the
%    Kronrod rule.
%
%    Output, real RESULTG, the approximate value of the
%    integral based on the Gauss-Legendre rule.
%
%    Output, real RESULTK, the approximate value of the
%    integral based on the Kronrod rule.
%
%    Output, real ERROR_EST, an estimate of the approximation
%    error.  This is computed by taking the sum of the absolute values of
%    the differences between the Gauss-Legendre and Kronrod rules
%    over each subinterval.  This is usually a good estimate of
%    the error in the value RESULTG.  The error in the Kronrod
%    estimate RESULTK is usually much smaller.
%
  resultg = 0.0;
  resultk = 0.0;
  error_est = 0.0;

  if ( a == b )
    return
  end

  if ( nk ~= 2 * ng + 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SUM_SUB_GK - Fatal error!\n' );
    fprintf ( 1, '  NK must equal 2 * NG + 1.\n' );
    fprintf ( 1, '  The input value was NG = %d\n', ng );
    fprintf ( 1, '  The input value was NK = %d\n', nk );
    error (  'SUM_SUB_GK - Fatal error!' );
  end

  h = ( b - a ) / nsub;

  for j = 1 : nsub

    xmid = a + 0.5 * h * ( 2 * j - 1 );

    partg = 0.0;
    partk = 0.0;

    for i = 1 : nk

      xval = xmid + 0.5 * h * xk(i);
      fk = func ( xval );
      partk = partk + 0.5 * h * wk(i) * fk;

      if ( mod ( i, 2 ) == 0 )
        partg = partg + 0.5 * h * wg(i/2) * fk;
      end

    end

    resultg = resultg + partg;
    resultk = resultk + partk;
    error_est = error_est + abs ( partk - partg );

  end

  return
end


text:
function [ m, d ] = easter_ds ( y )

%*****************************************************************************80
%
%% EASTER_DS computes the month and day of Easter for a Gregorian year.
%
%  Example:
%
%    Input:
%
%      Y = 2000
%
%    Output:
%
%      M = 4
%      D = 23
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Peter Duffett-Smith,
%    Practical Astronomy With Your Calculator,
%    Third Edition,
%    Cambridge University Press, 1996,
%    ISBN: 0-521-35699-7,
%    LC: QB62.5.D83.
%
%  Parameters:
%
%    Input, integer Y, the year, which must be 1583 or greater.
%    (The formula is only valid for years after the Gregorian calendar
%    was adopted.)
%
%    Output, integer M, D, the month and day of Easter.
%
  if ( y <= 0 )
    m = -1;
    d = -1;
    return
  end

  a = year_to_golden_number ( y );

  a = a - 1;

  b = floor ( y / 100 );
  c = mod ( y, 100 );

  dd = floor ( b / 4 );
  e = mod ( b, 4 );

  f = floor ( ( b + 8 ) / 25 );
  g = floor ( ( b - f + 1 ) / 3 );
  h = mod ( 19 * a + b - dd - g + 15, 30 );

  i = floor ( c / 4 );
  k = mod ( c, 4 );

  l = mod ( 32 + 2 * e + 2 * i - h - k, 7 );
  mm = floor ( ( a + 11 * h + 22 * l ) / 451 );

  m = floor ( ( h + l - 7 * mm + 114 ) / 31 );
  d = mod ( h + l - 7 * mm + 114, 31 ) + 1;

  return
end


text:
function path = grFindGeodesicPath(nodes, edges, ind0, ind1, edgeWeights)
%GRFINDGEODESICPATH Find a geodesic path between two nodes in the graph
%
%   PATH = grFindGeodesicPath(NODES, EDGES, NODE1, NODE2, WEIGHTS)
%   NODES and EDGES defines the graph, NODE1 and NODE2 are indices of the
%   node extremities, and WEIGHTS is the set of weights associated to each
%   edge.
%   The function returns a set of edge indices.
%
%
%   See also
%   grFindMaximalLengthPath
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-05-22,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% ensure weights are defined
if ~exist('edgeWeights', 'var')
    edgeWeights = ones(size(edges, 1), 1);
end

% check indices limits
nNodes = size(nodes, 1);
if max(ind0) > nNodes
    error('Start index exceed number of nodes in the graph');
end
if max(ind1) > nNodes
    error('End index exceed number of nodes in the graph');
end

% find a vertex opposite to the first extremity
dists = grPropagateDistance(nodes, edges, ind0, edgeWeights);

% iterate on neighbors of current node: choose next neighbor with smallest
% cumulated weight, until we are back on source node
path = [];
while true
    % find neighbor with lowest cumulated distance
    neighs = grAdjacentNodes(edges, ind1);
    neighDists = dists(neighs);
    indN = find(neighDists == min(neighDists), 1);
    ind2 = neighs(indN);

    if isempty(ind2)
        warning('graphs:grFindGeodesicPath', ...
            'No neighbor node found for node %d, graph may be not connected', ind1);
        break;
    end

    % add edge index to the path
    indE = find(sum(ismember(edges, [ind1 ind2]), 2) == 2, 1);
    path = [path indE]; %#ok<AGROW>

    % test if path is finished or not
    if ind2 == ind0
        break;
    end
    ind1 = ind2;
end

% reverse path direction
path = path(end:-1:1);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function[answer] = isfieldRec(struct, varargin)
% [answer] = isfieldRec(struct, varargin)
%
% Check if a struct has all the fields and subfields specified by their name strings.
% I.e. to see if struct a has child b with element c, call isfieldRec(a, 'b', 'c')
%
% Copyright by Holger Caesar, 2014

answer = true;
nargs = numel(varargin);

for i = 1 : nargs,
    fieldName = varargin{i};
    if ~isfield(struct, fieldName),
        answer = false;
        return;
    elseif i < nargs,
        struct = getfield(struct, fieldName);
    end;
end;

author: nightrome

repo: matconvnet-calvin