text:
function t_opf_knitro(quiet)
%T_OPF_FMINCON  Tests for FMINCON-based optimal power flow.

%   MATPOWER
%   Copyright (c) 2004-2017, 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 nargin < 1
    quiet = 0;
end

%% current mismatch, cartesian V
options = {
    {0, 0},
    {0, 1},
    {1, 0},
    {1, 1},
};

num_tests = 224;

t_begin(length(options)*num_tests, quiet);

[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;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;

casefile = 't_case9_opf';
if quiet
    verbose = 0;
else
    verbose = 0;
end

mpopt = mpoption('opf.violation', 1e-9);
mpopt = mpoption(mpopt, 'out.all', 0, 'verbose', verbose, 'opf.ac.solver', 'KNITRO');
mpopt = mpoption(mpopt, 'knitro.tol_x', 1e-8, 'knitro.tol_f', 1e-8);

for k = 1:length(options)
    if options{k}{1}, bal = 'I';  else, bal = 'S'; end  %% nodal balance
    if options{k}{2}, crd = 'c';  else, crd = 'p'; end  %% V coordinates
    t0 = sprintf('Knitro OPF (%s,%s) : ', bal, crd);

    if ~have_feature('knitro')
        t_skip(num_tests, 'Artelys Knitro not available');
        continue;
    end

    mpopt = mpoption(mpopt, 'opf.current_balance',  options{k}{1}, ...
                            'opf.v_cartesian',      options{k}{2} );

    %% set up indices
    ib_data     = [1:BUS_AREA BASE_KV:VMIN];
    ib_voltage  = [VM VA];
    ib_lam      = [LAM_P LAM_Q];
    ib_mu       = [MU_VMAX MU_VMIN];
    ig_data     = [GEN_BUS QMAX QMIN MBASE:APF];
    ig_disp     = [PG QG VG];
    ig_mu       = (MU_PMAX:MU_QMIN);
    ibr_data    = (1:ANGMAX);
    ibr_flow    = (PF:QT);
    ibr_mu      = [MU_SF MU_ST];
    ibr_angmu   = [MU_ANGMIN MU_ANGMAX];

    %% get solved AC OPF case from MAT-file
    load soln9_opf;     %% defines bus_soln, gen_soln, branch_soln, f_soln

    %% run OPF
    for s = 0:3
        mpopt = mpoption(mpopt, 'opf.start', s);
        t = sprintf('%s(start=%d): ', t0, s);
        [baseMVA, bus, gen, gencost, branch, f, success, et] = runopf(casefile, mpopt);
        t_ok(success, [t 'success']);
        t_is(f, f_soln, 3, [t 'f']);
        t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
        t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
        t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
        t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
        t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
        t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
        t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
        t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
        t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
        t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);
    end
    mpopt = mpoption(mpopt, 'opf.start', 0);    %% set 'opf.start' back to default

    %% run with automatic conversion of single-block pwl to linear costs
    t = [t0 '(single-block PWL) : '];
    mpc = loadcase(casefile);
    mpc.gencost(2, NCOST) = 2;
    [r, success] = runopf(mpc, mpopt);
    [f, bus, gen, branch] = deal(r.f, r.bus, r.gen, r.branch);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);
    if mpopt.opf.v_cartesian
        xr = [r.var.val.Vr;r.var.val.Vi;r.var.val.Pg;r.var.val.Qg;0;r.var.val.y];
    else
        xr = [r.var.val.Va;r.var.val.Vm;r.var.val.Pg;r.var.val.Qg;0;r.var.val.y];
    end
    t_is(r.x, xr, 8, [t 'raw x returned from OPF']);

    %% get solved AC OPF case from MAT-file
    load soln9_opf_Plim;       %% defines bus_soln, gen_soln, branch_soln, f_soln

    %% run OPF with active power line limits
    t = [t0 '(P line lim) : '];
    mpopt1 = mpoption(mpopt, 'opf.flow_lim', 'P');
    [baseMVA, bus, gen, gencost, branch, f, success, et] = runopf(casefile, mpopt1);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    t = [t0 '(P^2 line lim) : '];
    mpopt1 = mpoption(mpopt, 'opf.flow_lim', '2');
    [baseMVA, bus, gen, gencost, branch, f, success, et] = runopf(casefile, mpopt1);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %%-----  test OPF with quadratic gen costs moved to generalized costs  -----
    mpc = loadcase(casefile);
    mpc.gencost = [
        2   1500    0   3   0.11    5   0;
        2   2000    0   3   0.085   1.2 0;
        2   3000    0   3   0.1225  1   0;
    ];
    [baseMVA, bus_soln, gen_soln, gencost, branch_soln, f_soln, success, et] = runopf(mpc, mpopt);
    branch_soln = branch_soln(:,1:MU_ST);

    A = sparse(0,0);
    l = [];
    u = [];
    nb = size(mpc.bus, 1);      % number of buses
    ng = size(mpc.gen, 1);      % number of gens
    thbas = 1;                thend    = thbas+nb-1;
    vbas     = thend+1;       vend     = vbas+nb-1;
    pgbas    = vend+1;        pgend    = pgbas+ng-1;
    qgbas    = pgend+1;       qgend    = qgbas+ng-1;
    nxyz = 2*nb + 2*ng;
    N = sparse((1:ng)', (pgbas:pgend)', mpc.baseMVA * ones(ng,1), ng, nxyz);
    fparm = [ 1    0   0 1;
              1 -100 100 1;
              1  -10  10 1 ];
    H = 2 * spdiags(mpc.gencost(:, 5), 0, ng, ng);
    Cw = mpc.gencost(:, 6);
    mpc.gencost(:, 5:7) = 0;

    %% run OPF with quadratic gen costs moved to generalized costs
    t = [t0 'w/quadratic generalized gen cost : '];
    [r, success] = opf(mpc, A, l, u, mpopt, N, fparm, H, Cw);
    [f, bus, gen, branch] = deal(r.f, r.bus, r.gen, r.branch);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);
    t_is(r.cost.usr, f, 12, [t 'user cost']);

    %%-----  run OPF with legacy costs and deadzone  -----
    if mpopt.opf.v_cartesian
        t_skip(17, 'legacy cost example n/a to cartesian V case')
    else
        load soln9_opf;
        mpc = loadcase(casefile);
        mpc.N = sparse((1:nb)', (vbas:vend)', ones(nb,1), nb, nxyz);
        mpc.fparm = ones(nb,1) * [ 2 1.08 0.02 1e8 ];
        mpc.Cw = ones(nb, 1);
        t = [t0 'w/legacy cost, in deadzone : '];
        r = runopf(mpc, mpopt);
        [f, bus, gen, branch] = deal(r.f, r.bus, r.gen, r.branch);
        t_ok(r.success, [t 'success']);
        t_is(f, f_soln, 3, [t 'f']);
        t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
        t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
        t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
        t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
        t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
        t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
        t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
        t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
        t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
        t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);
        t_is(r.cost.usr, 0, 12, [t 'user cost']);

        t = [t0 'w/legacy cost, not in deadzone : '];
        mpc.fparm = ones(nb,1) * [ 2 1.08 0.01 1e8 ];
        r = runopf(mpc, mpopt);
        [f, bus, gen, branch] = deal(r.f, r.bus, r.gen, r.branch);
        t_ok(r.success, [t 'success']);
        t_is(f, 9009.0890, 3, [t 'f']);
        t_is([min(bus(:, VM)) mean(bus(:, VM)) max(bus(:, VM))], ...
            [1.066624, 1.083980, 1.091698], 5, [t 'bus voltage']);
        t_is(r.cost.usr, 1673.065465, 5, [t 'user cost']);
    end

    %%-----  run OPF with extra linear user constraints & costs  -----
    %% single new z variable constrained to be greater than or equal to
    %% deviation from 1 pu voltage at bus 1, linear cost on this z
    %% get solved AC OPF case from MAT-file
    if mpopt.opf.v_cartesian
        t_skip(14, 'lin constraint/cost example n/a to cartesian V case')
    else
        load soln9_opf_extras1;   %% defines bus_soln, gen_soln, branch_soln, f_soln
        A = sparse([1;1;2;2],[10;25;10;25],[-1;1;1;1],2,25);
        u = [Inf; Inf];
        l = [-1; 1];

        N = sparse(1, 25, 1, 1, 25);    %% new z variable only
        fparm = [1 0 0 1];              %% w = r = z
        H = sparse(1,1);                %% no quadratic term
        Cw = 100;

        t = [t0 'w/extra constraints & costs 1 : '];
        [r, success] = opf(casefile, A, l, u, mpopt, N, fparm, H, Cw);
        [f, bus, gen, branch] = deal(r.f, r.bus, r.gen, r.branch);
        t_ok(success, [t 'success']);
        t_is(f, f_soln, 3, [t 'f']);
        t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
        t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
        t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
        t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
        t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
        t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
        t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
        t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
        t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
        t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);
        t_is(r.var.val.z, 0.025419, 6, [t 'user variable']);
        t_is(r.cost.usr, 2.5419, 4, [t 'user cost']);
    end

    %%-----  test OPF with capability curves  -----
    mpc = loadcase('t_case9_opfv2');
    %% remove angle diff limits
    mpc.branch(1, ANGMAX) = 360;
    mpc.branch(9, ANGMIN) = -360;

    %% get solved AC OPF case from MAT-file
    load soln9_opf_PQcap;   %% defines bus_soln, gen_soln, branch_soln, f_soln

    %% run OPF with capability curves
    t = [t0 'w/capability curves : '];
    [baseMVA, bus, gen, gencost, branch, f, success, et] = runopf(mpc, mpopt);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %%-----  test OPF with angle difference limits  -----
    mpc = loadcase('t_case9_opfv2');
    %% remove capability curves
    mpc.gen(2:3, [PC1, PC2, QC1MIN, QC1MAX, QC2MIN, QC2MAX]) = zeros(2,6);

    %% get solved AC OPF case from MAT-file
    load soln9_opf_ang;   %% defines bus_soln, gen_soln, branch_soln, f_soln

    %% run OPF with angle difference limits
    t = [t0 'w/angle difference limits : '];
    [baseMVA, bus, gen, gencost, branch, f, success, et] = runopf(mpc, mpopt);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  1, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);
    t_is(branch(:,ibr_angmu ), branch_soln(:,ibr_angmu ),  2, [t 'branch angle mu']);

    %%-----  test OPF with ignored angle difference limits  -----
    %% get solved AC OPF case from MAT-file
    load soln9_opf;   %% defines bus_soln, gen_soln, branch_soln, f_soln

    %% run OPF with ignored angle difference limits
    t = [t0 'w/ignored angle difference limits : '];
    mpopt1 = mpoption(mpopt, 'opf.ignore_angle_lim', 1);
    [baseMVA, bus, gen, gencost, branch, f, success, et] = runopf(mpc, mpopt1);
    %% ang limits are not in this solution data, so let's remove them
    branch(1, ANGMAX) = 360;
    branch(9, ANGMIN) = -360;
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %% angle bounded above by 0, unbounded below
    %% for issue/18
    t = [t0 'w/angle difference limit = 0 : '];
    mpc = loadcase(casefile);
    b = 5;
    mpc.branch(b, ANGMAX) = 0;
    r = runopf(mpc, mpopt);
    t_ok(success, [t 'success']);
    diff = r.bus(r.branch(b, F_BUS), VA) - r.bus(r.branch(b, T_BUS), VA);
    t_is(diff, 0, 5, [t 'angle diff']);

    %%-----  OPF with ref bus not = bus 1, ref angle not = 0  -----
    t = [t0 'ref bus ~= 1, ref ang ~= 0 : '];
    mpc = loadcase(casefile);
    mpc.bus([1;3], BUS_TYPE) = [PV; REF];   %% swap bus types
    bus_soln([1;3], BUS_TYPE) = bus_soln([3;1], BUS_TYPE);   %% swap bus types
    mpc.bus(3, VA) = 3.3014277;
    r = runopf(mpc, mpopt);
    [success, f, bus, gen, branch] = deal(r.success, r.f, r.bus, r.gen, r.branch);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %%-----  test OPF with opf.use_vg  -----
    %% get solved AC OPF case from MAT-file
    load soln9_opf_vg;  %% defines bus_soln, gen_soln, branch_soln, f_soln
                        %%    and bus_soln1, gen_soln1, branch_soln1, f_soln1

    %% run with opf.use_vg = 1
    t = [t0 'w/opf.use_vg = 1 : '];
    mpc = loadcase(casefile);
    mpc.gen = mpc.gen([1 2 1 3], :);
    mpc.gencost = mpc.gencost([1 2 1 3], :);
    mpc.gen([1 3], [PMAX PMIN]) = mpc.gen([1 3], [PMAX PMIN]) / 2;
    mpc.gen(3, [QMIN, QMAX]) = 0;   %% no reactive capability for gen 3
    mpc.gencost([1 3], COST:end) = mpc.gencost([1 3], COST:end) / 2;
    mpc.gen(1, VG) = 1.05;
    mpc.gen(3, VG) = 1.06;
    mpopt1 = mpoption(mpopt, 'opf.use_vg', 1);
    r = runopf(mpc, mpopt1);
    t_ok(r.success, [t 'success']);
    t_is(r.f, f_soln, 3, [t 'f']);
    t_is(   r.bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   r.bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   r.bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   r.bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   r.gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   r.gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   r.gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(r.branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(r.branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(r.branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %% run with opf.use_vg = 0.9
    t = [t0 'w/opf.use_vg = 0.9 : '];
    mpopt1 = mpoption(mpopt, 'opf.use_vg', 0.9);
    r = runopf(mpc, mpopt1);
    t_ok(r.success, [t 'success']);
    t_is(r.f, f_soln1, 3, [t 'f']);
    t_is(   r.bus(:,ib_data   ),    bus_soln1(:,ib_data   ), 10, [t 'bus data']);
    t_is(   r.bus(:,ib_voltage),    bus_soln1(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   r.bus(:,ib_lam    ),    bus_soln1(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   r.bus(:,ib_mu     ),    bus_soln1(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   r.gen(:,ig_data   ),    gen_soln1(:,ig_data   ), 10, [t 'gen data']);
    t_is(   r.gen(:,ig_disp   ),    gen_soln1(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   r.gen(:,ig_mu     ),    gen_soln1(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(r.branch(:,ibr_data  ), branch_soln1(:,ibr_data  ), 10, [t 'branch data']);
    t_is(r.branch(:,ibr_flow  ), branch_soln1(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(r.branch(:,ibr_mu    ), branch_soln1(:,ibr_mu    ),  2, [t 'branch mu']);

    t = [t0 'hi-deg polynomial Pg costs : '];
    mpc = loadcase(casefile);
    mpc.gencost = [
        2   1500    0   6   1e-6/5  0   0   0   0   0;
        2   3000    0   5   1e-4/4  0   0   0   0   0;
        2   2000    0   3   1/2     0   0   0   0   0;
    ];
    r = runopf(mpc, mpopt);
    [f, bus, gen, branch] = deal(r.f, r.bus, r.gen, r.branch);
    t_ok(r.success, [t 'success']);
    t_is(f, 11899.4652, 4, [t 'f']);
    t_is(gen(:, PG), [100.703628; 88.719864; 128.679485], 5, [t 'Pg']);
    t_is(gen(:, QG), [8.151114; -27.105856; 12.471647], 5, [t 'Qg']);
    t_is([min(bus(:, VM)) mean(bus(:, VM)) max(bus(:, VM))], ...
        [1.059191 1.079404 1.1], 5, [t 'bus voltage']);

    t = [t0 'hi-deg polynomial Qg costs : '];
    mpc = loadcase(casefile);
    mpc.gencost = [
        2   1500    0   6   1e-6/5  0   0   0   0   0;
        2   3000    0   5   1e-4/4  0   0   0   0   0;
        2   2000    0   3   1/2     0   0   0   0   0;
        2   0       0   6   1e-6/5  0   0   0   0   0;
        2   0       0   3   1/2     0   0   0   0   0;
        2   0       0   5   1e-4/4  0   0   0   0   0;
    ];
    r = runopf(mpc, mpopt);
    [f, bus, gen, branch] = deal(r.f, r.bus, r.gen, r.branch);
    t_ok(r.success, [t 'success']);
    t_is(f, 11998.8345, 4, [t 'f']);
    t_is(gen(:, PG), [100.853092; 89.258007; 128.518423], 5, [t 'Pg']);
    t_is(gen(:, QG), [-20.834352; -7.966093; 34.591995], 5, [t 'Qg']);
    t_is([min(bus(:, VM)) mean(bus(:, VM)) max(bus(:, VM))], ...
        [1.007729 1.050776 1.1], 5, [t 'bus voltage']);

    %% OPF with user-defined nonlinear constraints
    t = [t0 'w/nonlin eq constraint : '];
    mpc = loadcase('case30');
    mpc.user_constraints.nle = {
        {'Pg_usr', 1, 'opf_nle_fcn1', 'opf_nle_hess1', {'Pg'}, {}}
    };
    r = runopf(mpc, mpopt);
    t_ok(r.success, [t 'success']);
    t_is(r.gen(1, PG) * r.gen(2, PG) / 100, r.gen(6, PG), 8, [t 'Pg(1)*Pg(2)=Pg(6)']);
    t_is(r.gen(6, PG), 20.751163, 5, [t 'Pg(6)']);
    %% check for indexing error, issue #77
    t_is(r.raw.pimul(61), 0, 6, [t 'issue #77']);
    t_is(r.mu.nln.u(61), 0, 6, [t 'issue #77']);

    %% OPF with all buses isolated
    t = [t0 'all buses isolated : '];
    mpc.bus(:, BUS_TYPE) = NONE;
    try
        r = runopf(mpc, mpopt);
        t_is(r.success, 0, 12, [t 'success = 0']);
    catch
        t_ok(0, [t 'unexpected fatal error']);
    end

    %% OPF with no branch limits
    t = [t0 'w/no branch limits : '];
    mpc = loadcase(casefile);
    mpc.branch(:, RATE_A) = 0;
    r = runopf(mpc, mpopt);
    t_ok(r.success, [t 'success']);
    t_is(r.f, 5496.128635, 4, [t 'f']);
    t_is(r.gen(:, PG), [90; 220.463932; 10], 5, [t 'Pg']);
    t_is([min(r.bus(:, VM)) mean(r.bus(:, VM)) max(r.bus(:, VM))], ...
        [1.070692 1.090449 1.1], 5, [t 'bus voltage']);
end

if have_feature('octave')
    warning(s1.state, file_in_path_warn_id);
end

t_end;

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [Var_Rel,Var_Abs]=CoreVarn(C)
%COREVARN  Calculates the 'core-variance'
%
%[Var_Rel,Var_Abs]=CoreVarn(C)
%Calculates the 'core-variance' using not means of
%variables(columns) but the avg of all elements in the matrix
%
%Var_Rel : Variance of squares in percent
%Var_Abs : Absolute variance of squares

% Copyright (C) 1995-2006  Rasmus Bro & Claus Andersson
% Copenhagen University, DK-1958 Frederiksberg, Denmark, rb@life.ku.dk
%
% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 2 of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
% You should have received a copy of the GNU General Public License along with
% this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
% Street, Fifth Floor, Boston, MA  02110-1301, USA.

C=C.^2;
p=prod(size(C));
mn=sum(C(:))/p;

Var_Abs=sum(sum( (C-mn).^2) );

Var_Rel=100*Var_Abs/(p*(p-1)*mn^2);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function r = betarnd(a,b,m,n);
%BETARND Random matrices from beta distribution.
%   R = BETARND(A,B) returns a matrix of random numbers chosen
%   from the beta distribution with parameters A and B.
%   The size of R is the common size of A and B if both are matrices.
%   If either parameter is a scalar, the size of R is the size of the other
%   parameter. Alternatively, R = BETARND(A,B,M,N) returns an M by N matrix.

%   Reference:
%      [1]  L. Devroye, "Non-Uniform Random Variate Generation",
%      Springer-Verlag, 1986

%   Copyright (c) 1993-98 by The MathWorks, Inc.
%   $Revision$  $Date: 2004/02/10 18:58:56 $

if nargin < 2,
    error('Requires at least two input arguments');
end

if nargin == 2
    [errorcode rows columns] = rndcheck(2,2,a,b);
end

if nargin == 3
    [errorcode rows columns] = rndcheck(3,2,a,b,m);
end

if nargin == 4
    [errorcode rows columns] = rndcheck(4,2,a,b,m,n);
end

if errorcode > 0
    error('Size information is inconsistent.');
end

r = zeros(rows,columns);

% Use Theorem 4.1, case A (Devroye, page 430) to derive beta
%   random numbers as a ratio of gamma random numbers.
if prod(size(a)) == 1
    a1 = a(ones(rows,1),ones(columns,1));
    g1 = gamrnd(a1,1);
else
    g1 = gamrnd(a,1);
end
if prod(size(b)) == 1
    b1 = b(ones(rows,1),ones(columns,1));
    g2 = gamrnd(b1,1);
else
    g2 = gamrnd(b,1);
end
r = g1 ./ (g1 + g2);

% Return NaN if b is not positive.
if any(any(b <= 0));
    if prod(size(b) == 1)
        tmp = NaN;
        r = tmp(ones(rows,columns));
    else
        k = find(b <= 0);
        tmp = NaN;
        r(k) = tmp(ones(size(k)));
    end
end

% Return NaN if a is not positive.
if any(any(a <= 0));
    if prod(size(a) == 1)
        tmp = NaN;
        r = tmp(ones(rows,columns));
    else
        k = find(a <= 0);
        tmp = NaN;
        r(k) = tmp(ones(size(k)));
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function varargout = process_evt_read( varargin )
% PROCESS_EVT_READ: Read the values from one/several stim or response channels and detect triggers.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2012-2021
%          Raymundo Cassani, 2022
%          Marc Lalancette, 2022

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Read from channel';
    sProcess.Category    = 'File';
    sProcess.SubGroup    = 'Events';
    sProcess.Index       = 40;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/StimDelays?highlight=%28Read+events+from+channel%29#Detection_of_the_button_responses';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'raw'};
    sProcess.OutputTypes = {'data', 'raw'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Option: Stim channels
    sProcess.options.stimchan.Comment = 'Event channels: ';
    sProcess.options.stimchan.Type    = 'text';
    sProcess.options.stimchan.Value   = '';
    % Option: Value/Bit
    sProcess.options.trackmode.Comment = {'Value: detect the changes of channel value', ...
                                          'Bit: detect the changes for each bit independently', ...
                                          'TTL: detect peaks of 5V/12V on an analog channel (baseline=0V)', ...
                                          'RTTL: detect peaks of 0V on an analog channel (baseline!=0V)'};
    sProcess.options.trackmode.Type    = 'radio';
    sProcess.options.trackmode.Value   = 1;
    % Option: Accept zeros
    sProcess.options.zero.Comment = 'Accept zeros as trigger values';
    sProcess.options.zero.Type    = 'checkbox';
    sProcess.options.zero.Value   = 0;
    % Option: Min event duration
    sProcess.options.min_duration.Comment = 'Reject events shorter than: ';
    sProcess.options.min_duration.Type    = 'value';
    sProcess.options.min_duration.Value   = {0, 'samples', 0};
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInput) %#ok<DEFNU>
    % Return all the input files
    OutputFiles = {};

    % ===== GET OPTIONS =====
    % Stim channel not specified
    StimChan = sProcess.options.stimchan.Value;
    if isempty(StimChan)
        bst_report('Error', sProcess, [], 'Event channel was not specified.');
        return
    end
    % Event type
    switch (sProcess.options.trackmode.Value)
        case 1,  EventsTrackMode = 'value';
        case 2,  EventsTrackMode = 'bit';
        case 3,  EventsTrackMode = 'ttl';
        case 4,  EventsTrackMode = 'rttl';
    end
    % Other options
    isAcceptZero = sProcess.options.zero.Value;
    MinDuration = sProcess.options.min_duration.Value{1};

    % ===== GET FILE DESCRIPTOR =====
    % Load the raw file descriptor
    isRaw = strcmpi(sInput.FileType, 'raw');
    if isRaw
        DataMat = in_bst_data(sInput.FileName, 'F');
        sFile = DataMat.F;
    else
        sFile = in_fopen(sInput.FileName, 'BST-DATA');
    end
    % Load channel file
    ChannelMat = in_bst_channel(sInput.ChannelFile);
    % Check if specified channels are available
    iChannels = channel_find(ChannelMat.Channel, StimChan);
    if isempty(iChannels)
        bst_report('Error', sProcess, sInput, ['Channel name(s) "' StimChan '" does not exist.']);
        return
    end

    % ===== DETECTION =====
    % CTF: Read separately upper and lower bytes
    if ismember(sFile.format, {'CTF', 'CTF-CONTINUOUS'})
        % Detect separately events on the upper and lower bytes of the STIM channel
        eventsU = Compute(sFile, ChannelMat, [StimChan '__U'], EventsTrackMode, isAcceptZero, MinDuration);
        eventsL = Compute(sFile, ChannelMat, [StimChan '__L'], EventsTrackMode, isAcceptZero, MinDuration);
        % If there are events on both: add marker U/L
        if ~isempty(eventsU) && ~isempty(eventsL)
            for iEvt = 1:length(eventsU)
                eventsU(iEvt).label = ['U', eventsU(iEvt).label];
            end
            for iEvt = 1:length(eventsL)
                eventsL(iEvt).label = ['L', eventsL(iEvt).label];
            end
        end
        events = [eventsL, eventsU];
    else
        events = Compute(sFile, ChannelMat, StimChan, EventsTrackMode, isAcceptZero, MinDuration);
    end

    % ===== SAVE RESULT =====
    % Only save changes if something was change
    if ~isempty(events)
        % Import new events in file structure
        sFile = import_events(sFile, [], events);
        % Report changes in .mat structure
        if isRaw
            DataMat.F = sFile;
        else
            DataMat.Events = sFile.events;
        end
        % Save file definition
        bst_save(file_fullpath(sInput.FileName), DataMat, 'v6', 1);
        % Report number of detected events
        bst_report('Info', sProcess, sInput, sprintf('%s: %d events detected in %d categories', StimChan, size([events.times],2), length(events)));
    end
    % Return all the input files
    OutputFiles{end+1} = sInput.FileName;
end

%% ===== COMPUTE =====
function [events, EventsTrackMode, StimChan] = Compute(sFile, ChannelMat, StimChan, EventsTrackMode, isAcceptZero, MinDuration)
    % Parse inputs
    if (nargin < 6)
        MinDuration = 0;
    end
    if (nargin < 5)
        isAcceptZero = 0;
    end
    if (nargin < 4)
        EventsTrackMode = 'ask';
    end
    if (nargin < 3)
        StimChan = [];
    end
    % Get some information
    ch_names = {ChannelMat.Channel.Name};
    samplesBounds = round(sFile.prop.times .* sFile.prop.sfreq);
    events = [];

    % ===== GET STIM CHANNEL =====
    if isempty(StimChan) || isequal(StimChan, 'ask') || isequal(StimChan, {'ask'})
        % Find possible event channels
        switch (sFile.format)
            case 'FIF'
                % Channel name must contains 'STI'
                iStiChan = find(~cellfun(@isempty, strfind(ch_names, 'STI')));
            case '4D'
                % TRIGGER Channels
                iStiChan = channel_find(ChannelMat.Channel, 'Stim');
            case {'CTF', 'CTF-CONTINUOUS'}
                % TRIGGER Channels
                iStiChan = channel_find(ChannelMat.Channel, 'Stim');
            case 'KDF'
                % Status Channel
                iStiChan = channel_find(ChannelMat.Channel, 'KDF');
            case 'EEG-BDF'
                % Status Channel
                iStiChan = channel_find(ChannelMat.Channel, 'BDF');
        end
        % No valid events channel
        if isempty(iStiChan)
            bst_error('No valid events channel found in this file.', 'Read events', 0);
            return
        end
        % If only one choice: select it by default
        if (length(iStiChan) == 1)
            StimChan = ch_names(iStiChan);
        % Else: offer multiple choices to the user
        else
            StimChan = java_dialog('checkbox', ...
                ['You can try to rebuild the events list using one or more technical tracks, or ' 10 ...
                 'ignore this step and process the file as continuous recordings without events.' 10 10 ...
                 'Available technical tracks: '], ...
                 'Read events', [], ch_names(iStiChan));
            if isempty(StimChan)
                events = [];
                return
            end
        end
    end
    % CTF: Select only upper or lower bytes
    isCtfUp = 0;
    isCtfLow = 0;
    if ismember(sFile.format, {'CTF', 'CTF-CONTINUOUS'}) && (length(StimChan) > 3)
        if strcmp(StimChan(end-2:end), '__U')
            isCtfUp = 1;
            StimChan = StimChan(1:end-3);
        elseif strcmp(StimChan(end-2:end), '__L')
            isCtfLow = 1;
            StimChan = StimChan(1:end-3);
        end
    end

    % ===== ASK READ MODE =====
    if strcmpi(EventsTrackMode, 'ask')
        res = java_dialog('question', ['Please select the interpretation mode at each time sample:

' ...
                                       '- Value: detect the changes of value on the trigger channel
' ...
                                       '- Bit: detect the changes for each bit of the channel independently
' ...
                                       '- TTL: detect peaks of 5V/12V on an analog channel (baseline=0V)
' ...
                                       '- RTTL: detect peaks of 0V on an analog channel (baseline!=0V)
' ...
                                       '- Ignore: do not read trigger channel

'], ...
                                       'Type of events', [], {'Value','Bit','TTL','RTTL','Ignore','Cancel'}, 'value');
        if isempty(res) || strcmpi(res, 'Cancel')
            events = -1;
            return
        end
        EventsTrackMode = lower(res);
    end
    % Ignore trigger channel
    if strcmpi(EventsTrackMode, 'ignore')
        return;
    end

    % ===== READ STIM CHANNELS =====
    % Intialize returned events structure
    events = repmat(db_template('event'), 0);
    % Get channel indices
    iChannels = channel_find(ChannelMat.Channel, StimChan);
    if isempty(iChannels)
        disp(['EVENTS> Error: Channel name(s) "' StimChan '" does not exist.']);
    end
    StimChan = {ChannelMat.Channel(iChannels).Name};
    % Define optimal block size for reading
    blockLength = 6000;
    % Adapt block size to FIF block size
    if strcmpi(sFile.format, 'FIF') && isfield(sFile.header, 'raw') && isfield(sFile.header.raw, 'rawdir') && ~isempty(sFile.header.raw.rawdir)
        fifBlockSize = double(sFile.header.raw.rawdir(1).nsamp);
        blockLength = fifBlockSize * max(1, round(blockLength / fifBlockSize));
    end
    % Process by blocks (if not: out of memory)
    totalLength = (samplesBounds(2) - samplesBounds(1) + 1);
    nbBlocks = ceil(totalLength / blockLength);
    % Progress bar, channels
    isProgressBar = bst_progress('isVisible');
    bst_progress('start', 'Import events', 'Reading events channels...', 0, length(iChannels) * nbBlocks);

    tracks_name = {};
    tracks_vals = {};
    tracks_smps = {};
    % Process each channel
    for iChannel = 1 : length(iChannels)
        % Increment progress bar
        track_prev = [];
        track_vals = [];
        track_smps = [];
        % Progress bar, blocks
        bst_progress('text', ['Reading events in channel ' StimChan{iChannel} ' ...']);
        for iBlock = 1:nbBlocks
            % Increment progress bar
            bst_progress('inc', 1);
            % === READ BLOCK ===
            % Get samples indices for this block
            samplesBlock = samplesBounds(1) + [(iBlock - 1) * blockLength, iBlock * blockLength - 1];
            samplesBlock(2) = min(samplesBlock(2), samplesBounds(2));
            % Read block of data
            [track, ~] = in_fread(sFile, ChannelMat, 1, samplesBlock, iChannels(iChannel));
            % Round values
            track = fix(track);
            % Keep only track changes
            if iBlock == 1
                track_smps(1) = 1;
                track_vals(1) = track(1);
                track_prev = track(1);
            end
            % === KEEP ALL CHANGES ===
            ixDiffTrack = find(diff([track_prev, track]));
            track_vals = [track_vals, track(ixDiffTrack)];
            track_smps = [track_smps, samplesBlock(1) + ixDiffTrack];
            % Saving ending state for next block
            track_prev = track(end);
        end

        % CTF: Read separately Upper and Lower bytes
        if ismember(sFile.format, {'CTF', 'CTF-CONTINUOUS'})
            % Events are read as int32, while they are actually uint32: fix negative values
            track_vals(track_vals<0) = track_vals(track_vals<0) + 2^32;
            % Keep only upper or lower bytes
            if isCtfUp
                track_vals = fix(track_vals / 2^16);
            elseif isCtfLow
                track_vals = double(bitand(uint32(track_vals), 2^16-1));
            end
        % Other formats
        else
            % Old code: Might be useless and/or detrimental to the reading of the events
            % tracks = reshape(double(typecast(int16(tracks(:)), 'uint16')), size(tracks));
            % Use the same fix as for CTF files
            track_vals(track_vals<0) = track_vals(track_vals<0) + 2^32;
        end

        % === SEPARATE TRACKS ===
        % Each bit of each channel is interpreted as a track
        switch lower(EventsTrackMode)
            case 'bit'
                % Convert track in binary values
                track_bin = double(fliplr(dec2bin(track_vals))');
                track_bin(track_bin == '0') = 0;
                track_bin(track_bin == '1') = 1;
                nBit = size(track_bin, 1);
                tracks_bit_name = cell(1, nBit);
                tracks_bit_vals = cell(1, nBit);
                tracks_bit_smps = cell(1, nBit);
                % Data for each binary track
                for iBit = 1 : nBit
                    if (length(StimChan) > 1)
                        tracks_bit_name{iBit} = sprintf('%s_%d', StimChan{iChannel}, iBit);
                    else
                        tracks_bit_name{iBit} = sprintf('%d', iBit);
                    end
                    % Keep changes in each binary tracks
                    ixDiffTrackBit = find(diff([track_prev, track_bin(iBit, :)]));
                    tracks_bit_vals{iBit} = track_bin(iBit, ixDiffTrackBit);
                    tracks_bit_smps{iBit} = track_smps(ixDiffTrackBit);
                end
                % Append to other tracks
                tracks_name = [tracks_name, tracks_bit_name];
                tracks_vals = [tracks_vals, tracks_bit_vals];
                tracks_smps = [tracks_smps, tracks_bit_smps];
            case 'value'
                tracks_name{end+1} = StimChan{iChannel};
                tracks_vals{end+1} = track_vals;
                tracks_smps{end+1} = track_smps;
            case {'ttl', 'rttl'}
                tracks_name{end+1} = StimChan{iChannel};
                track_vals = abs(round(track_vals));
                track_vals(track_vals ~= 0) = 1;
                tracks_vals{end+1} = track_vals;
                tracks_smps{end+1} = track_smps;
        end
    end

    % === PROCESS EACH TRACK SEPARATELY ===
    nTooShort = 0;
    for iTrack = 1:length(tracks_name)
        % Get the indices where something happens
        if strcmpi(EventsTrackMode, 'rttl')
            ixs = find(tracks_vals{iTrack} == 0);
        elseif isAcceptZero
            ixs = 1 : length(tracks_vals{iTrack});
        else
            ixs = find(tracks_vals{iTrack} ~= 0);
        end
        % Process each change individually
        for i = 1:length(ixs)
            if i < length(ixs)
                % Skip if shorter than MinDuration
                duration_i = tracks_smps{iTrack}(ixs(i)+1) - tracks_smps{iTrack}(ixs(i));
                if (MinDuration > 0) && (duration_i < MinDuration)
                    nTooShort = nTooShort + 1;
                    % Pass start sample to next event longer than MinDuration
                    tracks_smps{iTrack}(ixs(i)+1) = tracks_smps{iTrack}(ixs(i));
                    continue;
                end
            end
            % Build event name
            switch lower(EventsTrackMode)
                case 'bit'
                    label = tracks_name{iTrack};
                    if (MinDuration == 0) && isAcceptZero
                        if tracks_vals{iTrack}(ixs(i)) == 1
                            label = [label '-set'];
                        else
                            label = [label '-reset'];
                        end
                    end
                case 'value'
                    value = tracks_vals{iTrack}(ixs(i));
                    if (length(StimChan) > 1)
                        label = sprintf('%s_%d', tracks_name{iTrack}, value);
                    else
                        label = sprintf('%d', value);
                    end
                case {'ttl', 'rttl'}
                    label = tracks_name{iTrack};
            end
            % Find this event in list of events
            if ~isempty(events)
                iEvent = find(strcmpi({events.label}, label));
            else
                iEvent = [];
            end
            % If event does not exist yet: add it
            if isempty(iEvent)
                iEvent = length(events) + 1;
                events(iEvent).label      = label;
                events(iEvent).epochs     = [];
                events(iEvent).times      = [];
                events(iEvent).reactTimes = [];
                events(iEvent).select     = 1;
                events(iEvent).channels   = [];
                events(iEvent).notes      = [];
            end
            % Add occurrence of this event
            iOcc = length(events(iEvent).times) + 1;
            events(iEvent).epochs(iOcc)   = 1;
            events(iEvent).times(iOcc)    = (tracks_smps{iTrack}(ixs(i))-1) ./ sFile.prop.sfreq;
        end
    end
    % Display warning with removed events
    if (nTooShort > 0)
        disp(sprintf('BST> %d events shorter than %d sample(s) removed.', nTooShort, MinDuration));
    end
    % Close progress bar
    if ~isProgressBar
        bst_progress('stop');
    end
end

author: brainstorm-tools


text:
function out = iss_enc(signal, text, L_min, alpha, beta, gama)
%ISS_ENC is the function to hide data within audio using Improved
%   Spread Spectrum technique in time domain.
%
%   INPUTS VARIABLES
%       signal : Cover signal
%       text   : Message to hide
%       L_min  : Minimum segment length
%       alpha  : Embedding parameter 1
%       beta   : Embedding parameter 2
%       gama   : Embedding parameter 3
%
%   OUTPUTS VARIABLES
%       out    : Stego signal
%
%   Kadir Tekeli (kadir.tekeli@outlook.com)

if nargin < 3
    L_min = 1024;  % Setting a minimum value for segment length
end

if nargin < 4
    alpha = 0.005; beta = 1; gama = 0.5;
end

[s.len,s.ch] = size(signal);

bit = getBits(text);
L2  = floor(s.len/length(bit));
L   = max(L_min, L2);              %Length of segments
n   = floor(s.len/L);
N   = n - mod(n, 8);               %Number of segments

if (length(bit) > N)
    warning('Message is too long, being cropped...');
    bits=bit(1:N);
else
    bits=[bit, num2str(zeros(N-length(bit),1))'];
end

r = ones(L,1);
%r  = prng('password123',L);
pr = reshape(r*ones(1,N),N*L,1);
u  = alpha * pr;

out = signal;
ss = mixer(L,bits,-1,1);
iss = iss_improve(signal, N, L, alpha, r);

%%%%%%%%%%%%%%%%%%%%%%% EMBEDDING MESSAGE... %%%%%%%%%%%%%%%%%%%%%%%%
stego = signal(1:N*L,1) + (beta*ss - gama*iss).*u;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
out(:,1) = [stego; signal(N*L+1:s.len,1)];
end

function out = prng( key, L )
    pass = sum(double(key).*(1:length(key)));
    rand('seed', pass);
    out = 2*(rand(L, 1)>0.5)-1;
end

author: ktekeli

repo: audio-steganography-algorithms


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  battistn@tcnj.edu
% Date Simulation Created: August 28, 2019
% Institution: TCNJ
%
% Date IB2d Was Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%   .
%   .
%   .
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific fiber model or example, please let Nick (battistn@tcnj.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: updates the spring attributes!
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function springs_info = update_Springs(dt,current_time,xLag,yLag,springs_info)

%springs_info: col 1: starting spring pt (by lag. discretization)
%              col 2: ending spring pt. (by lag. discretization)
%              col 3: spring stiffness
%              col 4: spring resting lengths

% Contraction Frequency
freq = 0.8;

% CHANGE RESTING LENGTH BTWN SIDES OF JELLYFISH BELL
starting_muscle = 318+1;
ending_muscle = 318+39;
springs_info(starting_muscle:ending_muscle,4) = abs( cos(freq*current_time*pi) );

%NOTE: not 2*pi*ft b/c of abs()

author: nickabattista


text:
 function [smap, sinit] = ir_mri_sensemap_admm_3D_linear(ykj, varargin)
%function [smap, sinit] = ir_mri_sensemap_admm_3D_linear(ykj, [options])
%|
%| Regularized estimation of sensitivity maps for parallel MRI
%| using augmented Lagrangian methods with variable splitting.
%| See M.J. Allison et al., IEEE TMI, 32(3), 556-564, Mar. 2013.
%| (This code is less optimized/more general than the one used for the paper.)
%|
%| Signal model for "surface coil" image: y_j = s_j f_j + noise_j
%|	s_j	sensitivity map (relative to "bodycoil" reference image)
%|	f_j	unknown underlying object
%|	j: voxel index
%| Signal model for "body coil" image: z_j = f_j + noise_j.
%|
%| Cost function:
%|	argmin 1/2 ||z - Ds||_W^2 + lam/2 ||R s||^2
%| where:
%|	D = diag{y_j}
%|	W = diag{w_j}
%|	R is finite differencing matrix with non-periodic boundary conditions.
%|
%| Augmented Lagrangian for ADMM with u1=s and u0=C*s:
%|	L(u,s) = 1/2 ||h - A u||^2 + 1/2 ||u - Ts - eta||_V^2
%|	so u = [u1; u0], T = [I; C], h = [sqrt(W)*z; 0],
%|	A = [sqrt(W)*D 0 ; 0 sqrt(lam)*B], eta = [eta1; eta0],
%|	and V = [v1*I 0; 0 v0*I], where R = B*C and B has boundary effects
%|
%| in
%|	ykj	[nx ny nz ncoil]	noisy complex images for each coil
%|
%| options
%|	bodycoil	[nx ny nz]	body coil (reference) image
%|				(def: sum-of-squares with phase of 1st coil)
%|	init			type of initialization
%|					 (see ir_mri_sensemap_init.m)
%|	l2b			log_2(beta), regularization parameter (def: -5)
%|	order			regularization order: 1|2 (default: 2)
%|	niter			# of iterations (default: 500)
%|	isave			which iterations to save (number
%|					list, 'last', 'all') (default: last)
%|	thresh			thresholding value for initialization (default: 0.05)
%|	quiet			hide iteration display if 1 (default:0)
%|	maskD			mask specifying pixels that include real data
%|	etabtw			update eta between each alternating
%|					minimization step (def: 1)
%| out
%|	smap	[nx*ny ncoil (isave)]	denoised sensitivity maps
%|	sinit	[nx*ny ncoil]		initialization images
%|
%| Copyright 2013, Michael Allison, University of Michigan
%| 2015-08 Lianli Liu  Extended to 3D using LADMM method

% todo; see fessler/l/tex/paper/submit/liu,bias/
% argsADMM = { 'bodycoil', bodycoil_s, 'maskD', maskObj_s,        'conds', ksADMM, 'condu0', ku0ADMM,'l2b', l2b,'niter',niter,...

if nargin == 1 && streq(ykj, 'test'), ir_mri_sensemap_admm_test, return, end
if nargin < 2, ir_usage, end

% defaults
arg.bodycoil = [];
arg.l2b = -5;
arg.order = 2;
arg.niter = 500;
arg.isave = [];
arg.quiet = 0;
arg.thresh = 0.12;
arg.init = [];
arg.maskD = true(size(ykj(:,:,1)));
arg.v0 = [];
arg.v1 = [];
arg.conds = [];
arg.condu0 = [];
arg.etabtw = 1;

arg = vararg_pair(arg, varargin);

%% Check for common cases of misuse.
if ~islogical(arg.maskD)
	error('arg.maskD must be logical datatype')
end

if ~isempty(arg.conds)
	if ~isempty(arg.v1)
		fail(['You are specifying the condition number for s, ' ...
			'so you cannot specify v1 directly.'])
	end
end
if ~isempty(arg.condu0)
	if ~isempty(arg.v0)
		fail(['You are specifying the condition number for u0, ' ...
			'so you cannot specify v0 directly.'])
	end
end

% create parameters
arg.isave = iter_saver(arg.isave, arg.niter);
[nx ny nz ncoil] = size(ykj);
N = nx * ny * nz;
M = sum(arg.maskD);

% default reference image
if isempty(arg.bodycoil)
	arg.bodycoil = sqrt(sum(abs(ykj).^2, 3)); % sum-of-squares
	tmp = angle(ykj(:,:,1)); % crude estimate of phase of image f_j
	arg.bodycoil = arg.bodycoil .* exp(1i * angle(ykj(:,:,1)));
end
arg.bodycoil = arg.bodycoil(:); % lexicographical

%% create regularizer
%tic;
if nz > 1
	[C,B,phi] = ir_circ_zeroed_reg_3D(arg.order, true([nx ny nz]), 0);
else
	[C,B,phi] = ir_circ_zeroed_reg(arg.order, true([nx ny nz]), 0);
end
%printm('regulizer initialization:');
%toc

%% create output matix
smap = zeros(N,ncoil,length(arg.isave));
sinit = zeros(N,ncoil);

%% determine v0/v1 from condition numbers if needed (untimed)
if ~isempty(arg.condu0)
	arg.v0 = 2^arg.l2b / (arg.condu0 - 1);
end
if ~isempty(arg.conds)
	arg.v1 = (arg.v0 * (max(abs(phi)) - arg.conds * min(abs(phi)))) ...
		/ (arg.conds - 1);
end
condu0Act = (2^arg.l2b + arg.v0) / arg.v0;
condsAct = (max(abs(phi)) + arg.v1/arg.v0) / ...
		(min(abs(phi)) + arg.v1/arg.v0);
condu1Act = (1 + arg.v1) / arg.v1;
printm(['Actual v0:' num2str(arg.v0) ' v1:' num2str(arg.v1)]);
printm(['Actual K(u0):' num2str(condu0Act) ' K(s):' num2str(condsAct) ...
	' K(u1):' num2str(condu1Act)]);

%% ADMM Algorithm
for i = 1:ncoil % it is better to do this in parallel as they are independent.
	%yj = ykj(:,:,i);
	yj = ykj(:);

	%% initialize AL algorithm

	% initialize smap
	sk = ir_mri_sensemap_init_3D(arg.init, yj, arg.bodycoil, [nx ny nz], ...
			arg.thresh, arg.maskD);
	sk = sk(:);
	sinit(:,i) = sk;

	u1 = sk;
	u0 = C * u1;
	eta0 = zeros(size(u0));
	eta1 = zeros(N,1);

	% apply maskD
	yj = yj(arg.maskD);
	arg.bodycoilM = arg.bodycoil(arg.maskD);

	% init convergence params
	v1 = arg.v1;
	v0 = arg.v0;
	v1ov0 = arg.v1 / arg.v0;
	lamov0 = (2^arg.l2b) / arg.v0;

	% precompute inverse terms that dont change with coil or iteration
	if ~exist('GinvS','var')
		Ginv1 = zeros(N,1);
		Ginv1(arg.maskD) = abs(arg.bodycoilM).^2;
		Ginv1 = Ginv1 + arg.v1;
		Ginv1 = 1 ./ Ginv1;

		GinvS = phi(:) + v1ov0;
		GinvS = 1 ./ GinvS;
		GinvS = reshape(GinvS, [nx ny nz]);

		Ginv0 = lamov0 * B + 1; % B is 0 or 1 and DIAGONAL
		Ginv0 = 1 ./ Ginv0;
	end

	DMz = zeros(N,1);
	DMz(arg.maskD) = conj(arg.bodycoilM) .* yj;

	if any(arg.isave == 0)
		smap(:,i,arg.isave == 0) = sinit(:,i);
	end

	printm('pre computation:');
	%toc

	% begin iterations
	for j = 1:arg.niter

		% update sk
%		sk = ifftn(GinvS .* fftn(reshape(C' * (u0 - eta0) + v1ov0 * ...
%			(u1 - eta1), [nx ny nz])));
		rou = 1./49;
		%sk = sk-rou.*(sk+v1ov0*(C'*C)*sk-C' * (u0 - eta0) - v1ov0 * (u1 - eta1));
		temp = v1ov0*sk+(C'*C)*sk-v1ov0*(u1-eta1)-C'*(u0-eta0);
		sk = sk-rou*temp;

		sk = sk(:);

		% this is used often
		Cs = C * sk;

		if arg.etabtw
			eta0 = eta0 - (u0 - Cs);
			eta1 = eta1 - (u1 - sk);
		end

		% update u1 and u0
		u1 = Ginv1 .* (DMz + arg.v1 * (sk + eta1));
		u0 = Ginv0 .* (Cs + eta0);

		% update etas
		eta0 = eta0 - (u0 - Cs);
		eta1 = eta1 - (u1 - sk);

		% save subset of iterations
		if any(arg.isave == j)
			smap(:,i,arg.isave == j) = sk; % eta1;
		end

		% display progress counter
		if (mod(j,500) == 1) && (~arg.quiet)
		%	printm('outputting eta1');
			printm([num2str(j) ' of ' num2str(arg.niter)])
		end
	end % for j
end % for i

end % ir_mri_sensemap_admm()

% create_mask()
function mask = create_mask(img, thresh, dial)
img = abs(img);
img = img ./ max(max(img));

mask = (img >= thresh);

if dial > 0
	SE = strel('disk', dial, 0);

	mask = imdilate(mask, SE);
	mask = logical(mask);
end

end % create_mask()

function ir_mri_sensemap_admm_test()
%| A test function for the toolbox version of ir_mri_sensemap_admm.
%| Note: This function is not optimized like the one in the paper
%| (uses Fatrix etc.) and thus the results of the test will not match those published.
%|
%| Copywrite, M Allison, U of Michigan, 2013.

%% set parameters
ns = [0.0275,0.005]; % noise level for [bcoil, surfcoil] images
threshD = 0.13;
sCoil = 1;
order = 2;
l2b = 5;
niter = 1000;
ksADMM = 650;
ku0ADMM = 255;

%% create senstivity data
% create data (based on M Allison file generateSENSEdata.m)
printm('Creating sensitivity data.')
f.dir = [path_find_dir('mri') '/../data/mri/'];
f.xtrue = [f.dir 'brainweb_t1.jpg'];
ftrue = double6(imread(f.xtrue)');
ftrue = ftrue(2:end-1,2:end-1); % make it 256^2

% add phase component to real image
[sX, sY] = size(ftrue);
tmp1 = [1:sX] ./ sX;
tmp2 = [1:sY] ./ sY;
[X,Y] = meshgrid(tmp1,tmp2);
ctrue = (Y .* X .* pi) - pi/2; ctrue = ctrue';
ftrue = ftrue.*exp(j*ctrue);
clear tmp1; clear tmp2; clear X; clear Y; clear ctrue

%normalize image
ftrue = ftrue ./ max(max(abs(ftrue)));

%create sense maps (uses mri_sensemap_sim.m from toolbox)
ncoil = 4;
smap = mri_sensemap_sim('nx', sX, 'ny', sY, 'dx', 192/sX, ...
		'ncoil', ncoil, 'orbit', 90*ncoil, 'rcoil', 100);
if ncoil == 4, smap = smap(:,:,[4 1 3 2]); end
smap = smap ./ max(max(max(abs(smap))));

bcoilI = ftrue + ns(1) * (randn(size(ftrue)) + 1i * randn(size(ftrue)));
lcoilI = smap .* repmat(ftrue, [1 1 ncoil]) + ns(2) * ...
	(randn(size(smap)) + 1i * randn(size(smap)));

% flip so it is in correct orientation for radiologists
ftrue = fliplr(ftrue);
bcoilI = fliplr(bcoilI);
lcoilI = flipdim(lcoilI,2);
smap = flipdim(smap,2);

%% create mask of real data
maskObj = create_mask(bcoilI, threshD, 0);

%% normalize coils and apply W mask outside of code (for Cholesky)
bcoilZ = bcoilI;
mB = max(max(abs(bcoilI)));
bcoilI = bcoilI ./ mB;
lcoilI = lcoilI ./ mB;
bcoilZ(~maskObj) = 0;

%% select a coil to work with
lcoilI = lcoilI(:,:,sCoil);
smap = smap(:,:,sCoil);
printm('-- Done')

%% Compute the Cholesky based estimate for this small problem
if 1
	printm('Starting Cholesky based method')
	printm('-- Computing C');

	[jnk,jnk,jnk,C] = ir_circ_zeroed_reg(order,true(size(bcoilI)),1); % same as above three lines

	% implement maskR
	bI = col(bcoilZ(:));
	sI = col(lcoilI(:));

	% compute using sparse Cholesky
	printm('-- Computing H');
	bI2 = conj(bI) .* bI;
	H = spdiag(bI2) + 2^(l2b) * C' * C;
	printm('-- Computing rhs');
	rhs = (conj(bI) .* sI);
	printm('-- Computing lhs');
	tic;
	lhs = H \ double(rhs);
	toc
	est_chol = lhs;
	est_chol = embed(est_chol,true(size(bcoilI)));
	printm('-- Done')
end

%% Run the ADMM Estimator
printm('Starting ADMM based method')
argsADMM = {'order', order, 'bodycoil', bcoilI, 'l2b', ...
		l2b, 'niter', niter, 'isave', 'last', 'thresh', threshD, ...
		'conds', ksADMM, 'condu0', ku0ADMM, 'maskD', maskObj, 'etabtw', ...
	1, 'init', 'avg'};
tADMM= tic;
[est_admm] = ir_mri_sensemap_admm(lcoilI, argsADMM{:});
toc(tADMM)
est_admm = reshape(est_admm,size(bcoilI));
printm('-- Done')

%% compare results to truth and Cholesky estimator
im pl 2 3
im(1, smap)
im(2, est_chol)
im(3, est_admm)
im(5, est_chol-smap)
im(6, est_admm-smap)

end % ir_mri_sensemap_admm_test()

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%% 2.5d grid-based datmo
% Inputs: Velodyne points and GPS/IMU localization
% Output: Detection and tracking moving objects (DATMO)
% Alireza Asvadi, 2015 Apr
%% clear memory & command window
clc
clear all
close all
%% setting
[st, stack]   = stt;                                          % set parameters
x.trc         = []; ncn = []; k = []; nzn = [];               % keep object locations for plot
%% datmo
for frame     = st.st.st : st.st.tn

%% observations
[mat, pts]    = obs(st, frame);                               % points (with trns) and grid
%% robust motion detection / moving object detection
[brm, stack]  = mdl(stack, pts, mat, st);                     % background model
[frm, pts, ocn, szn] = frg(brm, pts, mat, st);                % foreground model & centroid of objects in the world coordinate
%% multiple object tracking
k             = klmi(ocn, ncn, nzn, k, frame, st);            % initialize new kalman filter
[k, ncn, nzn] = asc(k, ocn, szn);                             % data association + new candidates
k             = klmf(k);                                      % kalman tracking
x             = plt(x, k, pts, frm, st, frame);               % record object locations and plot

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function p08_title ( )

%*****************************************************************************80
%
%% P08_TITLE prints a title for problem 08.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    None
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Problem 08\n' );
  fprintf ( 1, '  Name:       RC01\n' );
  fprintf ( 1, '              Crandall, page 49, #1\n' );
  fprintf ( 1, '  Region:     0 <= X(i) <= 1\n' );
  fprintf ( 1, '  Integrand:  F(X) = sin^2 ( pi/4 * sum ( X(i) ) )\n' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function evalRes = wrapperScene(typ, trSet, valSet, testSet, classMapping, sceneMapping, objectDir)
	paths = getPaths();

	modelFileName = trainSceneModel({trSet, valSet}, paths, typ, classMapping, sceneMapping, objectDir);
	outputFileName = testSceneModel(testSet, paths, modelFileName);
	evalRes = benchmarkScene(outputFileName, sceneMapping, testSet);
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
function pf=lpcff2pf(ff)
%LPCFF2PF Convert complex spectrum to power spectrum PF=(FF)

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: lpcff2pf.m,v 1.4 2007/05/04 07:01:38 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pf=abs(ff).^2;

author: decouples


text:
function robTwoDynamicsDisplay(model, spaceNum)

% ROBTWODYNAMICSDISPLAY Display the robot dynamics model.

% FGPLVM

if nargin > 1
  spacing = repmat(32, 1, spaceNum);
else
  spaceNum = 0;
  spacing = [];
end
spacing = char(spacing);
fprintf(spacing);
fprintf('Second Tailored dynamics model for robot:\n')
fprintf(spacing);
fprintf('  Length distribution value of a: %2.4f\n', model.a)
fprintf(spacing);
fprintf('  Length distribution value of b: %2.4f\n', model.b)
fprintf(spacing);
fprintf('  Length distribution mix coefficient: %2.4f\n', model.mixR);
fprintf(spacing);
fprintf('  Angle distribution component 1 sigma2: %2.4f\n', model.sigma2One);
fprintf(spacing);
fprintf('  Angle distribution component 1 mix coefficient: %2.4f\n', model.mixThetaOne);
fprintf(spacing);
fprintf('  Angle distribution component 2 sigma2: %2.4f\n', model.sigma2Two);
fprintf(spacing);
fprintf('  Angle distribution component 2 mix coefficient: %2.4f\n', model.mixThetaTwo);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ n_data, n, k, x, b ] = bernstein_poly_01_values ( n_data )

%*****************************************************************************80
%
%% BERNSTEIN_POLY_01_VALUES returns some values of the Bernstein polynomials.
%
%  Discussion:
%
%    The Bernstein polynomials are assumed to be based on [0,1].
%
%    The formula for the Bernstein polynomials is
%
%      B(N,I)(X) = [N!/(I!*(N-I)!)] * (1-X)^(N-I) * X^I
%
%    In Mathematica, the function can be evaluated by:
%
%      Binomial[n,i] * (1-x)^(n-i) * x^i
%
%  First values:
%
%    B(0,0)(X) = 1
%
%    B(1,0)(X) =      1-X
%    B(1,1)(X) =               X
%
%    B(2,0)(X) =     (1-X)^2
%    B(2,1)(X) = 2 * (1-X)   * X
%    B(2,2)(X) =               X^2
%
%    B(3,0)(X) =     (1-X)^3
%    B(3,1)(X) = 3 * (1-X)^2 * X
%    B(3,2)(X) = 3 * (1-X)   * X^2
%    B(3,3)(X) =               X^3
%
%    B(4,0)(X) =     (1-X)^4
%    B(4,1)(X) = 4 * (1-X)^3 * X
%    B(4,2)(X) = 6 * (1-X)^2 * X^2
%    B(4,3)(X) = 4 * (1-X)   * X^3
%    B(4,4)(X) =               X^4
%
%  Special values:
%
%    B(N,I)(X) has a unique maximum value at X = I/N.
%
%    B(N,I)(X) has an I-fold zero at 0 and and N-I fold zero at 1.
%
%    B(N,I)(1/2) = C(N,K) / 2^N
%
%    For a fixed X and N, the polynomials add up to 1:
%
%      Sum ( 0 <= I <= N ) B(N,I)(X) = 1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer N, the degree of the polynomial.
%
%    Output, integer K, the index of the polynomial.
%
%    Output, real X, the argument of the polynomial.
%
%    Output, real B, the value of the polynomial B(N,K)(X).
%
  n_max = 15;

  b_vec = [ ...
     0.1000000000000000E+01, ...
     0.7500000000000000E+00, ...
     0.2500000000000000E+00, ...
     0.5625000000000000E+00, ...
     0.3750000000000000E+00, ...
     0.6250000000000000E-01, ...
     0.4218750000000000E+00, ...
     0.4218750000000000E+00, ...
     0.1406250000000000E+00, ...
     0.1562500000000000E-01, ...
     0.3164062500000000E+00, ...
     0.4218750000000000E+00, ...
     0.2109375000000000E+00, ...
     0.4687500000000000E-01, ...
     0.3906250000000000E-02 ];

  k_vec = [ ...
    0, ...
    0, 1, ...
    0, 1, 2, ...
    0, 1, 2, 3, ...
    0, 1, 2, 3, 4 ];

  n_vec = [ ...
    0, ...
    1, 1, ...
    2, 2, 2, ...
    3, 3, 3, 3, ...
    4, 4, 4, 4, 4 ];

  x_vec = [ ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    k = 0;
    x = 0.0;
    b = 0.0;
  else
    n = n_vec(n_data);
    k = k_vec(n_data);
    x = x_vec(n_data);
    b = b_vec(n_data);
  end

  return
end


text:
%%%%%%%%%%%%%%%%%%%% SHOW EXTRINSIC RESULTS %%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('show_camera'),
    show_camera = 1;
end;

if ~exist('n_ima')|~exist('fc'),
   fprintf(1,'No calibration data available.\n');
   return;
end;

check_active_images;

if n_ima ~= 0,
if ~exist(['omc_' num2str(ind_active(1))]),
   fprintf(1,'No calibration data available.\n');
   return;
end;
end;

%if ~exist('no_grid'),
   no_grid = 0;
%end;

if n_ima ~= 0,
if ~exist(['n_sq_x_' num2str(ind_active(1))]),
   no_grid = 1;
end;
else
    no_grid = 1;
end;

if ~exist('alpha_c'),
   alpha_c = 0;
end;

if 0,

err_std = std(ex');

fprintf(1,'\n\nCalibration results without principal point estimation:\n\n');
fprintf(1,'Focal Length:     fc = [ %3.5f   %3.5f]\n',fc);
fprintf(1,'Principal point:  cc = [ %3.5f   %3.5f]\n',cc);
fprintf(1,'Distortion:       kc = [ %3.5f   %3.5f   %3.5f   %3.5f]\n',kc);
fprintf(1,'Pixel error:      err = [ %3.5f   %3.5f]\n\n',err_std);

end;

% Color code for each image:

colors = 'brgkcm';

%%% Show the extrinsic parameters

if n_ima ~= 0,
if ~exist('dX'),
   eval(['dX = norm(Tc_' num2str(ind_active(1)) ')/10;']);
   dY = dX;
end;
else
    dX = 1;
end;

IP = 5*dX*[1 -alpha_c 0;0 1 0;0 0 1]*[1/fc(1) 0 0;0 1/fc(2) 0;0 0 1]*[1 0 -cc(1);0 1 -cc(2);0 0 1]*[0 nx-1 nx-1 0 0 ; 0 0 ny-1 ny-1 0;1 1 1 1 1];
BASE = 5*dX*([0 1 0 0 0 0;0 0 0 1 0 0;0 0 0 0 0 1]);
IP = reshape([IP;BASE(:,1)*ones(1,5);IP],3,15);

if ishandle(4),
	figure(4);
   [a,b] = view;
else
   figure(4);
   a = 50;
   b = 20;
end;

if show_camera,
    figure(4);
    plot3(BASE(1,:),BASE(3,:),-BASE(2,:),'b-','linewidth',2);
    hold on;
    plot3(IP(1,:),IP(3,:),-IP(2,:),'r-','linewidth',2);
    text(6*dX,0,0,'X_c');
    text(-dX,5*dX,0,'Z_c');
    text(0,0,-6*dX,'Y_c');
    text(-dX,-dX,dX,'O_c');
else
    figure(4);
    clf;
    hold on;
end;

for kk = 1:n_ima,

   if active_images(kk);

      if exist(['X_' num2str(kk)]) & exist(['omc_' num2str(kk)]),

	 eval(['XX_kk = X_' num2str(kk) ';']);

	 if ~isnan(XX_kk(1,1))

	    eval(['omc_kk = omc_' num2str(kk) ';']);
	    eval(['Tc_kk = Tc_' num2str(kk) ';']);
	    N_kk = size(XX_kk,2);

	    if ~exist(['n_sq_x_' num2str(kk)]),
	       no_grid = 1;
	    else
	       eval(['n_sq_x = n_sq_x_' num2str(kk) ';']);
	       if isnan(n_sq_x(1)),
		  no_grid = 1;
	       end;
	    end;

	    if ~no_grid,
	       eval(['n_sq_x = n_sq_x_' num2str(kk) ';']);
	       eval(['n_sq_y = n_sq_y_' num2str(kk) ';']);
	       if (N_kk ~= ((n_sq_x+1)*(n_sq_y+1))),
		  no_grid = 1;
	       end;
	    end;

	    if ~isnan(omc_kk(1,1)),

	       R_kk = rodrigues(omc_kk);

	       YY_kk = R_kk * XX_kk + Tc_kk * ones(1,length(XX_kk));

	       uu = [-dX;-dY;0]/2;
	       uu = R_kk * uu + Tc_kk;

	       if ~no_grid,
		  YYx = zeros(n_sq_x+1,n_sq_y+1);
		  YYy = zeros(n_sq_x+1,n_sq_y+1);
		  YYz = zeros(n_sq_x+1,n_sq_y+1);

		  YYx(:) = YY_kk(1,:);
		  YYy(:) = YY_kk(2,:);
		  YYz(:) = YY_kk(3,:);

		  %keyboard;

		  figure(4);
		  hhh= mesh(YYx,YYz,-YYy);
		  set(hhh,'edgecolor',colors(rem(kk-1,6)+1),'linewidth',1); %,'facecolor','none');
		  %plot3(YY_kk(1,:),YY_kk(3,:),-YY_kk(2,:),['o' colors(rem(kk-1,6)+1)]);
		  text(uu(1),uu(3),-uu(2),num2str(kk),'fontsize',14,'color',colors(rem(kk-1,6)+1));
	       else

		  figure(4);
		  plot3(YY_kk(1,:),YY_kk(3,:),-YY_kk(2,:),['.' colors(rem(kk-1,6)+1)]);
		  text(uu(1),uu(3),-uu(2),num2str(kk),'fontsize',14,'color',colors(rem(kk-1,6)+1));

	       end;

	    end;

	 end;

      end;

   end;

end;

figure(4);rotate3d on;
axis('equal');
title('Extrinsic parameters (camera-centered)');
%view(60,30);
view(a,b);
grid on;
hold off;
axis vis3d;
axis tight;
set(4,'color',[1 1 1]);
if ~show_camera,
    xlabel('X_c');
    ylabel('Z_c');
    zlabel('<-- Y_c');
end;

set(4,'Name','3D','NumberTitle','off');

%fprintf(1,'To generate the complete movie associated to the optimization loop, try: check_convergence;\n');

if exist('h_switch2')==1,
    if ishandle(h_switch2),
        delete(h_switch2);
    end;
end;

if n_ima ~= 0,
    if show_camera,
        h_switch2 = uicontrol('Parent',4,'Units','normalized', 'Callback','show_camera=0;ext_calib;', 'Position',[1-.30 0.04  .30  .04],'String','Remove camera reference frame','fontsize',8,'fontname','clean','Tag','Pushbutton1');
    else
        h_switch2 = uicontrol('Parent',4,'Units','normalized', 'Callback','show_camera=1;ext_calib;', 'Position',[1-.30 0.04  .30  .04],'String','Add camera reference frame','fontsize',8,'fontname','clean','Tag','Pushbutton1');
    end;
end;

if exist('h_switch')==1,
    if ishandle(h_switch),
        delete(h_switch);
    end;
end;

if n_ima ~= 0,
h_switch = uicontrol('Parent',4,'Units','normalized', 'Callback','ext_calib2', 'Position',[1-.30 0  .30  .04],'String','Switch to world-centered view','fontsize',8,'fontname','clean','Tag','Pushbutton1');
end;

figure(4);
rotate3d on;

repo: ekfmonoslam


text:
clear all
clc

%% Prime Factors
% This code gets user input number, calculates and displays its prime factors.
% For this, first it determines prime numbers which are less than or equal to
% user input number. Then if the input is dividable by that prime number,
% it becomes one of input's prime factors.

%% Request user input
prompt = 'Input your number: ';
n = input(prompt);

%%
counter = 0; % initialize number of prime factors

if n <= 1
    disp('input must be positive integer greater than 1')
else if floor(n)~= n
       disp('input must be positive integer')
    else
        for i = 2:1:n
            if i == 2
               isprime = 1;
            else
            half_i = floor(i/2)+1;
            j = 2;
            while j <= half_i         %lines 16 to 30 check if i is prime or not.
                  residual = mod(i,j);
                 if residual == 0
                    isprime = 0;
                    break
                 else if j == half_i
                         isprime = 1;
                         break
                     else
                         j=j+1;
                     end
                 end
            end
            end
            if isprime == 1 && mod(n,i) == 0
                   counter=counter+1;
                   f(counter) = i;         % prime factors of n will be storing
            end
        end
    end
end

disp('Prime factors of input number are: ')
disp(f)

author: TheAlgorithms

repo: MATLAB-Octave