text:
function mpc = toggle_dcline(mpc, on_off)
%TOGGLE_DCLINE Enable, disable or check status of DC line modeling.
%   MPC = TOGGLE_DCLINE(MPC, 'on')
%   MPC = TOGGLE_DCLINE(MPC, 'off')
%   T_F = TOGGLE_DCLINE(MPC, 'status')
%
%   Enables, disables or checks the status of a set of OPF userfcn
%   callbacks to implement DC lines as a pair of linked generators.
%   While it uses the OPF extension mechanism, this implementation
%   works for simple power flow as well as OPF problems.
%
%   These callbacks expect to find a 'dcline' field in the input MPC,
%   where MPC.dcline is an ndc x 17 matrix with columns as defined
%   in IDX_DCLINE, where ndc is the number of DC lines.
%
%   The 'int2ext' callback also packages up flow results and stores them
%   in appropriate columns of MPC.dcline.
%
%   NOTE: Because of the way this extension modifies the number of
%   rows in the gen and gencost matrices, caution must be taken
%   when using it with other extensions that deal with generators.
%
%   Examples:
%       mpc = loadcase('t_case9_dcline');
%       mpc = toggle_dcline(mpc, 'on');
%       results1 = runpf(mpc);
%       results2 = runopf(mpc);
%
%   See also IDX_DCLINE, ADD_USERFCN, REMOVE_USERFCN, RUN_USERFCN.

%   MATPOWER
%   Copyright (c) 2011-2016, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

if strcmp(upper(on_off), 'ON')
    %% define named indices into data matrices
    c = idx_dcline;

    %% check for proper input data
    if ~isfield(mpc, 'dcline') || size(mpc.dcline, 2) < c.LOSS1
        error('toggle_dcline: case must contain a ''dcline'' field, an ndc x %d matrix.', c.LOSS1);
    end
    if isfield(mpc, 'dclinecost') && size(mpc.dcline, 1) ~= size(mpc.dclinecost, 1)
        error('toggle_dcline: number of rows in ''dcline'' field (%d) and ''dclinecost'' field (%d) do not match.', ...
            size(mpc.dcline, 1), size(mpc.dclinecost, 1));
    end
    l0 = mpc.dcline(:, c.LOSS0);
    l1 = mpc.dcline(:, c.LOSS1);
    k = find( (l0 + l1 .* mpc.dcline(:, c.PMIN) < 0) | ...
              (l0 + l1 .* mpc.dcline(:, c.PMAX) < 0) );
    if ~isempty(k)
        warning('toggle_dcline: loss can be negative for DC line from bus %d to %d\n', ...
            [mpc.dcline(k, c.F_BUS:c.T_BUS)]');
    end

    %% add callback functions
    %% note: assumes all necessary data included in 1st arg (mpc, om, results)
    %%       so, no additional explicit args are needed
    mpc = add_userfcn(mpc, 'ext2int', @userfcn_dcline_ext2int);
    mpc = add_userfcn(mpc, 'formulation', @userfcn_dcline_formulation);
    mpc = add_userfcn(mpc, 'int2ext', @userfcn_dcline_int2ext);
    mpc = add_userfcn(mpc, 'printpf', @userfcn_dcline_printpf);
    mpc = add_userfcn(mpc, 'savecase', @userfcn_dcline_savecase);
    mpc.userfcn.status.dcline = 1;
elseif strcmp(upper(on_off), 'OFF')
    mpc = remove_userfcn(mpc, 'savecase', @userfcn_dcline_savecase);
    mpc = remove_userfcn(mpc, 'printpf', @userfcn_dcline_printpf);
    mpc = remove_userfcn(mpc, 'int2ext', @userfcn_dcline_int2ext);
    mpc = remove_userfcn(mpc, 'formulation', @userfcn_dcline_formulation);
    mpc = remove_userfcn(mpc, 'ext2int', @userfcn_dcline_ext2int);
    mpc.userfcn.status.dcline = 0;
elseif strcmp(upper(on_off), 'STATUS')
    if isfield(mpc, 'userfcn') && isfield(mpc.userfcn, 'status') && ...
            isfield(mpc.userfcn.status, 'dcline')
        mpc = mpc.userfcn.status.dcline;
    else
        mpc = 0;
    end
else
    error('toggle_dcline: 2nd argument must be ''on'', ''off'' or ''status''');
end

%%-----  ext2int  ------------------------------------------------------
function mpc = userfcn_dcline_ext2int(mpc, mpopt, args)
%
%   mpc = userfcn_dcline_ext2int(mpc, mpopt, args)
%
%   This is the 'ext2int' stage userfcn callback that prepares the input
%   data for the formulation stage. It expects to find a 'dcline' field
%   in mpc as described above. The optional args are not currently used.
%   It adds two dummy generators for each in-service DC line, with the
%   appropriate upper and lower generation bounds and corresponding
%   entries in gencost. It also expands columns of any A and N matrices
%   accordingly, if present.

%% define named indices into data matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
c = idx_dcline;

%% initialize some things
if isfield(mpc, 'dclinecost')
    havecost = 1;
else
    havecost = 0;
end

%% save version with external indexing
mpc.order.ext.dcline = mpc.dcline;              %% external indexing
if havecost
    mpc.order.ext.dclinecost = mpc.dclinecost;  %% external indexing
end

%% work with only in-service DC lines
mpc.order.dcline.status.on  = find(mpc.dcline(:, c.BR_STATUS) >  0);
mpc.order.dcline.status.off = find(mpc.dcline(:, c.BR_STATUS) <= 0);

%% remove out-of-service DC lines
dc = mpc.dcline(mpc.order.dcline.status.on, :); %% only in-service DC lines
if havecost
    dcc = mpc.dclinecost(mpc.order.dcline.status.on, :);    %% only in-service DC lines
    mpc.dclinecost = dcc;
end
ndc = size(dc, 1);          %% number of in-service DC lines
o = mpc.order;

%%-----  convert stuff to internal indexing  -----
dc(:, c.F_BUS) = o.bus.e2i(dc(:, c.F_BUS));
dc(:, c.T_BUS) = o.bus.e2i(dc(:, c.T_BUS));
mpc.dcline = dc;

%%-----  create gens to represent DC line terminals  -----
%% ensure consistency of initial values of PF, PT and losses
%% (for simple power flow cases)
dc(:, c.PT) = dc(:, c.PF) - (dc(:, c.LOSS0) + dc(:, c.LOSS1) .* dc(:, c.PF));

%% create gens
fg = zeros(ndc, size(mpc.gen, 2));
fg(:, MBASE)        = 100;
fg(:, GEN_STATUS)   =  dc(:, c.BR_STATUS);  %% status (should be all 1's)
fg(:, PMIN)         = -Inf;
fg(:, PMAX)         =  Inf;
tg = fg;
fg(:, GEN_BUS)      =  dc(:, c.F_BUS);      %% from bus
tg(:, GEN_BUS)      =  dc(:, c.T_BUS);      %% to bus
fg(:, PG)           = -dc(:, c.PF);         %% flow (extracted at "from")
tg(:, PG)           =  dc(:, c.PT);         %% flow (injected at "to")
fg(:, QG)           =  dc(:, c.QF);         %% VAr injection at "from"
tg(:, QG)           =  dc(:, c.QT);         %% VAr injection at "to"
fg(:, VG)           =  dc(:, c.VF);         %% voltage set-point at "from"
tg(:, VG)           =  dc(:, c.VT);         %% voltage set-point at "to"
k = find(dc(:, c.PMIN) >= 0);           %% min positive direction flow
if ~isempty(k)                              %% contrain at "from" end
    fg(k, PMAX)     = -dc(k, c.PMIN);       %% "from" extraction lower lim
end
k = find(dc(:, c.PMAX) >= 0);           %% max positive direction flow
if ~isempty(k)                              %% contrain at "from" end
    fg(k, PMIN)     = -dc(k, c.PMAX);       %% "from" extraction upper lim
end
k = find(dc(:, c.PMIN) < 0);            %% max negative direction flow
if ~isempty(k)                              %% contrain at "to" end
    tg(k, PMIN)     =  dc(k, c.PMIN);       %% "to" injection lower lim
end
k = find(dc(:, c.PMAX) < 0);            %% min negative direction flow
if ~isempty(k)                              %% contrain at "to" end
    tg(k, PMAX)     =  dc(k, c.PMAX);       %% "to" injection upper lim
end
fg(:, QMIN)         =  dc(:, c.QMINF);      %% "from" VAr injection lower lim
fg(:, QMAX)         =  dc(:, c.QMAXF);      %% "from" VAr injection upper lim
tg(:, QMIN)         =  dc(:, c.QMINT);      %%  "to"  VAr injection lower lim
tg(:, QMAX)         =  dc(:, c.QMAXT);      %%  "to"  VAr injection upper lim

%% fudge PMAX a bit if necessary to avoid triggering
%% dispatchable load constant power factor constraints
fg(isload(fg), PMAX) = -1e-6;
tg(isload(tg), PMAX) = -1e-6;

%% set all terminal buses to PV (except ref bus)
refbus = find(mpc.bus(:, BUS_TYPE) == REF);
mpc.bus(dc(:, c.F_BUS), BUS_TYPE) = PV;
mpc.bus(dc(:, c.T_BUS), BUS_TYPE) = PV;
mpc.bus(refbus, BUS_TYPE) = REF;

%% expand A and N, if present
nb = size(mpc.bus, 1);
ng = size(mpc.gen, 1);
if isfield(mpc, 'A') && ~isempty(mpc.A)
    [mA, nA] = size(mpc.A);
    if nA >= 2*nb + 2*ng    %% assume AC dimensions
        mpc.A = [   mpc.A(:, 1:2*nb+ng)      sparse(mA, 2*ndc) ...
                    mpc.A(:, 2*nb+ng+(1:ng)) sparse(mA, 2*ndc) ...
                    mpc.A(:, (2*nb+2*ng+1:nA)) ];
    else                    %% assume DC dimensions
        mpc.A = [   mpc.A(:, 1:nb+ng)   sparse(mA, 2*ndc) ...
                    mpc.A(:, (nb+ng+1:nA)) ];
    end
end
if isfield(mpc, 'N') && ~isempty(mpc.N)
    [mN, nN] = size(mpc.N);
    if nN >= 2*nb + 2*ng    %% assume AC dimensions
        mpc.N = [   mpc.N(:, 1:2*nb+ng)      sparse(mN, 2*ndc) ...
                    mpc.N(:, 2*nb+ng+(1:ng)) sparse(mN, 2*ndc) ...
                    mpc.N(:, (2*nb+2*ng+1:nN)) ];
    else                    %% assume DC dimensions
        mpc.N = [   mpc.N(:, 1:nb+ng)   sparse(mN, 2*ndc) ...
                    mpc.N(:, (nb+ng+1:nN)) ];
    end
end

%% append dummy gens
mpc.gen = [mpc.gen; fg; tg];

%% gencost
if isfield(mpc, 'gencost') && ~isempty(mpc.gencost)
    ngcc = size(mpc.gencost, 2);    %% dimensions of gencost
    [pc, qc] = pqcost(mpc.gencost, ng); %% split out re/active costs
    if havecost                 %% user has provided costs
        ndccc = size(dcc, 2);       %% number of dclinecost columns
        ccc = max([ngcc; ndccc]);   %% number of columns in new gencost
        if ccc > ngcc               %% right zero-pad gencost
            pc = [pc zeros(ng, ccc-ngcc)];
            if ~isempty(qc)         %% pad Qg costs, too
                qc = [qc zeros(ng, ccc-ngcc)];
            end
            ngcc = ccc;
        end

        %% flip function across vertical axis and append to gencost
        %% (PF for DC line = -PG for dummy gen at "from" bus)
        for k = 1:ndc
            if dcc(k, MODEL) == POLYNOMIAL
                nc = dcc(k, NCOST);
                temp = dcc(k, NCOST+(1:nc));
                %% flip sign on coefficients of odd terms
                %% (every other starting with linear term,
                %%  that is, the next to last one)
                temp((nc-1):-2:1) = -temp((nc-1):-2:1);
            else  %% dcc(k, MODEL) == PW_LINEAR
                nc = dcc(k, NCOST);
                temp = dcc(k, NCOST+(1:2*nc));
                %% switch sign on horizontal coordinate
                xx = -temp(1:2:2*nc);
                yy =  temp(2:2:2*nc);
                temp(1:2:2*nc) = xx(end:-1:1);
                temp(2:2:2*nc) = yy(end:-1:1);
            end
            padding = zeros(1, ngcc-NCOST-length(temp));
            gck = [dcc(k, 1:NCOST) temp padding];

            %% append to gencost
            pc = [pc; gck];
        end
        %% define zero cost for "to" end gen
        dcgc = ones(ndc, 1) * [2 0 0 2 zeros(1, ngcc-4)];
    else
        %% use zero cost as default on "from" end gen
        dcgc = ones(ndc, 1) * [2 0 0 2 zeros(1, ngcc-4)];
        pc = [pc; dcgc];
    end
    %% always use zero cost on "to" end gen
    pc = [pc; dcgc];

    %% reassemble gencost
    if ~isempty(qc)
        qc = [qc; dcgc; dcgc];  %% always use zero cost Qg, both ends
        mpc.gencost = [pc; qc];
    else
        mpc.gencost = pc;
    end
end

%%-----  formulation  --------------------------------------------------
function om = userfcn_dcline_formulation(om, mpopt, args)
%
%   om = userfcn_dcline_formulation(om, mpopt, args)
%
%   This is the 'formulation' stage userfcn callback that defines the
%   user constraints for the dummy generators representing DC lines.
%   It expects to find a 'dcline' field in the mpc stored in om, as
%   described above. By the time it is passed to this callback,
%   MPC.dcline should contain only in-service lines and the from and
%   two bus columns should be converted to internal indexing. The
%   optional args are not currently used.
%
%   If Pf, Pt and Ploss are the flow at the "from" end, flow at the
%   "to" end and loss respectively, and L0 and L1 are the linear loss
%   coefficients, the the relationships between them is given by:
%       Pf - Ploss = Pt
%       Ploss = L0 + L1 * Pf
%   If Pgf and Pgt represent the injections of the dummy generators
%   representing the DC line injections into the network, then
%   Pgf = -Pf and Pgt = Pt, and we can combine all of the above to
%   get the following constraint on Pgf ang Pgt:
%       -Pgf - (L0 - L1 * Pgf) = Pgt
%   which can be written:
%       -L0 <= (1 - L1) * Pgf + Pgt <= -L0

%% define named indices into data matrices
c = idx_dcline;

%% initialize some things
mpc = om.get_mpc();
dc = mpc.dcline;
ndc = size(dc, 1);              %% number of in-service DC lines
ng  = size(mpc.gen, 1) - 2*ndc; %% number of original gens/disp loads

%% constraints
nL0 = -dc(:, c.LOSS0) / mpc.baseMVA;
L1  =  dc(:, c.LOSS1);
Adc = [sparse(ndc, ng) spdiags(1-L1, 0, ndc, ndc) speye(ndc, ndc)];

%% add them to the model
om.add_lin_constraint('dcline', Adc, nL0, nL0, {'Pg'});

%%-----  int2ext  ------------------------------------------------------
function results = userfcn_dcline_int2ext(results, mpopt, args)
%
%   results = userfcn_dcline_int2ext(results, mpopt, args)
%
%   This is the 'int2ext' stage userfcn callback that converts everything
%   back to external indexing and packages up the results. It expects to
%   find a 'dcline' field in the results struct as described for mpc
%   above. It also expects that the last 2*ndc entries in the gen and
%   gencost matrices correspond to the in-service DC lines (where ndc is
%   the number of rows in MPC.dcline. These extra rows are removed from
%   gen and gencost and the flow is taken from the PG of these gens and
%   placed in the flow column of the appropiate dcline row. Corresponding
%   columns are also removed from any A and N matrices, if present. The
%   optional args are not currently used.

%% define named indices into data matrices
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
c = idx_dcline;

%% initialize some things
o = results.order;
k = find(o.ext.dcline(:, c.BR_STATUS));
ndc = length(k);                    %% number of in-service DC lines
ng  = size(results.gen, 1) - 2*ndc; %% number of original gens/disp loads
nb  = size(results.bus, 1);

%% extract dummy gens
fg = results.gen(ng    +(1:ndc), :);
tg = results.gen(ng+ndc+(1:ndc), :);

%% remove dummy gens
results.gen     = results.gen(1:ng, :);
if isfield(results, 'gencost') && ~isempty(results.gencost)
    results.gencost = results.gencost(1:ng, :);
end

%% delete corresponding columns from A and N, if present
if isfield(results, 'A') && ~isempty(results.A)
    [mA, nA] = size(results.A);
    if nA >= 2*nb + 2*ng + 4*ndc    %% assume AC dimensions
        results.A = results.A(:, [1:2*nb+ng 2*nb+ng+2*ndc+(1:ng) 2*nb+2*ng+4*ndc+1:nA]);
    else                            %% assume DC dimensions
        results.A = results.A(:, [1:nb+ng nb+ng+2*ndc+1:nA]);
    end
end
if isfield(results, 'N') && ~isempty(results.N)
    [mN, nN] = size(results.N);
    if nN >= 2*nb + 2*ng + 4*ndc    %% assume AC dimensions
        results.N = results.N(:, [1:2*nb+ng 2*nb+ng+2*ndc+(1:ng) 2*nb+2*ng+4*ndc+1:nN]);
    else                            %% assume DC dimensions
        results.N = results.N(:, [1:nb+ng nb+ng+2*ndc+1:nN]);
    end
end

%% get the solved flows
results.dcline(:, c.PF) = -fg(:, PG);
results.dcline(:, c.PT) =  tg(:, PG);
results.dcline(:, c.QF) =  fg(:, QG);
results.dcline(:, c.QT) =  tg(:, QG);
results.dcline(:, c.VF) =  fg(:, VG);
results.dcline(:, c.VT) =  tg(:, VG);
if size(fg, 2) >= MU_QMIN
    results.dcline(:, c.MU_PMIN ) = fg(:, MU_PMAX) + tg(:, MU_PMIN);
    results.dcline(:, c.MU_PMAX ) = fg(:, MU_PMIN) + tg(:, MU_PMAX);
    results.dcline(:, c.MU_QMINF) = fg(:, MU_QMIN);
    results.dcline(:, c.MU_QMAXF) = fg(:, MU_QMAX);
    results.dcline(:, c.MU_QMINT) = tg(:, MU_QMIN);
    results.dcline(:, c.MU_QMAXT) = tg(:, MU_QMAX);
end

%%-----  convert stuff back to external indexing  -----
results.order.int.dcline = results.dcline;  %% save internal version
%% copy results to external version
o.ext.dcline(k, c.PF:c.VT) = results.dcline(:, c.PF:c.VT);
if size(results.dcline, 2) == c.MU_QMAXT
    o.ext.dcline(k, c.MU_PMIN:c.MU_QMAXT) = results.dcline(:, c.MU_PMIN:c.MU_QMAXT);
end
results.dcline = o.ext.dcline;              %% use external version

%%-----  printpf  ------------------------------------------------------
function results = userfcn_dcline_printpf(results, fd, mpopt, args)
%
%   results = userfcn_dcline_printpf(results, fd, mpopt, args)
%
%   This is the 'printpf' stage userfcn callback that pretty-prints the
%   results. It expects a results struct, a file descriptor and a MATPOWER
%   options struct. The optional args are not currently used.

%% define named indices into data matrices
c = idx_dcline;

%% options
SUPPRESS        = mpopt.out.suppress_detail;
if SUPPRESS == -1
    if size(results.bus, 1) > 500
        SUPPRESS = 1;
    else
        SUPPRESS = 0;
    end
end
OUT_ALL = mpopt.out.all;
OUT_BRANCH      = OUT_ALL == 1 || (OUT_ALL == -1 && ~SUPPRESS && mpopt.out.branch);
if OUT_ALL == -1
    OUT_ALL_LIM = ~SUPPRESS * mpopt.out.lim.all;
elseif OUT_ALL == 1
    OUT_ALL_LIM = 2;
else
    OUT_ALL_LIM = 0;
end
if OUT_ALL_LIM == -1
    OUT_LINE_LIM    = ~SUPPRESS * mpopt.out.lim.line;
else
    OUT_LINE_LIM    = OUT_ALL_LIM;
end
ctol = mpopt.opf.violation; %% constraint violation tolerance
ptol = 1e-4;                %% tolerance for displaying shadow prices

%%-----  print results  -----
dc = results.dcline;
ndc = size(dc, 1);
kk = find(dc(:, c.BR_STATUS) ~= 0);
if OUT_BRANCH
    fprintf(fd, '\n================================================================================');
    fprintf(fd, '\n|     DC Line Data                                                             |');
    fprintf(fd, '\n================================================================================');
    fprintf(fd, '\n Line    From     To        Power Flow           Loss     Reactive Inj (MVAr)');
    fprintf(fd, '\n   #      Bus     Bus   From (MW)   To (MW)      (MW)       From        To   ');
    fprintf(fd, '\n------  ------  ------  ---------  ---------  ---------  ---------  ---------');
    loss = 0;
    for k = 1:ndc
        if dc(k, c.BR_STATUS)   %% status on
            fprintf(fd, '\n%5d%8d%8d%11.2f%11.2f%11.2f%11.2f%11.2f', ...
                        k, dc(k, c.F_BUS:c.T_BUS), dc(k, c.PF:c.PT), ...
                        dc(k, c.PF) - dc(k, c.PT), dc(k, c.QF:c.QT) );
            loss = loss + dc(k, c.PF) - dc(k, c.PT);
        else
            fprintf(fd, '\n%5d%8d%8d%11s%11s%11s%11s%11s', ...
                        k, dc(k, c.F_BUS:c.T_BUS), '-  ', '-  ', '-  ', '-  ', '-  ');
        end
    end
    fprintf(fd, '\n                                              ---------');
    fprintf(fd, '\n                                     Total:%11.2f\n', loss);
end

if OUT_LINE_LIM == 2 || (OUT_LINE_LIM == 1 && ...
        (any(dc(kk, c.PF) > dc(kk, c.PMAX) - ctol) || ...
         any(dc(kk, c.MU_PMIN) > ptol) || ...
         any(dc(kk, c.MU_PMAX) > ptol)))
    fprintf(fd, '\n================================================================================');
    fprintf(fd, '\n|     DC Line Constraints                                                      |');
    fprintf(fd, '\n================================================================================');
    fprintf(fd, '\n Line    From     To          Minimum        Actual Flow       Maximum');
    fprintf(fd, '\n   #      Bus     Bus    Pmin mu     Pmin       (MW)       Pmax      Pmax mu ');
    fprintf(fd, '\n------  ------  ------  ---------  ---------  ---------  ---------  ---------');
    for k = 1:ndc
        if OUT_LINE_LIM == 2 || (OUT_LINE_LIM == 1 && ...
                (dc(k, c.PF) > dc(k, c.PMAX) - ctol || ...
                 dc(k, c.MU_PMIN) > ptol || ...
                 dc(k, c.MU_PMAX) > ptol))
            if dc(k, c.BR_STATUS)   %% status on
                fprintf(fd, '\n%5d%8d%8d', k, dc(k, c.F_BUS:c.T_BUS) );
                if dc(k, c.MU_PMIN) > ptol
                    fprintf(fd, '%11.3f', dc(k, c.MU_PMIN) );
                else
                    fprintf(fd, '%11s', '-  ' );
                end
                fprintf(fd, '%11.2f%11.2f%11.2f', ...
                            dc(k, c.PMIN), dc(k, c.PF), dc(k, c.PMAX) );
                if dc(k, c.MU_PMAX) > ptol
                    fprintf(fd, '%11.3f', dc(k, c.MU_PMAX) );
                else
                    fprintf(fd, '%11s', '-  ' );
                end
            else
                fprintf(fd, '\n%5d%8d%8d%11s%11s%11s%11s%11s', ...
                            k, dc(k, c.F_BUS:c.T_BUS), '-  ', '-  ', '-  ', '-  ', '-  ');
            end
        end
    end
    fprintf(fd, '\n');
end

%%-----  savecase  -----------------------------------------------------
function mpc = userfcn_dcline_savecase(mpc, fd, prefix, args)
%
%   mpc = userfcn_dcline_savecase(mpc, fd, prefix, args)
%
%   This is the 'savecase' stage userfcn callback that prints the M-file
%   code to save the 'dcline' field in the case file. It expects a
%   MATPOWER case struct (mpc), a file descriptor and variable prefix
%   (usually 'mpc.'). The optional args are not currently used.

%% define named indices into data matrices
c = idx_dcline;

%% save it
ncols = size(mpc.dcline, 2);
fprintf(fd, '\n%%%%-----  DC Line Data  -----%%%%\n');
if ncols < c.MU_QMAXT
    fprintf(fd, '%%\tfbus\ttbus\tstatus\tPf\tPt\tQf\tQt\tVf\tVt\tPmin\tPmax\tQminF\tQmaxF\tQminT\tQmaxT\tloss0\tloss1\n');
else
    fprintf(fd, '%%\tfbus\ttbus\tstatus\tPf\tPt\tQf\tQt\tVf\tVt\tPmin\tPmax\tQminF\tQmaxF\tQminT\tQmaxT\tloss0\tloss1\tmuPmin\tmuPmax\tmuQminF\tmuQmaxF\tmuQminT\tmuQmaxT\n');
end
template = '\t%d\t%d\t%d\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g';
if ncols == c.MU_QMAXT
    template = [template, '\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f'];
end
template = [template, ';\n'];
fprintf(fd, '%sdcline = [\n', prefix);
fprintf(fd, template, mpc.dcline.');
fprintf(fd, '];\n');

%% to do, make it save mpc.dclinecost too (somehow I forgot this)
%% have a look at the saving of gencost in savecase for template

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function Q = StructureMeasure(prediction,GT)
% StructureMeasure computes the similarity between the foreground map and
% ground truth(as proposed in "Structure-measure: A new way to evaluate
% foreground maps" [Deng-Ping Fan et. al - ICCV 2017])
% Usage:
%   Q = StructureMeasure(prediction,GT)
% Input:
%   prediction - Binary/Non binary foreground map with values in the range
%                [0 1]. Type: double.
%   GT - Binary ground truth. Type: logical.
% Output:
%   Q - The computed similarity score
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Check input
if (~isa(prediction,'double'))
    error('The prediction should be double type...');
end
if ((max(prediction(:))>1) || min(prediction(:))<0)
    error('The prediction should be in the range of [0 1]...');
end
if (~islogical(GT))
    error('GT should be logical type...');
end

y = mean2(GT);

if (y==0)% if the GT is completely black
    x = mean2(prediction);
    Q = 1.0 - x; %only calculate the area of intersection
elseif(y==1)%if the GT is completely white
    x = mean2(prediction);
    Q = x; %only calcualte the area of intersection
else
    alpha = 0.5;
    Q = alpha*S_object(prediction,GT)+(1-alpha)*S_region(prediction,GT);
    if (Q<0)
      Q=0;
    end
end

end

author: ArcherFMY


text:
function [oa, pa, K, CM] = USFE_PCA(HSI, Tr, Te, dim, Trees)
[nx,ny,nz]=size(HSI);
data=reshape(HSI,nx*ny,nz);
clear HSI
%% matlab PCA from the stats toolbox (C:\Program Files\MATLAB\R2018b\toolbox\stats\stats\pca.m)
[code] = pca(data');
FE_Mpca=reshape(code(:,1:dim),nx,ny,dim);
[acc_Mean,acc_std,CM]=RF_ntimes_overal(FE_Mpca,Tr,Te,Trees);
pa=acc_Mean(1:dim,1);
oa=acc_Mean(dim+2,1);
K=acc_Mean(dim+3,1);

author: BehnoodRasti


text:
  classdef MyProgressMonitorNew < matlab.net.http.ProgressMonitor
      properties
          ProgHandle
          Direction matlab.net.http.MessageType
          Value uint64
          NewDir matlab.net.http.MessageType = matlab.net.http.MessageType.Request
      end

      methods
          function obj = MyProgressMonitorNew
              obj.Interval = .01;
          end

          function done(obj)
              obj.closeit();
          end

          function delete(obj)
              obj.closeit();
          end

          function set.Direction(obj, dir)
              obj.Direction = dir;
              obj.changeDir();
          end

          function set.Value(obj, value)
              obj.Value = value;
              obj.update();
          end
      end

      methods (Access = private)
          function update(obj,~)
              % called when Value is set
              import matlab.net.http.*
              if ~isempty(obj.Value)
                  if isempty(obj.Max)
                      % no maximum means we don't know length, so message changes on
                      % every call
                      value = 0;
                      if obj.Direction == MessageType.Request
                          msg = sprintf('Sent %d bytes...', obj.Value);
                      else
                          msg = sprintf('Received %d bytes...', obj.Value);
                      end
                  else
                      % maximum known; update proportional value
                      value = double(obj.Value)/double(obj.Max);
                      if obj.NewDir == MessageType.Request
                          % message changes only on change of direction
                          if obj.Direction == MessageType.Request
                              msg = 'Sending...';
                          else
                              msg = 'Receiving...';
                          end
                      end
                  end
                  if isempty(obj.ProgHandle)
                      % if we don't have a progress bar, display it for first time
                      obj.ProgHandle = ...
                          waitbar(value, msg, 'CreateCancelBtn', @(~,~)cancelAndClose(obj));

                      obj.NewDir = MessageType.Response;
                  elseif obj.NewDir == MessageType.Request || isempty(obj.Max)
                      % on change of direction or if no maximum known, change message
                      waitbar(value, obj.ProgHandle, msg);
                      obj.NewDir = MessageType.Response;
                  else
                      % no direction change else just update proportional value
                      waitbar(value, obj.ProgHandle);
                  end
              end

              function cancelAndClose(obj)
                  % Call the required CancelFcn and then close our progress bar. This is
                  % called when user clicks cancel or closes the window.
                  obj.CancelFcn();
                  obj.closeit();
              end
          end

          function changeDir(obj,~)
              % Called when Direction is set or changed. Leave the progress bar displayed.
              obj.NewDir = matlab.net.http.MessageType.Request;
          end
      end

      methods (Access=private)
          function closeit(obj)
              % Close the progress bar by deleting the handle so CloseRequestFcn isn't
              % called, because waitbar calls our cancelAndClose(), which would cause
              % recursion.
              if ~isempty(obj.ProgHandle)
                  delete(obj.ProgHandle);
                  obj.ProgHandle = [];
              end
          end
      end
  end

author: CelsoReyes


text:
function L = obslikebernoulli(X,hmm)
%
% Evaluate likelihood of data given observation model, for one continuous trial
%
% INPUT
% X          N by ndim data matrix
% hmm        hmm data structure
%
% OUTPUT
% B          Likelihood of N data points
%
% Author: Cam Higgins, OHBA, University of Oxford

K = hmm.K;
[T,ndim] = size(X);

L = zeros(T,K);

for k=1:K
    % expectation of log(p) is psi(a) - psi(a+b):
    pterm = X .* repmat( psi(hmm.state(k).W.a) - psi(hmm.state(k).W.a + hmm.state(k).W.b) ,T,1);
    qterm = (~X) .* repmat(psi(hmm.state(k).W.b) - psi(hmm.state(k).W.a + hmm.state(k).W.b) ,T,1);
    L(:,k) = sum(pterm+qterm,2);
end

L = exp(L);
end

author: OHBA-analysis

repo: HMM-MAR


text:
clc;

close all;
clear all;

%discremeinateanalysis
linearRegression
%logisticregression
naivbayes
regresstree
ccaregression
svm

author: faridani

repo: MatlabNLP

save_path: github-repos/MATLAB/faridani-MatlabNLP


text:
% sgwt_demo1 : SGWT for swiss roll data set
%
% This demo builds the SGWT for the swiss roll synthetic data set. It
% computes a set of scales adapted to the computed upper bound on the
% spectrum of the graph Laplacian, and displays the scaling function and
% the scaled wavelet kernels, as well as the corresponding frame bounds. It
% then computes the wavelets centered on a single vertex, and displays
% them. This essentally reproduces figure 3 from
% Hammond,Vangergheynst, Gribonval 2010.

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

% Create swiss roll point cloud
function sgwt_demo1
close all
fprintf('Welcome to SGWT demo #1\n');

npoints=500;
fprintf('Creating Swiss Roll point cloud with %g points\n',npoints);
dataparams=struct('n',npoints,'dataset',-1','noise',0,'state',0);
r=create_synthetic_dataset(dataparams);
x=rescale_center(r.x);

fprintf('Computing edge weights and graph Laplacian\n');
% Compute Weighted graph adjacency matrix, and graph Laplacian
d=distanz(x);
s=.1;
A=exp(-d.^2/(2*s^2));
L=full(sgwt_laplacian(A));

fprintf('Measuring largest eigenvalue, lmax = ');

%% Design filters for transform
Nscales=4;
lmax=sgwt_rough_lmax(L);
fprintf('%g\n',lmax);

fprintf('Designing transform in spectral domain\n');

designtype='abspline3';
%designtype='mexican_hat';
%designtype='meyer';
%designtype='simple_tf';

[g,t]=sgwt_filter_design(lmax,Nscales,'designtype',designtype);

arange=[0 lmax];
%% Display filter design in spectral domain
figure
sgwt_view_design(g,t,arange);
ylim([0 3])
set(gcf,'position',[0 780,600,300])
%% Chebyshev polynomial approximation
m=50; % order of polynomial approximation
fprintf('Computing Chebyshev polynomials of order %g for fast transform \n',m);
for k=1:numel(g)
  c{k}=sgwt_cheby_coeff(g{k},m,m+1,arange);
end

%% compute transform of delta at one vertex
jcenter=32; % vertex to center wavelets to be shown
fprintf('Computing forward transform of delta at vertex %g\n',jcenter);
N=size(L,1);
d=sgwt_delta(N,jcenter);
% forward transform, using chebyshev approximation
wpall=sgwt_cheby_op(d,L,c,arange);

fprintf('Displaying wavelets\n');
if (exist('OCTAVE_VERSION','builtin')~=0)
    % settings for octave
    msize=5;
    plotchar='s';
else
    % settings for native matlab
    msize=100;
    plotchar='.';
end

cp=[-1.4,-16.9,3.4]; % camera position
%% Visualize result

% show original point
ws=300;
figure;
xp=0; yp=ws+100;
set(gcf,'position',[xp,yp,ws-10,ws+10]);
scatter3(x(1,:),x(2,:),x(3,:),msize,d,plotchar);
set(gcf,'Colormap',[.5 .5 .5;1 0 0]);
clean_axes(cp);
title(sprintf('Vertex %g',jcenter));

% show wavelets
for n=1:Nscales+1
    wp=wpall{n};
    figure
    xp=mod(n,3)*(ws+10);
    yp=(1-floor((n)/3))*(ws+100);
    set(gcf,'position',[xp,yp,ws-10,ws+10]);
    scatter3(x(1,:),x(2,:),x(3,:),msize,wp,plotchar);
    colormap jet
    caxis([-1 1]*max(abs(wp)));
    clean_axes(cp);

    hcb=colorbar('location','north');
    cxt=get(hcb,'Xtick');
    cxt=[cxt(1),0,cxt(end)];
    set(hcb,'Xtick',cxt);
    cpos=get(hcb,'Position');
    cpos(4)=.02; % make colorbar thinner
    set(hcb,'Position',cpos);
    set(hcb,'Position',[.25 .91 .6 .02]);

    if n==1
      title('Scaling function');
    else
      title(sprintf('Wavelet at scale j=%g, t_j = %0.2f',n-1,t((n-1))));

    end
end

function clean_axes(cp)
xlim([-1 1]);ylim([-1 1]);zlim([-1 1]);
set(gca,'Xtick',[-1 0 1]);
set(gca,'Ytick',[-1 0 1]);
set(gca,'Ztick',[-1 0 1]);
axis square
set(gca,'CameraPosition',cp);

% rescale_center
% center input data at origin, then rescale so that all coordinates
% are between -1 and 1
%
% x should be dxN
function r=rescale_center(x)
N=size(x,2);
d=size(x,1);
x=x-repmat(mean(x,2),[1,N]);
c=max(abs(x(:)));
r=x/c;

repo: gspbox


text:
function a=bin2m(m)

%   a=bin2m(m)
%   converts from 0 and 1s to mseq in terms of -1 and 1

a=m*2-1;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function y = loggausspdf2(X, sigma)
% log pdf of Gaussian with zero mena
% Based on code written by Mo Chen (mochen@ie.cuhk.edu.hk). March 2009.
d = size(X,1);

[R,p]= chol(sigma);
if p ~= 0
    error('ERROR: sigma is not SPD.');
end
q = sum((R'\X).^2,1);  % quadratic term (M distance)
c = d*log(2*pi)+2*sum(log(diag(R)));   % normalization constant
y = -(c+q)/2;

author: lbasek

repo: image-denoising-benchmark


text:
function model = model_add_parts(model, lhs, ruleind, partner, ...
                                 filterind, numparts, psize, ...
                                 scale, coef_scale)
% Add a deformable parts to a rule.
%   model = model_add_parts(model, lhs, ruleind, partner, ...
%                           filterind, numparts, psize, ...
%                           scale, coef_scale)
%
% Return value
%   model       Object model
%
% Arguments
%   model       Object model
%   lhs         Parts are added to:
%   ruleind         model.rules{lhs}(ruleind)
%   partner     Partner ruleind:
%                   model.rules{lhs}(ruleind) and
%                   model.rules{lhs}(partner) are l/r mirror images of each other
%               Or, if length(partner) == 2:
%                   model.rules{lhs}(ruleind) and
%                   model.rules{partner(1)}(partner(2)) are l/r mirror images
%   filterind   Filter that parts are initialize from
%   numparts    Number of parts to add
%   psize       Size of each part
%   scale       Number of octaves down from lhs to place parts
%               (only scale = 0,1 have been tested)
%   coef_scale  Part filter coeficients are scaled by this value

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

% if the filter is mirrored, find its partner so mirrored
% parts can be added to it as well
if ~isempty(partner)
  if length(partner) == 2
    partner_lhs = partner(1);
    partner = partner(2);
  else
    partner_lhs = lhs;
  end
end

if nargin < 9
  coef_scale = 1;
end

source = model_get_block(model, model.filters(filterind));
pfilters = mkpartfilters(source, psize, numparts, scale);

for i = 1:numparts
  [model, symbolf] = model_add_terminal(model, 'w', coef_scale*pfilters(i).w);
  [model, N1] = model_add_nonterminal(model);

  % add deformation rule
  defoffset = 0;
  defparams = pfilters(i).alpha*[0.1 0 0.1 0];

  [model, rule] = model_add_def_rule(model, N1, symbolf, 'def_w', defparams);

  % add deformation symbols to rhs of rule
  anchor1 = pfilters(i).anchor;
  model.rules{lhs}(ruleind).rhs = [model.rules{lhs}(ruleind).rhs N1];
  model.rules{lhs}(ruleind).anchor = [model.rules{lhs}(ruleind).anchor anchor1];

  if ~isempty(partner)
    [model, symbolfp] = model_add_terminal(model, 'mirror_terminal', symbolf);
    [model, N2] = model_add_nonterminal(model);

    % add mirrored deformation rule
    model = model_add_def_rule(model, N2, symbolfp, 'mirror_rule', rule);

    x = pfilters(i).anchor(1);
    y = pfilters(i).anchor(2);
    % add deformation symbols to rhs of rule
    x2 = 2^scale*size(source, 2) - x - psize(2);
    anchor2 = [x2 y scale];
    model.rules{partner_lhs}(partner).rhs = ...
      [model.rules{partner_lhs}(partner).rhs N2];
    model.rules{partner_lhs}(partner).anchor = ...
      [model.rules{partner_lhs}(partner).anchor anchor2];
  end
end

author: rbgirshick

repo: voc-dpm

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


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

%%begin

% In the field of robotics there are many possible ways of representing
% orientations of which the most common are:
% - orthonormal rotation matrices (3x3),
% - three angles (1x3), and
% - quaternions.

% A rotation of pi/2 about the x-axis can be represented as an orthonormal rotation
% matrix

R = rotx(pi/2)
% which we can see is a 3x3 matrix.

% Such a matrix has the property that it's columns (and rows) are sets of orthogonal
% unit vectors.  The determinant of such a matrix is always 1

det(R)

% Let's create a more complex rotation

R = rotx(30, 'deg') * roty(50, 'deg') * rotz(10, 'deg')
% where this time we have specified the rotation angle in degrees.

% Any rotation can be expressed in terms of a single rotation about some axis
% in space

[theta,vec] = tr2angvec(R)
% where theta is the angle (in radians) and vec is unit vector representing the
% direction of the rotation axis.

% Commonly rotations are represented by Euler angles

eul = tr2eul(R)
% which are three angles such that R = rotz(a)*roty(b)*rotz(c), ie. the rotations
% required about the Z, then then the Y, then the Z axis.

% Rotations are also commonly represented by roll-pitch-yaw  angles

rpy = tr2rpy(R)
% which are three angles such that R = rotx(r)*roty(p)*rotz(y), ie. the rotations
% required about the X, then then the Y, then the Z axis.

% We can investigate the effects of rotations about different axes
% using this GUI based demonstration.  The menu buttons allow the rotation
% axes to be varied
%  *** close the window when you are done.
tripleangle('rpy', 'wait')

% The final useful form is the quaternion which comprises 4 numbers.  We can create
% a quaternion from an orthonormal matrix

q = Quaternion(R)
% where we can see that it comprises a scalar part and a vector part.  To convert back

q.R
% which is the same of the value of R we determined above.

% Quaternions are a class and the orientations they represent can be compounded, just
% as we do with rotation matrices by multiplication.

% First we create two quaternions

q1 = Quaternion( rotx(pi/2) )
q2 = Quaternion( roty(pi/2) )

% then the rotation of q1 followed by q2 is simply

q1 * q2

% We can also take the inverse of a Quaternion

q1 * inv(q1)
% which results in a null rotation (zero vector part)

author: petercorke

repo: robotics-toolbox-matlab


text:
function view = meshLighting(view, msh, L, distance);
% Set the lighting parameters for a surface mesh.
%
% view = meshLighting(<view>, <msh>, <L, distance=dialog>);
%
% view: mrVista view w/ mesh loaded. Defaults to selected gray.
%
% msh: mesh structure. Defaults to current mesh in view.
%
% L: structure with two fields:
%	ambient: [r g b] ambient lighting parameter (not sure details myself?)
%	diffuse: [r g b] diffuse lighting parameter
%
% distance: flag to move lights closer to, or further from, the mesh origin
% (and the mesh). If 1, don't move; <1, move closer; >1, move further away.
%
% If either L or distance are omitted, pops up a dialog.
%
% ras, 04/2008
if notDefined('view'), view = getSelectedGray; end
if notDefined('msh'), msh = viewGet(view,'currentmesh'); end

if isnumeric(msh), msh = view.mesh{msh};		end

lights = meshGet(msh,'lights');

% get params from dialog, if they're not provided
if notDefined('L') | notDefined('distance')
	%% put up a dialog for editing lighting levels --
	% for this function, we set the coefficients for all
	% lights to be the same. Build up an 'L' light struct
	% which will eventually be handed to mrMesh:
	if iscell(lights)
		L.diffuse = lights{1}.diffuse;
		L.ambient = lights{1}.ambient;
	else
		L.diffuse = lights(1).diffuse;
		L.ambient = lights(1).ambient;
	end

	% create the dialog
	dlg(1).fieldName = 'ambient';
	dlg(1).style	 = 'edit';
	dlg(1).string	 = 'Ambient Light Level';
	dlg(1).value	 = num2str(L.ambient);

	dlg(2).fieldName = 'diffuse';
	dlg(2).style	 = 'edit';
	dlg(2).string	 = 'Light Diffusion Coeffcient';
	dlg(2).value	 = num2str(L.diffuse);

	dlg(3).fieldName = 'distance';
	dlg(3).style	 = 'edit';
	dlg(3).string	 = ['Move lights closer (<1) or further (>1) ' ...
						'from current distance?'];
	dlg(3).value	 = '1';

	% get the response
	resp = generalDialog(dlg, 'Set Mesh Lighting');

	% parse the response
	L.diffuse = str2num(resp.diffuse);
	L.ambient = str2num(resp.ambient);
	distance  = str2num(resp.distance);
end

% apply the user settings to each light
host = meshGet(msh, 'host');
windowID = meshGet(msh, 'windowID');
for n = 1:length(lights)
	if iscell(lights)
		L.actor = lights{n}.actor;
		if distance ~= 1
			L.origin = distance .* lights{n}.origin;
		end

		lights{n} = mergeStructures(lights{n}, L);
	else
		L.actor = lights(n).actor;
		if distance ~= 1
			L.origin = distance .* lights(n).origin;
		end

		lights(n) = mergeStructures(lights(n), L);
	end
     mrMesh(host, windowID, 'set', L);
end

% update mesh, lights in view
msh = meshSet(msh, 'lights', lights);
view = viewSet(view, 'mesh', msh);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function image = read_mrtrix (filename)

% function: image = read_mrtrix (filename)
%
% returns a structure containing the header information and data for the MRtrix
% format image 'filename' (i.e. files with the extension '.mif' or '.mih').

f = fopen (filename, 'r');
assert(f ~= -1, 'error opening %s', filename);
L = fgetl(f);
if ~strncmp(L, 'mrtrix image', 12)
  fclose(f);
  error('%s is not in MRtrix format', filename);
end

transform = [];
dw_scheme = [];

while 1
  L = fgetl(f);
  if ~ischar(L), break, end;
  L = strtrim(L);
  if strcmp(L, 'END'), break, end;
  d = strfind (L,':');
  if isempty(d)
    disp (['invalid line in header: ''' L ''' - ignored']);
  else
    key = lower(strtrim(L(1:d(1)-1)));
    key = strrep(key, '-', '_'); % replace unpermitted '-' character in matlab fieldnames with '_'
    key = strrep(key, ' ', '_'); % replace unpermitted ' '
    value = strtrim(L(d(1)+1:end));
    if strcmp(key, 'dim')
      image.dim = str2num(char(split_strings (value, ',')))';
    elseif strcmp(key, 'vox')
      image.vox = str2num(char(split_strings (value, ',')))';
    elseif strcmp(key, 'layout')
      image.layout = value;
    elseif strcmp(key, 'datatype')
      image.datatype = value;
    elseif strcmp(key, 'transform')
      transform(end+1,:) = str2num(char(split_strings (value, ',')))';
    elseif strcmp(key, 'file')
      file = value;
    elseif strcmp(key, 'dw_scheme')
      dw_scheme(end+1,:) = str2num(char(split_strings (value, ',')))';
    else
      image = add_field (image, key, value);
    end
  end
end
fclose(f);

if ~isempty(transform)
  image.transform = transform;
  image.transform(4,:) = [ 0 0 0 1 ];
end

if ~isempty(dw_scheme)
  image.dw_scheme = dw_scheme;
end

if ~isfield (image, 'dim') || ~exist ('file') || ...
  ~isfield (image, 'layout') || ~isfield (image, 'datatype')
  error('critical entries missing in header - not reading data');
end

layout = split_strings(image.layout, ',');
order = (abs(str2num (char(layout)))+1)';

[ file, offset ] = strtok(file);
if isempty(offset), offset = 0; else offset = str2num(char(offset)); end
[f,g] = fileparts(filename);
if strcmp(file,'.'), file = filename; else file = fullfile (f, file); end

datatype = lower(image.datatype);
byteorder = datatype(end-1:end);

if strcmp(byteorder, 'le')
  f = fopen (file, 'r', 'l');
  datatype = datatype(1:end-2);
elseif strcmp(byteorder, 'be')
  f = fopen (file, 'r', 'b');
  datatype = datatype(1:end-2);
else
  if strcmp(datatype, 'bit')
    datatype = 'bit1';
    f = fopen(file, 'r', 'b');
  else
    f = fopen(file, 'r');
  end
end

assert(f ~= -1, 'error opening %s', filename);

fseek (f, offset, -1);
image.data = fread (f, inf, datatype);
fclose (f);

order(order)= 1:size(order,2);
image.data = reshape (image.data, image.dim(order));
image.data = ipermute (image.data, order);
for i=1:size(order,2)
  if layout{i}(1) == '-'
    image.data = flipdim(image.data, i);
  end
end

function S = split_strings (V, delim)
  S = {};
  while size(V,2) > 0
    [R, V] = strtok(V,delim);
    S{end+1} = R;
  end

function image = add_field (image, key, value)
  if isfield (image, key)
    previous = getfield (image, key);
    if iscell (previous)
      image = setfield (image, key, [ previous {value} ]);
    else
      image = setfield (image, key, { previous, value });
    end
  else
    image = setfield (image, key, value);
  end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [x, infos] = smu_nmf(V, rank, in_options)
% Stochastic multiplicative update for non-negative matrix factorization (SMU-NMF) algorithm.
%
% Inputs:
%       matrix      V
%       rank        rank
%       options     options
% Output:
%       w           solution of w
%       infos       information
%
%
% This file is part of NMFLibrary.
%
% Created by H.Kasai on Mar. 28, 2017
%
% Change log:
%
%       Mar. 15, 2018 (Hiroyuki Kasai): Fixed algorithm.
%
%       May. 20, 2019 (Hiroyuki Kasai): Added initialization module.
%
%       Jul. 12, 2022 (Hiroyuki Kasai): Modified code structures.
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.W_sub_mode    = 'Precon';
    local_options.H_sub_mode    = 'STD';
    local_options.accel         = false;
    local_options.ls            = false;
    local_options.precon        = false;
    local_options.h_repeat      = 1;
    local_options.rep_mode      = 'fix';
    local_options.robust        = false;

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % initialize factors
    init_options = options;
    [init_factors, ~] = generate_init_factors(V, rank, init_options);
    Wt = init_factors.W;
    H = init_factors.H;
    R = init_factors.R;

    % determine sub_mode
    if options.accel
        options.H_sub_mode = 'ACC';
    else
        if options.ls
            options.H_sub_mode = 'LS';
        else
            options.H_sub_mode = 'STD';
        end
        options.h_repeat = 1;
    end

    if options.precon
        options.W_sub_mode = 'Precon';
    else
        options.W_sub_mode = 'STD';
    end

    % permute samples
    if options.permute_on
        perm_idx = randperm(n);
    else
        perm_idx = 1:n;
    end
    V = V(:,perm_idx);
    H = H(:,perm_idx);

    % initialize
    method_name = sprintf('SMU-NMF (%s,%s)', options.W_sub_mode, options.H_sub_mode);
    epoch = 0;
    grad_calc_count = 0;

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    if strcmp(options.rep_mode, 'adaptive')
        rhoh = 1+(m+m*rank)/(1*(rank+1));
        alpha = 2;
        delta = 0.01;
    end

    % store initial info
    clear infos;
    [infos, f_val, optgap] = store_nmf_info(V, Wt, H, R, options, [], epoch, grad_calc_count, 0);

    if options.verbose > 1
        fprintf('SMU-NMF (%s,%s): Epoch = 0000, cost = %.16e, optgap = %.4e\n', options.W_sub_mode, options.H_sub_mode, f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main outer loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag] = check_stop_condition(epoch, infos, options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        cnt = 0;

        % main inner loop
        for t = 1 : options.batch_size : n - 1
            cnt = cnt + 1;

            % retrieve vt and ht
            vt = V(:,t:t+options.batch_size-1);
            ht = H(:,t:t+options.batch_size-1);

            % uddate ht
            Wtv = Wt.' * vt;
            WtW = Wt.' * Wt;
            if strcmp(options.H_sub_mode, 'ACC')
                if strcmp(options.rep_mode, 'adaptive')
                    gamma = 1;
                    eps0 = 1;
                    j = 1;
                    rhoh_alpha = rhoh*alpha;
                    %while j <= floor(1+rhoh*alpha) &&  gamma >= delta*eps0
                    while j <= rhoh_alpha && gamma >= delta*eps0
                        ht0 = ht;
                        ht = ht .* (Wtv) ./ (WtW * ht);
                        ht = ht + (ht<eps) .* eps;
                        if j == 1
                            eps0 = norm(ht0-ht);
                        end
                        gamma = norm(ht0-ht);
                        j = j+1;
                    end
                else
                    for ii=1:options.h_repeat
                        ht = ht .* (Wtv) ./ (WtW * ht);
                        ht = ht + (ht<eps) .* eps;
                    end
                end
            elseif strcmp(options.H_sub_mode, 'LS')
                ht = calc_nls_nmf(vt, Wt, 1e-16);
                ht = ht + (ht<eps) * eps;
            else
                ht = ht .* (Wtv) ./ (WtW * ht);
                ht = ht + (ht<eps) .* eps;
            end

            % update W
            delta_minus = vt * ht';
            delta_plus = Wt * (ht * ht');
            if strcmp(options.W_sub_mode, 'Precon')
                invhht = inv(ht * ht');
                invhht = max(invhht,0);
                Wt = Wt .* ((delta_minus * invhht) ./ (delta_plus * invhht));
            else
                Wt = Wt .* (delta_minus ./delta_plus);
            end
            Wt = Wt + (Wt<eps) .* eps;

            % store new h
            H(:,t:t+options.batch_size-1) = ht;

            grad_calc_count = grad_calc_count + m * options.batch_size;
        end

        % measure elapsed time
        elapsed_time = toc(start_time);

        % update epoch
        epoch = epoch + 1;

        % store info
        infos = store_nmf_info(V, Wt, H, R, options, infos, epoch, grad_calc_count, elapsed_time);

        % display info
        display_info(method_name, epoch, infos, options);

    end

    x.W = Wt;
    x.H(:,perm_idx) = H;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function img=minimizeTV(img,varargin)
%MINIMIZETV MATLAB wrapper for the CUDA stepest descend minimization of TV
% norm. Note that this does not minimize the TV noise, using the ROF model,
% this minimizes the TV alone, using gradient desced.
% Infinite iterations of this code will lead to a flat image.
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Ander Biguri
%--------------------------------------------------------------------------
if nargin >= 5
    [gpuids] = parse_inputs(varargin{3:length(varargin)});
else
    gpuids = GpuIds();
end
if nargin==1
    dtvg=1;
    ng=30;
else
    if nargin == 5
        dtvg=varargin{1};
        ng=varargin{2};
    else
        error('Wrong amount of inputs, 1 or 5 expected');
    end
end

img=minTV(img,dtvg,ng,gpuids.devices);

end

function [gpuids]=parse_inputs(varargin)
    %fprintf('parse_inputs0(varargin (%d))\n', length(varargin));
    if isempty(varargin)
        gpuids = GpuIds();
    else
        % create input parser
        p=inputParser;
        % add optional parameters
        addParameter(p,'gpuids', GpuIds());
        %execute
        parse(p,varargin{:});
        %extract
        gpuids=p.Results.gpuids;
    end
end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function GSS = parseScores(inputFile, predictor)
% parseScores
%	Parse the output from a predictor to generate the GSS
%
%	Input:
%	inputFile	a file with the output from the predictor
%	predictor	the predictor that was used. 'wolf' for WoLF PSORT, 'cello'
%               for CELLO, 'deeploc' for DeepLoc (opt, default 'wolf')
%
%	Output:
%	GSS         a gene scoring structure to be used in predictLocalization
%
%	The function normalizes the scores so that the best score for each gene
%	is 1.0.
%
%	Usage: GSS = parseScores(inputFile, predictor)

if nargin<2
    predictor='wolf';
else
    predictor=char(predictor);
end

fid=fopen(inputFile,'r');

if fid<1
    EM='Could not open file';
    dispEM(EM);
end

if strcmpi(predictor,'wolf')
    A=textscan(fid,'%s','Delimiter','\n','CommentStyle','#');

    %Each element should be for one gene, but some of them are on the form
    %"Pc20g11350: treating 9 X's as Glycines". Those should be removed
    I=~cellfun(@any,strfind(A{1},'treating'));

    B=regexp(A{1}(I),' ','split');

    %Reserve space for stuff
    GSS.compartments={};
    GSS.scores=[]; %Do not know number of comps yet
    GSS.genes=cell(numel(B),1);

    %Parsing is a bit cumbersome as ', ' is used as a delimiter in some
    %cases and ' ' in others. Use strrep to get rid of ','
    for i=1:numel(B)
        b=strrep(B{i},',','');
        GSS.genes{i}=b{1};

        %Then go through the compartments and add new ones as they are
        %found
        for j=2:2:numel(b)-1
            [~, J]=ismember(b(j),GSS.compartments);

            %Add new compartment if it does not exist
            if J==0
                GSS.compartments=[GSS.compartments;b(j)];
                J=numel(GSS.compartments);
                GSS.scores=[GSS.scores zeros(numel(B),1)];
            end

            GSS.scores(i,J)=str2double(b(j+1));
        end
    end
elseif strcmpi(predictor,'cello')
    fid=fopen(inputFile,'r');
    %Read the title line and fetch the list of compartments
    tline = fgetl(fid);
    tline=regexprep(tline,'^.+#Combined:\t','');
    tline=regexprep(tline,'\t#Most-likely-Location.+','');
    GSS.compartments=transpose(regexp(tline,'\t','split'));

    %Now iterate through the following lines in the file. Each row
    %corresponds to one gene and it consists of the scores for
    %compartments. Gene name is in the end of each line
    row=0;
    while 1
        row=row+1;
        tline = fgetl(fid);
        if ~ischar(tline)
            break;
        end
        tline=regexprep(tline,'^.+:\t','');
        tline=regexprep(tline,' .+','');
        tline=regexp(tline,'\t','split');
        GSS.scores(row,:)=str2double(tline(1:numel(GSS.compartments)));
        GSS.genes{row,1}=tline{1,end};
    end
elseif strcmpi(predictor,'deeploc')
    fid=fopen(inputFile,'r');
    %Read the title line and fetch the list of compartments
    tline = fgetl(fid);
    GSS.compartments=regexp(tline,'\t','split');
    GSS.compartments=GSS.compartments(3:end);

    %Now iterate through the following lines in the file. Each row
    %corresponds to one gene and it consists of the scores for
    %compartments. Gene name is in the end of each line
    row=0;
    while 1
        row=row+1;
        tline = fgetl(fid);
        if ~ischar(tline)
            break;
        end
        tline=regexp(tline,'\t','split');
        GSS.scores(row,:)=str2double(tline(3:end));
        GSS.genes{row,1}=tline{1,1};
    end
end

%Check if there are duplicate genes
[~, J, K]=unique(GSS.genes);

if numel(J)~=numel(K)
    EM='There are duplicate genes in the input file';
    dispEM(EM,false);
    GSS.genes=GSS.genes(J);
    GSS.scores=GSS.scores(J,:);
end

%Normalize
I=max(GSS.scores,[],2);
GSS.scores=bsxfun(@times, GSS.scores, 1./I);

fclose(fid);
end

author: SysBioChalmers

repo: RAVEN


text:
function [imdepthMin, imdepthMax, imdepthCol] = getDepthRange(bndinfo, glabels, dt, v0)
% getMinimumDepth(bndinfo, glabels, v0)
%
% Minimum depth is the depth of the foremost occluder.  Maximum depth is
% the depth of the current region, assuming that it is touching ground at
% its lowest point.  If a region is in the foreground, it's minimum depth
% may be a function, and the minimum depth is equal to the maximum depth.
% Otherwise, the minimum and maximum depths are single (different) values.
%
% log(Z) = log(f*y_c) - log(v_0-v_i)

global DO_DISPLAY;

%% Set some parameters

MIN_HV_DIST = 0.01;  % minimum distance between horizon and object contact
SKY_DEPTH = 2*(1 ./ MIN_HV_DIST);

%% Get pixel images for ground, vertical, and sky
labim = glabels(bndinfo.wseg);
vim = labim==2;
gim = labim==1;
sim = labim==3;

[imh, imw] = size(vim);

%% Make sure that horizon is above ground and below sky
scol = sum(sim, 2)>0;
gcol = sum(gim, 2)>0;
v0 = (1-v0) * imh;

try
minv = find(scol);  minv = minv(end-1);
maxv = find(gcol);  maxv = maxv(2);

if v0 < minv || v0 > maxv
    if minv > maxv
        v0 = maxv-1;
    else
        v0 = maxv-imh/20;
    end
end
%disp(num2str([v0 maxv minv]))
catch
end

v0 = 1 - v0/imh;
%disp(num2str(v0))

%figure(1), imagesc(cat(3, vim, gim, sim)), axis image

%% Get possible ground-vertical boundary pixels

stats = regionprops(bndinfo.wseg, 'BoundingBox', 'PixelIdxList');
bbox = vertcat(stats.BoundingBox);
idx = {stats.PixelIdxList};

% bbox = [x1 y1 x2 y2]
bbox = [bbox(:, 1) bbox(:, 2) bbox(:, 1)+bbox(:,3) bbox(:,2)+bbox(:,4)];

vfilt = [ones(3, 1) ; zeros(3, 1)];
gfilt = 1-vfilt;
boundaryim = (imfilter(double(vim), vfilt)==sum(vfilt(:))) & ...
    (imfilter(double(gim), gfilt)==sum(gfilt(:)));
boundaryim(end, :) = vim(end, :);
boundaryim = imdilate(boundaryim, ones(3, 1));

%figure(3), imagesc(boundaryim), axis image

bpts = cell(bndinfo.nseg, 1);
edges = bndinfo.edges.indices;
spLR = bndinfo.edges.spLR;
eim = zeros(bndinfo.imsize);
for k = 1:numel(edges)
    sp1 = spLR(k, 1);
    sp2 = spLR(k, 2);
    bndind = edges{k}(boundaryim(edges{k}));
    if ~isempty(bndind)
        if glabels(sp1)==2 % && glabels(sp2)==1
            bpts{sp1} = [bpts{sp1} ; bndind];
        elseif glabels(sp2)==2 % && glabels(sp1)==1
            bpts{sp2} = [bpts{sp2} ; bndind];
        end
    end
    eim(edges{k}) = 1;
end

%% Get a simple estimate of depth (trace columns upwards)

imdepthCol = zeros(bndinfo.imsize(1:2));
maxg = imh*ones(1, imw);
for y = imh:-1:1
    maxg(gim(y, :)) = y;
    imdepthCol(y, :) = maxg;
end
%imy = repmat((1:imh)', [1 imw]);
% imdepthCol = log(1 ./ max(v0 - (imh-imdepthCol)/imh, MIN_HV_DIST));
% imdepthCol(sim==1) = log(SKY_DEPTH);
for k = 1:numel(idx)
    imdepthCol(k) = (imh-median(imdepthCol(idx{k})))/imh;
end
% if DO_DISPLAY
%     figure(3), imagesc(imdepthCol), axis image
% end

%% Get the depth of each object that has a ground boundary

bx = zeros(numel(bpts), 2);
by = zeros(numel(bpts), 2);
bz = zeros(numel(bpts), 2);
bzmin = zeros(numel(bpts), 2);
bzmax = zeros(numel(bpts), 2);
isvalid = false(numel(bpts), 1);
for k = 1:numel(bpts)
    if numel(bpts{k})>5 || bbox(k, 4)+1>=imh

        [py, px] = getPerimeter(bndinfo, k, 50);
        bx(k, :) = bbox(k, [1 3]);

        cu = [];  cv = [];
        if  1 %bx(k, 2)-bx(k, 1) > imw/10  % only get line for larger objects

            try
                [cu, cv] = poly2contacts(dt, px, py, bndinfo.imsize, v0, 1.7, [], 0.5);
            catch
            end
            %plot(cu, cv, 'r*')
            cu = round(cu);  cv = round(cv);
            ind = (boundaryim(cv + (cu-1)*imh)); % | (cv==max(cv));
            cu = cu(ind);  cv = cv(ind);
        end

        [cu, ind] = unique(cu);
        cv = cv(ind);

        if ~isempty(cu)
            if numel(cu) == 1
                ypos = max(py);
                imslope = 0;
            elseif numel(cu)==2
                imslope = (cv(2)-cv(1)) / (cu(2)-cu(1));
                ypos = cv(1) - imslope*cu(1);
            else
                rline = robustfit(cu, cv);
                ypos = rline(1);
                imslope = rline(2);
            end
            by(k, :) = imslope*bx(k, :) + ypos;
            bz(k, :) = (1 ./ max(v0 - (imh-by(k, :))/imh, MIN_HV_DIST));
            isvalid(k) = true;
        end
    end
end

bzmin(isvalid, :) = bz(isvalid, :);
bzmax(isvalid, :) = bz(isvalid, :);
unknown = (~isvalid) & (glabels==2);

%% Assign a maximum depth to object of unknown depth and sky
missing = unknown;
for k = find(missing)'
    bx(k, :) = bbox(k, [1 3]);
    by(k, :) = bbox(k, 4);
    bzmax(k, :) = 1 ./ max(v0 - (imh-by(k, :))/imh, MIN_HV_DIST);
end

%% Assign a minimum depth to objects of unknown depth;

missing = unknown;

changed = 1;
while changed
    changed = 0;
    isvalid2 = isvalid;
    for k = find(missing)'
        bx(k, :) = bbox(k, [1 3]);

        bzmin(k, :) = bzmax(k, :);

        ind = find(spLR(:, 1)==k);
        neighbors = spLR(ind, 2);
        validn = isvalid(neighbors) & (bbox(neighbors, 4) > bbox(k, 4));
        neighbors = neighbors(validn);
        ind = ind(validn);

        maxov = 0;
        for tk2 = 1:numel(neighbors)
            k2 = neighbors(tk2);

            ov = min(bbox([k k2], 3)) - max(bbox([k k2], 1));
            %ov = numel(bndinfo.edges.indices{mod(ind(tk2)-1,bndinfo.ne)+1});

            if (ov > maxov) || (ov > 0.25*bbox(k, 3)-bbox(k,1))
                imslope = (by(k2, 2)-by(k2, 1))./(bx(k2, 2)-bx(k2, 1));
                ypos = by(k2, 1) - imslope*bx(k2, 1);
                by(k, :) = imslope*bx(k, :) + ypos;
                tmpz = 1 ./ max(v0 - (imh-by(k, :))/imh, MIN_HV_DIST);
                if all(tmpz < bzmin(k, :))
                    bzmin(k, :) = tmpz; %(1 ./ max(v0 - (imh-by(k, :))/imh, 0.01));
                end
                maxov = ov;
                missing(k) = false;
                isvalid2(k) = true;
                changed = 1;
            end

        end
    end

    isvalid = isvalid2;

end

depthim = zeros([imh imw]);
for k = 1:bndinfo.nseg
    [y, x] = ind2sub(bndinfo.imsize, idx{k});
    switch glabels(k)
        case 0
            z = bzmin(k, 1) + (x-bx(k, 1)) ./ ...
                (bx(k, 2)-bx(k, 1)) .* (bzmin(k,2)-bzmin(k,1));
        case 1
            z = 1 ./ max(v0 - (imh-y)/imh, MIN_HV_DIST);
        case 2
            z = bzmin(k, 1) + (x-bx(k, 1)) ./ ...
                (bx(k, 2)-bx(k, 1)) .* (bzmin(k,2)-bzmin(k,1));
        case 3
            z = SKY_DEPTH;
    end
    depthim(idx{k}) = log(z);
end

imdepthMin = depthim;

if DO_DISPLAY
    figure(1), imagesc(depthim), axis image, colormap jet
end
%print -f1 -djpeg99 ../tmp/gt_mindepth.jpg

%depthim = zeros([imh imw]);
for k = find(glabels==2)'
    [y, x] = ind2sub(bndinfo.imsize, idx{k});
    zk = bzmax(k, 1) + (x-bx(k, 1)) ./ ...
        (bx(k, 2)-bx(k, 1)) .* (bzmax(k,2)-bzmax(k,1));
    depthim(idx{k}) = log(zk);
end

if DO_DISPLAY
    figure(2), imagesc(depthim), axis image, colormap jet
end
%print -f2 -djpeg99 ../tmp/gt_maxdepth.jpg

imdepthMax = depthim;

author: Cloud-CV

repo: object-proposals


text:
function [analMu, analCov, diffMu, diffCov] = gpPosteriorMeanCovarTest(model, X)
% GPPOSTERIORMEANCOVARTEST Test the gradients of the mean and covariance.
% FORMAT
% DESC tests the gpPosteriorMeanCovar and gpPosteriorGradMeanCovar
% functions.
% ARG model : the model to test the gradients for.
% ARG X : the input locations to test the gradients for.
% RETURN analMu : the analytical gradients of the mean with respect
% to X.
% RETURN analCov : the analytical gradients of the covariance with respect
% to X.
% RETURN diffMu : the numerical gradients of the mean with respect
% to X.
% RETURN diffCov : the numerical gradients of the covariance with respect
% to X.
% RETURN delta : erros between the numerical and analytical
% gradients.
%
% SEEALSO : gpPosteriorMeanCovar, gpPosteriorGradMeanCovar
%
% COPYRIGHT : Neil D. Lawrence, 2006

% GP

[analMu, analCov] = gpPosteriorGradMeanCovar(model, X);
origX = X;
change = 1e-6;
for i = 1:size(X, 1)
  for j = 1:size(X, 2)
    X(i, j) = origX(i, j) - change;
    [muMinus, covMinus] = gpPosteriorMeanCovar(model, X);
    X(i, j) = origX(i, j) + change;
    [muPlus, covPlus] = gpPosteriorMeanCovar(model, X);
    X(i, j) = origX(i, j);
    diffMu{j}(i, :) = (muPlus(i, :) - muMinus(i, :))/(2*change);
    for k = 1:model.d
      % Not sure why the transpose is need here ... hope it isn't
      % two wrongs making a right ...
      diffCov{j, k}(:, i) = (covPlus{k}(i, :) - covMinus{k}(i, :))'/(2*change);
    end

  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% chrelerr.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function flag = chrelerr(fbest,stop)
% checks whether the required tolerance for a test function with known
% global minimum has already been achieved
% Input:
% fbest		function value to be checked
% stop(1)	relative error with which a global minimum with not too
%		small absolute value should be reached
% stop(2)	global minimum function value of a test function
% stop(3)	if abs(fglob) is very small, we stop if the function
%		value is less than stop(3)
% Output:
% flag          = 0 the required tolerance has been achieved
% 		= 1 otherwise

function flag = chrelerr(fbest,stop)
fglob = stop(2);
if fbest - fglob <= max(stop(1)*abs(fglob),stop(3))
  flag = 0;
else
  flag = 1;
end

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function M = stiefelstackedfactory(m, d, k)
% Stiefel(k, d)^m, represented as matrices of size m*d-by-k.
%
% function M = stiefelstackedfactory(m, d, k)
%
% Points on this manifold are matrices Y of size n x k, with n = m*d.
% Y is thought of as m matrices of size d x k each, stacked on top of each
% other. Call them Y1, ..., Ym. Each Yi is an orthonormal matrix, that is,
% its d rows are unit norm and are orthogonal to each other. Thus, this
% geometry is a product of Stiefel manifolds.
%
% To easily transform matrices Y to 3D arrays Y3 of size d x k x m such
% that each slice Y3(:, :, i) corresponds to one of the matrices Yi, use
% the functions
%
%    Y3 = M.to3D(Y)   and   Y = M.to2D(Y3).
%
% The ambient space R^(nxk) is endowed with the usual inner product
% <A, B> = trace(A'*B). This inner product is restricted to the tangent
% spaces of the present manifold, thus making it a Riemannian submanifold
% of the Euclidean space R^(nxk). Tangent vectors are represented as
% matrices of the same size as Y, and can likewise be converted to 3D
% arrays and back using to3D() and to2D().
%
% In dealing with this geometry, especially when dealing with the 3D array
% representations of points and tangent vectors, the tools multiprod,
% multitransp, multitrace, multiscale etc. available in Manopt are often
% useful.
%
% See also: stiefelfactory obliquefactory multiprod multitransp

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, May 4, 2015.
% Contributors:
% Change log:

    assert(k >= d, 'k must be at least as large as d.');

    n = m*d;

    M.name = @() sprintf('Manifold of %d orthonormal matrices of size %dx%d, stacked', m, d, k);

    M.dim = @() m*(k*d - .5*d*(d+1));

    M.size = @() [m, d, k];

    M.inner = @(x, d1, d2) d1(:).'*d2(:);

    M.norm = @(x, d) norm(d(:));

    M.dist = @(x, y) error('stiefelstackedfactory.dist not implemented yet.');

    M.typicaldist = @() sqrt(M.dim());

    % Convert a dxkxm matrix to an nxk matrix
    M.to2D = @to2D;
    function A2 = to2D(A3)
        A2 = reshape(multitransp(A3), [k, m*d])';
    end

    % Convert an nxk matrix to a dxkxm matrix
    M.to3D = @to3D;
    function A3 = to3D(A2)
        A3 = multitransp(reshape(A2', [k, d, m]));
    end

    % Given 2 3D matrices A and B of size dxkxm, returns a 3D matrix C of
    % size dxdxm such that each slice C(:, :, i) is the symmetric part of
    % the product A(:, :, i) * B(:, :, i)'. The name is short for
    % "symmetric-block-diagonal", because if A and B were transformed to
    % their 2D equivalents via to2D, then the output would contain the
    % symmetric parts of the diagonal blocks of A*B'.
    M.symbdiag = @symbdiag;
    function C = symbdiag(A, B)
        C = multisym(multiprod(A, multitransp(B)));
    end

    % Orthogonal projection from the ambient space R^(nxk) to the tangent
    % space at X.
    M.proj = @projection;
    function Zt = projection(Y, Z)
        Y3 = to3D(Y);
        Z3 = to3D(Z);
        Lambda = symbdiag(Y3, Z3);
        Zt3 = Z3 - multiprod(Lambda, Y3);
        Zt = to2D(Zt3);
    end

    M.tangent = M.proj;

    M.egrad2rgrad = M.proj;

    M.ehess2rhess = @ehess2rhess;
    function rhess = ehess2rhess(Y, egrad, ehess, Ydot)
        Y3 = to3D(Y);
        Ydot3 = to3D(Ydot);
        egrad3 = to3D(egrad);
        C = symbdiag(Y3, egrad3);
        CYdot = to2D(multiprod(C, Ydot3));
        rhess = projection(Y, ehess - CYdot);
    end

    M.retr = @retraction;
    function Y = retraction(Y, U, t)
        if nargin < 3
            t = 1.0;
        end
        Y = Y + t*U;
        Y3 = to3D(Y);
        for i = 1 : m
            % Orthonormalize the rows of Y3(:, :, i):
            [u, s, v] = svd(Y3(:, :, i), 'econ'); %#ok<ASGLU>
            Y3(:, :, i) = u*v';
            % Alternatively, one could also use qr_unique as retraction.
        end
        Y = to2D(Y3);
    end

    M.exp = @exponential;
    function Y = exponential(Y, U, t)
        if nargin == 2
            t = 1;
        end
        tU3 = multitransp(to3D(t*U));
        Y3 = multitransp(to3D(Y));
        % From a formula by Ross Lippert, Example 5.4.2 in AMS08.
        for i = 1 : m
            X = Y3(:, :, i);
            Z = tU3(:, :, i);
            Y3(:, :, i) = [X, Z] * ...
                          expm([  X'*Z , -Z'*Z ; eye(d) , X'*Z]) * ...
                          [ expm(-X'*Z) ; zeros(d) ];
            % We may loose orthonormality here. Just to be sure:
            [u, s, v] = svd(Y3(:, :, i), 'econ'); %#ok<ASGLU>
            Y3(:, :, i) = u*v';
        end
        Y = to2D(multitransp(Y3));
    end

    M.hash = @(Y) ['z' hashmd5(Y(:))];

    M.rand = @random;
    function Y = random()
        Y3 = zeros(d, k, m);
        for i = 1 : m
            [Q, unused] = qr(randn(k, d), 0); %#ok<ASGLU>
            Y3(:, :, i) = Q';
        end
        Y = to2D(Y3);
    end

    M.randvec = @randomvec;
    function U = randomvec(Y)
        U = projection(Y, randn(n, k));
        U = U / M.norm(Y, U);
    end

    M.lincomb = @matrixlincomb;

    M.zerovec = @(x) zeros(n, k);

    M.transp = @(x1, x2, u) projection(x2, u);

    M.vec = @(x, u_mat) u_mat(:);
    M.mat = @(x, u_vec) reshape(u_vec, [n, k]);
    M.vecmatareisometries = @() true;

end

author: NicolasBoumal

repo: manopt


text:
function rs=surrogate_test(s1,ntests,method,func)

%tstoolbox/@signal/surrogate_test
%   Syntax:
%     * rs=surrogate_test(s, ntests, method,func)
%
%   Input Arguments:
%     * s - has to be a real, scalar signal
%     * ntests - is the number of surrogate data sets to create
%     * method - method to generate surrogate data sets:
%          + 1: surrogate1
%          + 2: surrogate2
%          + 3: surrogate3
%     * func - string with matlab-code, have to return a signal object
%       with a scalar time series. The data to process is a signal object
%       referred by the qualifier s (see example).
%
%   Output Arguments:
%     * rs is a signal object with a three dimensional time series. The
%       first component is the result of the func function applied to the
%       original data set s. The second component is the mean of the
%       result of the func function applied to the ntests surrogate data
%       sets. The third component is the standard deviation. There is a
%       special plothint ('surrerrorbar') for the view function to show
%       this result in the common way.
%
%   surrogate_test runs an automatic surrogate data test task. It
%   generates ntests surrogate data sets an performs the func function to
%   each set. func is a string with matlab-code who returns a signal s
%   with a scalar time series.
%
%   Example:
%st = surrogate_test(s, 10, 1, 1, 'largelyap(embed(s,3,1,1), 128,20,10);');
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

s=s1;
s=eval(func);
x=[data(s)];
for tests=1:ntests
  switch(method)
   case 1
    s=surrogate1(s1);
   case 2
    s=surrogate2(s1);
   case 3
    s=surrogate3(s1);
   otherwise
    s=surrogate1(s1);
  end
  s=eval(func);
  x=[x data(s)];
end

a=[];
for test=1:length(x(:,1))
  a=[a; mean(x(test,2:end)) std(x(test,2:end))];
end

a=[x(:,1) a];

rs = signal(core(a),s1);	% special constructor calling syntax for working routines
%rs=setaxis(rs,1,a);
rs = addhistory(rs, ['Computed ' num2str(ntests) ' surrogate data function values']);
rs = setplothint(rs, 'surrerrorbar');
rs = addcommandlines(rs, 's = surrogate_test(s', ntests,method,func);

return

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [X,A,row,col] = reshape_hsi(I,A)
%RESHAPE_HSI Summary of this function goes here
%   Detailed explanation goes here
if nargin < 2
    A = [];
end

if ndims(I) == 3
    [row,col,B] = size(I);
    N = row*col;
    if ~isempty(A)
        [row,col,M] = size(A);
        A = reshape(A, N, M);
    end
    X = reshape(I, N, B);
else
    [N,B] = size(I);
    row = N;
    col = 1;
    X = I;
end

author: zhouyuanzxcv

repo: Hyperspectral


text:
% imagescloglog() - make an imagesc(0) plot with log y-axis and
%                   x-axis values
%
% Usage:  >> imagescloglog(times,freqs,data);
% Usage:  >> imagescloglog(times,freqs,data,clim,xticks,yticks,'key','val',...);
%
% Inputs:
%   times = vector of x-axis values (LOG spaced)
%   freqs = vector of y-axis values (LOG spaced)
%   data  = matrix of size (freqs,times)
%
% Optional inputs:
%   clim   = optional color limit
%   xticks = graduation for x axis
%   yticks = graduation for y axis
%   ...    = 'key', 'val' properties for figure
%
% Author: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 4/2003

% Copyright (C) 4/2003 Arnaud Delorme, SCCN/INC/UCSD, scott@sccn.ucsd.edu
%
% 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

function imagescloglog(times,freqs,data,clim, xticks, yticks, varargin)

  if size(data,1) ~= length(freqs)
      fprintf('logfreq(): data matrix must have %d rows!\n',length(freqs));
      return
  end
  if size(data,2) ~= length(times)
      fprintf('logfreq(): data matrix must have %d columns!\n',length(times));
      return
  end
  if min(freqs)<= 0
      fprintf('logfreq(): frequencies must be > 0!\n');
      return
  end
  try, icadefs; catch, warning('Using MATLAB default colormap'); end

  steplog = log(times(2))-log(times(1)); % same for all points
  realborders = [exp(log(times(1))-steplog/2) exp(log(times(end))+steplog/2)];
  newtimes    = linspace(realborders(1), realborders(2), length(times));

  % regressing 3 times
  border  = mean(newtimes(2:end)-newtimes(1:end-1))/2; % automatically added to the borders in imagesc
  newtimes = linspace(realborders(1)+border, realborders(2)-border, length(times));
  border  = mean(newtimes(2:end)-newtimes(1:end-1))/2; % automatically added to the borders in imagesc
  newtimes = linspace(realborders(1)+border, realborders(2)-border, length(times));
  border  = mean(newtimes(2:end)-newtimes(1:end-1))/2; % automatically added to the borders in imagesc
  newtimes = linspace(realborders(1)+border, realborders(2)-border, length(times));

  % problem with log images in Matlab: border are automatically added
  % to account for half of the width of a line: but they are added as
  % if the data was linear. The commands below compensate for this effect

  steplog = log(freqs(2))-log(freqs(1)); % same for all points
  realborders = [exp(log(freqs(1))-steplog/2) exp(log(freqs(end))+steplog/2)];
  newfreqs    = linspace(realborders(1), realborders(2), length(freqs));

  % regressing 3 times
  border  = mean(newfreqs(2:end)-newfreqs(1:end-1))/2; % automatically added to the borders in imagesc
  newfreqs = linspace(realborders(1)+border, realborders(2)-border, length(freqs));
  border  = mean(newfreqs(2:end)-newfreqs(1:end-1))/2; % automatically added to the borders in imagesc
  newfreqs = linspace(realborders(1)+border, realborders(2)-border, length(freqs));
  border  = mean(newfreqs(2:end)-newfreqs(1:end-1))/2; % automatically added to the borders in imagesc
  newfreqs = linspace(realborders(1)+border, realborders(2)-border, length(freqs));

  if nargin == 4 & ~isempty(clim)
      imagesc(newtimes,newfreqs,data,clim);
  else
      imagesc(newtimes,newfreqs,data);
  end;

  set(gca, 'yscale', 'log', 'xscale', 'log');
  try colormap(DEFAULT_COLORMAP); catch, end;

  % puting ticks
  % ------------
  if nargin >= 5
      divs = xticks;
  else
      divs = linspace(log(times(1)), log(times(end)), 10);
      divs = ceil(exp(divs)); divs = unique_bc(divs); % ceil is critical here, round might misalign
                                               % out-of border label with within border ticks
  end;
  set(gca, 'xtickmode', 'manual');
  set(gca, 'xtick', divs);
  if nargin >= 6
      divs = yticks;
  else
      divs = linspace(log(freqs(1)), log(freqs(end)), 10);
      divs = ceil(exp(divs)); divs = unique_bc(divs); % ceil is critical here, round might misalign
                                               % out-of border label with within border ticks
  end;
  set(gca, 'ytickmode', 'manual');
  set(gca, 'ytick', divs);

  % additional properties
  % ---------------------
  set(gca, 'yminortick', 'off', 'xaxislocation', 'bottom', 'box', 'off', 'ticklength', [0.03 0], 'tickdir','out', 'color', 'none');
  if ~isempty(varargin)
      set(gca, varargin{:});
  end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [DietFormulation] = readInDietFromVMH(fileNameDiet)
% This function reads in the diet that has been created and downloaded from
% the https://www.vmh.life/#nutrition and converts it into the whole-body
% metabolic model consistent format.
%
% [DietFormulation] = readInDietFromVMH(fileNameDiet)
%
% INPUT
% fileNameDiet      File name
%
% OUTPUT
% DietFormulation   Diet definition
%
% Ines Thiele 2016-2019

[Numbers, Strings] = xlsread(fileNameDiet{1});

ColFlux = 1;% assumes that fluxValues are given in 2nd col

DietNames = Strings(2:end,6); % assumes that Rxn names are given in 6th column

DietNames = regexprep(DietNames,'EX_','Diet_EX_');
DietNames = regexprep(DietNames,'\(e\)','\[d\]');
% Diet exchanges for all individuals
Diets = cellstr(num2str((Numbers(:,ColFlux))));

DietFormulation = [DietNames  Diets];

author: opencobra

repo: cobratoolbox


text:
function Q = randrot(n, N)
% Generates uniformly random rotation matrices.
%
% function Q = randrot(n, N)
%
% Q is an n-by-n-by-N array such that each slice Q(:, :, i) is a random
% orthogonal matrix of size n of determinant +1 (i.e., a matrix in SO(n)),
% sampled from the Haar measure (uniform distribution).
%
% By default, N = 1.
%
% Complexity: N times O(n^3).
% Theory in Diaconis and Shahshahani 1987 for the uniformity on O(n);
% With details in Mezzadri 2007,
% "How to generate random matrices from the classical compact groups."
%
% To ensure matrices in SO(n), we permute the two first columns when
% the determinant is -1.
%
% See also: randskew qr_unique randunitary

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Sept. 25, 2012.
% Contributors:
% Change log:
%   June 18, 2019 (NB)
%       Now generating all initial random matrices in one shot (which
%       should be more efficient) and calling qr_unique.

    if nargin < 2
        N = 1;
    end

    if n == 1
        Q = ones(1, 1, N);
        return;
    end

    % Generated as such, Q is uniformly distributed over O(n): the group
    % of orthogonal matrices; see Mezzadri 2007.
    Q = qr_unique(randn(n, n, N));

    for k = 1 : N

        % If a slice of Q is in O(n) but not in SO(n), we permute its two
        % first columns to negate its determinant. This ensures the new
        % slice is in SO(n), uniformly distributed.
        if det(Q(:, :, k)) < 0
            Q(:, [1 2], k) = Q(:, [2 1], k);
        end

    end

end


text:
function l = exp_length(x)
% Get the number of elements (parameters) in any expression.
% Length = exp_length(Expression)
%
% In:
%   Expression  : any expression
%
% Out:
%   Length      : the expression's length; highest allowed index allowed for exp_part.
%
% See also:
%   exp_head, exp_part, exp_parts
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-04-26

if isfield(x,{'head','parts'})
    % struct with fields .head & .parts
    l = length(x.parts);
elseif isfield(x,'tracking') && isfield(x.tracking,'expression')
    % impure expressions viewed as expressions
    l = exp_length(x.tracking.expression);
else
    % most MATLAB data structures are not further broken down when reflected into the expression
    % system
    l = 1;
end

repo: REST


text:
function [E]=EnergyOfLabelField(segmentation,potential,width,height,class_number)
n=size(segmentation,1);
segmentation=reshape(segmentation,[width height]);
Nei8=imstack2vectors(NeiX(segmentation));
E=zeros(n,class_number);
for i=1:class_number
    E(:,i)=sum(Nei8~=i,2);
end
E=E*potential;
end

author: Sable

repo: mcbench-benchmarks


text:
% LOTKA-VOLTERRA system
% System identification: DelayDMDc

clear all, close all, clc
figpath = '../FIGURES/';
datapath = '../DATA/';
addpath('../utils');

%% Generate Data
InputSignalType = 'sine2';%prbs; chirp; noise; sine2; sphs; mixed
Ndelay = 1;
getTrainingData

%% DMDc: B = unknown  and with time delay coordinates
if Ndelay == 1
    ModelName = 'DMDc';
elseif Ndelay>1
    ModelName = 'DelayDMDc';
end

numOutputs = size(Hx,1); numInputs = size(Hu,1); numVar = 2;
r1 = size(Hx,1); r2 = size(Hx,1);
[sysmodel_DMDc,U,Up] = DelayDMDc_MV(Hx,Hu,size(Hx,1),size(Hx,1),dt,size(Hx,1),size(Hu,1),2);

%% Prediction over training phase
[xDMDc,~] = lsim(sysmodel_DMDc,Hu,tspan(1:Nt),Hx(:,1));
xDMDc = xDMDc(:,end-1:end);
xDMDc = xDMDc + repmat(xmean,[Nt 1]);

%% Show validation
clear ph
figure,box on,
ccolors = get(gca,'colororder');
ph(1) = plot(tspan,x(:,1),'-','Color',ccolors(1,:),'LineWidth',1); hold on
ph(2) = plot(tspan,x(:,2),'-','Color',ccolors(2,:),'LineWidth',1);
ph(3) = plot(tspan(Ndelay:Nt+Ndelay-1),xDMDc(:,1),'--','Color',ccolors(1,:)-[0 0.2 0.2],'LineWidth',2);
ph(4) = plot(tspan(Ndelay:Nt+Ndelay-1),xDMDc(:,2),'--','Color',ccolors(2,:)-[0.1 0.2 0.09],'LineWidth',2);
xlim([0 100])
xlabel('Time')
ylabel('Population size')
legend(ph([1,3]),'True',ModelName)
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
print('-depsc2', '-loose', '-cmyk', [figpath,'EX_LOTKA_SI_',ModelName,'_',InputSignalType,'.eps']);

%% Prediction
% Reference
tspanV   = [100:dt:200];
xA      = xv;
tA      = tv;

% Model
if Ndelay == 1
    x0      = [x(end,1:2)];
    Hunew   = [u(end),uv(1:end)];
    [xB,tB] = lsim(sysmodel_DMDc,Hunew,tspanV,[x0-[xmean]]');
elseif Ndelay > 1
    x0      = [x(end-Ndelay+1,1:2),x(end,1:2)];
    Hunew   = [ u(end-Ndelay+1:end),uv(1:end-Ndelay);
        u(end),uv(1:end-1)];
    [xB,tB] = lsim(sysmodel_DMDc,Hunew,tspanV,[x0-[xmean]]');
    xB = xB(:,3:4); xB = xB + repmat(xmean,[size(xB,1) 1]);
end

xB = xB + repmat(xmean,[length(tB) 1]);

%% Show training and prediction
VIZ_SI_Validation

%% Save Data
Model.name = 'DelayDMDc';
Model.sys = sysmodel_DMDc;
Model.Ndelay = Ndelay;
Model.dt = dt;
save(fullfile(datapath,['EX_LOTKA_SI_',ModelName,'_',InputSignalType,'.mat']),'Model')

author: eurika-kaiser

repo: SINDY-MPC


text:
function populateStationsCombo(handles, hRefStationsCombo)
    maData = getappdata(handles.ma_MainGUI,'ma_data');

    if(isfield(maData.spacecraft,'stations'))
        stations = maData.spacecraft.stations;
        if(isempty(stations))
            set(hRefStationsCombo, 'Enable','off');
            set(hRefStationsCombo,'String', 'None Available');
        else
            set(hRefStationsCombo, 'Enable','on');

            stnStr = cell(0,1);
            for(i=1:length(stations)) %#ok<*NO4LP>
                stnStr{end+1} = stations{i}.name; %#ok<AGROW>
            end
            set(hRefStationsCombo,'String',stnStr);
        end
    else
        set(hRefStationsCombo, 'Enable','off');
        set(hRefStationsCombo,'String', 'None Available');
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% Replace zero values in a vector with arbitrary small values
%
% Description
%  This is a convenient function to avoid methods crashing because
%  of zero values in a signal.
%
% Input
%  x    : A vector containing possibly zero values.
%
% Output
%  y    : The same vector as x with zero values replaced.
%
% References
%  https://github.com/covarep/covarep/issues/75
%
% Copyright (c) 2011 University of Crete - Computer Science Department
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <gad27@cam.ac.uk>

function y = replacezeros(x)

    y = x;
    idx = find(x==0);
    y(idx) = eps*rand(1,length(idx));

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function [volume_pix_size] = mrGetVolPixSize(voldr,subject)
%function [volume_pix_size] = mrGetVolPixSize(voldr,subject)
%
% PURPOSE: Prompt the user for the volume pixel size
% AUTHOR:  Poirson
% DATE:    03.16.98
% HISTORY: Based on routine by Geoff Boynton
% NOTES:   The default values [256/240,256/240,1/1] mean that
%          24 cm of brain image are interpolated on the 256 pixels
%          in the  x and y dimensions (units are pixels/mm).
%          And that the volume plane thickness is 1 pixel per 1 mm.
%
% 09.17.99 RFD: cleaned up the code a bit.
% 04.01.01 ARW : Shouldn't that default 1mm  be a 124/150mm?
% See if we can just read the values
if nargin<2 % new way- full path is passed in
	paramsFile = [voldr,filesep,'UnfoldParams'];
else % old method- pass two variable that must be combined (why?)
  	paramsFile = [voldr,filesep,subject,filesep,'UnfoldParams'];
end

if ~exist([paramsFile '.mat'], 'file')
  disp ([paramsFile ' not found.']);
  disp ('Please enter values.');
  disp('Enter size of volume anatomy pixels/mm in x,y and z directions');
  volume_pix_size=input('Default is [256/240,256/240,124/150]: ');

  if isempty(volume_pix_size)
    volume_pix_size=[256/240,256/240,124/150];
  end

  %Create the file 'UnfoldParams.mat' in the subject's volume anatomy directory
  estr=(['save ' paramsFile ' volume_pix_size']);
  eval(estr);
else 	% Load in the infolding parameters (volme_pix_sizes)
  estr=(['load ' paramsFile]);
  eval(estr);
  volume_pix_size = volume_pix_size;
end

return

author: vistalab

repo: vistasoft


text:
%% Example of dop2Dbox Applied to Illumination Correction
%
% (C) 2013 Matthew harker and Paul O'Leary
% Institute for Automation
% University of Leoben
% A-8700 Leoben
% Austria
%
% email: office@harkeroleary.org
%
% This program demonstrates the application of discrete orthogonal
% polynomials to illumination correction in an image.
%
% This tool requires the discrete orthogonal toolbox: dopBox:
%
%   http://www.mathworks.com/matlabcentral/fileexchange/41250
%
%%
close all;
clear;
%
%% Load the Test Data
% \cellName{loadData}
%
% This is an image of a metallurgical sample taken with a microscope
%
load metallurgy;
%
fig1 = figure;
imagesc( D );
colormap(gray);
axis image;
%
%\caption{Example metallurgical image, note the inhomogeneous illumination}
%
%% Generate the Basis Functions for the X and Y Directions
% \cellName{genBases}
%
% Determine the size of the image
%
[ny, nx] = size( D );
%
% Define the degree of the basis functions for the x and y directions
%
degreeX = 2;
degreeY = 2;
%
% The function call dop uses the number of basis functions not degree
%
noBfsX = degreeX + 1;
noBfsY = degreeY + 1;
%
% Generate the discrete orthogonal basis functions
%
Bx = dop( nx, noBfsX );
By = dop( ny, noBfsY );
%
%% Compute the 2D Polynomial Spectrum
% \cellName{computeSpec}
%
% Only a few basis functions have been used. Consequently the 2D polynomial
% is of low degree. The fit will model the global intensity of the image.
%
Sp = By' * D * Bx;
%
% generate scale vectors for the x and y directions
%
xScale = 0:degreeX;
yScale = 0:degreeY;
%
fig2 = figure;
imagesc( xScale, yScale, Sp );
xlabel('Degree in $$x$$');
ylabel('Degree in $$y$$');
colorbar;
%
% \caption{2D polynomial spectrum}
%
%% Compute the Illumination Approximation
% \cellName{approx}
%
% Reconstructing the 2D data set corresponds to the global illumination
% approximation.
%
Z = By * Sp * Bx';
%
fig3 = figure;
imagesc( Z );
colorbar;
colormap(gray);
%
% \caption{2D Polynomial approximation to the image. The low degree of
% approximation ensures that the details in the image are not modelled.
% In this manner the approximation corresponds to the global illumination}
%
%% Correct the Image Intensity
% \cellName{correct}
%
C = D - Z;
%
fig4 = figure;
imagesc( C );
colorbar;
colormap(gray);
%
% \caption{Metallurgical image corrected for the non-uniform illumination.}


text:
% Test file for chebtech2 constructor.
% Here, we check populate().  (This function is not user-facing.)

function pass = test_constructor(pref)

% Get preferences:
if ( nargin < 1 )
    pref = chebtech.techPref();
end
% Set the tolerance:
tol = 100*pref.chebfuneps;

% Initialize with default data:
data = chebtech.parseDataInputs(struct());

%%
% Test on a scalar-valued function:
pref.extrapolate = 0;
pref.refinementFunction = 'nested';
f = @(x) sin(x);
g = populate(chebtech2, f, data, pref);
x = chebtech2.chebpts(length(g.coeffs));
values = g.coeffs2vals(g.coeffs);
pass(1) = norm(f(x) - values, inf) < tol;
pass(2) = abs(vscale(g) - sin(1)) < eps && g.ishappy && eps < tol;

pref.extrapolate = 1;
pref.refinementFunction = 'nested';
g = populate(chebtech2, f, data, pref);
x = chebtech2.chebpts(length(g.coeffs));
values = g.coeffs2vals(g.coeffs);
pass(3) = norm(f(x) - values, inf) < tol;
pass(4) = norm(vscale(g) - sin(1), inf) < tol && logical(eps);

pref.extrapolate = 0;
pref.refinementFunction = 'resampling';
g = populate(chebtech2, f, data, pref);
x = chebtech2.chebpts(length(g.coeffs));
values = g.coeffs2vals(g.coeffs);
pass(5) = norm(f(x) - values, inf) < tol;
pass(6) = abs(vscale(g) - sin(1)) < eps && logical(eps);

pref.extrapolate = 1;
pref.refinementFunction = 'resampling';
g = populate(chebtech2, f, data, pref);
x = chebtech2.chebpts(length(g.coeffs));
values = g.coeffs2vals(g.coeffs);
pass(7) = norm(f(x) - values, inf) < tol;
pass(8) = norm(vscale(g) - sin(1), inf) < tol && logical(eps);

%%
% Test on an array-valued function:
pref.extrapolate = 0;
pref.refinementFunction = 'nested';
f = @(x) [sin(x) cos(x) exp(x)];
g = populate(chebtech2, f, data, pref);
x = chebtech2.chebpts(length(g.coeffs));
values = g.coeffs2vals(g.coeffs);
pass(9) = norm(f(x) - values, inf) < tol;

pref.extrapolate = 1;
pref.refinementFunction = 'nested';
g = populate(chebtech2, f, data, pref);
x = chebtech2.chebpts(length(g.coeffs));
values = g.coeffs2vals(g.coeffs);
pass(10) = norm(f(x) - values, inf) < tol;

pref.extrapolate = 0;
pref.refinementFunction = 'resampling';
g = populate(chebtech2, f, data, pref);
x = chebtech2.chebpts(length(g.coeffs));
values = g.coeffs2vals(g.coeffs);
pass(11) = norm(f(x) - values, inf) < tol;

pref.extrapolate = 1;
pref.refinementFunction = 'resampling';
g = populate(chebtech2, f, data, pref);
x = chebtech2.chebpts(length(g.coeffs));
values = g.coeffs2vals(g.coeffs);
pass(12) = norm(f(x) - values, inf) < tol;

%%
% Some other tests:

% This should fail with an error:
try
    f = @(x) x + NaN;
    populate(chebtech2, f, data, pref);
    pass(13) = false;
catch ME
    pass(13) = strcmp(ME.message, 'Too many NaNs/Infs to handle.');
end

% As should this:
try
    f = @(x) x + Inf;
    populate(chebtech2, f, data, pref);
    pass(14) = false;
catch ME
    pass(14) = strcmp(ME.message, 'Too many NaNs/Infs to handle.');
end

% Test that the extrapolation option avoids endpoint evaluations.
pref.extrapolate = 1;
try
    populate(chebtech2, @(x) [F(x) F(x)], data, pref);
    pass(15) = true;
catch ME %#ok<NASGU>
    pass(15) = false;
end
pref.extrapolate = 0;

    function y = F(x)
        if ( any(abs(x) == 1) )
            error('Extrapolate should prevent endpoint evaluation.');
        end
        y = sin(x);
    end

% Check that things don't crash if pref.minSamples and pref.maxLength are equal.
try
    pref.minSamples = 8;
    pref.maxLength = 8;
    populate(chebtech2, @sin, data, pref);
    pass(16) = true;
catch
    pass(16) = false;
end

%%
% Test logical-valued functions:
f = chebtech2(@(x) x > -2);
g = chebtech2(1);
pass(17) = normest(f - g) < eps;

f = chebtech2(@(x) x < -2);
g = chebtech2(0);
pass(18) = normest(f - g) < eps;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun