text:
function varargout = linearProgramming(maximizeFlag, f, A, b, lb, ub, constraintTypes, variableTypes, options)
% Provides a common interface to several linear programming solvers.
% Checks each solver for errors.
%
% That is it solves the problem:
%   max f'*x
%   subject to A*x~b, lb<=x<=ub
% where "~" can be equality, or double or single-sided inequality
% constraints if the chosen linear programming solver supports those
% types of constraints. Additionally the variables, x, can, if the chosen
% linear programming solver supports it, either continuous, integer, and
% boolean valued.
%
% Constraint Types: etiher a single value, or a column vector or values:
%   'F' Free (unbounded) variable (the constraint is ignored).
%   'U' Variable with upper bound ( A(i,:)*x <= b(i)).
%   'S' Fixed Variable (A(i,:)*x = b(i)).
%   'L' Variable with lower bound (A(i,:)*x >= b(i)).
%   'D' Double-bounded variable (A(i,:)*x >= -b(i) and A(i,:)*x <= b(i)).
%
% Variable Types: either a single value, or a column vector or values:
%   'C' Continuous variable.
%   'I' Integer variable
%   'B' Binary variable
%
% Output
%  x
%  lambda,
%  fopt
%  errorFlag
%  errorMsg
%  extra
%
% Requirements: at least one of the following linear programing solvers
% - glpk -- GNU Linear programming kit, to solve linear programming problems
%   http://sourceforge.net/projects/glpkmex/
% - lp_solve -- open source solve linear programming solver
%   http://sourceforge.net/projects/lpsolve
% - clp -- open source solve linear programming solver
%   http://control.ee.ethz.ch/~joloef/mexclp.zip
% - bpmpd -- open source solve linear programming solver
%   http://www.pserc.cornell.edu/bpmpd/
% - qsopt -- open source solve linear programming solver
%   http://control.ee.ethz.ch/~joloef/mexqsopt.msql
% - Optimization Toolbox - has linprog, MATLAB's linear programming
%   solver. Note: linprog is not vevy good. We suggest you use another
%   solver instead.
%   http://www.mathworks.com/products/optimization/
%
% Author: Jonathan Karr
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 10/31/2008

%- A is non-empty
%- b is real column vector of length ncols(A)
%- lb, ub are real column vectors of length nrows(A)
validateattributes(A, {'numeric'}, {'nonempty'});
validateattributes(b, {'numeric'}, {'real', 'size', [size(A, 1) 1]});
if isempty(lb)
    lb = -Inf(size(A, 2), 1);
else
    validateattributes(lb, {'numeric'}, {'real', 'size', [size(A, 2) 1]});
end
if isempty(ub)
    ub = -Inf(size(A, 2), 1);
else
    validateattributes(ub, {'numeric'}, {'real', 'size', [size(A, 2) 1]});
end

%error check length sizes of contraint and variable types, expand if
%necessary
if isscalar(constraintTypes)
    if ~ismembc(constraintTypes, 'DFLSUdflsu')
        throw(MException('ComputationUtil:LinearProgramming', 'linear programming: invalid constraints types'));
    end
    constraintTypes = constraintTypes(ones(length(b), 1), 1);
elseif length(constraintTypes) ~= length(b)
    throw(MException('ComputationUtil:LinearProgramming', 'linear programming: invalid constraints types'));
elseif ~all(ismembc(constraintTypes, 'DFLSUdflsu'))
    throw(MException('ComputationUtil:LinearProgramming', 'linear programming: invalid constraints types'));
end

if isscalar(variableTypes)
    if ~ismembc(variableTypes, 'BCIbci')
        throw(MException('ComputationUtil:LinearProgramming', 'linear programming: invalid variable types'));
    end
    variableTypes = variableTypes(ones(size(A, 2), 1), 1);
elseif length(variableTypes) ~= length(lb)
    throw(MException('ComputationUtil:LinearProgramming', 'linear programming: invalid variable types'));
elseif ~all(ismembc(variableTypes, 'BCIbci'))
    throw(MException('ComputationUtil:LinearProgramming', 'linear programming: invalid variable types'));
end

%compute sense from maximizeFlag
%if maximize='maximize'->sense=-1
%otherwise->sense=1
sense = 1 - 2 *strcmp(maximizeFlag, 'maximize');

%solver linear programming problem using chosen solver
varargout = cell(max(nargout, 1), 1);
switch options.solver(1:2)
    case 'gl'; [varargout{:}] = runglpk(   f, A, b, lb, ub, sense, constraintTypes, variableTypes, options);
    case 'li'; [varargout{:}] = runlinprog(f, A, b, lb, ub, sense, constraintTypes, variableTypes, options);
    case 'bp'; [varargout{:}] = runbpmpd(  f, A, b, lb, ub, sense, constraintTypes, variableTypes, options);
    case 'cl'; [varargout{:}] = runclp(    f, A, b, lb, ub, sense, constraintTypes, variableTypes, options);
    case 'lp'; [varargout{:}] = runlpsolve(f, A, b, lb, ub, sense, constraintTypes, variableTypes, options);
    case 'qs'; [varargout{:}] = runqsopt(  f, A, b, lb, ub, sense, constraintTypes, variableTypes, options);
    case 'to'; [varargout{:}] = runtomlab( f, A, b, lb, ub, sense, constraintTypes, variableTypes, options);
    case 'cp'; [varargout{:}] = runcplex(  f, A, b, lb, ub, sense, constraintTypes, variableTypes, options);
    case 'gu'; [varargout{:}] = rungurobi( f, A, b, lb, ub, sense, constraintTypes, variableTypes, options);
    case 'mo'; [varargout{:}] = runmosek(  f, A, b, lb, ub, sense, constraintTypes, variableTypes, options);
    otherwise; throw(MException('ComputationUtil:LinearProgramming', 'Support for %s not implemented', options.solver));
end

function [x, lambda, fopt, errorFlag, errorMsg, extra] = runlinprog(f, A, b, lb, ub, sense, constraintTypes, variableTypes, options)

%setup equality, inequality constraints based on constraint types
Aeq = [];
beq = [];
if sum(constraintTypes == 'S') == length(constraintTypes)
    Aeq = A;
    A = [];
    beq = b;
    b = [];
elseif sum(constraintTypes == 'U') == length(constraintTypes)
elseif sum(constraintTypes == 'L') == length(constraintTypes)
    A = -A;
    b = -b;
else
    throw(MException('ComputationUtil:LinearProgramming', 'linprog: constraint types invalid'));
end

%throw error if integer or boolean variable types requested, linprog
%doesn't have this capability
if sum(variableTypes == 'C') ~= length(variableTypes)
    throw(MException('ComputationUtil:LinearProgramming', 'linprog: variable types invalid'));
end

%linprog options
if isfield(options.solverOptions, 'linprog')
    linprogoptions = options.solverOptions.linprog;
else
    linprogoptions = struct;
end

try
    [x, fopt, exitflag, ~, lambda] = linprog(sense * f, A, b, Aeq, beq, lb, ub, [], linprogoptions);

    lambda = struct('reducedCosts', NaN(size(A, 2), 1), 'shadowPrices', lambda);
    fopt = fopt * sense;
    errorFlag = exitflag ~= 1;
    switch exitflag
        case  1; errorMsg = 'Function converged to a solution x.';
        case  0; errorMsg = 'Number of iterations exceeded options.MaxIter.';
        case -2; errorMsg = 'No feasible point was found.';
        case -3; errorMsg = 'Problem is unbounded.';
        case -4; errorMsg = 'NaN value was encountered during execution of the algorithm.';
        case -5; errorMsg = 'Both primal and dual problems are infeasible.';
        case -7; errorMsg = 'Search direction became too small. No further progress could be made.';
    end
catch exception
    errorFlag = 1;
    errorMsg = exception.message;
    x = NaN(size(lb));
    lambda = struct('reducedCosts', NaN(size(A, 2), 1), 'shadowPrices', NaN(size(A, 1), 1));
    fopt = NaN;
end

extra = struct;

function [x, lambda, fopt, errorFlag, errorMsg, extra] = runglpk(f, A, b, lb, ub, sense, constraintTypes, variableTypes, options)

%make structure of glpk options
if isfield(options.solverOptions, 'glpk')
    glpkoptions = options.solverOptions.glpk;
else
    glpkoptions = struct;
end

%call glpk
try
    [x, fopt, status, extra] = glpkcc(f, A, b, lb, ub, constraintTypes, variableTypes, sense, glpkoptions);
    lambda = struct('reducedCosts', extra.redcosts, 'shadowPrices', extra.lambda);
    errorFlag = ~(status == 5 || status == 2);
    switch status
        %General Errors
        case 5; errorMsg = 'Solution is optimal';
        case 2; errorMsg = 'Solution is feasible';
        case 1; errorMsg = 'Solution is undefined';
        case 3; errorMsg = 'Solution is infeasible';
        case 4; errorMsg = 'No feasible solution exists';
        case 6; errorMsg = 'Solution is unbounded';

        %Simplex method Errors:
        case 101; errorMsg = 'Invalid basis';
        case 102; errorMsg = 'Singular matrix';
        case 103; errorMsg = 'Ill-conditioned matrix';
        case 104; errorMsg = 'Invalid bounds';
        case 105; errorMsg = 'Solver failed';
        case 106; errorMsg = 'Objective lower limit reached';
        case 107; errorMsg = 'Objective upper limit reached';
        case 108; errorMsg = 'Iteration limit exceeded';
        case 109; errorMsg = 'Time limit exceeded';
        case 110; errorMsg = 'No primal feasible solution';

        %Interior point method, mixed integer problem Errors:
        case 204; errorMsg = 'Unable to start the search.';
        case 205; errorMsg = 'Objective function lower limit reached.';
        case 206; errorMsg = 'Objective function upper limit reached.';
        case 207; errorMsg = 'Iterations limit exhausted.';
        case 208; errorMsg = 'Time limit exhausted.';
        case 209; errorMsg = 'No feasible solution.';
        case 210; errorMsg = 'Numerical instability.';
        case 211; errorMsg = 'Problems with basis matrix.';
        case 212; errorMsg = 'No convergence (interior).';
        case 213; errorMsg = 'No primal feasible solution (LP presolver).';
        case 214; errorMsg = 'No dual feasible solution (LP presolver).';

        otherwise, errorMsg = sprintf('Invalid error code %d', status);
    end
catch exception
    errorFlag = 1;
    errorMsg = exception.message;
    x = NaN(size(lb));
    lambda = struct('reducedCosts', NaN(size(A, 2), 1), 'shadowPrices', NaN(size(A, 1), 1));
    fopt = NaN;
end

extra = struct;

function [x, lambda, fopt, errorFlag, errorMsg, extra] = runlpsolve(f, A, b, lb, ub, sense, constraintTypes, variableTypes, options)
%options
verbose = 3;
presolve = 0; %none
scaling = 4 + 64 + 128; %geometric + equilibrate + integers
if isfield(options.solverOptions, 'lp_solve')
    if isfield(options.solverOptions.lp_solve, 'verbose'), verbose = options.solverOptions.lp_solve.verbose; end
    if isfield(options.solverOptions.lp_solve, 'presolve'), presolve = options.solverOptions.lp_solve.presolve; end
    if isfield(options.solverOptions.lp_solve, 'scaling'), scaling = options.solverOptions.lp_solve.scaling; end
end

%convert variable types to lp_solve format'
if any(variableTypes == 'B'); throw(MException('ComputationUtil:LinearProgramming', 'lp_solve: boolean variables are not allowed')); end;
xint = NaN(size(variableTypes));
xint(variableTypes == 'I') = 1;
xint(variableTypes == 'C') = 0;

%convert constraints types to lp_solve format
if any(constraintTypes == 'F'); throw(MException('ComputationUtil:LinearProgramming', 'lp_solve: free constraints are not allowed')); end;
if any(constraintTypes == 'D'); throw(MException('ComputationUtil:LinearProgramming', 'lp_solve: double-sided constraints are not allowed')); end;
con_types = NaN(size(constraintTypes));
con_types(constraintTypes == 'L') = 2;
con_types(constraintTypes == 'S') = 3;
con_types(constraintTypes == 'U') = 1;

try
    lp = mxlpsolve('make_lp', size(A, 1), size(A, 2));
    mxlpsolve('set_verbose', lp, verbose);
    mxlpsolve('set_mat', lp, A);
    mxlpsolve('set_rh_vec', lp, b);
    mxlpsolve('set_obj_fn', lp, f);
    mxlpsolve('set_sense', lp, (sense == -1) + 0);
    mxlpsolve('set_constr_type', lp, con_types);
    mxlpsolve('set_bounds', lp, lb, ub);
    mxlpsolve('set_int', lp, xint);
    mxlpsolve('set_presolve', lp, presolve, mxlpsolve('get_presolveloops', lp));
    mxlpsolve('set_scaling', lp, scaling);

    status = mxlpsolve('solve', lp);
    errorFlag =~ status == 0 || status == 1 || status == 11 || status == 12;
    if ~errorFlag
        [fopt, x, duals] = mxlpsolve('get_solution', lp);
        reducedCosts = mxlpsolve('get_reduced_costs', lp);
    else
        fopt = NaN;
        x = NaN(size(lb));
        duals = NaN(size(A, 1), 1);
        reducedCosts = NaN(size(A, 2), 1);
    end

    extra = struct;
    if nargin >= 6
        tmpFileName = ['lpsolve.' datestr(now, 30) '.txt'];

        mxlpsolve('set_outputfile', lp, tmpFileName);
        mxlpsolve('print_scales', lp);

        fid = fopen(tmpFileName, 'r');
        fgetl(fid);
        fgetl(fid);
        line = fgetl(fid);
        extra.objScaling = str2double(line(32:end));
        extra.rowScaling = zeros(size(A, 1), 1);
        extra.colScaling = zeros(size(A, 2), 1);
        for i = 1:size(A, 1)
            line = fgetl(fid);
            extra.rowScaling(i) = str2double(line(32:end));
        end
        for i = 1:size(A, 2)
            line = fgetl(fid);
            extra.colScaling(i) = str2double(line(32:end));
        end
        fclose(fid);
    end

    mxlpsolve('delete_lp', lp);

    if nargin >= 6
        delete(tmpFileName);
    end

    lambda = struct('reducedCosts', reducedCosts, 'shadowPrices', duals);
    switch status
        case -2; errorMsg = 'Out of memory';
        case 0; errorMsg = 'An optimal solution was obtained';
        case 1; errorMsg = ['The model is sub-optimal. Only happens if there are integer variables and there is already an integer solution found. The solution is not guaranteed the most optimal one.\n'...
                '* A timeout occured (set via set_timeout or with the -timeout option in lp_solve)\n'...
                '* set_break_at_first was called so that the first found integer solution is found (-f option in lp_solve)\n'...
                '* set_break_at_value was called so that when integer solution is found that is better than the specified value that it stops (-o option in lp_solve)\n'...
                '* set_mip_gap was called (-g/-ga/-gr options in lp_solve) to specify a MIP gap\n'...
                '* An abort function is installed (put_abortfunc) and this function returned TRUE\n'...
                '* At some point not enough memory could not be allocated'];
        case 2; errorMsg = 'The model is infeasible';
        case 3; errorMsg = 'The model is unbounded';
        case 4; errorMsg = 'The model is degenerative';
        case 5; errorMsg = 'Numerical failure encountered';
        case 6; errorMsg = 'The abort routine returned TRUE. See put_abortfunc';
        case 7; errorMsg = 'A timeout occurred. A timeout was set via set_timeout';
        case 9; errorMsg = 'The model could be solved by presolve. This can only happen if presolve is active via set_presolve';
        case 10; errorMsg = 'The B&B routine failed';
        case 11; errorMsg = 'The B&B was stopped because of a break-at-first (see set_break_at_first) or a break-at-value (see set_break_at_value)';
        case 12; errorMsg = 'A feasible B&B solution was found';
        case 13; errorMsg = 'No feasible B&B solution found';
    end
catch exception
    errorFlag = 1;
    errorMsg = exception.message;
    x = NaN(size(lb));
    lambda = struct('reducedCosts', NaN(size(A, 2), 1), 'shadowPrices', NaN(size(A, 1), 1));
    fopt = NaN;
    extra = struct;
end

function [x, lambda, fopt, errorFlag, errorMsg, extra] = runqsopt(f, A, b, lb, ub, sense, constraintTypes, variableTypes, options)

%setup equality, inequality constraints based on constraint types
Aeq = [];
beq = [];
if sum(constraintTypes == 'S') == length(constraintTypes)
    Aeq = A;
    A = [];
    beq = b;
    b = [];
elseif sum(constraintTypes == 'U') == length(constraintTypes)
elseif sum(constraintTypes == 'L') == length(constraintTypes)
    A = -A;
    b = -b;
else
    throw(MException('ComputationUtil:LinearProgramming', 'qsopt: constraint types invalid'));
end

%throw error if integer or boolean variable types requested, qsopt
%doesn't have this capability
if sum(variableTypes == 'C') ~= length(variableTypes)
    throw(MException('ComputationUtil:LinearProgramming', 'qsopt: variable types invalid'));
end

%make structure of qsopt options
qsoptoptions = struct;
if isfield(options.solverOptions, 'qsopt')
    qsoptoptions = options.solverOptions.qsopt;
end

try
    [x, dual, status] = qsopt(sense * f, A, b, Aeq, beq, lb, ub, qsoptoptions);
    lambda = struct('reducedCosts', NaN(size(A, 2), 1), 'shadowPrices', dual);
    fopt = f' * x;
    errorFlag = status ~= 1;
    switch status
        case 1; errorMsg = 'optimal';
        case 2; errorMsg = 'infeasible';
        case 3; errorMsg = 'unbounded';
        case 4; errorMsg = 'iteration limit';
        case 5; errorMsg = 'time limit';
        case 6; errorMsg = 'other problem';
    end
catch exception
    errorFlag = 1;
    errorMsg = exception.message;
    x = NaN(size(lb));
    lambda = struct('reducedCosts', NaN(size(A, 2), 1), 'shadowPrices', NaN(size(A, 1), 1));
    fopt = NaN;
end

extra = struct;

function [x, lambda, fopt, errorFlag, errorMsg, extra] = runclp(f, A, b, lb, ub, sense, constraintTypes, variableTypes, options)

%setup equality, inequality constraints based on constraint types
Aeq = [];
beq = [];
if sum(constraintTypes == 'S') == length(constraintTypes)
    Aeq = A;
    A = [];
    beq = b;
    b = [];
elseif sum(constraintTypes == 'U') == length(constraintTypes)
elseif sum(constraintTypes == 'L') == length(constraintTypes)
    A = -A;
    b = -b;
else
    throw(MException('ComputationUtil:LinearProgramming', 'clp: constraint types invalid'));
end

%throw error if integer or boolean variable types requested, clp
%doesn't have this capability
if sum(variableTypes == 'C') ~= length(variableTypes)
    throw(MException('ComputationUtil:LinearProgramming', 'clp: variable types invalid'));
end

%make structure of clp options
clpoptions = struct;
if isfield(options.solverOptions, 'clp')
    clpoptions = options.solverOptions.clp;
end

try
    [x, dual, status] = clp(zeros(length(f)), sense * f, A, b, Aeq, beq, lb, ub, clpoptions);
    lambda = struct('reducedCosts', NaN(size(A, 2), 1), 'shadowPrices', dual);
    fopt = f' * x;
    errorFlag = status ~= 0;
    switch status
        case 0; errorMsg = 'optimal';
        case 1; errorMsg = 'infeasible';
        case 2; errorMsg = 'unbounded';
    end
catch exception
    errorFlag = 1;
    errorMsg = exception.message;
    x = NaN(size(lb));
    lambda = struct('reducedCosts', NaN(size(A, 2), 1), 'shadowPrices', NaN(size(A, 1), 1));
    fopt = NaN;
end

extra = struct;

function [x, lambda, fopt, errorFlag, errorMsg, extra] = runbpmpd(f, A, b, lb, ub, sense, constraintTypes, variableTypes, ~)

%throw error if integer or boolean variable types requested, bpmpd
%doesn't have this capability
if sum(variableTypes == 'C') ~= length(variableTypes)
    throw(MException('ComputationUtil:LinearProgramming', 'bpmpd: variable types invalid'));
end

%convert constraints types to bpmpd format
if any(constraintTypes == 'F'); throw(MException('ComputationUtil:LinearProgramming', 'bpmpd: free constraints are not allowed')); end;
if any(constraintTypes == 'D'); throw(MException('ComputationUtil:LinearProgramming', 'bpmpd: double-sided constraints are not allowed')); end;
e = NaN(size(constraintTypes));
e(constraintTypes == 'L') = 1;
e(constraintTypes == 'S') = 0;
e(constraintTypes == 'U') = -1;

%throw error if integer or boolean variable types requested, bpmpd
%doesn't have this capability
if sum(variableTypes == 'C') ~= length(variableTypes)
    throw(MException('ComputationUtil:LinearProgramming', 'bpmpd: variable types invalid'));
end

try
    llist = (1:length(lb))';
    ulist = (1:length(ub))';
    [x, dual, ~, w, errorMsg] = bp(zeros(length(f)), A, b, sense * f, e, llist, lb, ulist, ub, bpopt,0);
    lambda = struct('reducedCosts', w, 'shadowPrices', dual);
    fopt = f' * x;
    errorFlag =~ strcmp(errorMsg, 'optimal solution');
catch exception
    errorFlag = 1;
    errorMsg = exception.message;
    x = NaN(size(lb));
    lambda = struct('reducedCosts', NaN(size(A, 2), 1), 'shadowPrices', NaN(size(A, 1), 1));
    fopt = NaN;
end

extra = struct;

function [x, lambda, fopt, errorFlag, errorMsg, extra] = runtomlab( f, A, b, lb, ub, sense, constraintTypes, ~, options)

b_L = -Inf(size(b));
b_U =  Inf(size(b));
b_L(constraintTypes == 'L') = b(constraintTypes == 'L');
b_U(constraintTypes == 'U') = b(constraintTypes == 'U');
b_L(constraintTypes == 'S') = b(constraintTypes == 'S');
b_U(constraintTypes == 'S') = b(constraintTypes == 'S');
b_L(constraintTypes == 'D') = -b(constraintTypes == 'D');
b_U(constraintTypes == 'D') =  b(constraintTypes == 'D');

prob = lpAssign(sense * f, A, b_L, b_U, lb, ub, [], [], [], [], [], [], [], [], []);

solver = 'minos';
prob.optParam = struct();
PriLev = 0;
if isfield(options.solverOptions, 'tomlab')
    if isfield(options.solverOptions.tomlab, 'solver'), solver = options.solverOptions.tomlab.solver; end
    if isfield(options.solverOptions.tomlab, 'optParam'), prob.optParam = options.solverOptions.tomlab.optParam; end
    if isfield(options.solverOptions.tomlab, 'PriLev'), PriLev = options.solverOptions.tomlab.PriLev; end
end

try
    result = tomRun(solver, prob, PriLev);
    errorFlag = result.ExitFlag ~= 0;
    errorMsg = result.ExitText;
    fopt = sense * result.f_k;
    x = result.x_k;
    lambda = struct('reducedCosts', NaN(size(A, 2), 1), 'shadowPrices', NaN(size(A, 1), 1));
catch exception
    errorFlag = 1;
    errorMsg = exception.message;
    x = NaN(size(lb));
    lambda = struct('reducedCosts', NaN(size(A, 2), 1), 'shadowPrices', NaN(size(A, 1), 1));
    fopt = NaN;
end

extra = struct;

function [x, lambda, fopt, errorFlag, errorMsg, extra] = runcplex(f, A, b, lb, ub, sense, constraintTypes, variableTypes, options)
constraintTypes = upper(constraintTypes);
if ~ismembc(constraintTypes, 'LSU')
    throw(MException('ComputationUtil:LinearProgramming', 'linear programming: invalid constraints types'));
end

Aineq = [
    -A(constraintTypes == 'L', :)
    A(constraintTypes == 'U', :)
    ];
bineq = [
    -b(constraintTypes == 'L', :)
    b(constraintTypes == 'U', :)
    ];

Aeq = A(constraintTypes == 'S', :);
beq = b(constraintTypes == 'S', :);

try
    [x, fopt, exitFlag, output] =  cplexmilp(...
        sense * f, ...
        Aineq, bineq, ...
        Aeq, beq, ...
        [], [], [], lb, ub, variableTypes');
    errorFlag = exitFlag ~= 1;
    errorMsg = output.message;
catch exception
    errorFlag = 1;
    errorMsg = exception.message;
    x = NaN(size(lb));
    fopt = NaN;
end
lambda = struct('reducedCosts', NaN(size(A, 2), 1), 'shadowPrices', NaN(size(A, 1), 1));

extra = struct('status', output.cplexstatus);

function [x, lambda, fopt, errorFlag, errorMsg, extra] = rungurobi( f, A, b, lb, ub, sense, constraintTypes, variableTypes, options)
b_L = -Inf(size(b));
b_U =  Inf(size(b));
b_L(constraintTypes == 'L') = b(constraintTypes == 'L');
b_U(constraintTypes == 'U') = b(constraintTypes == 'U');
b_L(constraintTypes == 'S') = b(constraintTypes == 'S');
b_U(constraintTypes == 'S') = b(constraintTypes == 'S');
b_L(constraintTypes == 'D') = -b(constraintTypes == 'D');
b_U(constraintTypes == 'D') =  b(constraintTypes == 'D');

xint = NaN(size(variableTypes));
xint(variableTypes == 'I') = 1;
xint(variableTypes == 'C') = 0;

PriLev = 0;
grbControl = struct();
if isfield(options.solverOptions, 'gurobi')
    if isfield(options.solverOptions.gurobi, 'PriLev'), PriLev = options.solverOptions.gurobi.PriLev; end
    if isfield(options.solverOptions.gurobi, 'grbControl'), grbControl = options.solverOptions.gurobi.grbControl; end
end

try
    [x, ~, v, rc, fopt, ~, ~, status] = ...
        gurobi(sense * f, A, lb, ub, b_L, b_U, [], ...
        grbControl, PriLev, xint);
    fopt = fopt * sense;
    [errorMsg, exitFlag] = grbStatus(status);
    errorFlag = exitFlag ~= 0;
    lambda = struct('reducedCosts', rc, 'shadowPrices', v);
catch exception
    errorFlag = 1;
    errorMsg = exception.message;
    x = NaN(size(lb));
    lambda = struct('reducedCosts', NaN(size(A, 2), 1), 'shadowPrices', NaN(size(A, 1), 1));
    fopt = NaN;
end

extra = struct;

function [x, lambda, fopt, errorFlag, errorMsg, extra] = runmosek( f, A, b, lb, ub, sense, constraintTypes, ~, options)
if sense == 1
    cmd = 'minimize';
else
    cmd = 'maximize';
end

b_L = -Inf(size(b));
b_U =  Inf(size(b));
b_L(constraintTypes == 'L') = b(constraintTypes == 'L');
b_U(constraintTypes == 'U') = b(constraintTypes == 'U');
b_L(constraintTypes == 'S') = b(constraintTypes == 'S');
b_U(constraintTypes == 'S') = b(constraintTypes == 'S');
b_L(constraintTypes == 'D') = -b(constraintTypes == 'D');
b_U(constraintTypes == 'D') =  b(constraintTypes == 'D');

lb = max(lb, -1e7);
ub = min(ub,  1e7);

opts = struct(...
    'MSK_IPAR_LOG', 0, 'MSK_IPAR_MAX_NUM_WARNINGS', 1e5, ...
    'MSK_IPAR_INTPNT_SCALING', 0, 'MSK_IPAR_SIM_SCALING', 0, 'MSK_IPAR_SIM_SCALING_METHOD', 1);
if isfield(options.solverOptions, 'mosek')
    opts = options.solverOptions.mosek;
end

try
    result = msklpopt(f, A, b_L, b_U, lb, ub, opts, cmd);
    fopt = result.sol.itr.pobjval;
    x = result.sol.itr.xx;
    errorFlag = result.rcode ~= 0;
    errorMsg = result.rmsg;
    lambda = struct(...
        'reducedCosts', result.sol.itr.slx + result.sol.itr.sux, ...
        'shadowPrices', result.sol.itr.slc + result.sol.itr.suc);
catch exception
    errorFlag = 1;
    errorMsg = exception.message;
    x = NaN(size(lb));
    lambda = struct('reducedCosts', NaN(size(A, 2), 1), 'shadowPrices', NaN(size(A, 1), 1));
    fopt = NaN;
end

extra = struct;

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function varargout = freq_response_localwfs_vss(X,xs,src,conf)
%FREQ_RESPONSE_LOCALWFS_VSS frequency response for local WFS
%
%   Usage: [S,f] = freq_response_localwfs_vss(X,xs,src,conf)
%
%   Input parameters:
%       X           - listener position / m
%       xs          - position of virtual source / m
%       src         - source type of the virtual source
%                         'pw' -plane wave
%                         'ps' - point source
%                         'fs' - focused source
%       conf        - configuration struct (see SFS_config)
%
%   Output parameters:
%       S           - simulated frequency response
%       f           - corresponding frequency axis / Hz
%
%   FREQ_RESPONSE_LOCALWFS_VSS(X,xs,src,conf) simulates the frequency response
%   of a source synthesized by local WFS using virtual secondary sources (VSS)
%   at the given virtual microphone position X.
%   The length in samples of the frequency response is given by conf.N. The
%   actual calculation is done via sound_field_mono() and a loop over frequency
%   f.
%
%   See also: sound_field_mono_localwfs_vss, time_response_localwfs_vss

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 4;
nargmax = 4;
narginchk(nargmin,nargmax);
isargposition(X);
isargxs(xs);
isargchar(src);
isargstruct(conf);

%% ===== Configuration ==================================================
N = conf.N;
showprogress = conf.showprogress;
useplot = conf.plot.useplot;
% Check type of secondary sources to use
if strcmp('2D',conf.dimension)
    greens_function = 'ls';
else
    greens_function = 'ps';
end

%% ===== Computation ====================================================
% Disable progress bar and plotting for sound_field_imp()
conf.showprogress = false;
conf.plot.useplot = false;
% Get the position of the loudspeakers
x0_real = secondary_source_positions(conf);
% Generate frequencies (10^1-10^4.3)
f = logspace(0,4.3,N)';
S = zeros(size(f));
% Get the result for all frequencies
for ii = 1:length(f)
    if showprogress, progress_bar(ii,length(f)); end
    [D, x0] = driving_function_mono_localwfs_vss(x0_real,xs,src,f(ii),conf);
    % Calculate sound field at the listener position
    P = sound_field_mono(X(1),X(2),X(3),x0,greens_function,D,f(ii),conf);
    S(ii) = abs(P);
end

% Return parameter
if nargout>0, varargout{1}=S; end
if nargout>1, varargout{2}=f; end

%% ===== Plotting ========================================================
if nargout==0 || useplot
    figure;
    figsize(conf.plot.size(1),conf.plot.size(2),conf.plot.size_unit);
    semilogx(f,db(S));
    set(gca,'XTick',[10 100 250 1000 5000 20000]);
    ylabel('amplitude / dB');
    xlabel('frequency / Hz');
end

author: sfstoolbox

repo: sfs-matlab


text:
function spm_save(f,var,varargin)
% Save text and numeric data to file
% FORMAT spm_save(f,var,opts,...)
% f     - filename (can be gzipped) {csv,tsv,json,txt,mat,npy}
% var   - data array or structure
% opts  - optional inputs to be passed on to lower level function
%__________________________________________________________________________
% Copyright (C) 2018-2019 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_save.m 7668 2019-09-27 10:44:45Z guillaume $

ext = lower(spm_file(f,'ext'));
if isempty(ext), ext = ''; end
switch ext

    case 'gz'
        fgz = spm_file(f,'basename');
        spm_save(fgz,var,varargin{:});
        gzip(fgz);
        delete(fgz);

    case {'csv','tsv'}
        if strcmp(ext,'csv')
            delim = ',';
        else
            delim = sprintf('\t');
        end
        if isstruct(var) || iscell(var) || isnumeric(var) || islogical(var)
            if isstruct(var)
                fn = fieldnames(var);
                var = struct2cell(var)';
                for i=1:numel(var)
                    var{i} = var{i}(:);
                    if ~iscell(var{i})
                        var{i} = cellstr(num2str(var{i},16));
                        var{i}(cellfun(@(x) strcmp(x,'NaN'),var{i})) = {'n/a'};
                    end
                end
                var = [fn'; var{:}];
            elseif iscell(var)
                var = cellfun(@(x) num2str(x,16), var, 'UniformOutput',false);
            elseif isnumeric(var) || islogical(var)
                var = num2cell(var);
                var = cellfun(@(x) num2str(x,16), var, 'UniformOutput',false);
            end
            try, var = strtrim(var); end

            fid = fopen(f,'Wt');
            if fid == -1
                error('Unble to write file %s.', f);
            end
            for i=1:size(var,1)
                for j=1:size(var,2)
                    if isempty(var{i,j})
                        var{i,j} = 'n/a';
                    elseif any(var{i,j} == delim)
                        var{i,j} = ['"' var{i,j} '"'];
                    end
                    fprintf(fid,'%s',var{i,j});
                    if j < size(var,2), fprintf(fid,delim); end
                end
                fprintf(fid,'\n');
            end
            fclose(fid);
        elseif isa(var,'table')
            writetable(var,f,'FileType','text','Delimiter',delim);
        else
            error('Unknown data type.');
        end

    case 'json'
        if nargin < 3, varargin = {struct([])}; end
        spm_jsonwrite(f,var,varargin{:});

    case {'','txt','md'}
        var = cellstr(var);
        fid = fopen(f,'Wt');
        if fid == -1
            error('Unable to write file %s.', f);
        end
        for i=1:numel(var)
            fprintf(fid,'%s\n',var{i});
        end
        fclose(fid);

    case 'mat'
        if nargin < 3, varargin = {spm_get_defaults('mat.format')}; end
        if isstruct(var)
            save(f,'-struct','var',varargin{:});
        else
            error('Unable to write file %s.', f);
        end

    case 'npy'
        fid = fopen(f,'W');
        fwrite(fid,[147 'NUMPY'],'uint8');
        fwrite(fid,[2 0],'uint8');
        dt  = containers.Map(...
            {'logical','uint8','uint16','uint32','uint64','int8','int16','int32','int64','single','double'},...
            {'b1','u1','u2','u4','u8','i1','i2','i4','i8','f4','f8'});
        hdr = ['{''descr'': ''<' dt(class(var)) ''', ' ...
            '''fortran_order'': True, ' ...
            '''shape'': (' sprintf('%d, ',size(var)) '), }' sprintf('\n')];
        len = 6+2+4+numel(hdr);
        hdr = [hdr blanks(ceil(len/8)*8 - len)];
        fwrite(fid,numel(hdr),'uint32');
        fwrite(fid,hdr,'uint8');
        fwrite(fid,var,class(var));
        fclose(fid);

    otherwise
        error('Unknown file format.');

end

author: spm


text:
function nvar = p19_nvar ( option )

%*****************************************************************************80
%
%% P19_NVAR sets the number of variables for problem 19.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Output, integer NVAR, the number of variables.
%
  nvar = 13;

  return
end

author: johannesgerer

repo: jburkardt-m


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% function [yc,His] = GaussNewton(fctn,yc,varargin)
%
% Gauss-Newton scheme with variable line search for minimizing J = fctn(yc)
%
% Input:
%   fctn        function handle
%   yc          starting guess
%   varargin    optional parameter, see below
%
% Output:
%   yc          numerical optimizer (current iterate)
%   his         iteration history
%
% Note: the linear systems are solved using solveLinearSystem.m
%==============================================================================

function [yc,His] = GaussNewton(fctn,yc,varargin)

if nargin==0
    help(mfilename)
    E9_Hands_NPIR
    yc = 'endOfMinimalExample';
    return;
end;

% parameter initialization -----------------------------------------------
maxIter      = 10;              % maximum number of iterations
tolJ         = 1e-3;            % for stopping, objective function
tolY         = 1e-2;            %   - " -     , current value
tolG         = 1e-2;            %   - " -     , norm of gradient
lineSearch   = @Armijo;         % linesearch scheme
LSmaxIter    = 10;              % maximum number of line search iterations
LSreduction  = 1e-4;            % minimal reduction in line search
vecNorm      = @norm;           % norm to be used for dJ and dy
solver       = [];              % linear solver
yStop        = [];              % used for stopping in multi-level framework
Jstop        = [];              %
Plots        = @(iter,para) []; % for plots;

for k=1:2:length(varargin),     % overwrites default parameter
  eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

if ~isa(Plots,'function_handle') && (Plots == 0 || strcmp(Plots,'off')),
  Plots        = @(iter,para) []; % for plots;
end;

if isempty(yStop), yStop  = yc; end; % yStop used for stopping only
% -- end parameter setup   ----------------------------------------------

% some output
if isstring(solver) || ischar(solver)
    solverName = solver;
elseif isa(solver,'function_handle')
    solverName = func2str(solver);
else
    error(1);
end
FAIRmessage(sprintf('JM 2011/01/02 : %s / %s / %s',...
  mfilename,func2str(lineSearch),solverName))
fprintf('[ maxIter=%s / tolJ=%s / tolY=%s / tolG=%s / numel(yc)=%d]\n',...
  num2str(maxIter),num2str(tolJ),num2str(tolY),num2str(tolG),length(yc));

% -- initialize  ---------------------------------------------------------
STOP = zeros(5,1);

if isempty(Jstop),
  % evaluate objective function for stopping values and plots
  [Jstop,para] = fctn(yStop); Jstop = abs(Jstop) + (Jstop == 0);
  Plots('stop',para);
end;

% evaluate objective function for starting values and plots
[Jc,para,dJ,H] = fctn(yc);
Plots('start',para);
iter = 0; yOld = 0*yc; Jold = Jc; y0 = yc;

hisStr    = {'iter','J','Jold-J','|\nabla J|','|dy|','LS'};
his        = zeros(maxIter+2,6);
his(1,1:3) = [-1,Jstop,Jstop-Jc];
his(2,:)   = [0,Jc,Jstop-Jc,vecNorm(dJ),vecNorm(yc-yStop),0];

% some output
fprintf('%4s %-12s %-12s %-12s %-12s %4s\n%s\n',...
  hisStr{:},char(ones(1,64)*'-'));
dispHis = @(var) ...
  fprintf('%4d %-12.4e %-12.3e %-12.3e %-12.3e %4d\n',var);
dispHis(his(1,:));
dispHis(his(2,:));
% -- end initialization   ------------------------------------------------

%-- start the iteration --------------------------------------------------
while 1,
  % check stopping rules
  STOP(1) = (iter>0) && abs(Jold-Jc)   <= tolJ*(1+abs(Jstop));
  STOP(2) = (iter>0) && (norm(yc-yOld) <= tolY*(1+norm(yStop)));
  STOP(3) = norm(dJ)                   <= tolG*(1+abs(Jstop));
  STOP(4) = norm(dJ)                   <= 1e6*eps;
  STOP(5) = (iter >= maxIter);
  if all(STOP(1:3)) || any(STOP(4:5)), break;  end;

  iter = iter + 1;

  % solve the Gauss-Newton System
  [dy,solver] = solveLinearSystem(-dJ',H,solver);

  % check descent direction
  % note: descent is not granted if using an iterative solver
  descent =   dJ * dy;
  if descent > 0,
    warning('no descent direction, switch to -dy!')
    dy      = -dy;
  end;

  % perform Armijo line-search
  [t,yt,LSiter] = lineSearch(fctn,yc,dy,Jc,dJ,...
      'para',para,'LSmaxIter',LSmaxIter,'LSreduction',LSreduction);
  if (t == 0),
      break;
  end; % break if line-search fails

  % save old values and update
  yOld = yc; Jold = Jc; yc = yt;
  [Jc,para,dJ,H] = fctn(yc); % evalute objective function

  % some output
  his(iter+2,:) = [iter,Jc,Jold-Jc,vecNorm(dJ),vecNorm(yc-yOld),LSiter];
  dispHis(his(iter+2,:));
  para.normdY = vecNorm(yc - yOld);
  Plots(iter,para);
% pause
end;%while; % end of iteration loop
%-------------------------------------------------------------------------
Plots(iter,para);

% clean up
His.str = hisStr;
His.his = his(1:iter+2,:);
fprintf('STOPPING:\n');
fprintf('%d[ %-10s=%16.8e <= %-25s=%16.8e]\n',STOP(1),...
  '(Jold-Jc)',(Jold-Jc),'tolJ*(1+|Jstop|)',tolJ*(1+abs(Jstop)));
fprintf('%d[ %-10s=%16.8e <= %-25s=%16.8e]\n',STOP(2),...
  '|yc-yOld|',norm(yc-yOld),'tolY*(1+norm(yc)) ',tolY*(1+norm(yStop)));
fprintf('%d[ %-10s=%16.8e <= %-25s=%16.8e]\n',STOP(3),...
  '|dJ|',norm(dJ),'tolG*(1+abs(Jstop))',tolG*(1+abs(Jstop)));
fprintf('%d[ %-10s=%16.8e <= %-25s=%16.8e]\n',STOP(4),...
  'norm(dJ)',norm(dJ),'eps',1e3*eps);
fprintf('%d[ %-10s=  %-14d >= %-25s=  %-14d]\n',STOP(5),...
  'iter',iter,'maxIter',maxIter);

FAIRmessage([mfilename,' : done !'],'=');

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


text:
function Z = projectData(X, U, K)
%PROJECTDATA Computes the reduced data representation when projecting only
%on to the top k eigenvectors
%   Z = projectData(X, U, K) computes the projection of
%   the normalized inputs X into the reduced dimensional space spanned by
%   the first K columns of U. It returns the projected examples in Z.
%

% You need to return the following variables correctly.
Z = zeros(size(X, 1), K);

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the projection of the data using only the top K
%               eigenvectors in U (first K columns).
%               For the i-th example X(i,:), the projection on to the k-th
%               eigenvector is given as follows:
%                    x = X(i, :)';
%                    projection_k = x' * U(:, k);
%

subsetU = U(:, 1:K);
Z = X * subsetU;

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

end

author: vugsus

repo: coursera-machine-learning


text:
function [xSmooth,PSmooth,xUpd,PUpd]=EKalmanBatchSmoother(xInit,PInit,z,h,HJacob,f,FJacob,R,Q,kD,numIter)
%%EKALMANBATCHSMOOTHER Run the extended forward-backward Kalman smoother
%                for nonlinear dynamic and measurement models on a batch of
%                measurements. The smoothed result at one time step or
%                along the entire batch is available. The initial predicted
%                states cannot be uninformative.
%
%INPUTS: xInit The predicted state at the time of the initial measurement
%              in z.
%        PInit The covariance matrix associated with the predicted state at
%              the time of the initial measurement in z.
%            z The zDim X N matrix of measurements for the whole batch.
%            h A NX1 cell array of function handles for the measurement
%              function that transform the state into the measurement
%              domain at each step. If the same measurement function is
%              used for all steps in the batch, then h can just be the
%              single function handle used.
%       HJacob A NX1 cell array of function handles for the measurement
%              Jacobian matrix that each takes the target state as a
%              parameter. If a single measurement Jacobian matrix is used
%              for all steps of the batch, then HJacob can just be the
%              single function handle used. If an empty matrix is passed or
%              the parameter is omitted, then HJacob will be found using
%              numerical differentiation via the numDiff function with
%              default parameters.
%            f A NX1 cell array of function handles for the state
%              transition function that transform the state into the
%              measurement domain at each step. If the same measurement
%              function is used for all steps in the batch, then h can just
%              be the single function handle used.
%       FJacob A NX1 cell array of function handles for the state
%              transition Jacobian matrix that each takes the target state
%              as a parameter. If a single measurement Jacobian matrix is
%              used for all steps of the batch, then HJacob can just be the
%              single function handle used. If an empty matrix is passed or
%              the parameter is omitted, then HJacob will be found using
%              numerical differentiation via the numDiff function with
%              default parameters.
%            R The zDim X zDim X N hypermatrix of measurement covariance
%              matrices. Alternatively, if all of the measurement
%              covariance matrices are the same, one can just pass a single
%              zDim X zDim matrix.
%            Q The xDim X xDim X (N-1) hypermatrix of process noise
%              covariance matrices. Alternatively, if all of the process
%              noise covariance matrices are the same, one can just pass a
%              single xDim X xDim matrix.
%           kD The discrete time-step at which the smoothed state estimate
%              is desired, where z(:,1) is at discrete time-step 1 (not 0).
%              If kD is omitted or an empty matrix is passed, then
%              results along the entire batch are obtained.
%      numIter The number of iterations to perform if an iterated EKF is
%              desired. The default is zero. That is, just use the
%              standard update without any additional iterations.
%
%OUTPUTS: xEst The xDimXN smoothed state estimates at all steps if kD is
%              not provided or the xDimX1 smoothed state estimate at step
%              kD if kD is provided.
%         PEst The covariance matrices associated with the smoothed
%              state estimates. This is xDimXxDimXN for the whole batch
%              if kD is not provided and is xDimXxDim if kD is provided.
%         xUpd The xDimXN state estimates of the forward filter (not
%              smoothed) at all times.
%         PUpd The xDimXxDimXN covariance matrices of the forward filter
%              (not smoothed) at all times.
%
%This function implements an extended Kalman smoother modified from Chapter
%8.6 of [1].
%
%The smoothing iteration algorithm is taken from Section III of [2]. This
%provides an algorithm to iterate the smoothed state, but not the smoothed
%covariance matrix.
%
%REFERENCES:
%[1] Y. Bar-Shalom, X. R. Li, and T. Kirubarajan, Estimation with
%    Applications to Tracking and Navigation. New York: John Wiley and
%    Sons, Inc, 2001.
%[2] L. A. Johnston and V. Krishnamurthy, "Derivation of a sawtooth
%    iterated extended Kalman smoother via the AECM algorithm," IEEE
%    Transactions on Signal Processing, vol. 49, no. 9, pp. 1899-1909,
%    2001.
%
%March 2015 David Karnick, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if nargin<10
    kD=[];
end
if nargin<11
    numIter=0;
end

xDim=size(xInit,1);
zDim=size(z,1);
N=size(z,2);

if(isempty(FJacob))
    FJacob=@(x)numDiff(x,f,xDim);
end
if(isempty(HJacob))
    HJacob=@(x)numDiff(x,h,zDim);
end

if(isa(HJacob,'function_handle'))
    HJacob=repmat({HJacob},N,1);
end
if(isa(h,'function_handle'))
    h=repmat({h},N,1);
end
if(isa(FJacob,'function_handle'))
    FJacob=repmat({FJacob},N,1);
end
if(isa(f,'function_handle'))
    f=repmat({f},N,1);
end
if(size(R,3)==1)
    R=repmat(R,1,1,N);
end
if(size(Q,3)==1)
    Q=repmat(Q,1,1,N-1);
end

%Run the Kalman filter forward
xPred=zeros(xDim,N);
PPred=zeros(xDim,xDim,N);
xUpd=zeros(xDim,N);
PUpd=zeros(xDim,xDim,N);

%The first step uses the priors
xPred(:,1)=xInit;
PPred(:,:,1)=PInit;

%The rest of the steps
for curStep=1:(N-1)
    [xUpd(:,curStep),PUpd(:,:,curStep)]=EKFUpdate(xPred(:,curStep),PPred(:,:,curStep),z(:,curStep),R(:,:,curStep),h{curStep},HJacob{curStep},numIter);
    [xPred(:,curStep+1),PPred(:,:,curStep+1)]=DiscEKFPred(xUpd(:,curStep),PUpd(:,:,curStep),f{curStep},FJacob{curStep},Q(:,:,curStep));
end
[xUpd(:,end),PUpd(:,:,end)]=EKFUpdate(xPred(:,N),PPred(:,:,N),z(:,N),R(:,:,N),h{N},HJacob{N},numIter);

%Run the backwards Kalman smoother
xSmooth=zeros(xDim,N);
PSmooth=zeros(xDim,xDim,N);

xSmooth(:,N)=xUpd(:,N);
PSmooth(:,:,N)=PUpd(:,:,N);
for curStep=(N-1):-1:1
    F=FJacob{curStep}(xUpd(:,curStep));
    C=PUpd(:,:,curStep)*F'/PPred(:,:,curStep+1);
    xSmooth(:,curStep)=xUpd(:,curStep)+C*(xSmooth(:,curStep+1)-xPred(:,curStep+1));
    PSmooth(:,:,curStep)=PUpd(:,:,curStep)+C*(PSmooth(:,:,curStep+1)-PPred(:,:,curStep+1))*C';

    for curIter=1:numIter
        F=FJacob{curStep}(xSmooth(:,curStep));
        H=HJacob{curStep}(xSmooth(:,curStep));
        B=pinv(pinv(PUpd(:,:,curStep))+H'*inv(R(:,:,curStep))*H+F'*inv(Q(:,:,curStep))*F);
        xSmooth(:,curStep)=xUpd(:,curStep)+B*(...
            H'*inv(R(:,:,curStep))*(z(:,curStep)-h{curStep}(xPred(:,curStep))-H*(xPred(:,curStep)-xSmooth(:,curStep)))...
            +F'*inv(Q(:,:,curStep))*(xSmooth(:,curStep+1)-xPred(:,curStep+1)));
    end
end

if(~isempty(kD))
    xSmooth=xSmooth(:,kD);
    PSmooth=PSmooth(:,:,kD);
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:
%BITWISE_AND  Calculates the per-element bit-wise conjunction of two arrays or an array and a scalar
%
%     dst = cv.bitwise_and(src1, src2)
%     dst = cv.bitwise_and(src1, src2, 'OptionName', optionValue, ...)
%
% ## Input
% * __src1__ first input array or a scalar.
% * __src2__ second input array or a scalar. In case both are array, they must
%   have the same size and type.
%
% ## Output
% * __dst__ output array that has the same size and type as the input arrays.
%
% ## Options
% * __Mask__ optional operation mask, 8-bit single channel array, that
%   specifies elements of the output array to be changed. Not set by default.
% * __Dest__ Used to initialize the output `dst` when a mask is used. Not set
%   by default.
%
% Computes bitwise conjunction of the two arrays (`dst = src1 & src2`).
%
% The function calculates the per-element bit-wise logical conjunction for:
%
% * Two arrays when `src1` and `src2` have the same size:
%
%       dst(I) = src1(I) AND src2(I) if mask(I) != 0
%
% * An array and a scalar when `src2` is constructed from Scalar or has the
%   same number of elements as `size(src1,3)`:
%
%       dst(I) = src1(I) AND src2 if mask(I) != 0
%
% * A scalar and an array when `src1` is constructed from Scalar or has the
%   same number of elements as `size(src2,3)`:
%
%       dst(I) = src1 AND src2(I) if mask(I) != 0
%
% In case of floating-point arrays, their machine-specific bit representations
% (usually IEEE754-compliant) are used for the operation. In case of
% multi-channel arrays, each channel is processed independently. In the second
% and third cases above, the scalar is first converted to the array type.
%
% See also: cv.bitwise_or, cv.bitwise_xor, cv.bitwise_not, bitand
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% Fig. 5.47  Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
%  Figure 5.47 Root locus for the heat exchanger, with and without delay

clf
numG=1;
denG=conv([10 1],[60 1]);
sysG=tf(numG, denG);
rlocus(sysG);
axis([-.8 .8 -.6 .6]);
%pause;
hold on
gtext('locus without delay')
Td=5;
numD1 = 1;
denD1=[Td 1];
numGD1=conv(numG,numD1);
denGD1=conv(denG,denD1);
sysGD1= tf(numGD1, denGD1);
rlocus(sysGD1);

gtext('locus with (0,1) lag Pade')
title('Fig. 5.47 Root locus for heat exchanger with and without delay')
%pause;
numDe2=[Td^2/12 -Td/2 1];
denDe2=[Td^2/12 +Td/2 1];
numGD=conv(numG,numDe2);
denGD=conv(denG,denDe2);
sysGD= tf(numGD, denGD);
rlocus(sysGD);

gtext('locus with(2,2) Pade')
grid on
% title('Fig. 5.47 Root locus for heat exchanger with and without delay')
% z=0:.1:.9;
% wn= .1:.1:.6;
% sgrid(z, wn)
hold off

author: Sable

repo: mcbench-benchmarks


text:
function [Fx,Fy] = snakeForce(emap,mode,mu,niter)
%SNAKEFORCE Components of external force for use in the snake algorithm.
%   [Fx,Fy] = SNAKEFORCE(EMAP,MODE,MU,NITER) uses input edge map, EMAP,
%   to compute the force component images Fx and Fy. These images are of
%   the same size as EMAP and contain the values of Fx and Fy at all
%   points of EMAP. For example, Fx(i,j) is the x-component of the
%   force at coordinates (i,j) of EMAP. If MODE = 'MOG' (the default),
%   force components Fx and Fy of the gradient of EMAP are computed
%   (remember, the gradient is a 2D vector). This mode does not require
%   MU nor NITER. If MODE = 'gvf' then the force is based on the
%   gradient vector flow of the edge map. Option 'gvf' requires that MU
%   and the number of iterations, NITER, be provided.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% PRELIMINARIES.
% Default.
if nargin == 1
	mode = 'MOG';
end
% Work with lower case.
mode = lower(mode);

% COMPUTE THE FORCE COMPONENTS.
% Note that MATLAB function gradient works with (c,r) intead of (r,c)
% coordinates, as we do in the book.
switch mode
	case 'mog'
      % Unnormalized gradient of the edge map.
      [Fy,Fx] = gradient(emap);
	case 'gvf'
      % Gradient vector flow of the edge map.
      % Compute the magnitude of the gradient squared.
      [egy,egx] = gradient(emap);
      gradMagSq = (egx.^2 + egy.^2); % To be used later.
      % Initialize GVF to the gradient of the edge map.
      vx = egx;
      vy = egy;
      % Iterate to find vx and vy.
      for I = 1:niter
         % The 4 in the following expressions cancels the 1/4 in the
         % kernel used by the MATLAB Laplacian function del2.
         vx = vx + mu*4*del2(vx) - gradMagSq.*(vx - egx);
         vy = vy + mu*4*del2(vy) - gradMagSq.*(vy - egy);
      end
      % Set the forces equal to the components of the gradient vector
      % flow [see Eq. (12-17)].
      Fx = vx;
      Fy = vy;
end

author: dipum

repo: dipum-toolbox

save_path: github-repos/MATLAB/dipum-dipum-toolbox


text:
%% Copyright (C) 2008 Eric Chassande-Mottin
%% Copyright (C) 2011 Carnë Draug
%% Copyright (C) 2016, 2018, 2022 Colin B. Macdonald
%%
%% This program is free software; you can redistribute it and/or modify it under
%% the terms of the GNU General Public License as published by the Free Software
%% Foundation; either version 3 of the License, or (at your option) any later
%% version.
%%
%% This program is distributed in the hope that it will be useful, but WITHOUT
%% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
%% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
%% details.
%%
%% You should have received a copy of the GNU General Public License along with
%% this program; if not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defun laguerreL (@var{n}, @var{x})
%% Evaluate Laguerre polynomials.
%%
%% Compute the value of the Laguerre polynomial of order @var{n}
%% for each element of @var{x}.
%% For example, the Laguerre polynomial of order 14 evaluated at
%% the point 6 is
%% @example
%% @group
%% @c doctest: +SKIP_IF(compare_versions (OCTAVE_VERSION(), '6.0.0', '<'))
%% laguerreL (14, 6)
%%   @result{} 0.9765
%% @end group
%% @end example
%%
%% This implementation uses a three-term recurrence directly on the values
%% of @var{x}.  The result is numerically stable, as opposed to evaluating
%% the polynomial using the monomial coefficients.  For example, we can
%% compare the above result to a symbolic construction:
%% @example
%% @group
%% syms x
%% L = laguerreL (14, x);
%% exact = subs (L, x, 6)
%%   @result{} exact = (sym)
%%       34213
%%       ─────
%%       35035
%% @end group
%% @end example
%% If we extract the monomial coefficients and numerically evaluate the
%% polynomial at a point, the result is rather poor:
%% @example
%% @group
%% coeffs = sym2poly (L);
%% @c doctest: +XFAIL_IF(compare_versions (OCTAVE_VERSION(), '6.0.0', '<'))
%% polyval (coeffs, 6)
%%   @result{} 0.9765
%% err = ans - double (exact);
%% num2str (err, '%.3g')
%%   @result{} -1.68e-11
%% @end group
%% @end example
%% So please don't do that!  The numerical @code{laguerreL} function
%% does much better:
%% @example
%% @group
%% err = laguerreL (14, 6) - double (exact)
%%   @result{} err = 9.9920e-16
%% @end group
%% @end example
%%
%% @seealso{@@sym/laguerreL}
%% @end defun

function L = laguerreL(n, x)

  if (nargin ~= 2)
    print_usage ();
  end

  if (any (n < 0) || any (mod (n, 1) ~= 0))
    error('second argument "n" must consist of positive integers');
  end

  if (~isscalar (n) && isscalar (x))
    x = x*ones (size (n));
  elseif (~isscalar (n) && ~isscalar (x) && ~isequal (size (n), size (x)))
    error ('inputs must be same size or scalar')
  end

  L0 = ones (size (x), class(x));
  L1 = 1 - x;

  if (isscalar (n))
    if (n == 0)
      L = L0;
    elseif (n == 1)
      L = L1;
    else
      for k = 2:n
        L = (2*k-1-x)/k .* L1 - (k-1)/k * L0;
        L0 = L1;
        L1 = L;
      end
    end
  else
    L = L0;
    L(n >= 1) = L1(n >= 1);
    maxn = max (n(:));
    for k = 2:maxn
      I = (n >= k);  % mask for entries still to be updated
      L(I) = (2*k - 1 - x(I))/k .* L1(I) - (k - 1)/k * L0(I);
      L0 = L1;
      L1 = L;
    end
  end

end

%!error laguerreL (1)
%!error laguerreL (1, 2, 3)

%!assert (isequal (laguerreL (0, rand), 1))

%!test
%! x = rand;
%! assert (isequal (laguerreL (1, x), 1 - x))

%!test
%! x=rand;
%! y1=laguerreL(2, x);
%! p2=[.5 -2 1];
%! y2=polyval(p2,x);
%! assert(y1 - y2, 0, 10*eps);

%!test
%! x=rand;
%! y1=laguerreL(3, x);
%! p3=[-1/6 9/6 -18/6 1];
%! y2=polyval(p3,x);
%! assert(y1 - y2, 0, 20*eps);

%!test
%! x=rand;
%! y1=laguerreL(4, x);
%! p4=[1/24 -16/24 72/24 -96/24 1];
%! y2=polyval(p4,x);
%! assert(y1 - y2, 0, 30*eps)

%!error <positive integer> laguerreL(1.5, 10)
%!error <same size or scalar> laguerreL([0 1], [1 2 3])
%!error <same size or scalar> laguerreL([0 1], [1; 2])

%!test
%! % numerically stable implementation (in n)
%! L = laguerreL (10, 10);
%! Lex = 1763/63;
%! assert (L, Lex, -eps)
%! L = laguerreL (20, 10);
%! Lex = -177616901779/14849255421;  % e.g., laguerreL(sym(20),10)
%! assert (L, Lex, -eps)

%!test
%! % vectorized x
%! L = laguerreL (2, [5 6 7]);
%! Lex = [3.5 7 11.5];
%! assert (L, Lex, eps)

%!test
%! L = laguerreL (0, [4 5]);
%! assert (L, [1 1], eps)

%!test
%! % vector n
%! L = laguerreL ([0 1 2 3], [4 5 6 9]);
%! assert (L, [1 -4 7 -26], eps)

%!test
%! % vector n, scalar x
%! L = laguerreL ([0 1 2 3], 6);
%! assert (L, [1 -5 7 1], eps)

%!assert (isa (laguerreL (0, single (1)), 'single'))
%!assert (isa (laguerreL (1, single ([1 2])), 'single'))
%!assert (isa (laguerreL ([1 2], single ([1 2])), 'single'))

repo: octsympy


text:
function vn = normalizeVector(v)
%NORMALIZEVECTOR Normalize a vector to have norm equal to 1
%
%   V2 = normalizeVector(V);
%   Returns the normalization of vector V, such that ||V|| = 1. V can be
%   either a row or a column vector.
%
%   When V is a MxN array, normalization is performed for each row of the
%   array.
%
%   Example:
%   vn = normalizeVector([3 4])
%   vn =
%       0.6000   0.8000
%   vectorNorm(vn)
%   ans =
%       1
%
%   See Also:
%   vectors2d, vectorNorm
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 29/11/2004.
%

%   HISTORY
%   2005-01-14 correct bug
%   2009-05-22 rename as normalizeVector
%   2011-01-20 use bsxfun

dim = size(v);

if dim(1)==1 || dim(2)==1
    % in case of one vector, the norm is a scalar
    vn = v / sqrt(sum(v.^2));
else
    % for several vectors, need to adapt size of norm
    vn = bsxfun(@rdivide, v, sqrt(sum(v.^2, 2)));
    %same as: vn = v./repmat(sqrt(sum(v.*v, 2)), [1 dim(2)]);
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function map = black2whiteColorMap(varargin)
%Linear white to black gray-scale color map

map = makeColorMap([0 0 0],[0.7 0.7 0.7],[1 1 1]);

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
%            [Lden, Ld] = deninfac(symLden, L,dense,DAt, d, absd, qblkstart,pars)
% DENINFAC
% Uses pars.maxuden as max. allowable |L(i,k)|. Otherwise num. reordering.
%
% ******************** INTERNAL FUNCTION OF SEDUMI ********************
%
% See also sedumi

function [Lden, Ld] = deninfac(symLden, L,dense,DAt, d, absd, qblkstart,pars)
%
% 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
%

% ------------------------------------------------------------
% INITIALIZE: check the numerical tolerance parameters in "pars",
%   and supply with defaults if necessary.
% ------------------------------------------------------------
if nargin < 7
    pars.maxuden = 10;                  % Stability of product-form factors
else
    if ~isfield(pars,'maxuden')
        pars.maxuden = 10;                % Stability of product-form factors
    end
end
% --------------------------------------------------
% Scaling: let Ad = Aden * Dden, as [LP, QBLK, QNORM, QTR]
% with QBLK=DAt.denq, and the others dense.A(:,j)*dj,
% with dj = sqrt(d.l(k)) or dj=sqrt(d.det(k)).
% NOTE that A*P(d)*A' = ADA + Ad * Ad'.
% --------------------------------------------------
if ~isempty(dense.cols)
    i1 = dense.l+1; i2 = i1 + length(dense.q);
    Ad = [dense.A(:,1:i1-1), DAt.denq, dense.A(:,i2:end), dense.A(:,i1:i2-1)];
    smult = [d.l(dense.cols(1:dense.l)); ones(length(dense.q),1); ...
        adenscale(dense,d,qblkstart); -d.det(dense.q)];
    % --------------------------------------------------
    % Correct for sparse part, i.e. let LAD = L\Ad
    % (this also handles L.perm)
    % --------------------------------------------------
    LAD = sparfwslv(L,Ad, symLden.LAD);
    clear Ad;
    % --------------------------------------------------
    % Factor the dense columns as "diag+rank-1"=diag(L.d)+LAD(:,k)*LAD(:,k)'
    % --------------------------------------------------'
    [Lden, Ld] = dpr1fact(LAD, L.d, symLden, smult, pars.maxuden);
    Lden.dz = symLden.dz;
    Lden.first = symLden.first;
    Lden.perm = symLden.perm;
    clear LAD;
else
    % --------------------------------------------------
    % If no dense columns
    % --------------------------------------------------
    Lden.betajc = 0;
    Ld = L.d;
    Lden.rowperm = 1:length(L.d);
end
% ------------------------------------------------------------
% If still Ld(i) = 0, then set to something positive. Otherwise
% we cannot do PCG.
% ------------------------------------------------------------
skip = find(L.skip);
if ~isempty(skip)
    dtol = pars.canceltol * absd(L.perm(skip));
    dtol = max(dtol, pars.abstol);
    Ld(skip(Ld(skip) <= dtol)) = 1;
end

repo: REST


text:
function [ortirfmatrix]=irfsim_new(beta,D,n,m,p,k,horizon)

% [irfmatrix ortirfmatrix]=bear.irfsim(beta,D,n,m,p,k,horizon)
% computes IRF matrices and orthogonalised IRF matrices
% inputs:  - vector 'beta': vectorised form of VAR coefficients (defined in 1.1.12)
%          - matrix 'D': structural matrix for the OLS model (defined in 2.3.3)
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'm': number of exogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'p': number of lags included in the model (defined p 7 of technical guide)
%          - integer 'k': number of coefficients to estimate for each equation in the BVAR model (defined p 7 of technical guide)
%          - integer 'horizon': number of IRF periods
% outputs: - matrix 'irfmatrix': record of the series of irf matrices
%          - matrix 'ortirfmatrix': record of the series of orthogonalised irf matrices

% first reshape beta to obtain B
B=reshape(beta,k,n);

% deal with shocks in turn
for ii=1:n

% create a matrix of zeros of dimension p*n
Y=zeros(p,n);
% set the value of the last row, column i, equal to 1
Y(p,ii)=1;

   % repeat the algorithm from period T+1 to period T+h
   for jj=1:horizon-1

   % step 1
   % use the function lagx to obtain the matrix temp, containing the endogenous regressors
   temp=bear.lagx(Y,p-1);

   % step 2
   % define the vector X
   X=[temp(end,:) zeros(1,m)];

   % step 3
   % obtain the predicted value for T+jj
   yp=X*B;

   % step 4
   % concatenate yp at the top of Y
   Y=[Y;yp];

   % repeat until values are obtained for T+h
   end

% consider 'Y' and trim the (p-1) initial periods: what remains is the series of IRFs for period T to period T+h-1
Y=Y(p:end,:);

% record the results in the matrix irfmatrix

   % loop over periods
   for jj=1:horizon

      % loop over variables
      for kk=1:n
      irfmatrix(kk,ii,jj)=Y(jj,kk);
      end

   end

% conduct the same process with shocks in other variables
end

% obtain now orthogonalised IRFs
% loop over periods
for ii=1:horizon
ortirfmatrix(:,:,ii)=irfmatrix(:,:,ii)*D;
end

author: european-central-bank

repo: BEAR-toolbox


text:
%clear all

%fname = 'D:\MatData\umeda data2\modified_motion01_03.mat';
%%load(fname)
%
%ydata = ydata(2:3,:,:);
%
%xdata(:,:,ind) = [];
%ydata(:,:,ind) = [];

%fsave = 'test_spike.mat';
%
%load(fsave)

ypred = predict_output(xtest, Model, parm);
ydev  = predict_variance(xtest, Model, parm);

ymin = min(ydev(:));
ydev = ydev - ymin;

yerr = abs(ypred - ytest);

[err,rcor] = mean_sq_error(yerr,ydev)

NY = 2; NX = 1;
nfig = 1;
subplot(NY,NX,nfig)
%plot(yerr(1,:,1))
plot(yerr(1,:,1)/max(yerr(1,:,1)))
hold on

%nfig = nfig+1;
subplot(NY,NX,nfig)
%plot(ydev(1,:,1),'-r')
hold on
plot(ydev(1,:,1)/max(ydev(1,:,1)),'-r')

return

[M,T,Ntr] = size(xdata);
[N,T,Ntr] = size(ydata);

x = reshape(xdata,[M,T*Ntr]);
y = reshape(ydata,[N,T*Ntr]);

NY = 5; NX = 2;
nfig = 0;

for n=1:N
	nfig = nfig+1;
	if nfig > NX*NY, figure; nfig=1; end
	subplot(NY,NX,nfig)
	plot(y(n,:))
end

figure

NY = 6; NX = 2;

for n=1:M
	subplot(NY,NX,n)
	plot(x(n,:))
end

%save(fsave,'xdata','ydata')

author: KamitaniLab

repo: GenericObjectDecoding


text:
function [M]=Fd_mtx(d, a, bound)
%Finite difference approximation of scalar diffusion equation in QTT
%   [M]=FD_MTX(D, A, BOUND)
%   Generate finite difference matrix from diffusion coefficient A(n,n) or A(n,n,n)
%   M - n^D-by-n^D sparse matrix
%   D - dimensions,
%   BOUND:
%       0 - Dirichlet,
%       1 - Neuman

%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

n = size(a, 1)-2;
if (d==1)
  bar_a = zeros(n,3);
  for i=1:n

     bar_a(i,2)=(a(i+2)+2*a(i+1)+a(i))*0.25;
     bar_a(i,1)=(a(i)+a(i+1))*0.5;
     bar_a(i,3)=(a(i+1)+a(i+2))*0.5;
  end
  D1=2*eye(n)-diag(ones(n-1,1),1)-diag(ones(n-1,1),-1);
  %M=D1*(n+1)^2;
  M=D1;
  bar_a_M=spdiags(bar_a,[-1,0,1],n,n);
  M = M.*bar_a_M;
  return;
end

if (d==2)
    if (bound==0)
        bar_a = zeros(n^2, 5);
        for i=1:n
            for j=1:n
                bar_a((i-1)*n+j, 3) = (a(i+1,j+1)+a(i+2,j+1)+a(i+1,j+2)+a(i+2,j+2))*0.25;
                if (j>1) bar_a((i-1)*n+j-1, 2) = (a(i+1,j+1)+a(i+2,j+1))*0.5; end;
                if (j<n) bar_a((i-1)*n+j+1, 4) = (a(i+1,j+2)+a(i+2,j+2))*0.5; end;
                if (i>1) bar_a((i-1-1)*n+j, 1) = (a(i+1,j+1)+a(i+1,j+2))*0.5; end;
                if (i<n) bar_a((i+1-1)*n+j, 5) = (a(i+2,j+1)+a(i+2,j+2))*0.5; end;
            end;
        end;
        D1 = spdiags(-1*ones(n,1), [-1], n, n)+spdiags(2*ones(n,1), [0], n, n)+spdiags(-1*ones(n,1), [1], n, n);
        M = kron(D1*(n+1)^2, speye(n))+kron(speye(n), D1*(n+1)^2);
        bar_a_M = spdiags(bar_a, [-n, -1, 0, 1, n], n^2, n^2);
        M = M.*bar_a_M;
    else
        wa = zeros(n+3,n+3);
        wa(1:n+2,1:n+2)=a;
        wa(n+3,1:n+2)=a(n+2,:);
        wa(1:n+2,n+3)=a(:,n+2);
        wa(n+3,n+3)=a(n+2,n+2);
        n = n+2;
        bar_a = zeros(n^2, 5);
        for i=1:n
            for j=1:n
                bar_a((i-1)*n+j, 3) = (wa(i,j)+wa(i+1,j)+wa(i,j+1)+wa(i+1,j+1))*0.25;
                if (j>1) bar_a((i-1)*n+j-1, 2) = (wa(i,j)+wa(i+1,j))*0.5; end;
                if (j<n) bar_a((i-1)*n+j+1, 4) = (wa(i,j+1)+wa(i+1,j+1))*0.5; end;
                if (i>1) bar_a((i-1-1)*n+j, 1) = (wa(i,j)+wa(i,j+1))*0.5; end;
                if (i<n) bar_a((i+1-1)*n+j, 5) = (wa(i+1,j)+wa(i+1,j+1))*0.5; end;
            end;
        end;
        D1 = spdiags(-1*ones(n,1), [-1], n, n)+spdiags(2*ones(n,1), [0], n, n)+spdiags(-1*ones(n,1), [1], n, n);
        D1(1,1)=1;
        D1(n,n)=1;
        M = kron(D1*(n-1)^2, speye(n))+kron(speye(n), D1*(n-1)^2);
        bar_a_M = spdiags(bar_a, [-n, -1, 0, 1, n], n^2, n^2);
        M = M.*bar_a_M;
    end;
end;

if (d==3)
    if (bound==0)
        bar_a = zeros(n^3, 7);
        wa=a;
        for i=1:n
            for j=1:n
                for k=1:n
                    bar_a((i-1)*n^2+(j-1)*n+k, 4) = (wa(i+1,j+1,k+1)+wa(i+2,j+1,k+1)+wa(i+1,j+2,k+1)+wa(i+2,j+2,k+1)+wa(i+1,j+1,k+2)+wa(i+2,j+1,k+2)+wa(i+1,j+2,k+2)+wa(i+2,j+2,k+2))*0.125;
                    if (k>1) bar_a((i-1)*n^2+(j-1)*n+k-1, 3) = (wa(i+1,j+1,k+1)+wa(i+2,j+1,k+1)+wa(i+1,j+2,k+1)+wa(i+2,j+2,k+1))*0.25; end;
                    if (k<n) bar_a((i-1)*n^2+(j-1)*n+k+1, 5) = (wa(i+1,j+1,k+2)+wa(i+2,j+1,k+2)+wa(i+1,j+2,k+2)+wa(i+2,j+2,k+2))*0.25; end;

                    if (j>1) bar_a((i-1)*n^2+(j-1-1)*n+k, 2) = (wa(i+1,j+1,k+1)+wa(i+2,j+1,k+1)+wa(i+1,j+1,k+2)+wa(i+2,j+1,k+2))*0.25; end;
                    if (j<n) bar_a((i-1)*n^2+(j+1-1)*n+k, 6) = (wa(i+1,j+2,k+1)+wa(i+2,j+2,k+1)+wa(i+1,j+2,k+2)+wa(i+2,j+2,k+2))*0.25; end;

                    if (i>1) bar_a((i-1-1)*n^2+(j-1)*n+k, 1) = (wa(i+1,j+1,k+1)+wa(i+1,j+2,k+1)+wa(i+1,j+1,k+2)+wa(i+1,j+2,k+2))*0.25; end;
                    if (i<n) bar_a((i+1-1)*n^2+(j-1)*n+k, 7) = (wa(i+2,j+1,k+1)+wa(i+2,j+2,k+1)+wa(i+2,j+1,k+2)+wa(i+2,j+2,k+2))*0.25; end;
                end;
            end;
        end;
        D1 = spdiags(-1*ones(n,1), [-1], n, n)+spdiags(2*ones(n,1), [0], n, n)+spdiags(-1*ones(n,1), [1], n, n);
        M = kron(kron(D1*(n+1)^2, speye(n)), speye(n))+kron(kron(speye(n), D1*(n+1)^2), speye(n)) + kron(kron(speye(n), speye(n)), D1*(n+1)^2);
        bar_a_M = spdiags(bar_a, [-n^2, -n, -1, 0, 1, n, n^2], n^3, n^3);
        M = M.*bar_a_M;
    else
        wa = zeros(n+3,n+3,n+3);
        wa(1:n+2, 1:n+2, 1:n+2) = a;
        wa(n+3, 1:n+2, 1:n+2) = a(n+2, :, :);
        wa(1:n+2, n+3, 1:n+2) = a(:, n+2, :);
        wa(1:n+2, 1:n+2, n+3) = a(:, :, n+2);
        wa(1:n+2, n+3, n+3) = a(:, n+2, n+2);
        wa(n+3, 1:n+2, n+3) = a(n+2, :, n+2);
        wa(n+3, n+3, 1:n+2) = a(n+2, n+2, :);
        wa(n+3,n+3,n+3)=a(n+2,n+2,n+2);
        n = n+2;
        bar_a = zeros(n^3, 7);
        for i=1:n
            for j=1:n
                for k=1:n
                    bar_a((i-1)*n^2+(j-1)*n+k, 4) = (wa(i,j,k)+wa(i+1,j,k)+wa(i,j+1,k)+wa(i+1,j+1,k)+wa(i,j,k+1)+wa(i+1,j,k+1)+wa(i,j+1,k+1)+wa(i+1,j+1,k+1))*0.125;
                    if (k>1) bar_a((i-1)*n^2+(j-1)*n+k-1, 3) = (wa(i,j,k)+wa(i+1,j,k)+wa(i,j+1,k)+wa(i+1,j+1,k))*0.25; end;
                    if (k<n) bar_a((i-1)*n^2+(j-1)*n+k+1, 5) = (wa(i,j,k+1)+wa(i+1,j,k+1)+wa(i,j+1,k+1)+wa(i+1,j+1,k+1))*0.25; end;

                    if (j>1) bar_a((i-1)*n^2+(j-1-1)*n+k, 2) = (wa(i,j,k)+wa(i+1,j,k)+wa(i,j,k+1)+wa(i+1,j,k+1))*0.25; end;
                    if (j<n) bar_a((i-1)*n^2+(j+1-1)*n+k, 6) = (wa(i,j+1,k)+wa(i+1,j+1,k)+wa(i,j+1,k+1)+wa(i+1,j+1,k+1))*0.25; end;

                    if (i>1) bar_a((i-1-1)*n^2+(j-1)*n+k, 1) = (wa(i,j,k)+wa(i,j+1,k)+wa(i,j,k+1)+wa(i,j+1,k+1))*0.25; end;
                    if (i<n) bar_a((i+1-1)*n^2+(j-1)*n+k, 7) = (wa(i+1,j,k)+wa(i+1,j+1,k)+wa(i+1,j,k+1)+wa(i+1,j+1,k+1))*0.25; end;
                end;
            end;
        end;
        D1 = spdiags(-1*ones(n,1), [-1], n, n)+spdiags(2*ones(n,1), [0], n, n)+spdiags(-1*ones(n,1), [1], n, n);
        D1(1,1)=1;
        D1(n,n)=1;
        M = kron(kron(D1*(n-1)^2, speye(n)), speye(n))+kron(kron(speye(n), D1*(n-1)^2), speye(n)) + kron(kron(speye(n), speye(n)), D1*(n-1)^2);
        bar_a_M = spdiags(bar_a, [-n^2, -n, -1, 0, 1, n, n^2], n^3, n^3);
        M = M.*bar_a_M;
    end;
end;

end

author: oseledets

repo: TT-Toolbox


text:
function draw_full_net(net,name)
%input: net (DAGNN format)  output_filename output: jpg
%This code is modified from http://stackoverflow.com/questions/5518200/automatically-generating-a-diagram-of-function-calls-in-matlab
OutputName = 'net';
if ~isempty(name)
    OutputName = name;
end
dotFile = [OutputName '.dot'];
fid = fopen(dotFile, 'w');
fprintf(fid, 'digraph G {\n');
fprintf(fid,'rankdir=LR;\n');
for i = 1:numel(net.layers)
    for j = 1:numel(net.layers(i).inputs)  % multi-input to one output
        fprintf( fid,'%s -> %s', char(net.layers(i).inputs(j)), char(net.layers(i).outputs));
        block = net.layers(i).block;
        block_class = class(block);
        if (isequal(block_class ,'dagnn.Conv'))
            fh = block.size(1);
            fw = block.size(2);
            fcount=block.size(4);
            step=block.stride(1);
            fprintf (fid,'[label = %s%dx%dx%dStride%d];\n',block_class(7:end),fh,fw,fcount,step);
        elseif isequal(block_class ,'dagnn.Pooling')
            fh = block.poolSize(1);
            fw = block.poolSize(2);
            step=block.stride(1);
            fprintf (fid,'[label = %s%dx%dStride%d];\n',block_class(7:end),fh,fw,step);
       % elseif isequal(block_class ,'dagnn.DropOut')
           % rate = block.rate;
            %fprintf (fid,'[label = %s%.1f];\n',block_class(7:end),rate);
        else
            fprintf (fid,'[label = %s];\n',block_class(7:end));
        end
    end
end
fprintf(fid, '}\n');
fclose(fid);

% Render to image
imageFile = [OutputName '.png'];
% Assumes the GraphViz bin dir is on the path; if not, use full path to dot.exe
cmd = sprintf('dot -Tpng -Gsize="2048,2048" "%s" -o"%s"', dotFile, imageFile);  % for better view, you can use number bigger than 32
system(cmd);
fprintf('Wrote to %s\n', imageFile);
%im = imread(imageFile);
%imshow(im);

author: layumi

repo: Image-Text-Embedding


text:
% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.
function c = colorize2(im, color)

    if nargin < 2
        color = [1 1 1];
    end
    c = [];
    for i=1:numel(color)
        c = cat(3, c, im*color(i));
    end

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [J, grad] = costFunction(theta, X, y)
%COSTFUNCTION Compute cost and gradient for logistic regression
%   J = COSTFUNCTION(theta, X, y) computes the cost of using theta as the
%   parameter for logistic regression and the gradient of the cost
%   w.r.t. to the parameters.

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

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

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Note: grad should have the same dimensions as theta
%

% X.shape = [m, n], y.shape = [m, 1], theta.shape = [n, 1]
% y_pred = sigmoid(X*theta), y_pred.shape = [m, 1]
y_pred = sigmoid(X*theta);
J = -1.0 / m * sum( y.*log(y_pred) + (1-y).*log(1-y_pred) );
grad = 1.0 / m .* (X' * (y_pred - y));

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

end

author: imLogM


text:
function U = randunitary(n, N)
% Generates uniformly random unitary matrices.
%
% function U = randunitary(n, N)
%
% U is a n-by-n-by-N array such that each slice U(:, :, i) is a random
% unitary matrix of size n (i.e., a matrix in the unitary group U(n)),
% sampled from the Haar measure (uniform distribution).
%
% By default, N = 1.
%
% Complexity: N times O(n^3).
% For details on the algorithm, see Mezzadri 2007,
% "How to generate random matrices from the classical compact groups."
%
% See also: randrot qr_unique

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, June 18, 2019.
% Contributors:
% Change log:

    if nargin < 2
        N = 1;
    end

    if n == 1
        U = sign(randn(1, 1, N) + 1i*randn(1, 1, N));
        return;
    end

    % Generated as such, the slides of U are uniformly distributed over
    % U(n), the set of unitary matrices: see Mezzadri 2007, p597.
    U = qr_unique(randn(n, n, N) + 1i*randn(n, n, N));

end

author: NicolasBoumal

repo: manopt


text:
% new procedure: generate motor seed/map in trialRes space
% find top 100 cells correlating to motor, screen with Rh4+Rh5(+Rh6??) masks

% corr sweep (0.5-0.7) to get top %2 cells, then kmeans, save; rank by stim-lock, save
clear all; close all; clc

%% folder setup
isSaveFig = 1;
isPlotFig = 1;

outputDir = GetOutputDataDir;
saveDir = [];
saveDir{1} = fullfile(outputDir,'motor_map_notseed_lrRes_1_kmeans_030417');
saveDir{2} = fullfile(outputDir,'motor_map_notseed_lrRes_2_autoclus_030417');
if ~exist(saveDir{1}, 'dir'), mkdir(saveDir{1}), end;
if ~exist(saveDir{2}, 'dir'), mkdir(saveDir{2}), end;

%% init

hfig = figure;
InitializeAppData(hfig);
ResetDisplayParams(hfig);

setappdata(hfig,'isMotorseed',0);
setappdata(hfig,'isTrialRes',1);

%% run fish
range_fish = GetFishRange;%[1:3,5:18];
M_thres_reg = zeros(3,18);
M_numTopCorr = zeros(1,18);
M_motorseedRegs = cell(1,18);
M_compareMotorCellNumber = zeros(2,18);

for i_fish = range_fish
    ClusterIDs = [1,1];
    [~,~,~,stim,behavior,M_0] = LoadSingleFishDefault(i_fish,hfig,ClusterIDs);

    MASKs = getappdata(hfig,'MASKs');
    CellXYZ_norm = getappdata(hfig,'CellXYZ_norm');
    absIX = getappdata(hfig,'absIX');

    %% get motor regressors
    [~,~,regressor_m_raw] = GetMotorRegressor(behavior,i_fish);
    Reg = regressor_m_raw([1,3],:);

    %% regression (in tRes), thresholding by % of cells (instead of corr thres)

    %     [~,Reg_tRes] = GetTrialAvrLongTrace(hfig,Reg);
    %     [~,M_0_tRes] = GetTrialAvrLongTrace(hfig,M_0);

%     Reg = FindClustermeans(gIX_seed,M);
    Reg_LR = Reg-repmat(mean(Reg),2,1);

    %     Corr = corr(Reg_tRes_LR',M_0_tRes');
    Corr = corr(Reg_LR',M_0');
    [Corr_sorted,IX_corr] = sort(Corr,2,'descend');

    %     [corr_max,IX] = max(Corr,[],1);
    %     [~,I] = sort(corr_max,'descend');

    M_target_prctcell= [1,2,3];
    for i_numcell = 1:2,
        nCells_total = size(M_0,1);
        prctcell = M_target_prctcell(i_numcell);
        nCells_target = round(prctcell/100 * nCells_total);

        %% target cell number: counting cells in hindbrain only

        Msk_IDs = 114; % mask for full hindbrain

        % isScreenMskFromAllCells
        cIX = (1:length(absIX))';
        gIX = ones(size(cIX));
        [cIX_hb,gIX_hb] = ScreenCellsWithMasks(Msk_IDs,cIX,gIX,MASKs,CellXYZ_norm,absIX);

        I_hb = ismember(IX_corr,cIX_hb);
        cum_I_hb1 = cumsum(I_hb(1,:));
        cum_I_hb2 = cumsum(I_hb(2,:));
        lastIX{1} = find(cum_I_hb1==nCells_target,1,'first');
        lastIX{2} = find(cum_I_hb2==nCells_target,1,'first');

        % cut off the desired number of cells to display
        for i_lr = 1:2
            cIX = IX_corr(i_lr,1:lastIX{i_lr})';
            gIX = ceil((1:length(cIX))'/(length(cIX)/20));
            M = UpdateIndices_Manual(hfig,cIX,gIX);
            [~,M_tRes] = GetTrialAvrLongTrace(hfig,M);
            setappdata(hfig,'M',M_tRes);

            M_thres_reg(i_lr,i_fish) = Corr_sorted(i_lr,lastIX{i_lr});
            M_compareMotorCellNumber(i_lr,i_fish) = length(cIX);

            %% clustering: trying Autoclus here!
            if i_numcell==1
                %%
                numK = 16; % manual
                gIX = Kmeans_Direct(M_tRes,numK);
            else
                cIX_reg = cIX;%(1:size(M_0,1))';
                [cIX,gIX] = AutoClustering(cIX,gIX,M_0,cIX_reg);
                UpdateIndices_Manual(hfig,cIX,gIX);
            end
            %% Plot anat
            if isPlotFig,

                I = LoadCurrentFishForAnatPlot(hfig,cIX,gIX);
                f{i_lr} = DrawCellsOnAnat(I);
%                 f{i_lr} = figure('Position',[50,100,1400,800]);
%                 % isCentroid,isPlotLines,isPlotBehavior,isPlotRegWithTS
%                 subplot(121)
%                 setappdata(hfig,'isPlotBehavior',1);
%                 setappdata(hfig,'isStimAvr',0);
%                 setappdata(hfig,'isPlotLines',0);
%                 %             UpdateTimeIndex(hfig);
%                 DrawTimeSeries(hfig,cIX,gIX);
%
%                 % right plot
%                 ax = subplot(122);
%                 I = LoadCurrentFishForAnatPlot(hfig,cIX,gIX);
%                 DrawCellsOnAnat(I,ax);

            end
        end

        combineFiguresLR(f{1},f{2});

        % Save plot
        if isSaveFig,
            filename = fullfile(saveDir{i_numcell}, ['Fish',num2str(i_fish),'_motormap_',num2str(prctcell),'%_each']);
            saveas(gcf, filename, 'png');
            close(gcf)
        end

    end
end

%%
% range_fish excludes Fish 4
% M_compareMotorCellNumber(:,4) = NaN;
% figure;bar(M_compareMotorCellNumber')

author: xiuyechen

repo: FishExplorer

path: github-repos/MATLAB/xiuyechen-FishExplorer/FishExplorer-c61392cf0835480d64fc03c15f1992935fdc7106/figure scripts/SensoryMotor/arc code/motormap, tRes versions, comparisons, w Clusters or Autoclus to follow/fig3_motormap_trialRes_lrRes_notseed.m


text:
%run this test case with the command
%results = runtests('hmmerTests.m')
function tests = hmmerTests
tests = functiontests(localfunctions);
end

function testHmmer(testCase)
%This unit test comprises the functionality test for HMMER in RAVEN:
% 1. Check of parsed HMMER results against the expected.
%
% NOTE: as hmm and HMMER results files are time-specific, no checks for
% these files existence are done. Also, due to the way HMMER is utilized in
% getKEGGModelForOrganism (HMMER result files can be parsed only once all
% required hmm files are generated), the code segment involving HMMER
% results parsing is pasted in this test function. Should the parsing problems
% occur in the results processing, the code modifications shall be done in
% this function and getKEGGModelForOrganism respectively.

%%
%Get the directory for RAVEN Toolbox
[ST, I]=dbstack('-completenames');
ravenPath=fileparts(fileparts(fileparts(ST(I).file)));

%Identify the operating system
if isunix
    if ismac
        binEnd='.mac';
    else
        binEnd='';
    end
elseif ispc
    binEnd='.exe';
else
    dispEM('Unknown OS, exiting.')
    return
end

%Create empty structures needed for HMMER results
actHmmResult.genes = {};
actHmmResult.scores = [];

%Create structures that contain expected HMMER results
expHmmResult.genes = {'sp|P41947|MEL6_YEASX','sp|P41946|MEL5_YEASX', 'sp|P41945|MEL2_YEASX', 'sp|P04824|MEL1_YEASX'};
expHmmResult.scores = [10^-250, 10^-250, 10^-250, 10^-250];

%Generate temporary names for working directory and outFile
tmpDIR=tempname;
outFile=tempname;

%Run HMMER multi-threaded to use all logical cores assigned to MATLAB
cores = evalc('feature(''numcores'')');
cores = strsplit(cores, 'MATLAB was assigned: ');
cores = regexp(cores{2},'^\d*','match');
cores = cores{1};

%Create a temporary folder and copy multi-FASTA file there
[~, ~]=system(['mkdir "' tmpDIR '"']);

sourceDir = fileparts(which(mfilename));
copyfile(fullfile(sourceDir,'test_data','yeast_galactosidases.fa'),tmpDIR);
copyfile(fullfile(sourceDir,'test_data','human_galactosidases.fa'),tmpDIR);

%%
%Train a hidden Markov model
[~, ~]=system(['"' fullfile(ravenPath,'software','hmmer',['hmmbuild' binEnd]) '" --cpu "' num2str(cores) '" "' fullfile(tmpDIR,'human_galactosidases.hmm') '" "' fullfile(tmpDIR,'yeast_galactosidases.fa') '"']);

%Run a homology search against the newly-trained HMM
[~, output]=system(['"' fullfile(ravenPath,'software','hmmer',['hmmsearch' binEnd]) '" --cpu "' num2str(cores) '" "' fullfile(tmpDIR,'human_galactosidases.hmm') '" "' fullfile(tmpDIR,'yeast_galactosidases.fa') '"']);

%Save the output to a file
fid=fopen(outFile,'w');
fwrite(fid,output);
fclose(fid);

%%
%Parse the results
geneCounter=0;
fid=fopen(outFile,'r');
beginMatches=false;
while 1
    %Get the next line
    tline = fgetl(fid);

    %Abort at end of file
    if ~ischar(tline)
        break;
    end

    if and(beginMatches,strcmp(tline,'  ------ inclusion threshold ------'))
        break;
    end

    if beginMatches==false
        %This is how the listing of matches begins
        if any(strfind(tline,'E-value '))
            %Read one more line that is only padding
            tline = fgetl(fid);
            beginMatches=true;
        end
    else
        %If matches should be read
        if ~strcmp(tline,'   [No hits detected that satisfy reporting thresholds]') && ~isempty(tline)
            elements=regexp(tline,' ','split');
            elements=elements(cellfun(@any,elements));

            %Check if the match is below the treshhold
            score=str2double(elements{1});
            gene=elements{9};
            if score<=10^-50
                %If the score is exactly 0, change it to a very
                %small value to avoid NaN
                if score==0
                    score=10^-250;
                end
                %Check if the gene is added already and, is so, get
                %the best score for it
                geneCounter=geneCounter+1;
                actHmmResult.genes{geneCounter}=gene;
                actHmmResult.scores(geneCounter)=score;
            end
        else
            break;
        end
    end
end
fclose(fid);

%Remove the old tempfiles
delete([outFile '*']);

%Remove temporary folder, since testing is finished
[~, ~]=system(['rm "' tmpDIR '" -r']);

%%
%Test 1a: Check if HMMER results match the expected ones
verifyEqual(testCase,actHmmResult,expHmmResult);

%Test 1b: Modify actual HMMER results structure and check if test fails
actHmmResult.score(2)=1;
verifyNotEqual(testCase,actHmmResult,expHmmResult);
end

author: SysBioChalmers

repo: RAVEN


text:
% ir_mri_partial_fourier_3d_test2.m
% test 2D partial Fourier MRI method ir_mri_partial_fourier_3d.m
% to examine behavior of xu:01:pfi
% 2014-06-? Mai Le
% 2014-06-22 JF mods

%% --- Parameters that control recon quality

percentage = 0.65; % partial k-space fraction, must be greater than 0.5

window_step3 = 12; % controls transition band width for apodization
window_step8 = 6;

%% --- Chose test case

nz = 1;
% uncomment line below to test 3D
% even number just big enough to pass condition in PF_3D that requires apodization area of 16:
%nz = 2*ceil((2*max(window_step3,window_step8)+1) / (2*(percentage - 0.5)) / 2);

% Shepp-Logan
%xtrue = phantom('Modified Shepp-Logan');
if nz > 1
	xtrue = ellipse_im(64); % need even dims
	xtrue = xtrue(5:64,:); % stress test non-square
	xtrue = repmat(xtrue, [1 1 nz]);
im(xtrue) % fails in octave - todo try in matlab
return
	xtrue(:,:,1) = zeros(size(xtrue(:,:,1))); xtrue(10, 50) = 1; % impulse
else
	xtrue = ellipse_im(128);
	xtrue = xtrue(5:124,:); % stress test non-square
%	xtrue = zeros(size(xtrue)); xtrue(10, 80) = 1; % impulse
end
[nx ny nz] = size(xtrue)

% make it complex
xtrue = xtrue .* exp(-1i*0.05); % constant phase - simple case
xtrue = xtrue .* exp(-1i*2*pi*(([1:nx]'*[1:ny])/(nx*ny)).^2); % nonlinear phase

% Tough Case
%xtrue = 1 + rand(nx, ny, nz) .* exp(1i*0.05); % not quite constant phase :(

kspace = fftshift(fftn(xtrue));

if 0 % Artificial k-space for visualizing effect of algo
	kspace = rand(size(xtrue));
	xtrue = ifftn(ifftshift(k));
end

%% --- Select partial k-space and do PF recon

pf_location = [1 0 0]; % stress test
nx_pf = round(percentage*nx);
ny_pf = round(percentage*ny);
nz_pf = round(percentage*nz);
kx = (1:nx_pf) + pf_location(1) * (nx - nx_pf);
ky = (1:ny_pf) + pf_location(2) * (ny - ny_pf);
kz = (1:nz_pf) + pf_location(3) * (nz - nz_pf);

sampling = false(size(xtrue));
sampling(kx, ky, kz) = true;
% partial_kspace = sampling .* k;
partial_kspace = kspace(kx, ky, kz);

% log scale over 6 digits (single precision)
im_log = @(i, x, t) ...
	im(i, log10(abs(x) / max(abs(kspace(:))) + 1e-20), [-6 0], t);

ir_fontsize im_axes 5
ir_fontsize title 10
ir_fontsize label 10
im plc 3 5
im(1, abs(xtrue), 'true'), cbar
im(2, sampling), cbar
tmp = @(k, n) unique([1 n minmax(k)']);
xtick(tmp(kx, nx))
ytick(tmp(ky, ny))
im_log(6, kspace, 'full k-space'), cbar
im_log(7, embed(partial_kspace(:), sampling), 'partial'), cbar

if nz > 1
	full_dims = [nx ny nz];
else
	full_dims = [nx ny];
end

if ~isvar('img_est')
	[img_est, full_kspace] = ir_mri_partial_fourier_3d(...
		partial_kspace, full_dims, 'pf_location', pf_location, ...
		'chat', 1, 'show', 1, 'niter', 10, ...
		'window_step3', window_step3, 'window_step8', window_step8);
%		'init', xtrue, ... % of course it works
end

im(5, abs(img_est)), title '|img_est|', cbar
im(10, angle(img_est)), title '\angle < img_est', cbar

diff_img = img_est - xtrue;
%im(9, abs(img), '|Recon|'), cbar
im(11, abs(diff_img), 'Error'), cbar
xlabelf('max = %.3g', max(abs(diff_img(:))))
im_log(15, full_kspace, 'kspace est'), cbar
im_log(12, full_kspace - kspace, 'kspace err'), cbar

if 0
	im plc 1 1
	im_log(1, full_kspace - kspace, 'kspace err')
end

%im_log(11, full_kspace, 'estimated')

% conclusions:
% - perfect recovery not possible
% - missing corners of k-space filled in by convolution of fft of phase difference between current iterate and low frequency estimate

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function p = psnr( I, J )

% Compute PSNR between images
% Peak signal-to-noise ratio; built-in function not used since it
% assumes pixel values in interval [0, 255].
%
% Syntax:
%	p = psnr( I, J )
%
%
% Input:
%	I : Original image.
%
%	J : Denoised image.
%
%
% Output:
%	p : PSNR between I and J.

if ~isequal( size(I), size(J) )
	error('Images must be of the same size');
end

pixel_range = max(I(:)) - min(I(:));
MSE = norm( I - J, 'fro' );

p = 20*log10( pixel_range * sqrt(numel(I)) / MSE );

end


text:
classdef WBDetector < handle
    %WBDETECTOR  WaldBoost detector - Object Detection using Boosted Features
    %
    % Class for object detection using WaldBoost from [Sochman05].
    % Uses a Waldboost cascade and local binary patterns computed as integral
    % features for 2D object detection.
    %
    % ## Example
    % The basic usage is the following:
    %
    %     % train
    %     detector = cv.WBDetector();
    %     detector.train('/path/to/pos/', '/path/to/neg/');
    %     detector.write('model.xml');
    %
    %     % detect
    %     detector = cv.WBDetector();
    %     detector.read('model.xml');
    %     img = cv.imread('image.png', 'Grayscale',true);
    %     [bboxes,conf] = detector.detect(img);
    %     for i=1:numel(bboxes)
    %         img = cv.rectangle(img, bboxes{i}, 'Color',[0 255 0]);
    %     end
    %     imshow(img)
    %
    % ## References
    % [Sochman05]:
    % > J. Sochman and J. Matas. "WaldBoost - Learning for Time Constrained
    % > Sequential Detection", IEEE Conference on CVPR 2005, Vol 2, p. 150-156
    % > [PDF](https://dspace.cvut.cz/bitstream/handle/10467/9494/2005-Waldboost-learning-for-time-constrained-sequential-detection.pdf?sequence=1)
    %
    % See also: cv.WBDetector.WBDetector, cv.WBDetector.detect,
    %  detectPeopleACF, peopleDetectorACF, trainACFObjectDetector
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    methods
        function this = WBDetector()
            %WBDETECTOR  Create instance of WBDetector
            %
            %     detector = cv.WBDetector()
            %
            % See also: cv.WBDetector, cv.WBDetector.train, cv.WBDetector.read
            %
            this.id = WBDetector_(0, 'new');
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     detector.delete()
            %
            % See also: cv.WBDetector
            %
            if isempty(this.id), return; end
            WBDetector_(this.id, 'delete');
        end

        function read(this, fname_or_str, varargin)
            %READ  Read detector from file
            %
            %     detector.read(filename)
            %     detector.read(str, 'FromString',true)
            %     detector.read(..., 'OptionName', optionValue, ...)
            %
            % ## Input
            % * __filename__ Name of the file to read from.
            % * __str__ String containing serialized detector you want to load.
            %
            % ## Options
            % * __FromString__ Logical flag to indicate whether the input is a
            %   a filename or a string containing the serialized detector.
            %   default false
            %
            % See also: cv.WBDetector.write
            %
            WBDetector_(this.id, 'read', fname_or_str, varargin{:});
        end

        function varargout = write(this, filename)
            %WRITE  Write detector to file
            %
            %     detector.write(filename)
            %     str = detector.write(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % ## Output
            % * __str__ optional output. If requested, the object is persisted
            %   to a string in memory instead of writing to disk.
            %
            % See also: cv.WBDetector.read
            %
            [varargout{1:nargout}] = WBDetector_(this.id, 'write', filename);
        end

        function train(this, pos_samples, neg_imgs)
            %TRAIN  Train WaldBoost detector
            %
            %     detector.train(pos_samples, neg_imgs)
            %
            % ## Input
            % * **pos_samples** Path to directory with cropped positive
            %   samples.
            % * **neg_imgs** Path to directory with negative (background)
            %   images.
            %
            % See also: cv.WBDetector.detect
            %
            WBDetector_(this.id, 'train', pos_samples, neg_imgs);
        end

        function [bboxes,confidences] = detect(this, img)
            %DETECT  Detect objects on image using WaldBoost detector
            %
            %     [bboxes,confidences] = detector.detect(img)
            %
            % ## Input
            % * __img__ Input image for detection.
            %
            % ## Output
            % * __bboxes__ Bounding boxes coordinates.
            % * __confidences__ Confidence values for bounding boxes.
            %
            % See also: cv.WBDetector.WBDetector
            %
            [bboxes,confidences] = WBDetector_(this.id, 'detect', img);
        end
    end

end


text:
close all;
clear all;
clc;

addpath([cd '/utilies']);
load(['AR_EigenFace']);

%%%%%%%%%%%%%%%%%%%%%%%%
%FDDL parameter
%%%%%%%%%%%%%%%%%%%%%%%%
opts.nClass        =   100;
opts.wayInit       =   'PCA';
opts.lambda1       =   0.005;
opts.lambda2       =   0.05;
opts.nIter         =   15;
opts.show          =   true;
[Dict,Drls,CoefM,CMlabel] = FDDL(tr_dat,trls,opts);

%%%%%%%%%%%%%%%%%%%%%%%%
% Sparse Classification
%%%%%%%%%%%%%%%%%%%%%%%%
lambda   =   0.005;
nClass   =   opts.nClass;
weight   =   0.5;

td1_ipts.D    =   Dict;
td1_ipts.tau1 =   lambda;
if size(td1_ipts.D,1)>=size(td1_ipts.D,2)
   td1_par.eigenv = eigs(td1_ipts.D'*td1_ipts.D,1);
else
   td1_par.eigenv = eigs(td1_ipts.D*td1_ipts.D',1);
end

ID   =   [];
for indTest = 1:size(tt_dat,2)
    fprintf(['Totalnum:' num2str(size(tt_dat,2)) 'Nowprocess:' num2str(indTest) '\n']);
    td1_ipts.y          =      tt_dat(:,indTest);
    [opts]              =      IPM_SC(td1_ipts,td1_par);
    s                   =      opts.x;

    for indClass  =  1:nClass
        temp_s            =  zeros(size(s));
        temp_s(indClass==Drls) = s(indClass==Drls);
        zz                =  tt_dat(:,indTest)-td1_ipts.D*temp_s;
        gap(indClass)     =  zz(:)'*zz(:);

        mean_coef_c         =   CoefM(:,indClass);
        gCoef3(indClass)    =  norm(s-mean_coef_c,2)^2;
    end

    wgap3  = gap + weight*gCoef3;
    index3 = find(wgap3==min(wgap3));
    id3    = index3(1);
    ID     = [ID id3];
end

fprintf('%s%8f\n','reco_rate  =  ',sum(ID==ttls)/(length(ttls)));

author: hiroyuki-kasai

repo: ClassifierToolbox