text:
function [mustUSet, posMustU] = findMustU(model, minFluxesW, maxFluxesW, varargin)
% This function runs the second step of `optForce`, that is to solve a
% bilevel mixed integer linear programming  problem to find a first order
% `MustU` set.
%
% USAGE:
%
%    [mustUSet, posMustU] = findMustU(model, minFluxesW, maxFluxesW, varargin)
%
% INPUTS:
%    model:                      (structure) COBRA metabolic model
%                                with at least the following fields:
%
%                                 * .rxns - Reaction IDs in the model
%                                 * .mets - Metabolite IDs in the model
%                                 * .S -    Stoichiometric matrix (sparse)
%                                 * .b -    RHS of `Sv = b` (usually zeros)
%                                 * .c -    Objective coefficients
%                                 * .lb -   Lower bounds for fluxes
%                                 * .ub -   Upper bounds for fluxes
%    minFluxesW:                (double array) of size `n_rxns x 1`
%                               Minimum fluxes for each
%                               reaction in the model for wild-type strain.
%                               This can be obtained by running the
%                               function `FVAOptForce`.
%                               E.g.: `minFluxesW = [-90; -56];`
%    maxFluxesW:                (double array) of size `n_rxns x 1`
%                               Maximum fluxes for each
%                               reaction in the model for wild-type strain.
%                               This can be obtained by running the
%                               function `FVAOptForce`.
%
% OPTIONAL INPUTS:
%    constrOpt:                 (structure) structure containing
%                               additional contraints. Include here only
%                               reactions whose flux is fixed, i.e.,
%                               reactions whose lower and upper bounds have
%                               the same value. Do not include here
%                               reactions whose lower and upper bounds have
%                               different values. Such contraints should be
%                               defined in the lower and upper bounds of
%                               the model. The structure has the following
%                               fields:
%
%                                 * .rxnList - Reaction list (cell array)
%                                 * .values -  Values for constrained
%                                   reactions (double array)
%                                   E.g.: `struct('rxnList', {{'EX_gluc', 'R75', 'EX_suc'}}, 'values', [-100, 0, 155.5]');`
%    runID:                     (string) ID for identifying this run.
%                               Default: ['run' date hour].
%    outputFolder:              (string) name for folder in which results
%                               will be stored.
%                               Default: 'OutputsFindMustU'.
%    outputFileName:            (string) name for files in which results
%                               will be stored.
%                               Default: 'MustUSet'.
%    printExcel:                (double) boolean to describe wheter data
%                               must be printed in an excel file or not.
%                               Default: 1
%    printText:                 (double) boolean to describe wheter data
%                               must be printed in an plaint text file or not.
%                               Default: 1
%    printReport:               (double) 1 to generate a report in a plain
%                               text file. 0 otherwise.
%                               Default: 1
%    keepInputs:                (double) 1 to save inputs to run
%                               `findMustU.m`, 0 otherwise.
%                               Default: 1
%    printLevel:                (double) 1 to print results in console.
%                               0 otherwise.
%                               Default: 0
%
% OUTPUTS:
%    mustUSet:                  (cell array) Size: number of reactions found X 1
%                               cell array containing the
%                               reactions ID which belong to the `Must_U` Set
%    posMustU:                  (double array) Size: number of reactions found X 1
%                               double array containing the
%                               positions of reactions in the model.
%    outputFileName.xls         File containing one column array
%                               with identifiers for reactions in `MustU`. This
%                               file will only be generated if the user entered
%                               `printExcel = 1`. Note that the user can choose
%                               the name of this file entering the input
%                               `outputFileName = 'PutYourOwnFileNameHere';`
%    outputFileName.txt         File containing one column array
%                               with identifiers for reactions in MustU. This
%                               file will only be generated if the user entered
%                               `printText = 1`. Note that the user can choose
%                               the name of this file entering the input
%                               `outputFileName = 'PutYourOwnFileNameHere';`
%    outputFileName_Info.xls    File containing five column arrays.
%
%                                 * C1: identifiers for reactions in `MustU`
%                                 * C2: min fluxes for reactions according to FVA
%                                 * C3: max fluxes for reactions according to FVA
%                                 * C4: min fluxes achieved for reactions,
%                                   according to `findMustU`
%                                 * C5: max fluxes achieved for reactions,
%                                   according to `findMustU`
%                               This file will only be generated if the user
%                               entered `printExcel = 1`. Note that the user can
%                               choose the name of this file entering the input
%                               `outputFileName = 'PutYourOwnFileNameHere';`
%    outputFileName_Info.txt    File containing five column arrays.
%
%                                 * C1: identifiers for reactions in `MustU`.
%                                 * C2: min fluxes for reactions according to FVA
%                                 * C3: max fluxes for reactions according to FVA
%                                 * C4: min fluxes achieved for reactions,
%                                   according to `findMustU`
%                                 * C5: max fluxes achieved for reactions,
%                                   according to `findMustU`
%                               This file will only be generated if the user
%                               entered `printText = 1`. Note that the user can
%                               choose the name of this file entering the input
%                               `outputFileName = 'PutYourOwnFileNameHere';`
% NOTE:
%
%    This function is based in the GAMS files written by Sridhar
%    Ranganathan which were provided by the research group of Costas D.
%    Maranas. For a detailed description of the `optForce` procedure, please
%    see: `Ranganathan S, Suthers PF, Maranas CD (2010) OptForce: An
%    Optimization Procedure for Identifying All Genetic Manipulations
%    Leading to Targeted Overproductions. PLOS Computational Biology 6(4):
%    e1000744`. https://doi.org/10.1371/journal.pcbi.1000744
%
% .. Author: - Sebastian Mendoza, May 30th 2017, Center for Mathematical Modeling, University of Chile, snmendoz@uc.cl

optionalParameters = {'constrOpt', 'runID', 'outputFolder', 'outputFileName',  ...
    'printExcel', 'printText', 'printReport', 'keepInputs', 'verbose', 'printLevel'};

if (numel(varargin) > 0 && (~ischar(varargin{1}) || ~any(ismember(varargin{1},optionalParameters))))

    tempargin = cell(1,2*(numel(varargin)));
    for i = 1:numel(varargin)

        tempargin{2*(i-1)+1} = optionalParameters{i};
        tempargin{2*(i-1)+2} = varargin{i};
    end
    varargin = tempargin;

end

parser = inputParser();
parser.addRequired('model',@(x) isstruct(x) && isfield(x, 'S') && isfield(model, 'rxns')...
    && isfield(model, 'mets') && isfield(model, 'lb') && isfield(model, 'ub') && isfield(model, 'b')...
    && isfield(model, 'c'))
parser.addRequired('minFluxesW',@isnumeric)
parser.addRequired('maxFluxesW',@isnumeric)
parser.addParamValue('constrOpt', struct('rxnList', {{}}, 'values', []),@(x) isstruct(x) && isfield(x, 'rxnList') && isfield(x, 'values') ...
    && length(x.rxnList) == length(x.values) && length(intersect(x.rxnList, model.rxns)) == length(x.rxnList))
hour = clock; defaultRunID = ['run-' date '-' num2str(hour(4)) 'h' '-' num2str(hour(5)) 'm'];
parser.addParamValue('runID', defaultRunID, @(x) ischar(x))
parser.addParamValue('outputFolder', 'OutputsFindMustU', @(x) ischar(x))
parser.addParamValue('outputFileName', 'MustUSet', @(x) ischar(x))
parser.addParamValue('printExcel', 1, @(x) isnumeric(x) || islogical(x));
parser.addParamValue('printText', 1, @(x) isnumeric(x) || islogical(x));
parser.addParamValue('printReport', 1, @(x) isnumeric(x) || islogical(x));
parser.addParamValue('keepInputs', 1, @(x) isnumeric(x) || islogical(x));
parser.addParamValue('verbose', 1, @(x) isnumeric(x) || islogical(x));
parser.addParamValue('printLevel', 1, @(x) isnumeric(x) || islogical(x));

parser.parse(model, minFluxesW, maxFluxesW, varargin{:})
model = parser.Results.model;
minFluxesW = parser.Results.minFluxesW;
maxFluxesW = parser.Results.maxFluxesW;
constrOpt= parser.Results.constrOpt;
runID = parser.Results.runID;
outputFolder = parser.Results.outputFolder;
outputFileName = parser.Results.outputFileName;
printExcel = parser.Results.printExcel;
printText = parser.Results.printText;
printReport = parser.Results.printReport;
keepInputs = parser.Results.keepInputs;

printFlags = {'printLevel','verbose'};
%get the printLevel.
if all(~ismember(printFlags,parser.UsingDefaults))
    error('Either supply printLevel or verbose optional parameter')
else
    if any(~ismember(printFlags,parser.UsingDefaults))
        selected = ~ismember(printFlags,parser.UsingDefaults);
        printLevel = parser.Results.(printFlags{selected});
    else
        printLevel = parser.Results.printLevel;
    end
end

% correct size of constrOpt
if ~isempty(constrOpt.rxnList)
    if size(constrOpt.rxnList, 1) > size(constrOpt.rxnList,2); constrOpt.rxnList = constrOpt.rxnList'; end;
    if size(constrOpt.values, 1) > size(constrOpt.values,2); constrOpt.values = constrOpt.values'; end;
end

%current path
workingPath = pwd;
runID = [workingPath filesep runID];
%go to the path associate to the ID for this run.
if exist(runID, 'dir')~=7
    mkdir(runID);
end
cd(runID);
outputFolder = [runID filesep outputFolder];

% if the user wants to generate a report.
if printReport
    %create name for file.
    hour = clock;
    reportFileName = ['report-' date '-' num2str(hour(4)) 'h' '-' num2str(hour(5)) 'm.txt'];
    freport = fopen(reportFileName, 'w');
    % print date of running.
    fprintf(freport, ['findMustU.m executed on ' date ' at ' num2str(hour(4)) ':' num2str(hour(5)) '\n\n']);
    % print matlab version.
    fprintf(freport, ['MATLAB: Release R' version('-release') '\n']);

    %print each of the inputs used in this running.
    fprintf(freport, '\nThe following inputs were used to run OptForce: \n');
    fprintf(freport, '\n------INPUTS------\n');
    %print model.
    fprintf(freport, '\nModel:\n');
    for i = 1:length(model.rxns)
        rxn = printRxnFormula(model, model.rxns{i}, false);
        fprintf(freport, [model.rxns{i} ': ' rxn{1} '\n']);
    end
    %print lower and upper bounds, minimum and maximum values for each of
    %the reactions in wild-type and mutant strain
    fprintf(freport, '\nLB\tUB\tMin_WT\tMax_WT\n');
    for i = 1:length(model.rxns)
        fprintf(freport, '%6.4f\t%6.4f\t%6.4f\t%6.4f\n', model.lb(i), model.ub(i), minFluxesW(i), maxFluxesW(i));
    end

    %print constraints
    fprintf(freport,'\nConstrained reactions:\n');
    for i = 1:length(constrOpt.rxnList)
        fprintf(freport,'%s: fixed in %6.4f\n', constrOpt.rxnList{i}, constrOpt.values(i));
    end

    fprintf(freport,'\nrunID(Main Folder): %s \n\noutputFolder: %s \n\noutputFileName: %s \n',...
        runID, outputFolder, outputFileName);

    fprintf(freport,'\nprintExcel: %1.0f \n\nprintText: %1.0f \n\nprintReport: %1.0f \n\nkeepInputs: %1.0f  \n\nverbose: %1.0f \n',...
        printExcel, printText, printReport, keepInputs, printLevel);

end

% export inputs for running the optimization problem to find the MustU Set
inputFolder = [runID filesep 'InputsMustU'];
saveInputsMustSetsFirstOrder(model, minFluxesW, maxFluxesW, constrOpt, inputFolder);

% create a directory to save results if this don't exist
if ~exist(outputFolder, 'dir')
    mkdir(outputFolder);
end

%number of reactions
n_rxns = length(model.rxns);

% initilize sets can, must and empty arrays to store values.
can = zeros(n_rxns, 1);
can(minFluxesW ~= 0) = 1;
can(maxFluxesW ~= 0) = 1;
must = zeros(n_rxns, 1);
mustU = zeros(n_rxns, 1);
vmin = zeros(n_rxns, 1);
vmax = zeros(n_rxns, 1);

found = 0;
%while a solution is still found
while 1
    % create bilevel problem
    bilevelMILPproblem = buildBilevelMILPproblemForFindMustU(model, can, must, maxFluxesW, constrOpt);
    % solve problem
    MustUSol = solveCobraMILP(bilevelMILPproblem, 'printLevel', printLevel);

    if MustUSol.stat~=1
        break;
    else
        % if there is a solution
        found = 1;
        %find which reaction was found
        pos_actives = find(MustUSol.int);
        %update must sets
        must(pos_actives) = 1;
        mustU(pos_actives) = 1;
        %find minimum value for the reaction found
        vmin(pos_actives) = MustUSol.cont(pos_actives);
        %find maximum value for the reaction found
        model2 = changeObjective(model, model.rxns(pos_actives));
        fba = optimizeCbModel(model2);
        vmax(pos_actives) = fba.f;
    end
end

%print the results now.
if printReport; fprintf(freport, '\n------RESULTS------\n'); end;

%if user decide not to show inputs files for findMustU.gms
if ~keepInputs; rmdir(inputFolder,'s'); end;

if found
    %if a solution is found
    if printReport; fprintf(freport, '\na MustU set was found\n'); end;
    if printLevel; fprintf('a MustU set was found\n'); end;
    %find mustU sets
    posMustU = find(mustU);
    mustUSet = model.rxns(posMustU);
else
    %if no solution is found
    if printReport; fprintf(freport, '\na MustU set was not found\n'); end;
    if printLevel; fprintf('a MustU set was not found\n'); end;
    %initilize empty arrays
    mustUSet = {};
    posMustU = {};
end

% print info into an excel text file if required by the user
if printExcel
    if found
        currentFolder = pwd;
        cd(outputFolder);
        Info = [{'Reactions'},{'Min Flux in Wild-type strain'},{'Max Flux in Wild-type strain'},{'Min Flux in Mutant strain'},{'Max Flux in Mutant strain'}];
        Info = [Info; [model.rxns(posMustU), num2cell(minFluxesW(posMustU)), num2cell(maxFluxesW(posMustU)), num2cell(vmin(posMustU)), num2cell(vmax(posMustU))]];
        setupxlwrite();
        xlwrite([outputFileName '_Info'],Info);
        xlwrite(outputFileName, mustUSet);
        cd(currentFolder);
        if printLevel; fprintf(['MustU set was printed in ' outputFileName '.xls  \n']); end;
        if printReport; fprintf(freport, ['\nMustU set was printed in ' outputFileName '.xls  \n']); end;
    else
        if printLevel; fprintf('No mustU set was not found. Therefore, no excel file was generated\n'); end;
        if printReport; fprintf(freport, '\nNo mustU set was not found. Therefore, no excel file was generated\n'); end;
    end
end

% print info into a plain text file if required by the user
if printText
    if found
        currentFolder = pwd;
        cd(outputFolder);
        f = fopen([outputFileName '_Info.txt'], 'w');
        fprintf(f, 'Reactions\tMin Flux in Wild-type strain\tMax Flux in Wild-type strain\tMin Flux in Mutant strain\tMax Flux in Mutant strain\n');
        for i=1:length(posMustU)
            fprintf(f, '%s\t%4.4f\t%4.4f\t%4.4f\t%4.4f\n', model.rxns{posMustU(i)}, minFluxesW(posMustU(i)), maxFluxesW(posMustU(i)), vmin(posMustU(i)), vmax(posMustU(i)));
        end
        fclose(f);

        f = fopen([outputFileName '.txt'], 'w');
        for i = 1:length(posMustU)
            fprintf(f, '%s\n', mustUSet{i});
        end
        fclose(f);

        cd(currentFolder);
        if printLevel; fprintf(['MustU set was printed in ' outputFileName '.txt  \n']); end;
        if printReport; fprintf(freport, ['\nMustU set was printed in ' outputFileName '.txt  \n']); end;
    else
        if printLevel; fprintf('No mustU set was found. Therefore, no excel file was generated\n'); end;
        if printReport; fprintf(freport, '\nNo mustU set was found. Therefore, no excel file was generated\n'); end;
    end
end

%close file for saving report
if printReport; fclose(freport); end;
if printReport; movefile(reportFileName, outputFolder); end;

%go back to the original path
cd(workingPath);

end

function bilevelMILPproblem = buildBilevelMILPproblemForFindMustU(model, can, must, maxFluxesW, constrOpt)

if nargin < 5 || isempty(constrOpt.rxnList)
    ind_ic = [];
    b_ic = [];
    sel_ic = zeros(length(model.rxns), 1);
    sel_ic_b = zeros(length(model.rxns), 1);
else
    %get indices of rxns
    [~, ind_a, ind_b] = intersect(model.rxns, constrOpt.rxnList);
    aux = constrOpt.values(ind_b);
    %sort for rxn index
    [sorted, ind_sorted] = sort(ind_a);
    ind_ic = sorted;
    b_ic = aux(ind_sorted); if size(b_ic, 1) > size(b_ic, 2); b_ic = b_ic'; end;
    sel_ic = zeros(length(model.rxns), 1);
    sel_ic(ind_ic) = 1;
    sel_ic_b = zeros(length(model.rxns), 1);
    sel_ic_b(ind_ic) = b_ic;
end

%convert inputs
S = model.S;
ub = model.ub;
lb = model.lb;
% Dimensions
[n_mets, n_rxns] = size(S);

% indices of not contrained variables
ind_nic = setdiff(1:n_rxns, ind_ic);

% boolean vector for not constrained variables
sel_nic = zeros(n_rxns, 1);
sel_nic(ind_nic) = 1;
% boolean vector for integer variables
selRxns = ones(size(model.rxns));
sel_int = selRxns;
% bolean vector for reactions in can set and not in must set and not in
% constrained set of reactions
sel_c_nm_nc = can & ~must & sel_nic;

% Number of integer variables
n_int = sum(sel_int);
% Number of inner  constraints
n_ic = length(ind_ic);
% Number of inner variables not constrained
n_nic = length(ind_nic);

% Iic
Iic=selMatrix(sel_ic);
% Inic
Inic=selMatrix(sel_nic);

% Set variable types
vartype_bl(1: 7 * n_rxns + n_int + n_mets + 1) = 'C';
vartype_bl(n_rxns + 1:n_rxns + n_int) = 'B';

H = 1000;
bigM = 1000;

%   v(j)      y(j)      mu(j)     w(j)    deltam(j) deltap(j)  theta(j) thetap(j) labmda(i)    z
%|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|
%   n         n_int      n         n         n         n          n         n         m        1

% Set upper/lower bounds
lb_bl = [lb; zeros(6 * n_rxns + n_int + n_mets + 1, 1)]; %v(j)
ub_bl = [ub; H * ones(6 * n_rxns + n_int + n_mets + 1, 1)]; %v(j)
lb_bl(n_rxns + n_int + 1:n_rxns + n_int + n_rxns) = -1000; %mu(j)
lb_bl(2 * n_rxns + n_int + 1:2 * n_rxns + n_int + n_rxns) = -1000; %w(j)
lb_bl(7 * n_rxns + n_int + 1:7 * n_rxns + n_int + n_mets) = -1000; %lambda(i)
lb_bl(7 * n_rxns + n_int + n_mets + 1) = -1000; %z

%PRIMAL PROBLEM
%1) primal1 (n_mets equations)
%   S*v=0
A_p = [S zeros(n_mets, n_rxns * 6 + n_int + n_mets + 1)];
b_p = zeros(n_mets,1);
csense_p(1:n_mets) = 'E';

%2) primal 2, 3 and 7 (n_ic equations)
%   v_ic = b_ic
if n_ic > 0
    A_p = [A_p; Iic zeros(n_ic, n_rxns * 6 + n_int + n_mets + 1)];
    b_p = [b_p; b_ic'];
    csense_p(end + 1:end + n_ic) = 'E';
end

%3) primal 5 (n_ic equations)
%   -v(j) >= -ub(j)
A_p = [A_p; -Inic zeros(n_nic, n_rxns * 6 + n_int + n_mets + 1)];
b_p = [b_p; -ub(ind_nic)];
csense_p(end + 1:end + n_nic) = 'G';
%
%4) primal 6 (n_ic equations)
%   v(j) >= lb(j)
A_p = [A_p; Inic zeros(n_nic, n_rxns * 6 + n_int + n_mets + 1)];
b_p = [b_p; lb(ind_nic)];
csense_p(end + 1:end + n_nic) = 'G';

%DUAL PROBLEM
%1) dualcon1 (n_nic equations)
%   sum_i(lambda(i) * S(i,j))  +  deltam(j) -deltap(j) - y(j) =0
A_d = [];
b_d = [];
A_d = [A_d; zeros(n_nic, n_rxns) -Inic zeros(n_nic, n_rxns) zeros(n_nic, n_rxns) Inic -Inic zeros(n_nic, n_rxns) zeros(n_nic, n_rxns) S(:, ind_nic)' zeros(n_nic, 1)];
b_d = [b_d; zeros(n_nic, 1)];
csense_d(1:n_nic) = 'E';

%2) dualcon2 (n_ic equations)
%   sum_i(lambda(i) * S(i,j))  +  mu(j) - y(j) =0
A_d = [A_d; zeros(n_ic, n_rxns) -Iic Iic zeros(n_ic, n_rxns) zeros(n_ic, n_rxns) zeros(n_ic, n_rxns) zeros(n_ic, n_rxns) zeros(n_ic, n_rxns) S(:, ind_ic)' zeros(n_ic, 1)];
b_d = [b_d; zeros(n_ic, 1)];
csense_d(end + 1:end + n_ic) = 'E';

%OUTER PROBLEM
% bilevel_obj_up (1 equation)
%z=sum(w(j)-basemax(j) * y(j)) -> z -sum(w(j)) + sum(basemax(j) * y(j)) = 0 for all j in can y not in must and not in
%contraint_flux
A_bl = [zeros(1, n_rxns) (maxFluxesW.*sel_c_nm_nc)' zeros(1, n_rxns) -sel_c_nm_nc' zeros(1, 4 * n_rxns + n_mets) 1];
b_bl = 0;
csense_bl = 'E';

%primal_dual_up (1 equation)
% sum(w(j))  +  sum(b(j) * mu(j))  sum(deltap(j) * UB(j) - deltam(j) * LB(j)) = 0
A_bl = [A_bl; zeros(1, n_rxns + n_int) -sel_ic_b'  ones(1, n_rxns) -lb'.* sel_nic' ub'.*sel_nic' zeros(1, 2 * n_rxns + n_mets + 1)];
b_bl = [b_bl; 0];
csense_bl = [csense_bl, 'E'];

% bilevelcon0_up (1 equation)
%sum(w(j)-basemax(j) * y(j))> = 0
A_bl = [A_bl; zeros(1, n_rxns) (-maxFluxesW.*sel_c_nm_nc)' zeros(1, n_rxns) sel_c_nm_nc' zeros(1, 4 * n_rxns + n_mets + 1)];
b_bl = [b_bl; 0];
csense_bl = [csense_bl, 'G'];

% bilevelcon1 (j equations)
%w(j) - bigM*y(j) <= 0
A_bl = [A_bl; zeros(n_rxns, n_rxns) -bigM * speye(n_int) zeros(n_rxns, n_rxns) speye(n_rxns) zeros(n_rxns, 4 * n_rxns + n_mets + 1) ];
b_bl = [b_bl;zeros(n_rxns, 1)];
csense_bl(end + 1:end + n_rxns) = 'L';

% bilevelcon2 (j equations)
%w(j) + bigM * y(j) >= 0
A_bl = [A_bl; zeros(n_rxns, n_rxns) bigM * speye(n_rxns) zeros(n_rxns, n_rxns) speye(n_rxns) zeros(n_rxns, 4 * n_rxns + n_mets + 1) ];
b_bl = [b_bl;zeros(n_rxns, 1)];
csense_bl(end + 1:end + n_rxns) = 'G';

% bilevelcon3 (j equations)
%w(j) <= v(j)  +  bigM * (1-y(j))   ->    w(j) - v(j)  +  bigM * y(j)) <= bigM
A_bl = [A_bl; -speye(n_rxns) bigM * speye(n_rxns) zeros(n_rxns, n_rxns) speye(n_rxns) zeros(n_rxns, 4 * n_rxns + n_mets + 1) ];
b_bl = [b_bl;bigM * ones(n_rxns, 1)];
csense_bl(end + 1:end + n_rxns)  =  'L';

% bilevelcon4 (j equations)
%w(j) >= v(j) - bigM * (1-y(j))   ->    w(j) - v(j) - bigM * y(j)) >= -bigM
A_bl = [A_bl; -speye(n_rxns) -bigM * speye(n_rxns) zeros(n_rxns, n_rxns) speye(n_rxns) zeros(n_rxns, 4 * n_rxns + n_mets + 1) ];
b_bl = [b_bl; -bigM * ones(n_rxns, 1)];
csense_bl(end + 1:end + n_rxns)  =  'G';

% bilevelcon5 (1 equation)
%z >= 0.5
A_bl = [A_bl; zeros(1, 7 * n_rxns + n_int + n_mets) 1];
b_bl = [b_bl; 0.5];
csense_bl(end + 1)  =  'G';

%must_bin (1 equation)
% sum(y(j))=1
A_bl = [A_bl; zeros(1, n_rxns) sel_c_nm_nc' zeros(1, 6 * n_rxns + n_mets + 1) ];
b_bl = [b_bl; 1];
csense_bl = [csense_bl, 'E'];

%blocked_bin (1 equation)
% sum(y(j))=1
A_bl = [A_bl; zeros(1, n_rxns) ones(1, n_int) zeros(1, 6 * n_rxns + n_mets + 1)];
b_bl = [b_bl; 1];
csense_bl = [csense_bl, 'E'];

%Build bilevel matrices and vectors
A_bl_up = [A_bl; A_d; A_p];
b_bl_up = [b_bl; b_d; b_p];
csense_bl_up = [csense_bl, csense_d, csense_p];
c_bl_up = zeros(7 * n_rxns + n_int + n_mets + 1, 1); c_bl_up(end) = 1;

% Helper arrays for extracting solutions
sel_cont_sol = 1:n_rxns;
sel_int_sol = n_rxns + 1:n_rxns + n_int;

% Construct problem structure
bilevelMILPproblem.A = A_bl_up;
bilevelMILPproblem.b = b_bl_up;
bilevelMILPproblem.c = c_bl_up;
bilevelMILPproblem.csense = csense_bl_up;
bilevelMILPproblem.lb = lb_bl;
bilevelMILPproblem.ub = ub_bl;
bilevelMILPproblem.vartype = vartype_bl;
bilevelMILPproblem.contSolInd = sel_cont_sol;
bilevelMILPproblem.intSolInd = sel_int_sol;

% Initialize initial solution x0
bilevelMILPproblem.x0 = [];

% Maximize
bilevelMILPproblem.osense = -1;

% Set model for MILP problem
bilevelMILPproblem.model = model;

end

author: opencobra

repo: cobratoolbox


text:
function R = SigmoidResponse(img, sr_n, sr_sigma, sr_B)
%
%       R = SigmoidResponse(img, sr_n, sr_sigma, sr_B)
%
%       This function computes sigmoid response
%
%       input:
%           -img: an image
%           -sr_n: power
%           -sr_sigma: saturation parameter
%           -sr_B:
%
%       output:
%           -R: is the response
%
%     Copyright (C) 2011-14  Francesco Banterle
%
%     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/>.
%

if(~exist('sr_n','var'))
    sr_n = 0.73;
end

if(~exist('sr_sigma','var'))
    sr_sigma = 1.0;
end

if(~exist('sr_B','var'))
    sr_B = 1.0;
end

img_n = img.^sr_n;
R = img_n ./ (img_n + sr_sigma.^sr_n);
R = R * sr_B;

end

author: banterle


text:
function [ imgList ] = parseImg( loc )
% PARSEIMG
% parse image paths from given location
%
% Hyeonseob Nam, 2015
%

% image extension :
ext = {'jpg', 'png'};

% parse image
tmpList = {};
for i=1:length(ext)
   extList = dir(fullfile(loc, ['/img/*.', ext{i}]));
   tmpList = {tmpList{:}, extList(:).name};
end

% put prefix path to imgList
for i=1:length(tmpList)
   tmpList{i} = fullfile(loc, '/img/', tmpList{i});
end

imgList = tmpList;

end

repo: CREST-Release


text:
function [] = test()
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
randn('state',1212412414424234324);
rand('state',1212412414424234324);
dim_ambient = 1000;
dim_r = 10;
M = randn(dim_ambient,dim_r);
N = randn(dim_r,dim_ambient);
D0 = M*N/dim_r;

E0 = sign(randn(dim_ambient,dim_ambient));
inds = rand(dim_ambient)<0.7;
E0(inds) = 0;

D = D0 + E0;

D_hat = as_rpca(D,0.05,5,1.3);

error = max(max(abs(D0 - D_hat)))./max(max(abs(D0)));
disp(['recover error=' num2str(error)]);

D_hat = as_rpca(D,0.05,10,1.3);

error = max(max(abs(D0 - D_hat)))./max(max(abs(D0)));
disp(['recover error=' num2str(error)]);

D_hat = as_rpca(D,0.05,15,1.3);

error = max(max(abs(D0 - D_hat)))./max(max(abs(D0)));
disp(['recover error=' num2str(error)]);
end

author: andrewssobral

repo: lrslibrary


text:
function [flag,feedBack]=nonBlockingCheck_WithFeedback(t_Kuka, param)
% This functions is used with nonBlocking motion functions to check whether
% the non-blcokking motion is finished execution, in addition this function
% returns a feedback about the state of the robot.

%% Syntax:
% [flag]=nonBlockingCheck_WithFeedback(t_Kuka, param)

%% Arguments
% t_Kuka: is the TCP/IP connection object
% param: an integer value used to specify the feedback required from the
% robot (EEF force, EEF moment, joints torques, joint angles), please check
% the "feedBack" variable in the "return value" section-below for more info
% about the supported values of (param).

%% return value
% flag is a numerical value:
% flag=1: signifies the goal is reached
% flag=0: signifies the goal is not reached
% flag=-1: signifies an error
% feedBack: a cell array representing the feedback from the robot, the
% dimention of the array and the information in its cells depend
% on the value of the "param" argument as in the following table:
%-------------------------------------------------
% param | feedback
%   1   | force at EEF described in frame of EEF
%   2   | external torques at the joints
%   3   | measured torques at the joints
%   4   | joints positions from internal encoders
%--------------------------------------------------
% Copyright: Mohammad SAFEEA 12th-May-2019

global paramName;
global paramVal;
paramArg_CarLevel={'DcSeCarEEfFrelEEF_','DcSeCarExT_','DcSeCarMT_','DcSeCarJP_'};
sleeping_time=0.005;
x='eefpos';
%% if motion at EEF level
if(strcmp(x,paramName))
    % check size of array
    if(max(size(paramVal))==6)
    else
        flag=-1;
        return;
    end
% check if parameters are equal
    paramArg=paramArg_CarLevel{param};
    [feedBack,tempFlag]=GetFeedback( t_Kuka,paramArg );
    if tempFlag==false
        flag=1;
        pause(sleeping_time);
    else
        flag=0;
        pause(sleeping_time);
    end
    return;
end
%% if motion at joint space level
x='jpos';
if(strcmp(x,paramName))
    % check size of array
    if(max(size(paramVal))==7)
    else
        flag=-1;
        return;
    end
    % check if parameters are equal
    paramArg=paramArg_CarLevel{param};
    [feedBack,tempFlag]=GetFeedback( t_Kuka,paramArg );
    if tempFlag==false
        flag=1;
        pause(sleeping_time);
    else
        flag=0;
        pause(sleeping_time);
    end
    return;
end
% otherwise
flag=-1;
end

function [ feeadback,flag ] = GetFeedback( t_Kuka, paramArg )
EEEFpos={0,0,0,0,0,0};
theCommand=paramArg;

for i=1:6
    x=EEEFpos{i};
    theCommand=[theCommand,sprintf('%0.5f',x),'_'];
end
% send the message through network
fprintf(t_Kuka, theCommand);
message=fgets(t_Kuka);
[feeadback,N]=getDoubleFromString(message);
if(N==0) % when motion is finished, server returns "done+stopCharacter", raising this case
    flag=false; % segnifies motion is done
    message=fgets(t_Kuka);
    [feeadback,N]=getDoubleFromString(message);
    pause(0.1); % wait for sometime
else
    flag=true;
end
end

function [jPos,j]=getDoubleFromString(message)
n=max(max(size(message)));
j=0;
numString=[];
for i=1:n
    if message(i)=='_'
        j=j+1;
        jPos{j}=str2num(numString);
        numString=[];
    else
        numString=[numString,message(i)];
    end
end
if j==0
    jPos={0};
end
end

repo: KST-Kuka-Sunrise-Toolbox


text:
%% cleanDir
% Below is a demonstration of the features of the |cleanDir| function

%%
clear; close all; clc;

%% Syntax
% |cleanDir(pathName,extCell);|

%% Description
% This function can clean a folder by removing either all files in a folder
% or all files with the extensions specified in the cell array extCell.
%
% E.g. if one is interested in removing all .txt and .xml files from the
% folder specified by pathName on could use:
% pathName='D:\MATLAB\GIBBON'
% extCell={'txt','xml'}; %Extensions of files to delete
% cleanDir(pathName,extCell);

%% Examples
%

%%
% Create an example folder containing files with 3 different extensions,
% i.e. txt, csv and doc.

% Create temporary example folder
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
pathName=fullfile(defaultFolder,'data','temp','cleanCheck');
mkdir(pathName)

% Create files
n=3; %Number of files for each type
for q=1:1:n %Create n files
    fileID=fopen(fullfile(pathName,['temp',num2str(q),'.txt']),'w');
    fprintf(fileID,'%d\n',pi);
    fclose(fileID);
end

for q=1:1:n %Create n files
    fileID=fopen(fullfile(pathName,['temp',num2str(q),'.csv']),'w');
    fprintf(fileID,'%d\n',pi);
    fclose(fileID);
end

for q=1:1:n %Create n files
    fileID=fopen(fullfile(pathName,['temp',num2str(q),'.xml']),'w');
    fprintf(fileID,'%d\n',pi);
    fclose(fileID);
end

%%
% Show current folder content
disp('Old folder content:')
ls(pathName)

%%
% Clean up directory by removing all files whose extensions are a member of
% the extension set provided.

extCell={'csv','txt'}; %File extensions for files to remove
cleanDir(pathName,extCell)

%%
% Show current folder content
disp('New folder content:')
ls(pathName)

%%
% remove the temporary folder created for this example
rmdir(pathName,'s')

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%
% This code calculates the interevent distances and the correlation integral
% of a given earthquake distribution.
% Calculation of the 3D distances between all possible pairs
% (combination of n epicenters taken 2 at a time) of earthquakes of
% the given dataset.
% Francesco Pacchiani 1/2000
%
%
% Variables
%
N = size(E,1);				% N= # of events in the catalogue; E= Earthquake catalogue
pairdist = []; 			% pairdist= Vector of interevent distances
j = nchoosek(N,2);			% j= # of interevent distances calculated
pairdist = zeros((N-1),N);
depth = zeros((N-1),N);
k = 1;
%E.Latitude= (max(Da(:,2))+min(Da(:,2)))/2;

%Ho_Wb = waitbar(0,'Calculating the fractal dimension D');
%Hf_Cfig = gcf;
%Hf_child = get(groot,'children');
%set(Hf_child,'pointer','watch','papertype','A4');
%
%
% Calculation of the interevent distances in 2D plus the depths differences.
%
%
for i = 1:(N-1)

    lon1 = repmat(E(i,1), [(N-1),1]);
    lat1 = repmat(E(i,2), [(N-1),1]);
    depth1 = repmat(E(i,7), [(N-1),1]);

    E(i,:) = [];

    lon2 = E.Longitude;
    lat2 = E.Latitude;
    depth2 = E.Depth;

    pairdist(1:(N-1),k) = distance(lat1,lon1,lat2,lon2);
    depth(1:(N-1),k) = depth1-depth2;

    k = k + 1;
    E = newt2;

    %Waitbar((0.75/(N-1))*i, Ho_Wb);

end

clear i j k;
%
%
% Conversion of the interevent distances from degrees to kilometers and calculates
% the interevent distances in three dimensions.
%
%
if dtokm == 1
    pairdist = pairdist.*111;
end

pairdist = (pairdist.^2 + depth.^2).^0.5;		% pairdist = Interevent distances (vector).
clear depth;
%
%
% Calculation of the correlation integral using as input the
% pair distances computed above.
%
%
% Variables
%
d = 3;						%d = the dimension of the embedding volume.
rmax1 = max(pairdist);
rmin1 = min(pairdist);
rmax = max(rmax1);
rmin = min(rmin1);

clear rmax1 rmin1;

if rmin == 0
    rmin = 0.01;
end

lrmin = log10(rmin);
lrmax = log10(rmax);

%u = (log10(rmin):0.15:log10(rmax))';
%
% Defining the distance vector r in order that on the
% log-log graph all the points plot at equal distances from one another.
%
r = (logspace(lrmin, lrmax, 50))';
%r = zeros(size(u,1),1);
%r = 10.^u;
%
%
corint = [];						% corint= Vector of ?cumulative? correlation integral values for increasing interevent radius
cor = zeros(size(r,1),N);
i = 1;

for k = 1:50

    l = [];
    l = pairdist < r(k);
    test(k,1) = sum(l);
    %cor(1:size(r,1),i) = (sum(l));
    %i=i+1;

    %Waitbar(0.75 + (0.25/size(r,1))*k,Ho_Wb);
end

sumcor = sum(cor,2);
sumcor1 = sumcor./N;

for j = 0:20

    corint1(1:50,1,(j+1)) = sumcor1(:,1,(j+1)).^(1/(j-1));
end

corint = corint1(1:50, 1:21);

clear i j j1 k l cor sumcor;
close(Ho_Wb);
Hf_child = get(groot,'children');
set(Hf_child,'pointer','arrow');
%
%
% Plotting of the correlation integral in function of the interevent
% distance r.
%
%
Hf_Fig = figure_w_normalized_uicontrolunits('Numbertitle','off','Name','Fractal Dimension');
Hl_gr1 = loglog(r, corint,'ko', 'MarkerSize',5);
set(Hl_gr1,'MarkerSize',5);
title(sprintf('Correlation Integral versus %.0fD Interevent Distances', d));

Hf_Fig = figure_w_normalized_uicontrolunits('Numbertitle','off','Name','Fractal Dimension');
Hl_gr1 = loglog(r, cor,'r+');
set(Hl_gr1,'MarkerSize',7);
title(sprintf('Correlation Integral versus %.0fD Interevent Distances', d));
%
%
dofd = 'fd';
dofdim;

author: CelsoReyes


text:
function varargout = interp1(varargin)
%INTERP1 (overloaded)

switch class(varargin{3})

    case 'sdpvar'

        if nargin < 4
            varargin{4} = 'linear';
        end

        if nargin < 5
            varargin{5} = [];
        end

        if isa(varargin{2},'function_handle')
            varargin{2} = varargin{2}(varargin{1});
        end

        if numel(varargin{3}) > 1
            d = size(varargin{3});

            if min(d) > 1
                % Matrix case, vectorize
                varargout{1} = interp1(varargin{1},varargin{2},reshape(varargin{3},[],1),varargin{4},varargin{5});
                varargout{1} = reshape(varargout{1},d);
                return
            end

            if any(any(isnan(varargin{1})))
                error('Interpolation grid contains NaNs');
            end
            if any(any(isnan(varargin{2})))
                error('Interpolation data contains NaNs');
            end

            % If only one data sequence, place in row vector
            if min(size(varargin{1})) == 1
                varargin{1} = repmat(varargin{1}(:)',max(d),1);
            end
            if min(size(varargin{2})) == 1
                varargin{2} = repmat(varargin{2}(:)',max(d),1);
            end

            out = [];
            varargin{3} = reshape(varargin{3},1,max(d));
            for i = 1:max(d)
                    Y.type = '()';
                    Y.subs{1} = 1;
                    Y.subs{2} = i;
                    xij = subsref(varargin{3},Y);
                    out = [out;interp1(varargin{1}(i,:),varargin{2}(i,:),xij,varargin{4},varargin{5})];
            end
            if d(1)==1
                out = out';
            end
            varargout{1} = out;
            return
        else
            if min(size(varargin{1})) > 1 || min(size(varargin{2})) > 1
                error('Data should be vectors for scalar interpolant');
            end
            % Column vector assumed format
            varargin{1} = reshape(varargin{1},[],1);
            varargin{2} = reshape(varargin{2},[],1);
        end

        if ~isa(varargin{1},'double') || ~isa(varargin{2},'double')
            error('First 2 arguments in interp1 should be approximation data');
        end

        %if any(diff(varargin{1})<0)
        %    error('First arguments in interp1 should be monotonically increasing');
        %end

        if isequal(varargin{4},'graph')
            if  ~isconvexdata(varargin{1},varargin{2}) && ~isconvexdata(varargin{1},-varargin{2})
                error('Data has to be convex or concave for graph approximant');
            end
        end

        % Reorder arguments to make sure sdpvar is first argument.
        % Use local version of interp to deal with this
        % Arguments xi yi x method -> x xi yi method
        varargout{1} = yalmip('define','interp1_internal',varargin{[3 1 2 4 5]});

    case 'double'
        if isa(varargin{1},'double') && isa(varargin{2},'sdpvar')
            % User is trying to do nonlinear programming where the knot
            % function values in an interpolant are free variables
            zi = varargin{3};
            [n,m] = size(zi);
            if n*m == 1
                varargout{1} = yalmip('define','interp1_nonlinear',varargin{1},varargin{2},varargin{3},varargin{4});
            else
                out = [];
                zi = reshape(zi,1,n*m);
                for i = 1:n*m
                    Y.type = '()';
                    Y.subs{1} = 1;
                    Y.subs{2} = i;
                    zij = subsref(zi,Y);
                    out = [out  yalmip('define','interp1_nonlinear',varargin{1},varargin{2},zij,varargin{4})];
                end
                varargout{1} = reshape(out,n,m);
            end

        elseif isa(varargin{1},'sdpvar') && isa(varargin{2},'sdpvar')
            % User is trying to do nonlinear programming where the knot
            % function values and locations in an interpolant are free variables

            zi = varargin{3};
            [n,m] = size(zi);
            if n*m == 1
                varargout{1} = yalmip('define','interp1_nonlinear',[varargin{1},varargin{2}],[],varargin{3},varargin{4});
            else
                out = [];
                zi = reshape(zi,1,n*m);
                for i = 1:n*m
                    Y.type = '()';
                    Y.subs{1} = 1;
                    Y.subs{2} = i;
                    zij = subsref(zi,Y);
                    out = [out   yalmip('define','interp1_nonlinear',[varargin{1},varargin{2}],[],zij,varargin{4})];
                end
                varargout{1} = reshape(out,n,m);
            end
        else
            error('SDPVAR/INTERP1 called with strange arguments!');
        end

    otherwise
        error('SDPVAR/INTERP1 called with strange argument!');
end

function isconvex = isconvexdata(xi,yi)
finterp = (yi(1:end-2) + yi(3:end))/2;
if all(finterp >= yi(2:end-1))
    isconvex = 1;
else
    isconvex = 0;
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function rtf(plot_frq,flag_save)

close all
evalin('base','stop=0;');

%=========SET THE BASIC FIGURE=================
fig = figure('Position',[500,500,800,600],...
        'NumberTitle','off',...
        'Name','Scope',...
        'doublebuffer','on',...
        'HandleVisibility','on',...
        'KeyPressFcn', @keypress, ...
        'Renderer', 'openGL');
%=============================================

%=============OPEN THE DEVICE FOR RECORD======
sample_frequency = 44100;
samples_per_frame = 1024;
%plot_frq=10;
record_time=600;
samples_to_acquire = record_time * sample_frequency;

%PREPARE    THE    DEVICE
ai = analoginput('winsound');
chan = addchannel( ai, 1 );
set( ai, 'SampleRate', sample_frequency )
set( ai, 'SamplesPerTrigger', samples_to_acquire )
set(ai, 'StopFcn', @stop_dev)

sample_frequency = get( ai, 'SampleRate' );

%SETTING CALL BACK FUNCTIONS:
%The first for capture   the
%second     for      display
set(ai, 'SamplesAcquiredFcnCount',samples_per_frame);
set(ai, 'SamplesAcquiredFcn',@flag);
set(ai, 'TimerPeriod',(1/plot_frq));
set(ai, 'TimerFcn',@disply);

%=============SAVE THE CONFIGURATION======
plot_ref=plot(zeros(10,1));
fid=-1;

%SAVE THE CURRENT PARAMETERS:
name_of_file=sprintf('%s-%d','real-anal',(round(sample_frequency/samples_per_frame)));
remark={1,...
        zeros(samples_per_frame*20,1)',...
        0,...
        plot_ref,...
        plot_frq,...
        cputime,...
        flag_save,...
        -1,...
        name_of_file
        };

set(ai, 'UserData',remark)

%=============START TO RECORD================
fprintf ('To stop the program set <stop=1> or press q in the figure window\n');
start (ai)

%=============================================
%=========THE MAIN PROGRAM====================
%=============================================
%                   *
%                  ***
%                 *****
%                  ***
%                  ***
%                  ***
%                  ***
%                  ***
%                 *****
%                  ***
%                   *
%=============================================
%==========CALLBACK FUNCTIONS=================
%=============================================

%=========Keypress callback===========
function keypress(src, e)
  keypressed=get(gcf,'CurrentCharacter');

  % ignore raw control, shift, alt keys
  if keypressed
    % Quit
    if strcmp( keypressed, 'q')
        evalin('base','stop=1;');
    end
  end
return

%============FLAG FUNCTION===================
%This function activated when we capture
%certain    amount  of          samples
function flag(obj,event)

  % CHECK FOR STOP SIGNAL
  if  evalin('base','stop')
      stop(obj)
  end

% GET THE OLD DATA
 remark=get(obj,'UserData');
 flag_write=remark{1};  %Do I have to
 buffer=remark{3};      %What is the old picture
 flag_save=remark{7};   %Are we in saving mode?
 fid=remark{8};         %What file descriptor to save

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %IN CASE - DELETE/SAVE THE OLD DATA

 if flag_write>20

     % IN CASE WE HAVE TO SAVE - CLOSE THE OLD FILE AND MAKE A NEW
     if flag_save>0
        fclose(fid);
        name_of_data=sprintf('%s-%d.dat','dat',(round(cputime*1000)));
        fid=fopen(name_of_data,'w');
    end

      %DELETE OLD DATA
     flag_write=1;
     buffer=[];
     remark{1}=flag_write; % SET THE POSITION OF THE READING SHIFT
 end

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % TAKE THE NEW DATA

 samples_per_frame=get(obj,'SamplesAcquiredFcnCount');
 data=(getdata(obj,samples_per_frame))';

 % IN CASE - WRITE THE DATA
 if flag_save>0
     if fid==-1
         name_of_data=sprintf('%s-%d.dat','dat',(round(cputime*1000)));
         fid=fopen(name_of_data,'w');
     end
     fwrite(fid,(data*10000),'short');
     remark{8}=fid;
 end

 % Add to buffer
 buffer=[buffer data];
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 remark{3}=buffer;
 set(obj,'UserData',remark);

return

function stop_dev(obj,event)

     remark=get(obj,'UserData');
     if (remark{8}>0)                 %FID>0 == There is open file
         fclose (remark{8});
     end
     close all

     fprintf('\n\nThanks for using Erlich Real-Time scope\n');
     save (remark{9},'remark');

     delete(obj)
     clear obj
return

function disply(obj,event)

  sample_frequency=get(obj,'SampleRate');
  remark=get(obj,'UserData');
  refresh_frq=remark{5};
  read_shift=remark{1};

  ring=remark{2};
  buffer=remark{3};

  end_shift=min((read_shift+round(sample_frequency/refresh_frq)),length(buffer));
  new_data=buffer(read_shift:end_shift);
  ring=[ring new_data];
  ring(1:length(new_data))=[];
  remark{1}=end_shift;
  remark{2}=ring;

  start_display(ring,remark{4})

  set(obj,'UserData',remark);
return

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function lagplot(c)
   %lagplot

   % Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   % PREP PLOT
   figure('Color','w','Position',[50 50 600 500]);
   set(gcf,'DefaultAxesFontSize',14);
   imagesc(c.lags);
   title('Lag time for maximum correlation (s)');

   % ADD DATES TO AXES
   n = length(c.trig);
   ticvals = 1:round(n/25):n;
   set(gca,'XTick',ticvals);
   set(gca,'YTick',ticvals);
   yt = get(gca,'YTick');
   set(gca,'YTickLabel',datestr(c.trig(yt),'yyyy-mm-dd HH:MM'),'FontSize',6);

   % DRESS UP THE FIGURE
   cmap = load('colormap_lag.txt');
   colormap(cmap);
   colorbar;
   xlabel('Event number');
   ylabel('Event date');

   %PRINT OUT FIGURE
   set(gcf, 'paperorientation', 'portrait');
   set(gcf, 'paperposition', [1.25 2.5 6 6] );
   %print(gcf, '-depsc2', 'FIG_tartan.ps')
end

author: geoscience-community-codes

repo: GISMO


text:
function [stoi_PercCorr] = stoi_d2percCorr(stoi_d)
% Converts the stoi measure, d, to percent words correct unit
%
% Syntax:	[stoi_PercCorr] = stoi_d2percCorr(stoi_d)
%
% Inputs:
% 	stoi_d - The STOI measure in the raw units
%
% Outputs:
% 	stoi_PercCorr - The corressponding STOI value in Percent of Words
% 	Correct as defined for the IEEE English Library
%
%
%    References:
%      C. H. Taal, R. C. Hendriks, R. Heusdens, and J. Jensen. A Short-Time
%      Objective Intelligibility Measure for Time-Frequency Weighted Noisy
%      Speech. In Acoustics Speech and Signal Processing (ICASSP), pages
%      4214-4217. IEEE, 2010.
%
%      C. H. Taal, R. C. Hendriks, R. Heusdens, and J. Jensen. An Algorithm
%      for Intelligibility Prediction of Time-Frequency Weighted Noisy Speech.
%      IEEE Transactions on Audio, Speech and Language Processing,
%      19(7):2125-2136, 2011.

% Author: Jacob Donley
% University of Wollongong
% Email: jrd089@uowmail.edu.au
% Date: 05 August 2015
% Revision: 0.1
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f_IEEE_a = -17.4906;
f_IEEE_b = 9.6921;

stoi_PercCorr = 100 / (1 + exp(f_IEEE_a * stoi_d + f_IEEE_b ) );

end

author: jtkim-kaist

repo: Speech-enhancement


text:
function H= gridutil_addScale(mnt, varargin)
%GRIDUTIL_ADDSCALE - Add a Scale for X- and Y-Axis Units to a a Grid Plot
%
%Synopsis:
% H= gridutil_addScale(MNT, <OPT>)
%
%Input:
% MNT: struct for electrode montage, see setElectrodeMontage
% OPT: struct or property/value list of optional properties:
%  .XUnit - Unit of X-Axis, default: 'ms'
%  .YUnit - Unit of Y-Axis, default: '\muV'
%
% Usually called by grid_plot

props = {'XUnit',                   'ms',                   'CHAR';
         'YUnit',                   '\muV',                 'CHAR';
         'YDir',                    'normal',               'CHAR';
         'ScaleHPos',               'zeroleft',             'CHAR';
         'ScaleVPos',               'zeromiddle',           'CHAR';
         'ScaleFontSize',           get(gca, 'FontSize'),   'DOUBLE';
         'ScaleShowOrientation',    1                       'BOOL'};

if nargin==0,
  H= props; return
end
opt= opt_proplistToStruct(varargin{:});
opt= opt_setDefaults(opt, props);
opt_checkProplist(opt, props);

xLim= get(gca, 'xLim');
yLim= get(gca, 'yLim');
xt= get(gca, 'xTick');
xx= median(diff(xt));
yt= get(gca, 'yTick');
yy= median(diff(yt));

x= [mnt.box(1,:) mnt.scale_box(1)];
y= [mnt.box(2,:) mnt.scale_box(2)];
w= [mnt.box_sz(1,:) mnt.scale_box_sz(1)];
h= [mnt.box_sz(2,:) mnt.scale_box_sz(2)];
bs= 0.005;
siz= (1+2*bs)*[max(x+w) - min(x) max(y+h) - min(y)];
pos= [mnt.scale_box(1)-min(x) mnt.scale_box(2)-min(y) ...
      mnt.scale_box_sz(1) mnt.scale_box_sz(2)]./[siz siz];
pos= pos + [bs bs 0 0];
H.ax= axes('position', pos);
ud= struct('type','ERP plus', 'chan','scale');
set(H.ax, 'xLim',xLim, 'yLim',yLim, 'userData',ud);

if strncmp('zero', opt.ScaleHPos, 4),
  if xLim(1)<0-0.1*diff(xLim) & xLim(2)>xx,
    opt.ScaleHPos= 'Zero';
  else
    opt.ScaleHPos= opt.ScaleHPos(5:end);
    if isempty(opt.ScaleHPos), opt.ScaleHPos= 'middle'; end
  end
end
if strncmpi('zero', opt.ScaleVPos, 4),
  if yLim(1)<0-0.1*diff(yLim) & yLim(2)>yy+0.1*diff(yLim),
    opt.ScaleVPos= 'Zero';
  else
    opt.ScaleVPos= opt.ScaleVPos(5:end);
    if isempty(opt.ScaleVPos), opt.ScaleVPos= 'middle'; end
  end
end

if isfield(opt, 'CLim')
   set(H.ax, 'CLim', opt.CLim')
   if isfield(opt, 'colormap')
      colormap(opt.Colormap);
   end
   cb = Colorbar;
   opt.ScaleHPos = 'left';
   ap = get(H.ax, 'position');
   set(cb, 'position', [ap(1)+0.8*ap(3) ap(2)+0.1*ap(4) 0.1*ap(3) 0.8*ap(4)])
  set(get(cb,'yLabel'),'String', opt.zUnit)
end

switch(opt.ScaleHPos),
 case 'Zero',
  x0= 0;
 case 'left',
  x0= xLim(1) + 0.05*diff(xLim);
 case 'middle',
  x0= mean(xLim) - xx/2;
 case 'right',
  x0= xLim(2) - xx - 0.05*diff(xLim);
 otherwise,
  error('unimplemented choice for opt.ScaleHPos');
end
switch(opt.ScaleVPos),
 case 'Zero',
  y0= 0;
% case 'top',
%  y0= yLim(2),
 case 'middle',
  y0= mean(yLim) - yy/2;
% case 'bottom',
%  y0= yLim(1);
 otherwise,
  error('unimplemented choice for opt.ScaleVPos');
end

H.vline= line([x0 x0], [y0 y0+yy]);
set(H.vline, 'LineWidth',2, 'Color','k');
H.text_YUnit= text(x0, y0+yy/2, sprintf(' %g %s', yy, opt.YUnit));
set(H.text_YUnit, 'verticalAli','middle', 'FontSize',opt.ScaleFontSize);
if opt.ScaleShowOrientation,
  if strcmpi(opt.YDir, 'reverse'),
    sgn= '-';
  else
    sgn= '+';
  end
  H.text_sgn= text(x0, y0+yy, [' ' sgn]);
  set(H.text_sgn, 'verticalAli','bottom', 'FontSize',opt.ScaleFontSize);
end

H.hline= line([x0 x0+xx], [y0 y0]);
set(H.hline, 'LineWidth',2, 'Color','k');
yt= y0-0.03*diff(yLim);
H.text_XUnit= text(x0+xx/2, yt, ...
                   sprintf('%g %s', xx, opt.XUnit));
set(H.text_XUnit, 'horizontalAli','center', 'verticalAli','top', ...
                  'FontSize',opt.ScaleFontSize);
if x0==xLim(1), %% leftmost position
  set(H.text_XUnit, 'position',[x0 yt 0], 'horizontalAli','left');
end

axis off;

if nargout==0,
  clear H;
end

author: bbci


text:
function [result, M] = ft_warp_optim(input, target, method)

% FT_WARP_OPTIM determine intermediate positions using warping (deformation)
% the input cloud of points is warped to match the target.
% The strategy is to start with simpelest linear warp, followed by a more
% elaborate linear warp, which then is followed by the nonlinear warps up
% to the desired order.
%
% [result, M] = ft_warp_pnt(input, target, method)
%     input          contains the Nx3 measured 3D positions
%     target         contains the Nx3 template 3D positions
%     method         should be any of
%                     'rigidbody'
%                     'globalrescale'
%                     'traditional' (default)
%                     'nonlin1'
%                     'nonlin2'
%                     'nonlin3'
%                     'nonlin4'
%                     'nonlin5'
%
% The default warping method is a traditional linear warp with individual
% rescaling in each dimension. Optionally you can select a nonlinear warp
% of the 1st (affine) up to the 5th order.
%
% When available, this function will use the MATLAB optimization toolbox.
%
% See also FT_WARP_APPLY, FT_WARP_ERRROR

% Copyright (C) 2000-2013, Robert Oostenveld
%
% 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

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

% this can be used for printing detailled user feedback
fb = false;

if nargin<3
  method='traditional';
end

pos1 = input;
pos2 = target;

% The ft_warp_error function might be located in the private subdirectory of
% FieldTrip, i.e. only available to functions in the FieldTrip toolbox.
% The following line ensures that the function can also be found by the
% feval that is executed by the optimalization toolbox.
errorfun = str2func('ft_warp_error');

% determine whether the MATLAB Optimization toolbox is available and can be used
if ft_hastoolbox('optim')
  optimfun = @fminunc;
else
  optimfun = @fminsearch;
end

% set the options for the minimalisation routine
if isequal(optimfun,  @fminunc)
  options  = optimset('fminunc');
  options  = optimset(options, 'Display', 'off');
  options  = optimset(options, 'MaxIter', 1500);
  % options  = optimset(options, 'MaxFunEvals', '1000*numberOfVariables');
  options  = optimset(options, 'TolFun', 1e-4);
  options  = optimset(options, 'LargeScale', 'off');
elseif isequal(optimfun, @fminsearch)
  options  = optimset('fminsearch');
  options  = optimset(options, 'Display', 'off');
  options  = optimset(options, 'MaxIter', 4500);
else
  ft_warning('unknown optimization function "%s", using default parameters', func2str(optimfun));
end

if fb; fprintf('distance = %f\n', errorfun([0 0 0 0 0 0], pos1, pos2, 'rigidbody')); end

if isempty(method)
  ft_error('incorrect warping method specified');
end

% the warp is done in steps, starting simple and progressively getting more complex
level = find(strcmp(method, {
  'rigidbody'         % 1
  'globalrescale'     % 2
  'traditional'       % 3
  'nonlin1'           % 4
  'nonlin2'           % 5
  'nonlin3'           % 6
  'nonlin4'           % 7
  'nonlin5'           % 8
  }));

if isempty(level)
  ft_error('incorrect warping method specified');
end

if level>=1
  % do a rigid-body transformation (6 parameters)
  if fb; disp('rigidbody...'); end
  ri = [0 0 0 0 0 0];
  rf = optimfun(errorfun, ri, options, pos1, pos2, 'rigidbody');
  if fb; fprintf('distance = %f\n', errorfun(rf, pos1, pos2, 'rigidbody')); end
end

if level>=2
  % do a rigid-body + global rescaling transformation (7 parameters)
  if fb; disp('rigidbody + global rescaling...'); end
  gi = [rf 1];
  gf = optimfun(errorfun, gi, options, pos1, pos2, 'globalrescale');
  if fb; fprintf('distance = %f\n', errorfun(gf, pos1, pos2, 'globalrescale')); end
end

if level>=3
  % do a rigid-body + individual rescaling transformation (9 parameters)
  if fb; disp('rigidbody + individual rescaling...'); end
  ti = [gf gf(7) gf(7)];
  tf = optimfun(errorfun, ti, options, pos1, pos2, 'traditional');
  if fb; fprintf('distance = %f\n', errorfun(tf, pos1, pos2, 'traditional')); end
end

if level>=4
  % do a first order nonlinear transformation,
  if fb; disp('1st order nonlinear...'); end
  e1i = traditional(tf);
  e1i = [e1i(1:3,4) e1i(1:3,1:3)];  % reshuffle from homogenous into nonlinear
  e1f = optimfun(errorfun, e1i, options, pos1, pos2);
  if fb; fprintf('distance = %f\n', errorfun(e1f, pos1, pos2, 'nonlinear')); end
end

if level>=5
  % do a second order nonlinear transformation,
  if fb; disp('2nd order nonlinear...'); end
  e2i = [e1f zeros(3,6)];
  e2f = optimfun(errorfun, e2i, options, pos1, pos2);
  if fb; fprintf('distance = %f\n', errorfun(e2f, pos1, pos2, 'nonlinear')); end
end

if level>=6
  % do a third order nonlinear transformation,
  if fb; disp('3rd order nonlinear...'); end
  e3i = [e2f zeros(3,10)];
  e3f = optimfun(errorfun, e3i, options, pos1, pos2);
  if fb; fprintf('distance = %f\n', errorfun(e3f, pos1, pos2, 'nonlinear')); end
end

if level>=7
  % do a fourth order nonlinear transformation,
  if fb; disp('4th order nonlinear...'); end
  e4i = [e3f zeros(3,15)];
  e4f = optimfun(errorfun, e4i, options, pos1, pos2);
  if fb; fprintf('distance = %f\n', errorfun(e4f, pos1, pos2, 'nonlinear')); end
end

if level>=8
  % do a fifth order nonlinear transformation,
  if fb; disp('5th order nonlinear...'); end
  e5i = [e4f zeros(3,21)];
  e5f = optimfun(errorfun, e5i, options, pos1, pos2);
  if fb; fprintf('distance = %f\n', errorfun(e5f, pos1, pos2, 'nonlinear')); end
end

% return the estimated parameters of the highest level warp
% and compute the warped points
switch level
  case 1
    M = rf;
  case 2
    M = gf;
  case 3
    M = tf;
  case 4
    M = e1f;
  case 5
    M = e2f;
  case 6
    M = e3f;
  case 7
    M = e4f;
  case 8
    M = e5f;
end
result = ft_warp_apply(M, input, method);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function res = nsamples(this)
% Method for getting the number of samples per trial
% FORMAT res = nsamples(this)
% _______________________________________________________________________
% Copyright (C) 2008-2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: nsamples.m 5025 2012-10-31 14:44:13Z vladimir $

res = this.Nsamples;

author: spm


text:
function [count] = countall3graphlets(L)

% Count all 3-node subgraphs in an undirected graph
% without node labels and with unweighted edges
% Author: Nino Shervashidze - nino.shervashidze@tuebingen.mpg.de
% Copyright 2012 Nino Shervashidze
% Input: L - 1xn cell array - adjacency list
% Output: count - 1x4 vector of integers. count(i)= number of graphlets with
%                 3-i+1 edges (see 3graphlets.pdf)

n=length(L); % number of nodes
count=zeros(1,4);
w=[1/6, 1/4, 1/2];
for v1=1:n
  for v2=L{v1}
    cardinalities=card_inter(L{v1}, L{v2}, length(L{v1}), length(L{v2}));
    count(1)=count(1)+w(1)*cardinalities(3);
    count(2)=count(2)+w(2)*(cardinalities(1)+cardinalities(2)-2);
    count(3)=count(3)+w(3)*(n-sum(cardinalities));
  end
end
count(4)=n*(n-1)*(n-2)/6-sum(count(1:3));
end

function [n] = card_inter(o_set1, o_set2, l1, l2)
% Find the cardinality of the intersection of two ordered sets of lengths l1 and l2 respectively
% n(1)=o_set1\o_set2, n(2)=o_set2\o_set1, n(3)=(o_set2 inter o_set1)
n=zeros(1,3);
i=1; j=1;

while i<=l1 && j <=l2
  if o_set1(i)<o_set2(j) n(1)=n(1)+1; i=i+1;
  else
    if o_set1(i)>o_set2(j) n(2)=n(2)+1; j=j+1;
    else i=i+1; j=j+1; n(3)=n(3)+1;
    end
  end
end
n(1)=n(1)+l1-i+1;
n(2)=n(2)+l2-j+1;
end

author: muhanzhang

repo: DGCNN

save_path: github-repos/MATLAB/muhanzhang-DGCNN


text:
%  Compute the power f^k of the multivariate polynomials f
%
%    Syntax:  >> h = mvPolynPower(f,k)
%
%    Input:   f -- multivariate polynomials in coeff. matrices
%             k -- integer exponent of the power
%
%   Output:   multivariate polynomial  f^k


text:
function res = nchannels(this)
% returns number of channels
% FORMAT res = nchannels(this)
% _______________________________________________________________________
% Copyright (C) 2008-2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: nchannels.m 5025 2012-10-31 14:44:13Z vladimir $

if this.montage.Mind==0
    res = length(this.channels);
else
    res = length(this.montage.M(this.montage.Mind).channels);
end


text:
function [curImage,phImage] = mrThreshPhVol(samp,sampSize,volco,volph,volume,sagSize,numSlices,x,y,dataRange)
%
%
% PURPOSE:
%   Compute an image whose pixel values show the best phase angle of
%   the time series at each image point.
%
%
% AUTHOR:  Engel
%
%

% Variable Declarations
thr = [];			% Vector of 1s and 0s.  1 means co > thresh
				% 0 means co <= thresh.
global interpflag volslicut volslimin1 volslimax1;

if isempty(volco)
   disp ('Correlation data is not available.');
   return
end

if (interpflag)
	sinIm = mrExtractImgVol(sin(volph),sagSize,dataRange(2)-dataRange(1)+1,samp,dataRange);
	cosIm = mrExtractImgVol(cos(volph),sagSize,dataRange(2)-dataRange(1)+1,samp,dataRange);
	phImage = atan2(sinIm,cosIm);
else
	phImage = mrExtractImgVol(volph,sagSize,dataRange(2)-dataRange(1)+1,samp,dataRange);
end
phImage = phImage;

curImage = mrExtractImgVol(volume,sagSize,numSlices,samp);

co = mrExtractImgVol(volco,sagSize,dataRange(2)-dataRange(1)+1,samp,dataRange);

cutoff = get(volslicut,'value');
disp(['Cutoff = ',num2str(cutoff)]);
thr = co > cutoff;
curImage(thr) = -((phImage(thr)+pi)/(2*pi)*110);	%Scale negatives for colors

myShowImageVol(curImage,sampSize,max(curImage)*get(volslimin1,'value'),max(curImage)*get(volslimax1,'value'),x,y);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function R = rand_rotation(varargin)
% generate random rotation matrix

params = inputParser;
params.CaseSensitive = false;

params.addParameter('RotationBound',2*pi,...
    @(x) 0.0<=x && x<=2*pi);

params.parse(varargin{:});

RotationBound = params.Results.RotationBound;

angle = RotationBound*rand - RotationBound/2;
axis  = randn(3,1);
axis  = axis / norm(axis);
R     = axang2rotm([axis' angle]);
end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
% ----------------------------------------------------------------------------
% function hfssExportToFileFC(fid, Name, PTS, Solution, Freq, [Phase = 0])
%
% Description :
% -------------
% Evaluates the top stack element at a set of points specified in an
% external file (PTS) and exports the data to a file (Name).
%
% Parameters :
% ------------
% fid      - file identifier of the HFSS script file.
% Name     - name of the field where the stack element will be exported.
% PTS      - name of the PTS file containing the grid (with extension).
% Solution - solution being resolved.
% Freq     - frequency of the solution.
% Phase    - phase of the field in degrees.
%
% Note :
% ------
% 1) This function is bad documented in scripting.pdf. It only indicates
%    two arguments while in reallity it requires more.
% 2) A grid can be generated by the function hfssCreatePTS.
%
% Example :
% ---------
% hfssCreatePTS(fid, 'MyGrid', [0 0 0], [1 1 1], [.1 .1 .1], 'mm');
% hfssEnterQty(fid, 'E');
% hfssExportToFileTC(fid, 'EField', 'MyGrid.pts', 'Solution', Freq);
%
% ----------------------------------------------------------------------------

% ----------------------------------------------------------------------------
% CHANGELOG
%
% 06-Oct-2012: *Initial release.
% ----------------------------------------------------------------------------

% ----------------------------------------------------------------------------
% Written by Daniel R. Prado
% danysan@gmail.com / drprado@tsc.uniovi.es
% 06 October 2012
% ----------------------------------------------------------------------------
function hfssExportToFileFC(fid, Name, PTS, Solution, Freq, Phase)

% Arguments processor.
if (nargin < 5)
	error('Insufficient # of arguments !');
elseif (nargin < 6)
    Phase = [];
end

% Setup default arguments
if isempty(Phase)
    Phase = 0;
end

% Preamble
fprintf(fid, '\n');
fprintf(fid, 'Set oModule = oDesign.GetModule("FieldsReporter")\n');

% Command
fprintf(fid, 'oModule.ExportToFile "%s.fld", _\n', Name);
fprintf(fid, '"%s", _\n', PTS);
fprintf(fid, '"%s : LastAdaptive", _\n', Solution);
fprintf(fid, 'Array("Freq:=", "%fGHz", _\n', Freq);
fprintf(fid, '"Phase:=", "%fdeg"), true\n', Phase);

author: yuip

repo: hfss-api

save_path: github-repos/MATLAB/yuip-hfss-api


text:
function k = whiteKernDiagCompute(kern, x)

% WHITEKERNDIAGCOMPUTE Compute diagonal of WHITE kernel.
% FORMAT
% DESC computes the diagonal of the kernel matrix for the white noise kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : a vector containing the diagonal of the kernel matrix
% computed at the given points.
%
% SEEALSO : whiteKernParamInit, kernDiagCompute, kernCreate, whiteKernCompute
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% KERN

k = repmat(kern.variance, size(x, 1), 1);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function signal = flt_beamformer(varargin)
% Recovers activity from the given ROIs via a beamformer.
% function Signal = flt_laplace(Signal,ROIs,)
%
% In:
%   Signal : EEGLAB data set, either continuous or epoched
%
%   ROIs   : Regions of interest from which to recover signals.
%
%   FieldType : type of field to recover ('axial' or 'normal')
%
%   ReferenceType : type of referencing used before this filter ('nasion' or 'common_average')
%
%   OverrideOriginal : whether to override the original signal
%
% Out:
%   Signal : activity from the given ROIs
%
% Notes:
%   This function currently does not perform activity renormalization.
%
% Examples:
%   % recover activity from four regions
%   eeg = flt_beamformer(eeg,{'Precentral_L', 'Precentral_R', 'Frontal_Sup_L', 'Frontal_Sup_R'},'normal')
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2014-02-11

if ~exp_beginfun('filter') return; end

declare_properties('name','Beamformer', 'follows',{'flt_selchans','flt_repair_bursts','flt_clean_settings'}, 'cannot_follow','flt_ica', 'independent_channels',false, 'independent_trials',true);

arg_define(varargin, ...
    arg_norep({'signal','Signal'}), ...
    arg({'roi_labels','ROIs','ROILabels'},{}, {'Precentral_L', 'Precentral_R', 'Frontal_Sup_L', 'Frontal_Sup_R', 'Frontal_Sup_Orb_L', 'Frontal_Sup_Orb_R', 'Frontal_Mid_L', 'Frontal_Mid_R', 'Frontal_Mid_Orb_L', 'Frontal_Mid_Orb_R', 'Frontal_Inf_Oper_L', 'Frontal_Inf_Oper_R', 'Frontal_Inf_Tri_L', 'Frontal_Inf_Tri_R', 'Frontal_Inf_Orb_L', 'Frontal_Inf_Orb_R', 'Rolandic_Oper_L', 'Rolandic_Oper_R', 'Supp_Motor_Area_L', 'Supp_Motor_Area_R', 'Frontal_Sup_Medial_L', 'Frontal_Sup_Medial_R', 'Frontal_Med_Orb_L', 'Frontal_Med_Orb_R', 'Insula_L', 'Insula_R', 'Cingulum_Ant_L', 'Cingulum_Ant_R', 'Cingulum_Mid_L', 'Cingulum_Mid_R', 'Cingulum_Post_L', 'Cingulum_Post_R', 'Hippocampus_L', 'Hippocampus_R', 'ParaHippocampal_L', 'ParaHippocampal_R', 'Calcarine_L', 'Calcarine_R', 'Cuneus_L', 'Cuneus_R', 'Lingual_L', 'Lingual_R', 'Occipital_Sup_L', 'Occipital_Sup_R', 'Occipital_Mid_L', 'Occipital_Mid_R', 'Occipital_Inf_L', 'Occipital_Inf_R', 'Fusiform_L', 'Fusiform_R', 'Postcentral_L', 'Postcentral_R', 'Parietal_Sup_L', 'Parietal_Sup_R', 'Parietal_Inf_L', 'Parietal_Inf_R', 'SupraMarginal_L', 'SupraMarginal_R', 'Angular_L', 'Angular_R', 'Precuneus_L', 'Precuneus_R', 'Paracentral_Lobule_L', 'Paracentral_Lobule_R', 'Temporal_Sup_L', 'Temporal_Sup_R', 'Temporal_Pole_Sup_L', 'Temporal_Pole_Sup_R', 'Temporal_Mid_L', 'Temporal_Mid_R', 'Temporal_Pole_Mid_L', 'Temporal_Pole_Mid_R', 'Temporal_Inf_L', 'Temporal_Inf_R', 'Olfactory_L', 'Olfactory_R', 'Rectus_L', 'Rectus_R', 'Amygdala_L', 'Amygdala_R', 'Caudate_L', 'Caudate_R', 'Thalamus_L', 'Thalamus_R', 'Heschl_L', 'Heschl_R'}, ...
        'Cortical anchor locations. List of locations to which components shall be constrained. The first k components are encouraged to lie close to the given locations, in the order of appearance. This is experimental and currently requires a) 10-20 locations and b) Guido Nolte''s source analysis toolbox (not included).','experimental',true), ...
    arg({'field_type','FieldType'},'normal',{'normal','axial'},'Regions of interest. These are the regions from which to recover signals.'), ...
    arg({'reference_type','ReferenceType'},'common_average',{'nasion','common_average'},'Referencing scheme. This is the type of re-referencing that was applied before flt_beamformer.'), ...
    arg({'override_original','OverrideOriginal'},true,[],'Override original data. If checked, the original signals will be replaced by the recovery.'), ...
    arg({'override_chanlabels','OverrideChanlabels'},{},[],'Override channel labels.'), ...
    arg({'cov_shrinkage','CovShrinkage'},0,[],'Covariance shrinkage parameter. For better conditioning.'), ...
    arg({'cov_robust','CovRobust'},false,[],'Use robust covariance estimate.'), ...
    arg_norep('M'), ...
    arg_norep('channel_mask'));

% calculate spatial filter matrix, if necessary
if ~exist('M','var')
    if isempty(override_chanlabels)
        chan_labels = {signal.chanlocs.labels};
    else
        chan_labels = override_chanlabels;
    end
    [dummy,dummy,channel_mask,dummy,dummy,dummy,normField,leadField] = hlp_diskcache('filterdesign',@calc_beamformer_constraints,chan_labels,roi_labels,eye(signal.nbchan),reference_type); %#ok<ASGLU>
    if strcmp(field_type,'normal')
        LF = normField;
    elseif strcmp(field_type,'axial')
        LF = leadField(:,:);
    else
        error('Unsupported field type requested: %s',hlp_tostring(field_type,1000));
    end
    if cov_robust
        C = cov_blockgeom(signal.data(channel_mask,:)');
    else
        C = cov(signal.data(channel_mask,:)');
    end
    C = (1-cov_shrinkage)*C + cov_shrinkage*mean(trace(C))*eye(length(C));
    M = LF'/C;
end

% apply M
signal = utl_register_field(signal,'timeseries','srcpot',reshape(M*signal.data(channel_mask,:,:),[],size(signal.data,2),size(signal.data,3)));
signal.etc.roi_labels = roi_labels;

if override_original
    % override signal.data & relabel channels
    signal.data = signal.srcpot;
    if size(signal.data,1) == signal.nbchan
        signal.chanlocs = struct('labels',roi_labels);
    else
        signal.nbchan = size(signal.data,1);
        signal.chanlocs = hlp_nanocache('cached_labels',10,@make_labels,roi_labels);
    end
end

% append the M and ok arguments to the online expression
exp_endfun('append_online',{'M',M,'channel_mask',channel_mask});

function chanlocs = make_labels(roi_labels)
chanlocs = struct('labels',[cellfun(@(f){sprintf('%s_X',f)},roi_labels) cellfun(@(f){sprintf('%s_Y',f)},roi_labels) cellfun(@(f){sprintf('%s_Z',f)},roi_labels)]);

repo: REST


text:
%% (Internal) Add the path only if not was already added.
%
%   paths_added = addpath_if_not_added(str_path, position)
%
% Arguments:
%
%      + str_path: paths to add
%
%      + position: where to add. See addpath doc.
%
% Output:
%
%      + paths_added: Only the paths added
%
% Example:
%
% See also addpath
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 14/5/2014
% Birthdate  : 21/4/2015
% Copyright 2008-2015
%
function paths_added = addpath_if_not_added(str_path, position)

paths_added = [];

if( nargin < 2 || isempty(position) )
    position = '-end';
end

if( ischar(str_path) )
    str_path = cellstr(str_path);
end

all_path = path;
for each_str = rowvec(str_path)
    aux_str = each_str{1};
    aux_val = strfind(all_path, aux_str);
    if(isempty(aux_val))
        paths_added = [paths_added; each_str];
        addpath(aux_str, position);
    end
end

author: marianux

repo: ecg-kit

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


text:
function gamblers_ruin_simulation ( a_stakes, b_stakes, game_num )

%*****************************************************************************80
%
%% GAMBLERS_RUIN_SIMULATION simulates the game of gambler's ruin.
%
%  Discussion:
%
%    Two players, A and B, repeatedly toss a coin.
%    For heads, A wins one dollar from B;
%    For tails, B wins one dollar from A.
%    Play continues until one player is bankrupt.
%
%    This program "plays" the game GAME_NUM times, always starting from
%    the same initial configuration.
%
%    it keeps track of the number of coin tosses required, the number
%    of times the lead changes, and the number of times each player wins.
%
%    At the end, it prints some statistics, and plots histograms of the
%    length of the game, and the number of lead changes.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer A_STAKES, B_STAKES, the number of dollars that A and
%    B have initially.
%
%    Input, integer GAME_NUM, the number of games to simulate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GAMBLERS_RUIN_SIMULATION\n' );
  fprintf ( 1, '  MATLAB version\n' );

  if ( nargin < 1 )
    a_stakes = input ( 'Enter A_STAKES: ' );
  else
    a_stakes = str2num ( a_stakes );
  end

  if ( nargin < 2 )
    b_stakes = input ( 'Enter B_STAKES: ' );
  else
    b_stakes = str2num ( b_stakes );
  end

  if ( nargin < 3 )
    game_num = input ( 'Enter GAME_NUM, the number of games to play: ' );
  else
    game_num = str2num ( game_num );
  end

  a_wins = 0;
  b_wins = 0;
%
%  Play GAME_NUM games.
%
  for game = 1 : game_num

    step_num = 0;
    leader = '0';
    flip_num = -1;
    a = a_stakes;
    b = b_stakes;

    while ( 0 < a & 0 < b )

      step_num = step_num + 1;

      r = rand ( );

      if ( r <= 0.5 )
        a = a + 1;
        b = b - 1;
      else
        a = a - 1;
        b = b + 1;
      end

      if ( a_stakes < a & leader ~= 'A' )
        leader = 'A';
        flip_num = flip_num + 1;
      elseif ( a < a_stakes & leader ~= 'B' )
        leader = 'B';
        flip_num = flip_num + 1;
      end

    end

    if ( b == 0 )
      a_wins = a_wins + 1;
    else
      b_wins = b_wins + 1;
    end

    flip(game) = flip_num;
    step(game) = step_num;

  end
%
%  Average over the number of games.
%
  step_ave = sum ( step(1:game_num) ) / game_num;
  prob_a_win = a_wins / game_num;
  prob_b_win = b_wins / game_num;
%
%  Print statistics.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of games in simulation = %d\n', game_num );
  fprintf ( 1, '  A starts game with %d\n', a_stakes );
  fprintf ( 1, '  B starts game with %d\n', b_stakes );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Average number of steps = %f\n', step_ave );
  fprintf ( 1, '  Expected number of steps = %d\n', a_stakes * b_stakes );
  fprintf ( 1, '  Maximum number of steps = %d\n', max ( step ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Average chance of A winning = %f\n', prob_a_win );
  fprintf ( 1, '  Expected chance of A winning = %f\n', a_stakes / ( a_stakes + b_stakes ) );
  fprintf ( 1, '  Average chance of B winning = %f\n', prob_b_win );
  fprintf ( 1, '  Expected chance of B winning = %f\n', b_stakes / ( a_stakes + b_stakes ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Average number of flips = %f\n', sum ( flip ) / game_num );
  fprintf ( 1, '  Maximum number of flips = %d\n', max ( flip ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Initial flip table:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number Freq Prob\n' );
  fprintf ( 1, '\n' );
  for flip_num = 0 : 10
    i = find ( flip == flip_num );
    k = length ( i );
    fprintf ( 1, '  %4d  %6d  %f\n', flip_num, k, k / game_num );
  end
%
%  Plot steps.
%
  figure ( 1 )

  hist ( step, 40 )

  title ( 'Gambler''s ruin, number of steps' )
  xlabel ( 'Steps' )
  ylabel ( 'Frequency' )

  figure ( 2 )

  hist ( flip, 40 )
  title ( 'Gambler''s ruin, number of changes in the lead (flips)' )
  xlabel ( 'Flips' )
  ylabel ( 'Frequency' )

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function fx2 = p07_fx2 ( x )

%*****************************************************************************80
%
%% P07_FX2 evaluates the second derivative of the function for problem 7.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the abscissa.
%
%    Output, real FX2, the second derivative of the function at X.
%
  fx2 = 6.0 * x;

  return
end


text:
function varargout = SaveFMFSettings(varargin)
% SAVEFMFSETTINGS M-file for SaveFMFSettings.fig
%      SAVEFMFSETTINGS, by itself, creates a new SAVEFMFSETTINGS or raises the existing
%      singleton*.
%
%      H = SAVEFMFSETTINGS returns the handle to a new SAVEFMFSETTINGS or the handle to
%      the existing singleton*.
%
%      SAVEFMFSETTINGS('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in SAVEFMFSETTINGS.M with the given input arguments.
%
%      SAVEFMFSETTINGS('Property','Value',...) creates a new SAVEFMFSETTINGS or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before SaveFMFSettings_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to SaveFMFSettings_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help SaveFMFSettings

% Last Modified by GUIDE v2.5 04-Aug-2010 23:40:57

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @SaveFMFSettings_OpeningFcn, ...
                   'gui_OutputFcn',  @SaveFMFSettings_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before SaveFMFSettings is made visible.
function SaveFMFSettings_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to SaveFMFSettings (see VARARGIN)

% Choose default command line output for SaveFMFSettings
handles.output = struct;
handles.archstr = computer('arch');
if length(varargin) >= 1,
  handles.previous_values = varargin{1};
else
  handles.previous_values = struct;
end

% Output FPS
if isfield(handles.previous_values,'OutputFPS'),
  handles.output.OutputFPS = handles.previous_values.OutputFPS;
else
  handles.output.OutputFPS = 30;
end
set(handles.edit_OutputFPS,'String',num2str(handles.output.OutputFPS));

% Compression
if length(varargin) >= 2 && ~isempty(varargin{2}),
  handles.Compressions = varargin{2};
else
  switch handles.archstr,
    case 'win32',
      handles.Compressions = {'None','MSVC','RLE','Cinepak','Indeo5'};
    case 'win64',
      handles.Compressions = {'None','MSVC','RLE','IYUV'};
    otherwise
      handles.Compressions = {'None'};
  end
end
if isfield(handles.previous_values,'Compression') && ...
    ismember(handles.previous_values.Compression,handles.Compressions),
  handles.output.Compression = handles.previous_values.Compression;
elseif ismember('None',handles.Compressions'),
  handles.output.Compression = 'None';
else
  handles.output.Compression = 'Uncompressed AVI';
end

set(handles.popupmenu_Compression,'String',handles.Compressions,...
  'Value',find(strcmp(handles.Compressions,handles.output.Compression),1));

% Quality
if isfield(handles.previous_values,'Quality') && ...
    ~strcmp(handles.output.Compression,'None')
  handles.output.Quality = handles.previous_values.Quality;
else
  handles.output.Quality = 100;
end
set(handles.edit_Quality,'String',num2str(handles.output.Quality));
if strcmp(handles.output.Compression,'None'),
  set(handles.edit_Quality,'Enable','off');
end

% Start Frame
if isfield(handles.previous_values,'StartFrame'),
  handles.output.StartFrame = handles.previous_values.StartFrame;
else
  handles.output.StartFrame = 1;
end
set(handles.edit_StartFrame,'String',num2str(handles.output.StartFrame));

% End Frame
if isfield(handles.previous_values,'EndFrame'),
  handles.output.EndFrame = handles.previous_values.EndFrame;
else
  handles.output.EndFrame = inf;
end
set(handles.edit_EndFrame,'String',num2str(handles.output.EndFrame));

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes SaveFMFSettings wait for user response (see UIRESUME)
uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = SaveFMFSettings_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;
delete(handles.figure1);

% --- Executes on selection change in popupmenu_Compression.
function popupmenu_Compression_Callback(hObject, eventdata, handles)
% hObject    handle to popupmenu_Compression (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns popupmenu_Compression contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popupmenu_Compression
handles.output.Compression = handles.Compressions{get(hObject,'Value')};
if strcmp(handles.output.Compression,'None'),
  handles.output.Quality = 100;
  set(handles.edit_Quality,'Enable','off','String',num2str(handles.output.Quality));
else
  set(handles.edit_Quality,'Enable','on');
end
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function popupmenu_Compression_CreateFcn(hObject, eventdata, handles)
% hObject    handle to popupmenu_Compression (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function edit_Quality_Callback(hObject, eventdata, handles)
% hObject    handle to edit_Quality (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit_Quality as text
%        str2double(get(hObject,'String')) returns contents of edit_Quality as a double

v = str2double(get(hObject,'String'));
if isnan(v) || v <= 0 || v > 100,
  set(hObject,'String',num2str(handles.output.Quality));
else
  handles.output.Quality = v;
  guidata(hObject,handles);
end

% --- Executes during object creation, after setting all properties.
function edit_Quality_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit_Quality (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in pushbutton_Cancel.
function pushbutton_Cancel_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_Cancel (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

handles.output = struct;
guidata(hObject,handles);
uiresume(handles.figure1);

% --- Executes on button press in pushbutton_OK.
function pushbutton_OK_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_OK (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

uiresume(handles.figure1);

function edit_OutputFPS_Callback(hObject, eventdata, handles)
% hObject    handle to edit_OutputFPS (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit_OutputFPS as text
%        str2double(get(hObject,'String')) returns contents of edit_OutputFPS as a double

v = str2double(get(hObject,'String'));
if isnan(v) || v <= 0,
  set(hObject,'String',num2str(handles.output.OutputFPS));
else
  handles.output.OutputFPS = v;
  guidata(hObject,handles);
end

function edit_StartFrame_Callback(hObject, eventdata, handles)
% hObject    handle to edit_StartFrame (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit_StartFrame as text
%        str2double(get(hObject,'String')) returns contents of edit_StartFrame as a double
s = get(hObject,'String');

v = round(str2double(s));
if v <= handles.output.EndFrame && v >= 1,
  handles.output.StartFrame = v;
else
  set(hObject,'String',num2str(handles.output.StartFrame));
end
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function edit_StartFrame_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit_StartFrame (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function edit_EndFrame_Callback(hObject, eventdata, handles)
% hObject    handle to edit_EndFrame (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit_EndFrame as text
%        str2double(get(hObject,'String')) returns contents of edit_EndFrame as a double
s = strtrim(get(hObject,'String'));
if strcmpi(s,'inf') || strcmpi(s,'infty'),
  handles.output.EndFrame = inf;
  set(hObject,'String','Inf');
else
  v = round(str2double(s));
  if v >= 1 && v >= handles.output.StartFrame,
    handles.output.EndFrame = v;
  else
    set(hObject,'String',num2str(handles.output.EndFrame));
  end
end
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function edit_EndFrame_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit_EndFrame (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes when user attempts to close figure1.
function figure1_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: delete(hObject) closes the figure
handles.output = struct;
guidata(hObject,handles);
uiresume(handles.figure1);

author: kristinbranson

repo: JAABA


text:
function [ha] = m22ha(m2)
% Convert area from square meters to hectares.
% Chad A. Greene 2012
ha = m2*0.0001 ;

author: Sable

repo: mcbench-benchmarks


text:
function X = metoffice_remove_bins(X,maskfile)

X(~metoffice_get_mask(maskfile),:) = [];

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function sF = regularisation(nodes,y,lambda,varargin)
% computes a regularisation
% Syntax
%   sF = S2FunHarmonic.regularisation(S2Grid,f,lambda)
%   sF = S2FunHarmonic.regularisation(S2Grid,f,lambda,'bandwidth',
%        bandwidth,'node_weights',W,'fourier_weights',What)
%
% Input
%  S2Grid - grid on the sphere
%  f      - function values on the grid (may be multidimensional)
%  lambda - parameter for regularisation
%
% Options
%  bandwidth  - maximum harmonic degree
%  W          - weight w_n for the node nodes (default: Voronoi weights)
%  What       - weight what_{m,l} for the Fourier space (default Sobolev weights for s = 2)
%

[nodes,idx] = unique(nodes(:));
y = reshape(y(idx),length(nodes),[]);
N = length(nodes);
W = get_option(varargin,'node_weights');
if isempty(W)
  W = nodes.calcVoronoiArea;
end
bw = get_option(varargin,'bandwidth',floor(sqrt(2*pi/mean(W))));

What = get_option(varargin,'fourier_weights');
if isempty(What)
  s = 2;
  What = (2*(0:bw)+1).^(2*s);
  What = repelem(What,1:2:(2*bw+1))';
end

% initialize nfsft
nfsftmex('precompute',bw,1000,1,0);
plan = nfsftmex('init_advanced',bw,N,1);
nfsftmex('set_x',plan,[nodes.rho';nodes.theta']);
nfsftmex('precompute_x',plan);

y = W.*y;

% adjoint nfsft
nfsftmex('set_f',plan, y);
nfsftmex('adjoint',plan);
fhat = nfsftmex('get_f_hat_linear',plan);

[fhat,~] = pcg(@(x) afun(x,plan,lambda,W,What),fhat);

sF = S2FunHarmonic(fhat);

% finalize nfsft
nfsftmex('finalize',plan);

end

function y = afun(x,plan,lambda,W,What)
  nfsftmex('set_f_hat_linear',plan,x);
  nfsftmex('trafo',plan);
  y = nfsftmex('get_f',plan);

  y = W.*y;

  nfsftmex('set_f',plan,y);
  nfsftmex('adjoint',plan);
  y = nfsftmex('get_f_hat_linear',plan);

  y = y+lambda*What.*x;
end

author: mtex-toolbox

repo: mtex

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


text:
function fwriten(fid,a,N)
%FWRITEN Write a value to a file in text in a fixed length field.
%
% If A is an integer, it is padded with zeros to the left.  If A is a char,
% it is padded with spaces to the right.  N is the length (in characters)
% of the field.
%
% Written by: Tom Krauss, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

if isa(a,'integer')
    % This is the reasonable way to zeropad integers using sprintf, but it
    % doesn't seem to work for very large values of 'a' for versions of
    % MATLAB prior to 2011b.  Must have been a bug in the Matlab sprintf
    % implementation.
    % x = sprintf('%0*d', N, a);
    % The following workaround seems to work for all sizes of integers even
    % on versions of MATLAB prior to the 2011b fix.
    x = sprintf('%0*.0f', N, a);
else
    x = sprintf('%-*s', N, a);
end
fwrite(fid,x(1:N),'uint8');

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint