text:
function t_scale_load(quiet)
%T_SCALE_LOAD  Tests for code in SCALE_LOAD.

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

n_tests = 846;

t_begin(n_tests, quiet);

%% 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;

mpc = loadcase('t_auction_case');
mpc.gen(8, GEN_BUS) = 2;    %% multiple d. loads per area, same bus as gen
mpc.gen(8, [QG QMIN QMAX]) = [ 3 0 3 ];
mpc.gencost(7, COST:end) = [-30 -600 -20 -300 -10 -100 0 0];    % 10, 20, 30
mpc.gencost(8, COST:end) = [-30  -60 -20  -30 -10  -10 0 0];    % 1, 2, 3
mpc.gencost(9, COST:end) = [-30 -850 -10 -250  -5  -50 0 0];    % 10, 20, 30
%% split a load in 2 half-size duplicates at same bus ...
mpc.gen(8, [PG PMIN PMAX QG QMIN QMAX]) = mpc.gen(8, [PG PMIN PMAX QG QMIN QMAX]) / 2;
mpc.gencost(8, COST:end) = mpc.gencost(8, COST:end) / 2;
%% ... with one before gens in matrix
mpc.gen = [mpc.gen(8, :); mpc.gen];
mpc.gencost = [mpc.gencost(8, :); mpc.gencost];
gg = find(~isload(mpc.gen));
ld = find( isload(mpc.gen));
for k = 1:3
    a{k} = find(mpc.bus(:, BUS_AREA) == k); %% buses in area k
    tmp = find(ismember(mpc.gen(ld, GEN_BUS), a{k}));
    lda{k} = ld(tmp);                       %% disp loads in area k
end
for k = 1:3
    area(k).fixed.p = sum(mpc.bus(a{k}, PD));
    area(k).fixed.q = sum(mpc.bus(a{k}, QD));
    area(k).disp.p = -sum(mpc.gen(lda{k}, PMIN));
    area(k).disp.qmin = -sum(mpc.gen(lda{k}, QMIN));
    area(k).disp.qmax = -sum(mpc.gen(lda{k}, QMAX));
    area(k).disp.q = area(k).disp.qmin + area(k).disp.qmax;
    area(k).both.p = area(k).fixed.p + area(k).disp.p;
    area(k).both.q = area(k).fixed.q + area(k).disp.q;
end
total.fixed.p = sum(mpc.bus(:, PD));
total.fixed.q = sum(mpc.bus(:, QD));
total.disp.p = -sum(mpc.gen(ld, PMIN));
total.disp.qmin = -sum(mpc.gen(ld, QMIN));
total.disp.qmax = -sum(mpc.gen(ld, QMAX));
total.disp.q = total.disp.qmin + total.disp.qmax;
total.both.p = total.fixed.p + total.disp.p;
total.both.q = total.fixed.q + total.disp.q;
orig_gc = mpc.gencost;

%%-----  single load zone, one scale factor  -----
dmd = 2;
t = 'all fixed loads (PQ) * 2 : ';
bus = scale_load(dmd, mpc.bus);
t_is(sum(bus(:, PD)), dmd*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), dmd*total.fixed.q, 8, [t 'total fixed Q']);
opt = struct('which', 'FIXED');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
t_is(sum(bus(:, PD)), dmd*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), dmd*total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);

t = 'all fixed loads (P) * 2 : ';
opt = struct('pq', 'P');
bus = scale_load(dmd, mpc.bus, [], [], opt);
t_is(sum(bus(:, PD)), dmd*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
opt = struct('pq', 'P', 'which', 'FIXED');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
t_is(sum(bus(:, PD)), dmd*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);

t = 'all loads (PQ) * 2 : ';
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen);
t_is(sum(bus(:, PD)), dmd*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), dmd*total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), dmd*total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), dmd*total.disp.qmax, 8, [t 'total disp Qmax']);

t = 'all loads/costs (PQ) * 2 : ';
[bus, gen, gencost] = scale_load(dmd, mpc.bus, mpc.gen, [], [], mpc.gencost);
t_is(sum(bus(:, PD)), dmd*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), dmd*total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), dmd*total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), dmd*total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end), 2*orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all loads/costs (PQ) * 2 : ';
opt = struct('cost', 1);
[bus, gen, gencost] = scale_load(dmd, mpc.bus, mpc.gen, [], opt, mpc.gencost);
t_is(sum(bus(:, PD)), dmd*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), dmd*total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), dmd*total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), dmd*total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end), 2*orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all loads (PQ) * 2 : ';
opt = struct('cost', 0);
[bus, gen, gencost] = scale_load(dmd, mpc.bus, mpc.gen, [], opt, mpc.gencost);
t_is(sum(bus(:, PD)), dmd*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), dmd*total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), dmd*total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), dmd*total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end),   orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all loads (P) * 2 : ';
opt = struct('pq', 'P');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
t_is(sum(bus(:, PD)), dmd*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);

t = 'all loads/costs (P) * 2 : ';
opt = struct('pq', 'P');
[bus, gen, gencost] = scale_load(dmd, mpc.bus, mpc.gen, [], opt, mpc.gencost);
t_is(sum(bus(:, PD)), dmd*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end), 2*orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all disp loads (PQ) * 2 : ';
opt = struct('which', 'DISPATCHABLE');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
t_is(sum(bus(:, PD)), total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), dmd*total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), dmd*total.disp.qmax, 8, [t 'total disp Qmax']);

t = 'all disp loads/costs (PQ) * 2 : ';
opt = struct('which', 'DISPATCHABLE');
[bus, gen, gencost] = scale_load(dmd, mpc.bus, mpc.gen, [], opt, mpc.gencost);
t_is(sum(bus(:, PD)), total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), dmd*total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), dmd*total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end), 2*orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all disp loads (P) * 2 : ';
opt = struct('pq', 'P', 'which', 'DISPATCHABLE');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
t_is(sum(bus(:, PD)), total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);

t = 'all disp loads/costs (P) * 2 : ';
opt = struct('pq', 'P', 'which', 'DISPATCHABLE');
[bus, gen, gencost] = scale_load(dmd, mpc.bus, mpc.gen, [], opt, mpc.gencost);
t_is(sum(bus(:, PD)), total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end), 2*orig_gc(ld, COST:end), 8, [t 'gencost loads']);

%%-----  single load zone, one scale quantity  -----
dmd = 200;
t = 'all fixed loads (PQ) => total = 200 : ';
opt = struct('scale', 'QUANTITY');
bus = scale_load(dmd, mpc.bus, [], [], opt);
t_is(sum(bus(:, PD)), dmd, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), dmd/total.fixed.p*total.fixed.q, 8, [t 'total fixed Q']);
opt = struct('scale', 'QUANTITY', 'which', 'FIXED');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
t_is(sum(bus(:, PD)), dmd-total.disp.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), (dmd-total.disp.p)/total.fixed.p*total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);

t = 'all fixed loads (P) => total = 200 : ';
opt = struct('scale', 'QUANTITY', 'pq', 'P');
bus = scale_load(dmd, mpc.bus, [], [], opt);
t_is(sum(bus(:, PD)), dmd, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
opt = struct('scale', 'QUANTITY', 'pq', 'P', 'which', 'FIXED');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
t_is(sum(bus(:, PD)), dmd-total.disp.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);

t = 'all loads (PQ) => total = 200 : ';
opt = struct('scale', 'QUANTITY');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
t_is(sum(bus(:, PD)), dmd/total.both.p*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), dmd/total.both.p*total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd/total.both.p*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), dmd/total.both.p*total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), dmd/total.both.p*total.disp.qmax, 8, [t 'total disp Qmax']);

t = 'all loads/costs (PQ) => total = 200 : ';
opt = struct('scale', 'QUANTITY');
[bus, gen, gencost] = scale_load(dmd, mpc.bus, mpc.gen, [], opt, mpc.gencost);
t_is(sum(bus(:, PD)), dmd/total.both.p*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), dmd/total.both.p*total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd/total.both.p*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), dmd/total.both.p*total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), dmd/total.both.p*total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end), dmd/total.both.p*orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all loads (P) => total = 200 : ';
opt = struct('scale', 'QUANTITY', 'pq', 'P');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
t_is(sum(bus(:, PD)), dmd/total.both.p*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd/total.both.p*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);

t = 'all loads/costs (P) => total = 200 : ';
opt = struct('scale', 'QUANTITY', 'pq', 'P');
[bus, gen, gencost] = scale_load(dmd, mpc.bus, mpc.gen, [], opt, mpc.gencost);
t_is(sum(bus(:, PD)), dmd/total.both.p*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd/total.both.p*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end), dmd/total.both.p*orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all disp loads (PQ) => total = 200 : ';
opt = struct('scale', 'QUANTITY', 'which', 'DISPATCHABLE');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
t_is(sum(bus(:, PD)), total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd-total.fixed.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), (dmd-total.fixed.p)/total.disp.p*total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), (dmd-total.fixed.p)/total.disp.p*total.disp.qmax, 8, [t 'total disp Qmax']);

t = 'all disp loads/costs (PQ) => total = 200 : ';
opt = struct('scale', 'QUANTITY', 'which', 'DISPATCHABLE');
[bus, gen, gencost] = scale_load(dmd, mpc.bus, mpc.gen, [], opt, mpc.gencost);
t_is(sum(bus(:, PD)), total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd-total.fixed.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), (dmd-total.fixed.p)/total.disp.p*total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), (dmd-total.fixed.p)/total.disp.p*total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end), (dmd-total.fixed.p)/total.disp.p*orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all disp loads (P) => total = 200 : ';
opt = struct('scale', 'QUANTITY', 'pq', 'P', 'which', 'DISPATCHABLE');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
t_is(sum(bus(:, PD)), total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd-total.fixed.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);

t = 'all disp loads/costs (P) => total = 200 : ';
opt = struct('scale', 'QUANTITY', 'pq', 'P', 'which', 'DISPATCHABLE');
[bus, gen, gencost] = scale_load(dmd, mpc.bus, mpc.gen, [], opt, mpc.gencost);
t_is(sum(bus(:, PD)), total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd-total.fixed.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end), (dmd-total.fixed.p)/total.disp.p*orig_gc(ld, COST:end), 8, [t 'gencost loads']);

%%-----  3 zones, area scale factors  -----
t = 'area fixed loads (PQ) * [3 2 1] : ';
dmd = [3 2 1];
bus = scale_load(dmd, mpc.bus);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)*area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), dmd(k)*area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
end
opt = struct('which', 'FIXED');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)*area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), dmd(k)*area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
end

t = 'area fixed loads (P) * [3 2 1] : ';
dmd = [3 2 1];
opt = struct('pq', 'P');
bus = scale_load(dmd, mpc.bus, [], [], opt);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)*area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
end
opt = struct('pq', 'P', 'which', 'FIXED');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)*area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
end

t = 'all area loads (PQ) * [3 2 1] : ';
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)*area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), dmd(k)*area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)*area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), dmd(k)*area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), dmd(k)*area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
end

t = 'all area loads/costs (PQ) * [3 2 1] : ';
[bus, gen, gencost] = scale_load(dmd, mpc.bus, mpc.gen, [], [], mpc.gencost);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)*area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), dmd(k)*area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)*area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), dmd(k)*area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), dmd(k)*area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    t_is(gencost(lda{k}, COST:end), dmd(k)*orig_gc(lda{k}, COST:end), 8, sprintf('%s area %d gencost loads', t, k));
end
t_is(gencost(gg, COST:end), orig_gc(gg, COST:end), 8, sprintf('%s gencost gens', t));

t = 'all area loads (P) * [3 2 1] : ';
opt = struct('pq', 'P');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)*area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)*area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
end

t = 'all area loads/costs (P) * [3 2 1] : ';
opt = struct('pq', 'P');
[bus, gen, gencost] = scale_load(dmd, mpc.bus, mpc.gen, [], opt, mpc.gencost);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)*area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)*area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    t_is(gencost(lda{k}, COST:end), dmd(k)*orig_gc(lda{k}, COST:end), 8, sprintf('%s area %d gencost loads', t, k));
end
t_is(gencost(gg, COST:end), orig_gc(gg, COST:end), 8, sprintf('%s gencost gens', t));

t = 'area disp loads (PQ) * [3 2 1] : ';
opt = struct('which', 'DISPATCHABLE');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)*area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), dmd(k)*area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), dmd(k)*area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
end

t = 'area disp loads/costs (PQ) * [3 2 1] : ';
opt = struct('which', 'DISPATCHABLE');
[bus, gen, gencost] = scale_load(dmd, mpc.bus, mpc.gen, [], opt, mpc.gencost);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)*area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), dmd(k)*area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), dmd(k)*area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    t_is(gencost(lda{k}, COST:end), dmd(k)*orig_gc(lda{k}, COST:end), 8, sprintf('%s area %d gencost loads', t, k));
end
t_is(gencost(gg, COST:end), orig_gc(gg, COST:end), 8, sprintf('%s gencost gens', t));

t = 'area disp loads (P) * [3 2 1] : ';
opt = struct('pq', 'P', 'which', 'DISPATCHABLE');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)*area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
end

t = 'area disp loads/costs (P) * [3 2 1] : ';
opt = struct('pq', 'P', 'which', 'DISPATCHABLE');
[bus, gen, gencost] = scale_load(dmd, mpc.bus, mpc.gen, [], opt, mpc.gencost);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)*area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    t_is(gencost(lda{k}, COST:end), dmd(k)*orig_gc(lda{k}, COST:end), 8, sprintf('%s area %d gencost loads', t, k));
end
t_is(gencost(gg, COST:end), orig_gc(gg, COST:end), 8, sprintf('%s gencost gens', t));

%%-----  3 zones, area scale quantities  -----
t = 'area fixed loads (PQ) => total = [100 80 60] : ';
dmd = [100 80 60];
opt = struct('scale', 'QUANTITY');
bus = scale_load(dmd, mpc.bus, [], [], opt);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k), 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), dmd(k)/area(k).fixed.p*area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
end
opt = struct('scale', 'QUANTITY', 'which', 'FIXED');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)-area(k).disp.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), (dmd(k)-area(k).disp.p)/area(k).fixed.p*area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
end

t = 'area fixed loads (P) => total = [100 80 60] : ';
dmd = [100 80 60];
opt = struct('scale', 'QUANTITY', 'pq', 'P');
bus = scale_load(dmd, mpc.bus, [], [], opt);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k), 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
end
opt = struct('scale', 'QUANTITY', 'pq', 'P', 'which', 'FIXED');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)-area(k).disp.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
end

t = 'all area loads (PQ) => total = [100 80 60] : ';
opt = struct('scale', 'QUANTITY');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)/area(k).both.p*area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), dmd(k)/area(k).both.p*area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)/area(k).both.p*area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), dmd(k)/area(k).both.p*area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), dmd(k)/area(k).both.p*area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
end

t = 'all area loads (P) => total = [100 80 60] : ';
opt = struct('scale', 'QUANTITY', 'pq', 'P');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)/area(k).both.p*area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)/area(k).both.p*area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
end

t = 'area disp loads (PQ) => total = [100 80 60] : throws expected exception';
dmd = [100 80 60];
opt = struct('scale', 'QUANTITY', 'which', 'DISPATCHABLE');
err = 0;
try
    [bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
catch
    [msg, id] = lasterr;
    expected = 'scale_load: impossible to make zone 2 load equal 80 by scaling non-existent dispatchable load';
    if ~isempty(strfind(msg, expected))
        err = 1;
    end
end
t_ok(err, t);

t = 'area disp loads (PQ) => total = [100 74.3941 60] : ';
dmd = [100 area(2).fixed.p 60];
opt = struct('scale', 'QUANTITY', 'which', 'DISPATCHABLE');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)-area(k).fixed.p, 8, sprintf('%s area %d disp P', t, k));
    if k == 2
        t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
        t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    else
        t_is(-sum(gen(lda{k}, QMIN)), (dmd(k)-area(k).fixed.p)/area(k).disp.p*area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
        t_is(-sum(gen(lda{k}, QMAX)), (dmd(k)-area(k).fixed.p)/area(k).disp.p*area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    end
end

t = 'area disp loads (P) => total = [100 74.3941 60] : ';
opt = struct('scale', 'QUANTITY', 'pq', 'P', 'which', 'DISPATCHABLE');
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, [], opt);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)-area(k).fixed.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
end

%%-----  explict single load zone  -----
t = 'explicit single load zone';
load_zone = zeros(1, size(mpc.bus, 1));
load_zone([3 4]) = 1;
dmd = 2;
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, load_zone);
Pd = mpc.bus(:, PD);
Pd([3 4]) = dmd * Pd([3 4]);
t_is( bus(:, PD), Pd, 8, t);

%%-----  explict multiple load zone  -----
t = 'explicit multiple load zone';
load_zone = zeros(1, size(mpc.bus, 1));
load_zone([3 4]) = 1;
load_zone([7 8]) = 2;
dmd = [2 0.5];
[bus, gen] = scale_load(dmd, mpc.bus, mpc.gen, load_zone);
Pd = mpc.bus(:, PD);
Pd([3 4]) = dmd(1) * Pd([3 4]);
Pd([7 8]) = dmd(2) * Pd([7 8]);
t_is( bus(:, PD), Pd, 8, t);

%%------------------------------------------------------------------
%% mostly same tests below, but with MPC input/output

%%-----  single load zone, one scale factor  -----
dmd = 2;
t = 'all loads (PQ) * 2 : ';
opt = struct('cost', 0);
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
t_is(sum(bus(:, PD)), dmd*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), dmd*total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), dmd*total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), dmd*total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end),   orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all loads/costs (PQ) * 2 : ';
mpc1 = scale_load(dmd, mpc);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
t_is(sum(bus(:, PD)), dmd*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), dmd*total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), dmd*total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), dmd*total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end), 2*orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all loads (P) * 2 : ';
opt = struct('pq', 'P', 'cost', 0);
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
t_is(sum(bus(:, PD)), dmd*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end),   orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all loads/costs (P) * 2 : ';
opt = struct('pq', 'P');
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
t_is(sum(bus(:, PD)), dmd*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end), 2*orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all disp loads (PQ) * 2 : ';
opt = struct('which', 'DISPATCHABLE', 'cost', 0);
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
t_is(sum(bus(:, PD)), total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), dmd*total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), dmd*total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end),   orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all disp loads/costs (PQ) * 2 : ';
opt = struct('which', 'DISPATCHABLE');
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
t_is(sum(bus(:, PD)), total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), dmd*total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), dmd*total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end), 2*orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all disp loads (P) * 2 : ';
opt = struct('pq', 'P', 'which', 'DISPATCHABLE', 'cost', 0);
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
t_is(sum(bus(:, PD)), total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end),   orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all disp loads/costs (P) * 2 : ';
opt = struct('pq', 'P', 'which', 'DISPATCHABLE', 'cost', 1);
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
t_is(sum(bus(:, PD)), total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end), 2*orig_gc(ld, COST:end), 8, [t 'gencost loads']);

%%-----  single load zone, one scale quantity  -----
dmd = 200;
t = 'all fixed loads (PQ) => total = 200 : ';
opt = struct('scale', 'QUANTITY', 'which', 'FIXED');
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
t_is(sum(bus(:, PD)), dmd-total.disp.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), (dmd-total.disp.p)/total.fixed.p*total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end),   orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all fixed loads (P) => total = 200 : ';
opt = struct('scale', 'QUANTITY', 'pq', 'P', 'which', 'FIXED', 'cost', 1);
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
t_is(sum(bus(:, PD)), dmd-total.disp.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end),   orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all loads (PQ) => total = 200 : ';
opt = struct('scale', 'QUANTITY', 'cost', 0);
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
t_is(sum(bus(:, PD)), dmd/total.both.p*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), dmd/total.both.p*total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd/total.both.p*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), dmd/total.both.p*total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), dmd/total.both.p*total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end),   orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all loads/costs (PQ) => total = 200 : ';
opt = struct('scale', 'QUANTITY');
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
t_is(sum(bus(:, PD)), dmd/total.both.p*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), dmd/total.both.p*total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd/total.both.p*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), dmd/total.both.p*total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), dmd/total.both.p*total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end), dmd/total.both.p*orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all loads (P) => total = 200 : ';
opt = struct('scale', 'QUANTITY', 'pq', 'P', 'cost', 0);
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
t_is(sum(bus(:, PD)), dmd/total.both.p*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd/total.both.p*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end),   orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all loads/costs (P) => total = 200 : ';
opt = struct('scale', 'QUANTITY', 'pq', 'P');
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
t_is(sum(bus(:, PD)), dmd/total.both.p*total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd/total.both.p*total.disp.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end), dmd/total.both.p*orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all disp loads (PQ) => total = 200 : ';
opt = struct('scale', 'QUANTITY', 'which', 'DISPATCHABLE', 'cost', 0);
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
t_is(sum(bus(:, PD)), total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd-total.fixed.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), (dmd-total.fixed.p)/total.disp.p*total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), (dmd-total.fixed.p)/total.disp.p*total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end),   orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all disp loads/costs (PQ) => total = 200 : ';
opt = struct('scale', 'QUANTITY', 'which', 'DISPATCHABLE');
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
t_is(sum(bus(:, PD)), total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd-total.fixed.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), (dmd-total.fixed.p)/total.disp.p*total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), (dmd-total.fixed.p)/total.disp.p*total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end), (dmd-total.fixed.p)/total.disp.p*orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all disp loads (P) => total = 200 : ';
opt = struct('scale', 'QUANTITY', 'pq', 'P', 'which', 'DISPATCHABLE', 'cost', 0);
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
t_is(sum(bus(:, PD)), total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd-total.fixed.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end),   orig_gc(ld, COST:end), 8, [t 'gencost loads']);

t = 'all disp loads/costs (P) => total = 200 : ';
opt = struct('scale', 'QUANTITY', 'pq', 'P', 'which', 'DISPATCHABLE');
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
t_is(sum(bus(:, PD)), total.fixed.p, 8, [t 'total fixed P']);
t_is(sum(bus(:, QD)), total.fixed.q, 8, [t 'total fixed Q']);
t_is(-sum(gen(ld, PMIN)), dmd-total.fixed.p, 8, [t 'total disp P']);
t_is(-sum(gen(ld, QMIN)), total.disp.qmin, 8, [t 'total disp Qmin']);
t_is(-sum(gen(ld, QMAX)), total.disp.qmax, 8, [t 'total disp Qmax']);
t_is(gencost(gg, COST:end),   orig_gc(gg, COST:end), 8, [t 'gencost gens']);
t_is(gencost(ld, COST:end), (dmd-total.fixed.p)/total.disp.p*orig_gc(ld, COST:end), 8, [t 'gencost loads']);

%%-----  3 zones, area scale factors  -----
t = 'area fixed loads (PQ) * [3 2 1] : ';
dmd = [3 2 1];
opt = struct('which', 'FIXED', 'cost', 1);
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)*area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), dmd(k)*area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    t_is(gencost(lda{k}, COST:end), orig_gc(lda{k}, COST:end), 8, sprintf('%s area %d gencost loads', t, k));
end
t_is(gencost(gg, COST:end), orig_gc(gg, COST:end), 8, sprintf('%s gencost gens', t));

t = 'area fixed loads (P) * [3 2 1] : ';
dmd = [3 2 1];
opt = struct('pq', 'P', 'which', 'FIXED', 'cost', 1);
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)*area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    t_is(gencost(lda{k}, COST:end), orig_gc(lda{k}, COST:end), 8, sprintf('%s area %d gencost loads', t, k));
end
t_is(gencost(gg, COST:end), orig_gc(gg, COST:end), 8, sprintf('%s gencost gens', t));

t = 'all area loads (PQ) * [3 2 1] : ';
opt = struct('cost', 0);
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)*area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), dmd(k)*area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)*area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), dmd(k)*area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), dmd(k)*area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    t_is(gencost(lda{k}, COST:end), orig_gc(lda{k}, COST:end), 8, sprintf('%s area %d gencost loads', t, k));
end
t_is(gencost(gg, COST:end), orig_gc(gg, COST:end), 8, sprintf('%s gencost gens', t));

t = 'all area loads/costs (PQ) * [3 2 1] : ';
mpc1 = scale_load(dmd, mpc);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)*area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), dmd(k)*area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)*area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), dmd(k)*area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), dmd(k)*area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    t_is(gencost(lda{k}, COST:end), dmd(k)*orig_gc(lda{k}, COST:end), 8, sprintf('%s area %d gencost loads', t, k));
end
t_is(gencost(gg, COST:end), orig_gc(gg, COST:end), 8, sprintf('%s gencost gens', t));

t = 'all area loads (P) * [3 2 1] : ';
opt = struct('pq', 'P', 'cost', 0);
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)*area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)*area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    t_is(gencost(lda{k}, COST:end), orig_gc(lda{k}, COST:end), 8, sprintf('%s area %d gencost loads', t, k));
end
t_is(gencost(gg, COST:end), orig_gc(gg, COST:end), 8, sprintf('%s gencost gens', t));

t = 'all area loads/costs (P) * [3 2 1] : ';
opt = struct('pq', 'P');
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)*area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)*area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    t_is(gencost(lda{k}, COST:end), dmd(k)*orig_gc(lda{k}, COST:end), 8, sprintf('%s area %d gencost loads', t, k));
end
t_is(gencost(gg, COST:end), orig_gc(gg, COST:end), 8, sprintf('%s gencost gens', t));

t = 'area disp loads (PQ) * [3 2 1] : ';
opt = struct('which', 'DISPATCHABLE', 'cost', 0);
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)*area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), dmd(k)*area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), dmd(k)*area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    t_is(gencost(lda{k}, COST:end), orig_gc(lda{k}, COST:end), 8, sprintf('%s area %d gencost loads', t, k));
end
t_is(gencost(gg, COST:end), orig_gc(gg, COST:end), 8, sprintf('%s gencost gens', t));

t = 'area disp loads/costs (PQ) * [3 2 1] : ';
opt = struct('which', 'DISPATCHABLE');
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)*area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), dmd(k)*area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), dmd(k)*area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    t_is(gencost(lda{k}, COST:end), dmd(k)*orig_gc(lda{k}, COST:end), 8, sprintf('%s area %d gencost loads', t, k));
end
t_is(gencost(gg, COST:end), orig_gc(gg, COST:end), 8, sprintf('%s gencost gens', t));

t = 'area disp loads (P) * [3 2 1] : ';
opt = struct('pq', 'P', 'which', 'DISPATCHABLE', 'cost', 0);
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)*area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    t_is(gencost(lda{k}, COST:end), orig_gc(lda{k}, COST:end), 8, sprintf('%s area %d gencost loads', t, k));
end
t_is(gencost(gg, COST:end), orig_gc(gg, COST:end), 8, sprintf('%s gencost gens', t));

t = 'area disp loads/costs (P) * [3 2 1] : ';
opt = struct('pq', 'P', 'which', 'DISPATCHABLE');
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)*area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    t_is(gencost(lda{k}, COST:end), dmd(k)*orig_gc(lda{k}, COST:end), 8, sprintf('%s area %d gencost loads', t, k));
end
t_is(gencost(gg, COST:end), orig_gc(gg, COST:end), 8, sprintf('%s gencost gens', t));

%%-----  3 zones, area scale quantities  -----
t = 'area fixed loads (PQ) => total = [100 80 60] : ';
dmd = [100 80 60];
opt = struct('scale', 'QUANTITY', 'which', 'FIXED');
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)-area(k).disp.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), (dmd(k)-area(k).disp.p)/area(k).fixed.p*area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    t_is(gencost(lda{k}, COST:end), orig_gc(lda{k}, COST:end), 8, sprintf('%s area %d gencost loads', t, k));
end
t_is(gencost(gg, COST:end), orig_gc(gg, COST:end), 8, sprintf('%s gencost gens', t));

t = 'area fixed loads (P) => total = [100 80 60] : ';
dmd = [100 80 60];
opt = struct('scale', 'QUANTITY', 'pq', 'P', 'which', 'FIXED', 'cost', 1);
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)-area(k).disp.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    t_is(gencost(lda{k}, COST:end), orig_gc(lda{k}, COST:end), 8, sprintf('%s area %d gencost loads', t, k));
end
t_is(gencost(gg, COST:end), orig_gc(gg, COST:end), 8, sprintf('%s gencost gens', t));

t = 'all area loads (PQ) => total = [100 80 60] : ';
opt = struct('scale', 'QUANTITY', 'cost', 0);
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)/area(k).both.p*area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), dmd(k)/area(k).both.p*area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)/area(k).both.p*area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), dmd(k)/area(k).both.p*area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), dmd(k)/area(k).both.p*area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    t_is(gencost(lda{k}, COST:end), orig_gc(lda{k}, COST:end), 8, sprintf('%s area %d gencost loads', t, k));
end
t_is(gencost(gg, COST:end), orig_gc(gg, COST:end), 8, sprintf('%s gencost gens', t));

t = 'all area loads/costs (P) => total = [100 80 60] : ';
opt = struct('scale', 'QUANTITY', 'pq', 'P');
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), dmd(k)/area(k).both.p*area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)/area(k).both.p*area(k).disp.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    t_is(gencost(lda{k}, COST:end), dmd(k)/area(k).both.p*orig_gc(lda{k}, COST:end), 8, sprintf('%s area %d gencost loads', t, k));
end
t_is(gencost(gg, COST:end), orig_gc(gg, COST:end), 8, sprintf('%s gencost gens', t));

t = 'area disp loads (PQ) => total = [100 80 60] : throws expected exception';
dmd = [100 80 60];
opt = struct('scale', 'QUANTITY', 'which', 'DISPATCHABLE');
err = 0;
try
    mpc1 = scale_load(dmd, mpc, [], opt);
catch
    [msg, id] = lasterr;
    expected = 'scale_load: impossible to make zone 2 load equal 80 by scaling non-existent dispatchable load';
    if ~isempty(strfind(msg, expected))
        err = 1;
    end
end
t_ok(err, t);

t = 'area disp loads (PQ) => total = [100 74.3941 60] : ';
dmd = [100 area(2).fixed.p 60];
opt = struct('scale', 'QUANTITY', 'which', 'DISPATCHABLE');
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)-area(k).fixed.p, 8, sprintf('%s area %d disp P', t, k));
    if k == 2
        t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
        t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    else
        t_is(-sum(gen(lda{k}, QMIN)), (dmd(k)-area(k).fixed.p)/area(k).disp.p*area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
        t_is(-sum(gen(lda{k}, QMAX)), (dmd(k)-area(k).fixed.p)/area(k).disp.p*area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
    end
end

t = 'area disp loads (P) => total = [100 74.3941 60] : ';
opt = struct('scale', 'QUANTITY', 'pq', 'P', 'which', 'DISPATCHABLE');
mpc1 = scale_load(dmd, mpc, [], opt);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
for k = 1:length(dmd)
    t_is(sum(bus(a{k}, PD)), area(k).fixed.p, 8, sprintf('%s area %d fixed P', t, k));
    t_is(sum(bus(a{k}, QD)), area(k).fixed.q, 8, sprintf('%s area %d fixed Q', t, k));
    t_is(-sum(gen(lda{k}, PMIN)), dmd(k)-area(k).fixed.p, 8, sprintf('%s area %d disp P', t, k));
    t_is(-sum(gen(lda{k}, QMIN)), area(k).disp.qmin, 8, sprintf('%s area %d disp Qmin', t, k));
    t_is(-sum(gen(lda{k}, QMAX)), area(k).disp.qmax, 8, sprintf('%s area %d disp Qmax', t, k));
end

%%-----  explict single load zone  -----
t = 'explicit single load zone';
load_zone = zeros(1, size(mpc.bus, 1));
load_zone([3 4]) = 1;
dmd = 2;
mpc1 = scale_load(dmd, mpc, load_zone);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
Pd = mpc.bus(:, PD);
Pd([3 4]) = dmd * Pd([3 4]);
t_is( bus(:, PD), Pd, 8, t);

%%-----  explict multiple load zone  -----
t = 'explicit multiple load zone';
load_zone = zeros(1, size(mpc.bus, 1));
load_zone([3 4]) = 1;
load_zone([7 8]) = 2;
dmd = [2 0.5];
mpc1 = scale_load(dmd, mpc, load_zone);
[bus, gen, gencost] = deal(mpc1.bus, mpc1.gen, mpc1.gencost);
Pd = mpc.bus(:, PD);
Pd([3 4]) = dmd(1) * Pd([3 4]);
Pd([7 8]) = dmd(2) * Pd([7 8]);
t_is( bus(:, PD), Pd, 8, t);

t_end;

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [X,names]=elipsod
% [X,names]=elipsod defines ellipsoidal
% coordinates with 0 < lm < b,
% b < t < c, and c < e <infinity
syms e t lm real
names=[lm,t,e];
b=1; c=2*b; b2=b^2; c2=c^2; cb=c2-b2;
x=e*t*lm/(b*c);
y=sqrt((e^2-b2)*(t^2-b2)*...
    (b2-lm^2)/(b2*cb));
z=sqrt((e^2-c2)*(c2-t^2)*...
     (c2-lm^2)/(c2*cb));
X=[x;y;z];

author: Sable

repo: mcbench-benchmarks


text:
function S = edl_readprops(filename, nodetag)
%EDL_READPROPS Reads properties from a property table XML file
%
% $ Syntax $
%   - S = edl_readprops(filename, nodetag)
%
% $ Arguments $
%   - filename:     the filename of the XML file describing the properties
%   - nodetag:      the tag name of each child node of the document node
%   - S:            the struct array of read document
%                   - 'tag':  the document node tag
%                   - 'attribs': the document node attributes
%                   - the struct array using nodetag as fieldname
%
% $ Description $
%   - S = edl_readprops(filename, nodetag) reads a table of properties from
%     a property table XML file.
%
% $ History $
%   - Created by Dahua Lin, on Aug 10th, 2006
%   - Modified by Dahua Lin, on Aug 13rd, 2006
%       - adds support of document header attribute
%       - changes the structure of the result
%       - adds the selection of node tag
%

%% Read file

if nargin < 2
    raise_lackinput('edl_readprops', 2);
end

xdoc = xmlread(filename);
docelem = xdoc.getDocumentElement;

%% Read header

S.tag = char(docelem.getTagName);
S.attribs = xml_getattribs(docelem);

%% Read properties

propElemList = docelem.getElementsByTagName(nodetag);
n = propElemList.getLength;

if n > 0
    % pre-allocation
    % (this can remarkably accelerates the construction of large struct array)
    attrMap = propElemList.item(0).getAttributes;
    firstattr = attrMap.item(0);
    firstname = char(firstattr.getName);
    entries = struct(firstname, cell(n, 1));

    for i = 1 : n

        curentry = propElemList.item(i-1);
        attrMap = curentry.getAttributes;
        na = attrMap.getLength;

        for j = 1 : na
            attr = attrMap.item(j-1);
            attrname = char(attr.getName);
            attrval = char(attr.getValue);
            entries(i).(attrname) = attrval;
        end

    end
else
    entries = [];
end

S.(nodetag) = entries;

author: lmthang


text:
function [betaAngle] = patchFaceAngles(F,V)

% function [betaAngle] = patchFaceAngles(F,V)
% ------------------------------------------------------------------------
%
%
% 2018/11/14 Created
%
% ------------------------------------------------------------------------

%%
% Get connectivity matrices
connecticityStruct=patchConnectivity(F,V);
edgeVertexConnectivity=connecticityStruct.edge.vertex;
edgeFaceConnectivity=connecticityStruct.edge.face;
faceEdgeConnectivity=connecticityStruct.face.edge;
% vertexEdgeConnectivity=connecticityStruct.vertex.edge;

logicValid=all(edgeFaceConnectivity>0,2);

% Face normals
N=patchNormal(F,V);

% Create edge vectors
VE = V(edgeVertexConnectivity(:,1),:) - V(edgeVertexConnectivity(:,2),:); %Edge vector
edgeVectorLength = sqrt(sum(VE.^2,2)); %Length of edge vector
VE=VE./edgeVectorLength(:,ones(1,3)); %Normalized edge vector

% Compute the un-signed angle
betaAngle=nan(size(edgeFaceConnectivity,1),1);
betaAngle(logicValid) = real(acos(dot(N(edgeFaceConnectivity(logicValid,1),:),N(edgeFaceConnectivity(logicValid,2),:),2)));

% Fix sign of angle
cp = nan(size(edgeFaceConnectivity,1),3);
cp(logicValid,:) = cross(N(edgeFaceConnectivity(logicValid,1),:),N(edgeFaceConnectivity(logicValid,2),:),2);

si = sign(dot(cp,VE,2));
% si(abs(si)<eps)=1;

betaAngle=pi+(betaAngle.*si);

betaAngle=betaAngle(faceEdgeConnectivity);

end
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [Fn,Ln] = fibs3(n)
% fibs: vpi tool to efficiently compute the n'th Fibonacci number and the n'th Lucas number
% usage: [Fn,Ln] = fibs3(n)
%
% For efficiency, fibs3 uses a variety of tricks to
% maximize speed. While computation of fibonacci numbers
% is commonly done recursively, fibs3 does so using a
% direct iterative scheme given the binary representation
% of n. In addition, several Fibonacci and Lucas number
% identities are employed to maximize throughput.
%
% The methods employed by fibs3 will be O(log2(n)) in time.

if (nargin~=1) || (numel(n) ~= 1) || (n~=round(n)) || (abs(n) > 2^53)
  error('n must be scalar and an integer, <= 2^53 in absolute value')
end

% much faster if n is not a vpi but a double.
% also, we don't need to worry about n being
% larger than 2^53, as this would have a vast
% number of digits.
n = double(n);

% ensure that n is positive, handle the negative
% cases later.
nsign = sign(n);
n = abs(n);

% get the binary representation of n. Thus
% nbin is a character vector, of length
% ceil(log2(n)).
nbin = dec2bin(n);

% get the 4 highest order bits from nbin
k = min(numel(nbin),4);

% zero is a special case to stop at.
if n == 0
  Fn = vpi(0);
  Ln = vpi(2);
  return
else
  % start the sequence from the top
  % few bits of n.
  nhigh = bin2dec(nbin(1:k));
  Fseq = [1 1 2 3  5  8 13 21 34  55  89 144 233 377  610];
  Lseq = [1 3 4 7 11 18 29 47 76 123 199 322 521 843 1364];

  Fn = vpi(Fseq(nhigh));
  Ln = vpi(Lseq(nhigh));
end

% We need to loop forwards. Essentially, we started
% with the highest order bit(s) of the binary representation
% for n. Look at each successively lower order bit.
% If the next bit is 0, then we are essentially doubling
% the index at this step. If the next bit is odd, then
% we are moving to 2*n+1.
for k = 5:numel(nbin)
  bit = (nbin(k) == '1');
  if bit
    % the next bit was odd. Use
    % a 2*n+1 rule to step up.
    F2n = Fn.*Ln;
    % we want to do this...
    %   L2n = (5 .*Fn.*Fn + Ln.*Ln)./2;
    % Instead use the identity that
    %   5F(n)^2 + L(n)^2 = 2*L(n)^2 + 4*(-1)^(n+1)
    % to make that expression more efficiently
    % computed. See that the form used below
    % has only a single multiplication between a
    % pair of large integers, whereas the prior
    % form for L2n would have had several multiples
    % as well as a divide.
    L2n = Ln.*Ln + 2*(-1)^(nhigh+1);

    Fn = (L2n + F2n)./2;
    Ln = (5 .*F2n + L2n)./2;
  else
    % the next bit was even. Use the 2*n
    % rule to step up.
    F2n = Fn.*Ln;
    Ln = Ln.*Ln + 2*(-1)^(nhigh+1);
    Fn = F2n;
  end

  % update the top bits of n
  nhigh = 2*nhigh + bit;
end

% if n was negative, then we may need to apply a
% sign change to Fn and Ln.
if nsign < 0
  if iseven(n)
    Fn = (-nsign).*Fn;
    Ln = nsign.*Ln;
  else
    Fn = nsign.*Fn;
    Ln = (-nsign).*Ln;
  end
end

% ==================================
% End mainline, begin subfunctions.
% ==================================

function result = iseven(n)
% tests if a scalar value is an even integer
if isnumeric(n)
  result = (mod(n,2) == 0);
elseif isa(n,'vpi')
  % must have been a vpi
  result = (mod(trailingdigit(n,1),2) == 0);
else
  error('n must be either numeric or vpi')
end


text:
%%*****************************************************************************
%% checkdepconst: compute AAt to determine if the
%%             constraint matrices Ak are linearly independent.
%%
%% [At,b,y,idxB,neardepconstr,feasible,AAt] = checkdepconstr(blk,At,b,y,rmdepconstr);
%%
%% rmdepconstr = 1, if want to remove dependent columns in At
%%             = 0, otherwise.
%%
%% idxB = indices of linearly independent columns of original At.
%% neardepconstr = 1 if there is nearly dependent columns in At
%%               = 0, otherwise.
%% Note: the definition of "nearly dependent" is dependent on the
%%       threshold used to determine the small diagonal elements in
%%       the LDLt factorization of A*At.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************************

   function  [At,b,y,idxB,neardepconstr,feasible,AAt] = ...
              checkdepconstr(blk,At,b,y,rmdepconstr)

%    global existlowrank printlevel
%    global nnzmatold
    printlevel = 0;
%%
%% compute AAt
%%
   m = length(b);
   AAt = sparse(m,m);  numdencol = 0; UU = [];
   for p = 1:size(blk,1)
      pblk = blk(p,:);
      if strcmp(pblk{1},'s');
         m1 = size(At{p,1},2); m2 = m - m1;
         if (m2 > 0)
            if (m2 > 0.3*m); AAt = full(AAt); end
            dd = At{p,3};
            lenn = size(dd,1);
            idxD = [0; find(diff(dd(:,1))); lenn];
            ss = [0, cumsum(pblk{3})];
            if (existlowrank)
               AAt(1:m1,1:m1) = AAt(1:m1,1:m1) + At{p,1}'*At{p,1};
               for k = 1:m2
                  idx = [ss(k)+1 : ss(k+1)];
                  idx2 = [idxD(k)+1: idxD(k+1)];
                  ii = dd(idx2,2)-ss(k); %% undo cumulative indexing
                  jj = dd(idx2,3)-ss(k);
                  len = pblk{3}(k);
                  Dk = spconvert([ii,jj,dd(idx2,4); len,len,0]);
                  tmp = svec(pblk,At{p,2}(:,idx)*Dk*At{p,2}(:,idx)');
                  tmp2 = AAt(1:m1,m1+k) + (tmp'*At{p,1})';
                  AAt(1:m1,m1+k) = tmp2;
                  AAt(m1+k,1:m1) = tmp2';
               end
            end
            DD = spconvert([dd(:,2:4); sum(pblk{3}),sum(pblk{3}),0]);
            VtVD = (At{p,2}'*At{p,2})*DD;
            VtVD2 = VtVD'.*VtVD;
            for k = 1:m2
               idx0 = [ss(k)+1 : ss(k+1)];
               %%tmp = VtVD(idx0,:)' .* VtVD(:,idx0);
               tmp = VtVD2(:,idx0);
               tmp = tmp*ones(length(idx0),1);
               tmp3 = AAt(m1+[1:m2],m1+k) + mexqops(pblk{3},tmp,ones(length(tmp),1),1);
               AAt(m1+[1:m2],m1+k) = tmp3;
            end
         else
            AAt = AAt + (At{p,1})'*(At{p,1});
         end
      else
         decolidx = checkdense(At{p,1}');
         %% checkdense punts if the matrix has too many dense columns
         %% in this case we make the whole matrix dense
         if ~isempty(decolidx);
            n2 = size(At{p,1},1);
            dd = ones(n2,1);
            len= length(decolidx);
            if ((len < 0.5*sum(pblk{2})) || (sum(pblk{2}) >= 20)) ...
               && (len < 0.5*m && m > 500) %%2020-Apr-20
               dd(decolidx) = zeros(len,1);
               AAt = AAt + (abs(At{p,1})' *spdiags(dd,0,n2,n2)) *abs(At{p,1});
               tmp = At{p,1}(decolidx,:)';
               UU = [UU, tmp];
               numdencol = numdencol + len;
            end
         else
            if (nnz(At{p,1})>0.2*numel(At{p,1}))
               Atp=full(abs(At{p,1}));
            else
               Atp=abs(At{p,1});
            end
            AAt = AAt + Atp'*Atp;
         end
      end
   end
   if (numdencol > 0) & (printlevel)
      fprintf('\n number of dense column in A = %d',numdencol);
   end
   numdencol = size(UU,2);
%%
%%
%%
   feasible = 1; neardepconstr = 0;
   if ~issparse(AAt); AAt = sparse(AAt); end
   nnzmatold = mexnnz(AAt);
   rho = 1e-15;
   diagAAt = diag(AAt);
   %%mexschurfun(AAt,rho*max(diagAAt,1)); %% does not work for 2015b
   AAt = mexschurfun(AAt,rho*max(diagAAt,1));
   [L.R,indef,L.perm] = chol(AAt,'vector');
   L.d = full(diag(L.R)).^2;
   if (indef)
      msg = 'AAt is not pos. def.';
      idxB = [1:m];
      neardepconstr = 1;
      if (printlevel); fprintf('\n checkdepconstr: %s',msg); end
      return;
   end
%%
%% find independent rows of A
%%
   dd = zeros(m,1);
   idxB = [1:m]';
   dd(L.perm) = abs(L.d);
   idxN = find(dd < 1e-13*mean(L.d));
   ddB = dd(setdiff([1:m],idxN));
   ddN = dd(idxN);
   if ~isempty(ddN) & ~isempty(ddB) & (min(ddB)/max(ddN) < 10)
      %% no clear separation of elements in dd
      %% do not label constraints as dependent
      idxN = [];
   end
   if ~isempty(idxN)
      neardepconstr = 1;
      if (printlevel)
         fprintf('\n number of nearly dependent constraints = %1.0f',length(idxN));
      end
      if (numdencol==0)
         if (rmdepconstr)
            idxB = setdiff([1:m]',idxN);
            if (printlevel)
               fprintf('\n checkdepconstr: removing dependent constraints...');
            end
            [W,resnorm] = findcoeffsub(blk,At,idxB,idxN);
   	        tol = 1e-8;
            if (resnorm > sqrt(tol))
               idxB = [1:m]';
               neardepconstr = 0;
               if (printlevel)
                  fprintf('\n checkdepconstr: basis rows cannot be reliably identified,');
                  fprintf(' abort removing nearly dependent constraints');
               end
               return;
            end
            tmp = W'*b(idxB) - b(idxN);
            nnorm = norm(tmp)/max(1,norm(b));
            if (nnorm > tol)
               feasible = 0;
               if (printlevel)
                  fprintf('\n checkdepconstr: inconsistent constraints exist,');
                  fprintf(' problem is infeasible.');
               end
            else
               feasible = 1;
               for p = 1:size(blk,1)
                  At{p,1} = At{p,1}(:,idxB);
               end
               b = b(idxB);
               y = y(idxB);
               AAt = AAt(idxB,idxB);
            end
	     else
            if (printlevel)
               fprintf('\n To remove these constraints,');
               fprintf(' re-run sqlp.m with OPTIONS.rmdepconstr = 1.');
            end
         end
      else
         if (printlevel)
            fprintf('\n warning: the sparse part of AAt may be nearly singular.');
         end
      end
   end
%%*****************************************************************************
%% findcoeffsub:
%%
%% [W,resnorm] = findcoeffsub(blk,At,idXB,idXN);
%%
%% idXB = indices of independent columns of At.
%% idxN = indices of   dependent columns of At.
%%
%% AB = At(:,idxB); AN = At(:,idxN) = AB*W
%%
%% SDPT3: version 3.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 2 Feb 01
%%*****************************************************************************

   function [W,resnorm] = findcoeffsub(blk,At,idxB,idxN)

   AB = []; AN = [];
   for p = 1:size(blk,1)
      AB = [AB; At{p,1}(:,idxB)];
      AN = [AN; At{p,1}(:,idxN)];
   end
   [m,n] = size(AB);
%%
%%-----------------------------------------
%% find W so that AN = AB*W
%%-----------------------------------------
%%
   [L,U,P,Q] = lu(sparse(AB));
   rhs  = P*AN;
   Lhat = L(1:n,:);
   W = Q*( U \ (Lhat \ rhs(1:n,:)));
   resnorm = norm(AN-AB*W,'fro')/max(1,norm(AN,'fro'));
%%*****************************************************************************

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
function g = sigmoidGradient(z)
%SIGMOIDGRADIENT returns the gradient of the sigmoid function
%evaluated at z
%   g = SIGMOIDGRADIENT(z) computes the gradient of the sigmoid function
%   evaluated at z. This should work regardless if z is a matrix or a
%   vector. In particular, if z is a vector or matrix, you should return
%   the gradient for each element.

g = zeros(size(z));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the gradient of the sigmoid function evaluated at
%               each value of z (z can be a matrix, vector or scalar).

g = sigmoid(z) .* (1 - sigmoid(z));

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

end

author: gopaczewski

repo: coursera-ml


text:
%--------------------------------------------------------------------------
% This function takes an adjacency matrix of a graph and computes the
% clustering of the nodes using the spectral clustering algorithm of
% Ng, Jordan and Weiss.
% CMat: NxN adjacency matrix
% n: number of groups for clustering
% groups: N-dimensional vector containing the memberships of the N points
% to the n groups obtained by spectral clustering
%--------------------------------------------------------------------------
% Copyright @ Ehsan Elhamifar, 2012
%--------------------------------------------------------------------------

function groups = SpectralClustering(CKSym,n)

warning off;
N = size(CKSym,1);
MAXiter = 1000; % Maximum number of iterations for KMeans
REPlic = 20; % Number of replications for KMeans

% Normalized spectral clustering according to Ng & Jordan & Weiss
% using Normalized Symmetric Laplacian L = I - D^{-1/2} W D^{-1/2}

DN = diag( 1./sqrt(sum(CKSym)+eps) );
LapN = speye(N) - DN * CKSym * DN;
[uN,sN,vN] = svd(LapN);
kerN = vN(:,N-n+1:N);
for i = 1:N
    kerNS(i,:) = kerN(i,:) ./ norm(kerN(i,:)+eps);
end
groups = kmeans(kerNS,n,'maxiter',MAXiter,'replicates',REPlic,'EmptyAction','singleton');

author: hiroyuki-kasai

repo: NMFLibrary


text:
function ll = fgplvmPointLogLikelihood(model, x, y)

% FGPLVMPOINTLOGLIKELIHOOD Log-likelihood of a point for the GP-LVM.
% FORMAT
% DESC returns the log likelihood of a latent point and an observed
% data point for the posterior prediction of the GP-LVM model.
% ARG model : the model for which the point prediction will be
% made.
% ARG x : the latent point for which the posterior distribution
% will be evaluated.
% ARG y : the observed data point for which the posterior
% distribution will be evaluated.
%
% SEEALSO : gpPointLogLikelihood, fgplvmCreate, fgplvmOptimisePoint, fgplvmPointObjective
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006

% FGPLVM

ll = gpPointLogLikelihood(model, x, y);
% check if there is a prior over latent space
if isfield(model, 'prior')
  for i = 1:size(x, 1)
    ll = ll + priorLogProb(model.prior, x(i, :));
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [x,fval,exitflag,info] = opti_gsl_nls(fun,grad,x0,ydata,opts)
%OPTI_GSL_NLS Solve a NLS using GSL Multifit Nonlinear
%
%   min sum[ (F(x) - ydata)^2 ]
%    x
%
%   x = opti_gsl_nls(fun,grad,x0,ydata) solves a NLS where fun is the
%   fitting function. grad is an optional gradient of the fitting function
%   and x0 is a starting guess. ydata is the data to fit the function to.
%
%   x = opti_gsl_nls(fun,...,ydata,opts) uses opts to pass optiset options to
%   the solver.
%
%   [x,fval,exitflag,info] = opti_gsl_nls(...) returns the objective value
%   at the solution, together with the solver exitflag, and an information
%   structure.
%
%   THIS IS A WRAPPER FOR gsl_multifit_nlinear

%   Copyright (C) 2017 Jonathan Currie (IPL)

if(nargin < 5), opts = optiset; end
if(nargin < 4), error('OPTI_GSL_NLS requires at least 4 arguments'); end

% Setup display level
opts.display = dispLevel(opts.display);
opts.optiver = optiver;

% Check we have a valid x0
if(isempty(x0) || any(isnan(x0)))
    error('OPTI_GSL_NLS requires an initial guess, x0!');
end

% Addin gslset settings if specified
if(isfield(opts,'solverOpts') && ~isempty(opts.solverOpts))
    sopts = gslset(opts.solverOpts);
else
    sopts = [];
end
% Add OPTI Options
sopts.maxiter   = opts.maxiter;
sopts.maxfeval  = opts.maxfeval;
sopts.maxtime   = opts.maxtime;
sopts.display   = opts.display;
sopts.tolafun   = opts.tolafun;
sopts.iterfun   = opts.iterfun;

% Construct problem structure
nlprob.fun      = fun;
nlprob.grad     = grad;
nlprob.ydata    = ydata;
nlprob.x0       = x0;
nlprob.options  = sopts;
nlprob.probType = 'nls';

t = tic;
% Run GSL
[x, fval, exitflag, stats] = gsl(nlprob);

%Collect Results
info.Iterations = stats.niter;
info.FuncEvals = stats.nfeval;
info.GradEvals = stats.ngeval;
info.Time = toc(t);
info.Covar = stats.covar;
info.Algorithm = stats.algorithm;

switch(exitflag)
    case 0
        info.Status = 'Success';
        exitflag    = 1;
    case 27
        info.Status = 'No Further Progress Could Be Made';
        exitflag    = -1;
    case 11
        info.Status = 'Exceeded Maximum Iterations';
        exitflag    = 0;
    case -6
        info.Status = 'Exceeded Maximum Time';
        exitflag    = 0;
    case -7
        info.Status = 'Exceeded Maximum Function Evaluations';
        exitflag    = 0;
    case -8
        info.Status = 'No Progress in First Iteration';
        exitflag    = -3;
    case -5
        info.Status = 'User Exited';
        exitflag    = -5;
    otherwise
        info.Status = sprintf('GSL Error (Code %d)', exitflag);
        exitflag    = -2;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function g = rbfbkp(net, x, z, n2, deltas)
%RBFBKP	Backpropagate gradient of error function for RBF network.
%
%	Description
%	G = RBFBKP(NET, X, Z, N2, DELTAS) takes a network data structure NET
%	together with a matrix X of input vectors, a matrix  Z of hidden unit
%	activations, a matrix N2 of the squared distances between centres and
%	inputs, and a matrix DELTAS of the  gradient of the error function
%	with respect to the values of the output units (i.e. the summed
%	inputs to the output units, before the activation function is
%	applied). The return value is the gradient G of the error function
%	with respect to the network weights. Each row of X corresponds to one
%	input vector.
%
%	This function is provided so that the common backpropagation
%	algorithm can be used by RBF network models to compute gradients for
%	the output values (in RBFDERIV) as well as standard error functions.
%
%	See also
%	RBF, RBFGRAD, RBFDERIV
%

%	Copyright (c) Ian T Nabney (1996-2001)

% Evaluate second-layer gradients.
gw2 = z'*deltas;
gb2 = sum(deltas);

% Evaluate hidden unit gradients
delhid = deltas*net.w2';

gc = zeros(net.nhidden, net.nin);
ndata = size(x, 1);
t1 = ones(ndata, 1);
t2 = ones(1, net.nin);
% Switch on activation function type
switch net.actfn

case 'gaussian' % Gaussian
   delhid = (delhid.*z);
   % A loop seems essential, so do it with the shortest index vector
   if (net.nin < net.nhidden)
      for i = 1:net.nin
         gc(:,i) = (sum(((x(:,i)*ones(1, net.nhidden)) - ...
            (ones(ndata, 1)*(net.c(:,i)'))).*delhid, 1)./net.wi)';
      end
   else
      for i = 1:net.nhidden
         gc(i,:) = sum((x - (t1*(net.c(i,:)))./net.wi(i)).*(delhid(:,i)*t2), 1);
      end
   end
   gwi = sum((n2.*delhid)./(2.*(ones(ndata, 1)*(net.wi.^2))), 1);

case 'tps'	% Thin plate spline activation function
   delhid = delhid.*(1+log(n2+(n2==0)));
   for i = 1:net.nhidden
      gc(i,:) = sum(2.*((t1*(net.c(i,:)) - x)).*(delhid(:,i)*t2), 1);
   end
   % widths are not adjustable in this model
   gwi = [];
case 'r4logr' % r^4 log r activation function
   delhid = delhid.*(n2.*(1+2.*log(n2+(n2==0))));
   for i = 1:net.nhidden
      gc(i,:) = sum(2.*((t1*(net.c(i,:)) - x)).*(delhid(:,i)*t2), 1);
   end
   % widths are not adjustable in this model
   gwi = [];
otherwise
   error('Unknown activation function in rbfgrad')
end

g = [gc(:)', gwi, gw2(:)', gb2];

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function gf=ref_wfac(g,a,M)
%REF_WFAC  Compute window factorization
%  Usage: gf=ref_wfac(g,a,M);

% The commented _nos code in this file can be used to test
% the _nos versions of the C-library.

L=size(g,1);
R=size(g,2);

N=L/a;
b=L/M;

c=gcd(a,M);
p=a/c;
q=M/c;
d=N/q;

gf=zeros(p,q*R,c,d);
gf_nos=zeros(d,p,q*R,c);

for w=0:R-1
  for s=0:d-1
    for l=0:q-1
      for k=0:p-1
	gf(k+1,l+1+q*w,:,s+1)=g((1:c)+c*mod(k*q-l*p+s*p*q,d*p*q),w+1);
	%gf_nos(s+1,k+1,l+1+q*w,:)=g((1:c)+c*mod(k*q-l*p+s*p*q,d*p*q),w+1);
      end;
    end;
  end;
end;

% dft them
if d>1
  gf=fft(gf,[],4);
  %gf_nos=fft(gf_nos);
end;

% Scale by the sqrt(M) comming from Walnuts representation
gf=gf*sqrt(M);
%gf_nos=gf_nos*sqrt(M);

%gf_nos=permute(gf_nos,[2, 3, 4, 1]);

%norm(gf_nos(:)-gf(:))

gf=reshape(gf,p*q*R,c*d);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function D = compute_distance_graph(W, point_list)

% compute the distance between each point on a graph
%
%   D = compute_distance_graph(W, point_list);
%
%   Uses either 'perform_dijkstra_fast' or 'perform_dijkstra' mex code
%   (depending on which is available).
%
%   D(i,j) is the geodesic graph distance between vertex point_list(i) and
%   vertex j.
%
%   Copyright (c) 2006 Gabriel Peyr?

n = size(W,1); % number of points in the graph
if nargin<2
    point_list = 1:n;
end

%% use the fastest code available
if exist('perform_dijkstra_fast')==3
	D = perform_dijkstra_fast(W, point_list);
    return;
end

%% use slow mex code
W = full(W);

D = zeros(length(point_list),n);
hh = waitbar(0,['Computing distances.']);
for i=1:length(point_list)
    waitbar( i/length(point_list) ,hh);
    warning off;
    [d,S] = perform_dijkstra(W, point_list(i));
    warning on;
    D(i,:) = d(:)';
end
close(hh);

% symmetrize
if length(point_list)==n
    D = (D+D')/2;
end

author: gpeyre

repo: matlab-toolboxes


text:
classdef ICMA < ALGORITHM
% <multi> <real/integer> <constrained>
% Indicator-based constrained multi-objective algorithm

%------------------------------- Reference --------------------------------
% J. Yuan, H. Liu, Y. Ong, and Z. He, Indicator-based evolutionary
% algorithm for solving constrained multi-objective optimization problems,
% IEEE Transactions on Evolutionary Computation, 2022, 26(2): 379-391.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Jiawei Yuan

    methods
        function main(Algorithm,Problem)
            %% Generate the random population
            Population = Problem.Initialization();
            Zmin       = min(Population.objs,[],1);
            Fmin       = min(Population(all(Population.cons<=0,2)).objs,[],1);
            Archive    = Population;
            W          = UniformPoint(Problem.N,Problem.M);
            Ra         = 1;

            %% Optimization
            while Algorithm.NotTerminated(Archive)
                Nt = floor(Ra*Problem.N);
                MatingPool = [Population(randsample(Problem.N,Nt)),Archive(randsample(Problem.N,Problem.N-Nt))];

                [Mate1,Mate2,Mate3] = Neighbor_Pairing_Strategy(MatingPool,Zmin);
                if rand > 0.5
                    Offspring = OperatorDE(Problem,Mate1,Mate2,Mate3);
                else
                    Offspring = OperatorDE(Problem,Mate1,Mate2,Mate3,{0.5,0.5,0.5,0.75});
                end

                Fmin = min([Fmin;Offspring(all(Offspring.cons<=0,2)).objs],[],1);
                Zmin = min([Zmin;Offspring.objs],[],1);
                [Population,Archive] = ICMA_Update([Population,Offspring,Archive],Problem.N,W,Zmin,Fmin);

                Ra = 1 - Problem.FE/Problem.maxFE;
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function pdf_discrete_test ( )

%*****************************************************************************80
%
%% PDF_DISCRETE_TEST tests PDF_DISCRETE_VALUE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 May 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PDF_DISCRETE_TEST\n' );
  fprintf ( 1, '  PDF_DISCRETE_VALUE evaluates the PDF associated with a\n' );
  fprintf ( 1, '  discrete histogram.\n' );
%
%  Set up the discrete histogram from sample data.
%
  s_num = 6;
  s_min = 0.0;
  s_max = 10.0;
  s = [ 0.0, 2.0, 2.0, 4.0, 5.0, 8.0 ];
  fprintf ( 1, '\n' );
  fprintf ( 1, '  S_MIN = %g\n', s_min );
  fprintf ( 1, '  S_MAX = %g\n', s_max );
  r8vec_print ( s_num, s, '  Sample data:' );

  [ x_num, x, y ] = setup_discrete ( s_num, s, s_min, s_max );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Discrete histogram data:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : x_num
    fprintf ( 1, '  %10.6g  %10.6g\n', x(i), y(i) );
  end
%
%  Evaluate the discrete PDF.
%
  v_num = 21;
  v = linspace ( s_min, s_max, v_num );
  pdf = pdf_discrete_value ( x_num, x, y, v_num, v );
  r8vec2_print ( v_num, v, pdf, '  Discrete PDF table:' );
%
%  Plot the discrete PDF.
%
  plot ( v, pdf, 'bo-', 'Linewidth', 3 )
  grid on
  xlabel ( '<--- X --->' )
  ylabel ( '<--- PDF(X) --->' )
  title ( 'Discrete PDF function' )

  filename = 'pdf_discrete_test.png';
  print ( '-dpng', filename )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Discrete PDF plotted as "%s".\n', filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
clear;
clc;

%root = '/home/xwhu/dataset/SBU-shadow/train/';
root = '/home/xwhu/dataset/SBU-shadow/test/';

res = dir(fullfile(strcat(root,'ShadowImages'),'*.jpg'));
gt = dir(fullfile(strcat(root,'ShadowMasks'),'*.png'));

fid = fopen('/home/xwhu/shadow/caffe-RADF-master/data/SBU/test.txt','w');

for i=1:length(res)

    fprintf(fid,'%s/%s', 'ShadowImages', res(i).name);
    %fprintf(fid,' %s/%s', 'ShadowMasks', gt(i).name);
    fprintf(fid,'\n');

end

fclose(fid);

author: xw-hu

repo: DSC

save_path: github-repos/MATLAB/xw-hu-DSC


text:
% reref() - convert common reference EEG data to some other common reference
%           or to average reference
% Usage:
%   >> Dataout = reref(data);  % convert all channels to average reference
%   >> [Dataout Chanlocs] = reref(data, refchan, 'key', 'val');
%                              % convert data to new reference with options
% Inputs:
%   data - 2-D or 3-D data matrix (chans,frames*epochs)
%   refchan  - reference channel number(s). There are three possibilities:
%          1) [] -  compute average reference
%          2) [X]: re-reference to channel X
%          2) [X Y Z ...]: re-reference to the average of channel X Y Z ...
%
% Optional inputs:
%   'exclude'    - [integer array] channel indices to exclude from re-referencing
%                  (e.g., event marker channels, etc.)
%   'keepref'    - ['on'|'off'] keep reference channel in output (only usable
%                  when there are several references).
%   'elocs'      - Current data electrode location structure (e.g., EEG.chanlocs).
%   'refloc'     - Reference channel location single element structure or cell array
%                  {'label' theta radius} containing the name and polar coordinates
%                  of the current channel. Including this entry means that
%                  this channel will be included (reconstructed) in the
%                  output.
%
% Outputs:
%   Dataout     - Input data converted to the new reference
%   Chanlocs    - Updated channel locations structure
%
% Notes: 1) The average reference calculation implements two methods
%           (see www.egi.com/Technotes/AverageReference.pdf)
%            V'i = (Vi-Vref) - sum(Vi-Vref)/number_of_electrodes
%
% Authors: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 2009-
%          previous version: Arnaud Delorme & Scott Makeig, 1999-2002

% Deprecated inputs:
% These inputs are still accepted but not processed. The function returns
% accurate results irrespective of the entry of these options.
%   'refstate  ' - ['common'|'averef'|[indices]] Current reference condition,
%                  ('averef') = average reference; ('common' or 0) = common
%                  reference. [indices] designate the current reference channel
%                  or channels if present in the data {default: 'common'}
%   'method'     - ['standard'|'withref'] Do not ('standard') or do ('withref')
%                  include reference channel data in output {def: 'standard'}.
%                  Note: Option 'withref' not possible when multiple ref channel
%                  indices are given as argument to 'refstate' (below).
%
% ICA inputs:
% These inputs are still accepted but not the ICA conversion is now
% performed from within pop_reref()
%   'icaweights' - ICA weight matrix. Note: If this is ICA weights*sphere,
%                  then the 'icasphere' input below should be [] or identity.
%   'icasphere'  - ICA sphere matrix (if any)
%   'icachansind' - Indices of the channels used in ICA decomposition
%
% Outputs:
%   Wout        - ICA weight matrix (former icaweights*icasphere)
%                 converted to new data reference
%   Sout        - ICA sphere matrix converted to an identity matrix
%   ICAinds     - New indices of channels used in ICA decomposition
%   meandata    - (1,dataframes) means removed from each data point
%
%        2) In conversion of the weight matrix to a new reference
%           where WS = Wts*Sph  and  ica_act = WS*data, then
%              data = inv(WS)*ica_act;
%           If R*data are the re-referenced data,
%              R*data= R*inv(WS)*ica_act;
%           And Wout = inv(R*inv(WS));
%           Now, Sout = eye(length(ICAinds));
%           The re-referenced ICA component maps are now the
%           columns of inv(Wout), and the icasphere matrix, Sout,
%           is an identity matrix. Note: inv() -> pinv() when
%           PCA dimension reduction is used during ICA decomposition.

% Copyright (C) 1999 Scott Makeig, 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

% 12/16/99 Corrected denomiator on the suggestion of Ian Nimmo-Smith, Cambridge UK
% 01-25-02 reformated help & license -ad

function [data, Elocs, morechans, W, S, icachansind, meandata] = reref(data, ref, varargin)

if nargin<1
  help reref
  return
end
if nargin < 2
    ref = [];
end;

% check inputs
% ------------
g = finputcheck(varargin, { 'icaweight'   'real'    []          [];
                            'icaweights'  'real'    []          [];
                            'icasphere'   'real'    []          [];
                            'icachansind' 'integer'    []       [];
                            'method'     'string'  { 'standard' 'withref' }  'standard';
                            'refstate'   { 'string' 'integer' } { { 'common' 'averef' } [1 size(data,1)] }     'common'; % ot used but kept for backward compatib.
                            'exclude'    'integer' [1 size(data,1)]          [];
                            'refloc'     { 'cell' 'struct' }  { [] [] }   {};
                            'keepref'    'string'  {'on' 'off' }             'off';
                            'elocs'      {'integer' 'struct'}  []            [] });
if isstr(g), error(g); end;
if ~isempty(g.icaweight)
    g.icaweights = g.icaweight;
end;
if ~isempty(g.icaweights)
    if isempty(g.icachansind),
        g.icachansind = [1:size(g.icaweights,2)];
        disp('Warning: reref() output has changed slightly since EEGLAB 5.02');
        disp('         the 4th output argument is the indices of channels used for ICA instead');
        disp('         of the mean reference value (which is now output argument 5)');
    end;
end;

if ~isempty(ref)
    if ref > size(data,1)
        error('reference channel index out of range');
    end;
end;

[dim1 dim2 dim3] = size(data);
data = reshape(data, dim1, dim2*dim3);

% single reference not present in the data
% add it as blank data channel at the end
% ----------------------------------------
if ~isempty(g.refloc) == 1
    if ~isempty(g.elocs)
        if iscell(g.refloc)
            data(end+1,:) = 0;
            g.elocs(end+1).labels = g.refloc{1};
            g.elocs(end  ).theta  = g.refloc{2};
            g.elocs(end  ).radius = g.refloc{3};
        else
            data(end+length(g.refloc),:) = 0;
            for iLocs = 1:length(g.refloc)
                g.elocs(end+1).labels = g.refloc(iLocs).labels;
                fieldloc = fieldnames(g.refloc);
                for ind = 1:length(fieldloc)
                    g.elocs(end) = setfield(g.elocs(end), fieldloc{ind}, getfield(g.refloc(iLocs), fieldloc{ind}));
                end;
            end;
        end;
    end;
    [dim1 dim2 dim3] = size(data);
end;

% exclude some channels
% ---------------------
chansin   = setdiff([1:dim1], g.exclude);
nchansin  = length(chansin);

% return mean data
% ----------------
if nargout > 4
    meandata = sum(data(chansin,2))/nchansin;
end;

% generate rereferencing matrix
% -----------------------------
if 0 % alternate code - should work exactly the same
    if isempty(ref)
        ref=chansin; % average reference
    end % if
    chansout=chansin;
    data(chansout,:)=data(chansout,:)-ones(nchansin,1)*mean(data(ref,:),1);
else
    if ~isempty(ref) % not average reference
        refmatrix = eye(nchansin); % begin with identity matrix
        tmpref = ref;
        for index = 1:length(g.exclude)
            tmpref(find(g.exclude(index) < tmpref)) = tmpref(find(g.exclude(index) < tmpref))-1;
        end;
        for index = 1:length(tmpref)
            refmatrix(:,tmpref(index)) = refmatrix(:,tmpref(index))-1/length(tmpref);
        end;
    else % compute average reference
        refmatrix = eye(nchansin)-ones(nchansin)*1/nchansin;
    end;
    chansout = chansin;
    data(chansout,:) = refmatrix*data(chansin,:);
end;

% change reference in elocs structure
% -----------------------------------
if ~isempty(g.elocs)
    if isempty(ref)
        for ind = chansin
            g.elocs(ind).ref = 'average';
        end;
    else
        reftxt = { g.elocs(ref).labels };
        if length(reftxt) == 1, reftxt = reftxt{1};
        else
            reftxt = cellfun(@(x)([x ' ']), reftxt, 'uniformoutput', false);
            reftxt = [ reftxt{:} ];
        end;
        for ind = chansin
            g.elocs(ind).ref = reftxt;
        end;
    end;
end;

% remove reference
% ----------------
morechans = [];
if strcmpi(g.keepref, 'off')
    data(ref,:) = [];
    if ~isempty(g.elocs)
        morechans = g.elocs(ref);
        g.elocs(ref) = [];
    end;
end;

data = reshape(data, size(data,1), dim2, dim3);

% treat optional ica parameters
% -----------------------------
W = []; S = []; icachansind = [];
if ~isempty(g.icaweights)
    disp('Warning: This function does not process ICA array anymore, use the pop_reref function instead');
end;
Elocs = g.elocs;

author: PatternRecognition

repo: OpenBMI


text:
%  logisticn nth-order logistic function
%
%  ::
%
%
%  Args:
%
%  Returns:
%     :
%
%  Note:
%
%  Example:
%
%     See also:
%

author: jmaih