text:
function t_mplinsolve(quiet)
%T_MPLINSOLVE  Tests of MIPS/MATPOWER linear solvers.

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

if nargin < 1
    quiet = 0;
end

t_begin(180, quiet);

isoctave = exist('OCTAVE_VERSION', 'builtin') == 5;
if isoctave
    lu_warning_id = 'Octave:lu:sparse_input';
    s = warning('query', lu_warning_id);
    warning('off', lu_warning_id);
end

%% non-sparse A matrix crashes some MATLAB lu() calls on certain
%% combinations of MATLAB and macOS versions
mlv = have_feature('matlab', 'vnum');
skipcrash = ~isempty(mlv) && mlv < 8.003 && mlv > 7.013 && ...
    strcmp(computer, 'MACI64');

ijs = [
    1 1 1205.63;
    4 1 -1205.63;
    25 1 17.3611;
    28 1 -17.3611;
    43 1 1;
    2 2 1024;
    8 2 -1024;
    26 2 16;
    32 2 -16;
    3 3 1164.84;
    6 3 -1164.84;
    27 3 17.0648;
    30 3 -17.0648;
    1 4 -1205.63;
    4 4 2201.34;
    5 4 -453.695;
    9 4 -542.009;
    13 4 0.776236;
    14 4 -0.451721;
    18 4 -0.324515;
    25 4 -17.3611;
    28 4 39.4759;
    29 4 -10.5107;
    33 4 -11.6041;
    37 4 -3.30738;
    38 4 1.94219;
    42 4 1.36519;
    4 5 -453.695;
    5 5 581.372;
    6 5 -127.677;
    13 5 -0.451721;
    14 5 0.568201;
    15 5 -0.11648;
    28 5 -10.5107;
    29 5 16.0989;
    30 5 -5.58824;
    37 5 1.94219;
    38 5 -3.2242;
    39 5 1.28201;
    3 6 -1164.84;
    5 6 -127.677;
    6 6 1676.74;
    7 6 -384.227;
    14 6 -0.11648;
    15 6 0.163059;
    16 6 -0.0465791;
    27 6 -17.0648;
    29 6 -5.58824;
    30 6 32.4374;
    31 6 -9.78427;
    38 6 1.28201;
    39 6 -2.4371;
    40 6 1.15509;
    6 7 -384.227;
    7 7 1141.16;
    8 7 -756.935;
    15 7 -0.0465791;
    16 7 0.371828;
    17 7 -0.325249;
    30 7 -9.78427;
    31 7 23.4822;
    32 7 -13.698;
    39 7 1.15509;
    40 7 -2.77221;
    41 7 1.61712;
    2 8 -1024;
    7 8 -756.935;
    8 8 1925.77;
    9 8 -144.836;
    16 8 -0.325249;
    17 8 0.844105;
    18 8 -0.518855;
    26 8 -16;
    31 8 -13.698;
    32 8 35.6731;
    33 8 -5.97513;
    40 8 1.61712;
    41 8 -2.80473;
    42 8 1.1876;
    4 9 -542.009;
    8 9 -144.836;
    9 9 686.845;
    13 9 -0.324515;
    17 9 -0.518855;
    18 9 0.84337;
    28 9 -11.6041;
    32 9 -5.97513;
    33 9 17.5792;
    37 9 1.36519;
    41 9 1.1876;
    42 9 -2.55279;
    10 10 1207.63;
    13 10 -1205.63;
    34 10 17.3611;
    37 10 -17.3611;
    11 11 1026;
    17 11 -1024;
    35 11 16;
    41 11 -16;
    12 12 1166.84;
    15 12 -1164.84;
    36 12 17.0648;
    39 12 -17.0648;
    4 13 0.776236;
    5 13 -0.451721;
    9 13 -0.324515;
    10 13 -1205.63;
    13 13 2190.83;
    14 13 -447.892;
    18 13 -535.137;
    28 13 3.30738;
    29 13 -1.94219;
    33 13 -1.36519;
    34 13 -17.3611;
    37 13 39.1419;
    38 13 -10.5107;
    42 13 -11.6041;
    4 14 -0.451721;
    5 14 0.568201;
    6 14 -0.11648;
    13 14 -447.892;
    14 14 573.221;
    15 14 -122.862;
    28 14 -1.94219;
    29 14 3.2242;
    30 14 -1.28201;
    37 14 -10.5107;
    38 14 15.5829;
    39 14 -5.58824;
    5 15 -0.11648;
    6 15 0.163059;
    7 15 -0.0465791;
    12 15 -1164.84;
    14 15 -122.862;
    15 15 1669.87;
    16 15 -379.651;
    29 15 -1.28201;
    30 15 2.4371;
    31 15 -1.15509;
    36 15 -17.0648;
    38 15 -5.58824;
    39 15 31.8704;
    40 15 -9.78427;
    6 16 -0.0465791;
    7 16 0.371828;
    8 16 -0.325249;
    15 16 -379.651;
    16 16 1131.92;
    17 16 -750.073;
    30 16 -1.15509;
    31 16 2.77221;
    32 16 -1.61712;
    39 16 -9.78427;
    40 16 23.1242;
    41 16 -13.698;
    7 17 -0.325249;
    8 17 0.844105;
    9 17 -0.518855;
    11 17 -1024;
    16 17 -750.073;
    17 17 1914.92;
    18 17 -138.499;
    31 17 -1.61712;
    32 17 2.80473;
    33 17 -1.1876;
    35 17 -16;
    40 17 -13.698;
    41 17 35.2181;
    42 17 -5.97513;
    4 18 -0.324515;
    8 18 -0.518855;
    9 18 0.84337;
    13 18 -535.137;
    17 18 -138.499;
    18 18 676.012;
    28 18 -1.36519;
    32 18 -1.1876;
    33 18 2.55279;
    37 18 -11.6041;
    41 18 -5.97513;
    42 18 17.0972;
    19 19 1.88667;
    25 19 -1;
    20 20 1.54931;
    26 20 -1;
    21 21 1.78346;
    27 21 -1;
    22 22 0.666667;
    34 22 -1;
    23 23 0.666667;
    35 23 -1;
    24 24 0.666667;
    36 24 -1;
    1 25 17.3611;
    4 25 -17.3611;
    19 25 -1;
    2 26 16;
    8 26 -16;
    20 26 -1;
    3 27 17.0648;
    6 27 -17.0648;
    21 27 -1;
    1 28 -17.3611;
    4 28 39.4759;
    5 28 -10.5107;
    9 28 -11.6041;
    13 28 3.30738;
    14 28 -1.94219;
    18 28 -1.36519;
    4 29 -10.5107;
    5 29 16.0989;
    6 29 -5.58824;
    13 29 -1.94219;
    14 29 3.2242;
    15 29 -1.28201;
    3 30 -17.0648;
    5 30 -5.58824;
    6 30 32.4374;
    7 30 -9.78427;
    14 30 -1.28201;
    15 30 2.4371;
    16 30 -1.15509;
    6 31 -9.78427;
    7 31 23.4822;
    8 31 -13.698;
    15 31 -1.15509;
    16 31 2.77221;
    17 31 -1.61712;
    2 32 -16;
    7 32 -13.698;
    8 32 35.6731;
    9 32 -5.97513;
    16 32 -1.61712;
    17 32 2.80473;
    18 32 -1.1876;
    4 33 -11.6041;
    8 33 -5.97513;
    9 33 17.5792;
    13 33 -1.36519;
    17 33 -1.1876;
    18 33 2.55279;
    10 34 17.3611;
    13 34 -17.3611;
    22 34 -1;
    11 35 16;
    17 35 -16;
    23 35 -1;
    12 36 17.0648;
    15 36 -17.0648;
    24 36 -1;
    4 37 -3.30738;
    5 37 1.94219;
    9 37 1.36519;
    10 37 -17.3611;
    13 37 39.1419;
    14 37 -10.5107;
    18 37 -11.6041;
    4 38 1.94219;
    5 38 -3.2242;
    6 38 1.28201;
    13 38 -10.5107;
    14 38 15.5829;
    15 38 -5.58824;
    5 39 1.28201;
    6 39 -2.4371;
    7 39 1.15509;
    12 39 -17.0648;
    14 39 -5.58824;
    15 39 31.8704;
    16 39 -9.78427;
    6 40 1.15509;
    7 40 -2.77221;
    8 40 1.61712;
    15 40 -9.78427;
    16 40 23.1242;
    17 40 -13.698;
    7 41 1.61712;
    8 41 -2.80473;
    9 41 1.1876;
    11 41 -16;
    16 41 -13.698;
    17 41 35.2181;
    18 41 -5.97513;
    4 42 1.36519;
    8 42 1.1876;
    9 42 -2.55279;
    13 42 -11.6041;
    17 42 -5.97513;
    18 42 17.0972;
    1 43 1;
];
A = sparse(ijs(:, 1), ijs(:, 2), ijs(:, 3));
b = [
    0;
    0;
    0;
    0;
    0;
    0;
    0;
    0;
    0;
    0;
    0;
    0;
    -0.00896054;
    -0.0617829;
    -0.0772931;
    -0.0230638;
    -0.0185934;
    -0.02;
    -0.336;
    -0.2755;
    -0.353;
    0;
    0;
    0;
    1.3;
    1.55;
    1.4;
    0;
    -0.9;
    0;
    -1;
    0;
    -1.25;
    0;
    0;
    0;
    0.167;
    -0.042;
    0.2835;
    -0.171;
    0.2275;
    -0.259;
    0;
];
% ex = [      %% for original symmetrical system
% 0;
% 0.04612219791119214;
% 0.0509334351882598;
% -0.05953571031927305;
% -0.09461814898578046;
% -0.008909854578010561;
% -0.05785829019394401;
% -0.02232729212460287;
% -0.1137760871247425;
% -0.03062777824802364;
% -0.01013282572376477;
% -0.005330939680091628;
% -0.02914165388753019;
% -0.03376073204420303;
% 4.021341450281111e-05;
% -0.01727289094763518;
% -0.008382063634320435;
% -0.04854008812629265;
% -0.2663945795760685;
% -0.4548081594272799;
% -0.3787862287965494;
% -0.02580075363496279;
% -0.02801219343110931;
% -0.09165765332863518;
% -0.1665986614487812;
% -0.4291388294822789;
% -0.3225500876094941;
% 3.967864472606993;
% 5.577372159790927;
% 3.762341481664236;
% 5.858342308599034;
% 3.951628532808602;
% 6.471657726723339;
% -0.01720051102355974;
% -0.01867480495813735;
% -0.06110513277164123;
% -0.1239317474796463;
% -0.1848327901217308;
% -0.4283638085291988;
% 0.07640167050820651;
% -0.1319901818980452;
% 0.4366406661538687;
% 0.0007894844305239289;
% ];
A(8, 2) = -500;     %% make non-symmetrical
ex = [
   0.000000000000000
   0.015818764480650
   0.028627515055857
  -0.066470494209706
  -0.106338553472755
  -0.029556247483918
  -0.080778493974592
  -0.046876059940052
  -0.126748712178182
  -0.032582479109501
  -0.010646589107665
  -0.006489850319503
  -0.030082640379903
  -0.034809123585434
  -0.001497869375252
  -0.018761742104354
  -0.009797350675461
  -0.049351378784847
  -0.145999102975876
  -0.546882809268771
  -0.407105729011254
  -0.043399950168419
  -0.013587814915265
  -0.085187156417456
   0.060547872388504
  -0.571791005228199
  -0.373056783462411
   4.199136375632448
   5.730774734397821
   3.598555773194044
   5.477255481510613
   3.440677757696717
   6.533433414753445
  -0.028933314578930
  -0.009058547806115
  -0.056791466007356
  -0.206286450370091
  -0.246526980342321
  -0.398303747647950
   0.132430861414235
  -0.064740631105635
   0.368923140091142
  -8.288373070379230
];

t = ''''' : ';
x = mplinsolve(A, b, '');
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
x = mplinsolve(full(A), b, '');
t_is(x, ex, 12, [t 'x (full A)']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);

t = '\ : ';
x = mplinsolve(A, b, '\');
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
x = mplinsolve(full(A), b, '\');
t_is(x, ex, 12, [t 'x (full A)']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);

t = '\ (transpose) : ';
opt = struct('tr', 1);
x = mplinsolve(A', b, '\', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
x = mplinsolve(full(A'), b, '\', opt);
t_is(x, ex, 12, [t 'x (full A)']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);

t = 'LU : ';
x = mplinsolve(A, b, 'LU');
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    x = mplinsolve(full(A), b, 'LU');
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end

t = 'LU (transpose) : ';
opt = struct('tr', 1);
x = mplinsolve(A', b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    x = mplinsolve(full(A'), b, 'LU', opt);
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end

t = 'LU3 : ';
x = mplinsolve(A, b, 'LU3');
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    x = mplinsolve(full(A), b, 'LU3');
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end
t = 'LU, nout = 3, vec = 1, thresh = 1 : ';
opt = struct('nout', 3, 'vec', 1, 'thresh', 1);
x = mplinsolve(A, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    x = mplinsolve(full(A), b, 'LU', opt);
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end

t = 'LU3a : ';
x = mplinsolve(A, b, 'LU3a');
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    x = mplinsolve(full(A), b, 'LU3a');
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end
t = 'LU, nout = 3, vec = 1 : ';
opt = struct('nout', 3, 'vec', 1);
x = mplinsolve(A, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    x = mplinsolve(full(A), b, 'LU', opt);
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end

t = 'LU4 : ';
x = mplinsolve(A, b, 'LU4');
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
t = 'LU, nout = 4, vec = 1 : ';
opt = struct('nout', 4, 'vec', 1);
x = mplinsolve(A, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);

t = 'LU5 : ';
x = mplinsolve(A, b, 'LU5');
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
t = 'LU, nout = 5, vec = 1 : ';
opt = struct('nout', 5, 'vec', 1);
x = mplinsolve(A, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);

t = 'LU3m : ';
x = mplinsolve(A, b, 'LU3m');
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    x = mplinsolve(full(A), b, 'LU3m');
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end
t = 'LU, nout = 3, vec = 0, thresh = 1 : ';
opt = struct('nout', 3, 'vec', 0, 'thresh', 1);
x = mplinsolve(A, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    x = mplinsolve(full(A), b, 'LU', opt);
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end

t = 'LU3am : ';
x = mplinsolve(A, b, 'LU3am');
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    x = mplinsolve(full(A), b, 'LU3am');
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end
t = 'LU, nout = 3, vec = 0 : ';
opt = struct('nout', 3, 'vec', 0);
x = mplinsolve(A, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    x = mplinsolve(full(A), b, 'LU', opt);
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end

t = 'LU4m : ';
x = mplinsolve(A, b, 'LU4m');
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
t = 'LU, nout = 4, vec = 0 : ';
opt = struct('nout', 4, 'vec', 0);
x = mplinsolve(A, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);

t = 'LU5m : ';
x = mplinsolve(A, b, 'LU5m');
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
t = 'LU, nout = 5, vec = 0 : ';
opt = struct('nout', 5, 'vec', 0);
x = mplinsolve(A, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);

%%-----  repeat all LU tests with saved factors  -----
t = 'pre-factored LU : ';
f = mplinsolve(A);
x = mplinsolve(f, b);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    f = mplinsolve(full(A));
    x = mplinsolve(f, b);
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end

t = 'pre-factored LU3 : ';
f = mplinsolve(A, [], 'LU3');
x = mplinsolve(f, b, 'LU3');
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    f = mplinsolve(full(A), [], 'LU3');
    x = mplinsolve(f, b, 'LU3');
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end
t = 'pre-factored LU, nout = 3, vec = 1, thresh = 1 : ';
opt = struct('nout', 3, 'vec', 1, 'thresh', 1);
[x, f] = mplinsolve(A, b, 'LU', opt);
x = mplinsolve(f, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    opt = struct('nout', 3, 'vec', 1, 'thresh', 1);
    [x, f] = mplinsolve(full(A), b, 'LU', opt);
    x = mplinsolve(f, b, 'LU', opt);
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end

t = 'pre-factored LU3a : ';
[x1, f] = mplinsolve(A, b, 'LU3a');
x = mplinsolve(f, b, 'LU3a');
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    [x1, f] = mplinsolve(full(A), b, 'LU3a');
    x = mplinsolve(f, b, 'LU3a');
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end
t = 'pre-factored LU, nout = 3, vec = 1 : ';
opt = struct('nout', 3, 'vec', 1);
f = mplinsolve(A, [], 'LU', opt);
x = mplinsolve(f, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    opt = struct('nout', 3, 'vec', 1);
    f = mplinsolve(full(A), [], 'LU', opt);
    x = mplinsolve(f, b, 'LU', opt);
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end

t = 'pre-factored LU4 : ';
[x, f] = mplinsolve(A, b, 'LU4');
x = mplinsolve(f, b, 'LU4');
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
t = 'pre-factored LU, nout = 4, vec = 1 : ';
opt = struct('nout', 4, 'vec', 1);
f = mplinsolve(A, [], 'LU', opt);
x = mplinsolve(f, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);

t = 'pre-factored LU5 : ';
f = mplinsolve(A, [], 'LU5');
x = mplinsolve(f, b, 'LU5');
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
t = 'pre-factored LU, nout = 5, vec = 1 : ';
opt = struct('nout', 5, 'vec', 1);
[x, f] = mplinsolve(A, b, 'LU', opt);
x = mplinsolve(f, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);

t = 'pre-factored LU3m : ';
f = mplinsolve(A, [], 'LU3m');
x = mplinsolve(f, b, 'LU3m');
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    f = mplinsolve(full(A), [], 'LU3m');
    x = mplinsolve(f, b, 'LU3m');
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end
t = 'pre-factored LU, nout = 3, vec = 0, thresh = 1 : ';
opt = struct('nout', 3, 'vec', 0, 'thresh', 1);
[x, f] = mplinsolve(A, b, 'LU', opt);
x = mplinsolve(f, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    opt = struct('nout', 3, 'vec', 0, 'thresh', 1);
    f = mplinsolve(full(A), [], 'LU', opt);
    x = mplinsolve(f, b, 'LU', opt);
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end

t = 'pre-factored LU3am : ';
f = mplinsolve(A, [], 'LU3am');
x = mplinsolve(f, b, 'LU3am');
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    f = mplinsolve(full(A), [], 'LU3am');
    x = mplinsolve(f, b, 'LU3am');
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end
t = 'pre-factored LU, nout = 3, vec = 0 : ';
opt = struct('nout', 3, 'vec', 0);
f = mplinsolve(A, [], 'LU', opt);
x = mplinsolve(f, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    opt = struct('nout', 3, 'vec', 0);
    [x, f] = mplinsolve(full(A), b, 'LU', opt);
    x = mplinsolve(f, b, 'LU', opt);
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end

t = 'pre-factored LU4m : ';
[x1, f] = mplinsolve(A, b, 'LU4m');
x = mplinsolve(f, b, 'LU4m');
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
t = 'pre-factored LU, nout = 4, vec = 0 : ';
opt = struct('nout', 4, 'vec', 0);
f = mplinsolve(A, [], 'LU', opt);
x = mplinsolve(f, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);

t = 'pre-factored LU5m : ';
f = mplinsolve(A, [], 'LU5m');
x = mplinsolve(f, b, 'LU5m');
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
t = 'pre-factored LU, nout = 5, vec = 0 : ';
opt = struct('nout', 5, 'vec', 0);
[x, f] = mplinsolve(A, b, 'LU', opt);
x = mplinsolve(f, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);

%%-----  repeat all LU tests with saved factors AND transpose  -----
t = 'transposed LU : ';
f = mplinsolve(A');
x = mplinsolve(f, b, 'LU', struct('tr', 1));
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    f = mplinsolve(full(A'));
    x = mplinsolve(f, b, 'LU', struct('tr', 1));
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end

t = 'transposed LU3 : ';
f = mplinsolve(A', [], 'LU3');
x = mplinsolve(f, b, 'LU3', struct('tr', 1));
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    [x1, f] = mplinsolve(full(A'), b, 'LU3');
    x = mplinsolve(f, b, 'LU3', struct('tr', 1));
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end
t = 'transposed LU, nout = 3, vec = 1, thresh = 1 : ';
opt = struct('nout', 3, 'vec', 1, 'thresh', 1);
[x, f] = mplinsolve(A', b, 'LU', opt);
opt.tr = 1;
x = mplinsolve(f, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    opt = struct('nout', 3, 'vec', 1, 'thresh', 1);
    [x, f] = mplinsolve(full(A'), b, 'LU', opt);
    opt.tr = 1;
    x = mplinsolve(f, b, 'LU', opt);
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end

t = 'transposed LU3a : ';
f = mplinsolve(A', [], 'LU3a');
x = mplinsolve(f, b, 'LU3a', struct('tr', 1));
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    f = mplinsolve(full(A'), [], 'LU3a');
    x = mplinsolve(f, b, 'LU3a', struct('tr', 1));
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end
t = 'transposed LU, nout = 3, vec = 1 : ';
opt = struct('nout', 3, 'vec', 1);
f = mplinsolve(A', [], 'LU', opt);
opt.tr = 1;
x = mplinsolve(f, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    opt = struct('nout', 3, 'vec', 1);
    [x, f] = mplinsolve(full(A'), b, 'LU', opt);
    opt.tr = 1;
    x = mplinsolve(f, b, 'LU', opt);
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end

t = 'transposed LU4 : ';
[x1, f] = mplinsolve(A', b, 'LU4');
x = mplinsolve(f, b, 'LU4', struct('tr', 1));
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
t = 'transposed LU, nout = 4, vec = 1 : ';
opt = struct('nout', 4, 'vec', 1);
f = mplinsolve(A', [], 'LU', opt);
opt.tr = 1;
x = mplinsolve(f, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);

t = 'transposed LU5 : ';
f = mplinsolve(A', [], 'LU5');
x = mplinsolve(f, b, 'LU5', struct('tr', 1));
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
t = 'transposed LU, nout = 5, vec = 1 : ';
opt = struct('nout', 5, 'vec', 1);
f = mplinsolve(A', [], 'LU', opt);
opt.tr = 1;
x = mplinsolve(f, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);

t = 'transposed LU3m : ';
[x1, f] = mplinsolve(A', b, 'LU3m');
x = mplinsolve(f, b, 'LU3m', struct('tr', 1));
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    f = mplinsolve(full(A'), [], 'LU3m');
    x = mplinsolve(f, b, 'LU3m', struct('tr', 1));
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end
t = 'transposed LU, nout = 3, vec = 0, thresh = 1 : ';
opt = struct('nout', 3, 'vec', 0, 'thresh', 1);
f = mplinsolve(A', [], 'LU', opt);
opt.tr = 1;
x = mplinsolve(f, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    opt = struct('nout', 3, 'vec', 0, 'thresh', 1);
    [x, f] = mplinsolve(full(A'), b, 'LU', opt);
    opt.tr = 1;
    x = mplinsolve(f, b, 'LU', opt);
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end

t = 'transposed LU3am : ';
[x1, f] = mplinsolve(A', b, 'LU3am');
x = mplinsolve(f, b, 'LU3am', struct('tr', 1));
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    f = mplinsolve(full(A'), [], 'LU3am');
    x = mplinsolve(f, b, 'LU3am', struct('tr', 1));
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end
t = 'transposed LU, nout = 3, vec = 0 : ';
opt = struct('nout', 3, 'vec', 0);
f = mplinsolve(A', [], 'LU', opt);
opt.tr = 1;
x = mplinsolve(f, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
if skipcrash
    t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
else
    opt = struct('nout', 3, 'vec', 0);
    f = mplinsolve(full(A'), [], 'LU', opt);
    opt.tr = 1;
    x = mplinsolve(f, b, 'LU', opt);
    t_is(x, ex, 12, [t 'x (full A)']);
    t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
end

t = 'transposed LU4m : ';
f = mplinsolve(A', [], 'LU4m');
x = mplinsolve(f, b, 'LU4m', struct('tr', 1));
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
t = 'transposed LU, nout = 4, vec = 0 : ';
opt = struct('nout', 4, 'vec', 0);
f = mplinsolve(A', [], 'LU', opt);
opt.tr = 1;
x = mplinsolve(f, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);

t = 'transposed LU5m : ';
[x1, f] = mplinsolve(A', b, 'LU5m');
x = mplinsolve(f, b, 'LU5m', struct('tr', 1));
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
t = 'transposed LU, nout = 5, vec = 0 : ';
opt = struct('nout', 5, 'vec', 0);
[x, f] = mplinsolve(A', b, 'LU', opt);
opt.tr = 1;
x = mplinsolve(f, b, 'LU', opt);
t_is(x, ex, 12, [t 'x']);
t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);

%% PARDISO
if have_feature('pardiso')
    if have_feature('pardiso_object')
        tols = [6 5 12 12 6 5];     %% tolerances for PARDISO v6
    else
        tols = [13 13 12 12 1 2];   %% tolerances for PARDISO v5
    end
    vb = false;

    t = 'PARDISO (direct) : ';
    opt = struct('pardiso', struct('solver', 0, 'verbose', vb));
    x = mplinsolve(A, b, 'PARDISO', opt);
    t_is(x, ex, tols(1), [t 'x']);
    t_is(norm(b - A*x), 0, tols(2), [t '||b - A*x||']);

    t = 'PARDISO (direct, transpose) : ';
    opt = struct('pardiso', struct('solver', 0, 'tr', 1, 'verbose', vb));
    x = mplinsolve(A', b, 'PARDISO', opt);
    t_is(x, ex, tols(1), [t 'x']);
    t_is(norm(b - A*x), 0, tols(2), [t '||b - A*x||']);

    t = 'PARDISO (direct, symmetric indefinite) : ';
    opt = struct('pardiso', struct('solver', 0, 'mtype', -2, 'verbose', vb));
    x = mplinsolve(A, b, 'PARDISO', opt);
    t_is(x, ex, tols(3), [t 'x']);
    t_is(norm(b - A*x), 0, tols(4), [t '||b - A*x||']);

    t = 'PARDISO (iterative) : ';
    opt = struct('pardiso', struct('solver', 1, 'verbose', vb));
    [x, info] = mplinsolve(A, b, 'PARDISO', opt);
    t_is(x, ex, tols(5), [t 'x']);
    t_is(norm(b - A*x), 0, tols(6), [t '||b - A*x||']);
else
    t_skip(8, ['PARDISO not available']);
end

if isoctave
    warning(s.state, lu_warning_id);
end

t_end;


text:
clear all;
close all;
clc;
for i=1:50
image_left = ['../../road/ir/',num2str(i),'.jpg'];
image_right = ['../../road/vi/',num2str(i),'.jpg'];
fused_path = ['result/',num2str(i),'.bmp'];
image_input1=imread(image_left);
image_input2=imread(image_right);
disp(i);
addpath(genpath('dtcwt_toolbox'));
addpath(genpath('fdct_wrapping_matlab'));
addpath(genpath('nsct_toolbox'));

% [imagename1 imagepath1]=uigetfile('source_images\*.jpg;*.bmp;*.png;*.tif;*.tiff;*.pgm;*.gif','Please choose the first input image');
% image_input1=imread(strcat(imagepath1,imagename1));
% [imagename2 imagepath2]=uigetfile('source_images\*.jpg;*.bmp;*.png;*.tif;*.tiff;*.pgm;*.gif','Please choose the second input image');
% image_input2=imread(strcat(imagepath2,imagename2));

% figure;imshow(image_input1);
% figure;imshow(image_input2);

% if size(image_input1)~=size(image_input2)
%     error('two images are not the same size.');
% end

A=double(image_input1);
B=double(image_input2);

level=4;

tic;
% F = lp_fuse(A, B, level, 3, 3);       %LP
% F = rp_fuse(A, B, level, 3, 3);      %RP
%F = dwt_fuse(A, B, level);           %DWT
% F = dtcwt_fuse(A,B,level);           %DTCWT
F = curvelet_fuse(A,B,level+1);      %CVT
%F = nsct_fuse(A,B,[2,3,3,4]);        %NSCT
toc;

% figure;imshow(uint8(F));
imwrite(uint8(F),fused_path);

end

author: thfylsty

path: github-repos/MATLAB/thfylsty-Classic-and-state-of-the-art-image-fusion-methods/Classic-and-state-of-the-art-image-fusion-methods-5d9457df396f1ea6921e1b9b3703995205940862/GTF/A general framework for image fusion based on multi-scale transform and sparse representation/MST_main.m


text:
function evalRes = benchmarkSSInternal(outp, imSet, param)
% 	outp can be a directory containing the images as mats...
%		or it can be a struct with superpixels, and scores
%	outp.imList is the list of images...
%	outp.superpixels
%	outp.scores
%	param.threshPR, param.threshIU, param.thresh, param.imSet
	fprintf('Benchmarking SS - %s', outp);

	param.imSet = imSet;
	infoFile = param.infoFile;
	thresh = param.thresh;
	threshPR = param.threshPR;
	threshIU = param.threshIU;
	ignoreBck = param.ignoreBck;

	dt = getMetadata(infoFile);
	nclasses = length(dt.className);
	imList = getImageSet(imSet);

	tp = zeros(length(threshPR), nclasses);
	fp = tp; fn = tp; tn = tp;
	rawcounts = zeros([nclasses+1, nclasses+1, length(threshIU)]);

	parfor i = 1:length(imList)
		prob = getProbCube(outp, imList{i}, nclasses);
		gtim = getGroundTruth(imList{i}, 'segmentation', infoFile);
		[TP(:,:,i), FP(:,:,i), TN(:,:,i), FN(:,:,i), RAWCOUNTS(:,:,:,i)] = segmentationMetrics(prob, gtim, threshPR, threshIU);
		fprintf('.');
		[~, ~, ~, ~, rawCountsAll{i}] = segmentationMetrics(prob, gtim, 0, thresh);
	end
	tp = sum(TP,3);
	fp = sum(FP,3);
	tn = sum(TN,3);
	fn = sum(FN,3);
	rawcounts = sum(RAWCOUNTS, 4);
	fprintf('\n');

	%Calculate the classwise APs.
	detailsAP = [];
	if(length(threshPR) > 0)
		P = max(tp,1)./max(tp+fp,1);
		R = max(tp,1)./max(tp+fn,1);
		for i = 1:nclasses,
			ap(i) = VOCap(R(end:-1:1,i), P(end:-1:1,i));
		end
		detailsAP = struct('ap', ap, 'P', P, 'R', R, 'threshPR', threshPR);
	end

	%Calculate the IU over all thresholds..
	detailsIU = [];
	if(length(threshIU) > 0)
		if(ignoreBck)
			rawcounts(1,:,:) = 0;
		end
		for i = 1:length(threshIU)
			rc = rawcounts(:,:,i);
			accuracies(:,i) = diag(rc)./max(sum(rc,2)+sum(rc,1)'-diag(rc),1);
			freq = sum(rc,2);
		end
		accuracies = accuracies(2:end, :);
		freq = freq(2:end, :);
		avacc = mean(accuracies);
		fwavacc = (freq'*accuracies)./sum(freq);
		detailsIU = struct('accuracies', accuracies, 'avacc', avacc, 'fwavacc', fwavacc, 'threshIU', threshIU, 'freq', freq);
	end

	% Calculate for the supplied single threshold
	rawcounts = zeros(size(rawCountsAll{1}));
	for i = 1:length(rawCountsAll),
		rawcounts = rawcounts + full(rawCountsAll{i});
	end

	if(ignoreBck)
		rawcounts(1,:) = 0;
	end
	accuracies = diag(rawcounts)./(max(1,sum(rawcounts,2)+sum(rawcounts,1)'-diag(rawcounts)));
	accuracies = accuracies(2:end);
	freq = sum(rawcounts,2);
	freq = freq(2:end);
	avacc = mean(accuracies);
	fwavacc = (freq'*accuracies)./sum(freq);

	conf = 100*bsxfun(@rdivide, rawcounts, sum(rawcounts,2)+(sum(rawcounts,2)==0));

	className = dt.className;

	evalRes.accuracies = accuracies;
	evalRes.avacc = avacc;
	evalRes.fwavacc = fwavacc;
	evalRes.conf = conf;
	evalRes.rawcounts = rawcounts;
	evalRes.rawCountsAll = rawCountsAll;
	evalRes.imList = imList;
	evalRes.outDir = outp;
	evalRes.className = className;
	evalRes.detailsAP = detailsAP;
	evalRes.detailsIU = detailsIU;
	evalRes.param = param;
	evalRes.thresh = param.thresh;
end

author: s-gupta

repo: rgbd

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


text:
function [inOctave] = is_octave()
% is_octave Test if in GNU/Octave or Matlab
%
% The function returns true only if it is run in GNU/Octave.
%
% Output:
% - inOctave (boolean): True if run in GUN/Octave, false if run in Matlab.
%

try
    OCTAVE_VERSION;
    inOctave = 1;
catch
    inOctave = 0;
end

author: votchallenge

repo: toolkit-legacy


text:
function [c, r, t, cpulse, acq_codes, verbose] = tapas_physio_read_physlogfiles_siemens_tics(...
    log_files, cardiac_modality, verbose, varargin)
% reads out physiological time series of Siemens logfiles with Tics
% The latest implementation of physiological logging in Siemens uses Tics,
% i.e. time stamps of 2.5 ms duration that are reset every day at midnight.
% These are used as a common time scale in all physiological logfiles -
% even though individual sampling times may vary - including cardiac,
% respiratory, pulse oximetry and acquisition time data itself
%
% [c, r, t, cpulse, acq_codes, verbose] = tapas_physio_read_physlogfiles_siemens_tics(...
%    log_files, cardiac_modality, verbose)
%
% IN    log_files
%       .log_cardiac        contains ECG or pulse oximeter time course
%                           for Siemens: *_PULS.log or _ECG[1-4].log.
%       .log_respiration    contains breathing belt amplitude time course
%                           for Siemens: *_RESP.log
%       .sampling_interval  1 entry: sampling interval (seconds)
%                           for both log files
%                           2 entries: 1st entry sampling interval (seconds)
%                           for cardiac logfile, 2nd entry for respiratory
%                           logfile
%                           default: 2.5 ms (1/400 Hz)
%       cardiac_modality    'ECG' or 'PULS'/'PPU'/'OXY' to determine
%                           which channel data to be returned
%                           if not given, will be read out from file name
%                           suffix
%       verbose
%       .level              debugging plots are created if level >=3
%       .fig_handles        appended by handle to output figure
%
%       varargin            propertyName/value pairs, as follows
%           'ecgChannel'    string indicating which channel (combination)
%                           shall be used
%                           'c1' (='PULS'), 'c2', 'c3', 'c4', ...
%                           single channel data is returned
%                           'v12'(='v1'), 'v23', .... v34 (='v2')
%                           voltage difference channel 1 - 2, 2-3 etc. is
%                           returned
%                           'mean'  = 1/2 (v12 + v34);
% OUT
%   cpulse              time events of R-wave peak in cardiac time series (seconds)
%   r                   respiratory time series
%   t                   vector of time points (in seconds)
%   c                   cardiac time series (ECG or pulse oximetry)
%   acq_codes           slice/volume start events marked by number <> 0
%                       for time points in t
%                       10/20 = scan start/end;
%                       1 = ECG pulse; 2 = OXY max; 4 = Resp trigger;
%                       8 = scan volume trigger
%
% EXAMPLE
%   [ons_secs.cpulse, ons_secs.rpulse, ons_secs.t, ons_secs.c] =
%       tapas_physio_read_physlogfiles_siemens_tics(logfiles);
%
%   See also tapas_physio_main_create_regressors

% Author: Lars Kasper
% Created: 2014-09-08
% Copyright (C) 2014 Institute for Biomedical Engineering, ETH/Uni Zurich.
%
% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

%% read out values
DEBUG = verbose.level >= 3;

% channel selection different in ECG data, multiple lines for multiple
% channel, labeled in log file
switch upper(cardiac_modality)
    case 'ECG'
        defaults.ecgChannel = 'v1'; % 'mean'; 'v1'; 'v2', 'v3', 'v4'
    case {'PPU', 'OXY', 'PMU', 'PULS'}
        defaults.ecgChannel = 'PULS';
end

args = tapas_physio_propval(varargin, defaults);
tapas_physio_strip_fields(args);

hasRespirationFile = ~isempty(log_files.respiration);
hasCardiacFile = ~isempty(log_files.cardiac);

% Cardiac and respiratory sampling intervals are ignored, since Tics are
% assumed to be counted in files
if ~isempty(log_files.sampling_interval)
    verbose = tapas_physio_log( ...
        'Ignoring given sampling intervals, using dt = 2.5 ms (tics) instead', ...
        verbose, 1);
end

dtTics = 2.5e-3;
dtCardiac = dtTics;
dtRespiration = dtTics;

%% Define outputs
t           = [];
acq_codes   = [];
tCardiac    = [];
cacq_codes  = [];
c           = [];
cpulse      = [];
tRespiration = [];
racq_codes  = [];
r           = [];
rpulse      = [];

%% Read resp. file
if hasRespirationFile

    C = tapas_physio_read_columnar_textfiles(log_files.respiration, 'RESP');
    nColumns = numel(C); % different file formats indicated by diff number of columns

    extTriggerSignals = [];
    switch nColumns
        case 3 % Cologne format
            r           = double(C{2});
            rSignals    = double(C{3});
        case {4,5}
            r           = double(C{3});
            rSignals    = ~cellfun(@isempty, C{4});
    end

    if nColumns == 5
        extTriggerSignals = ~cellfun(@isempty, C{5});
    end

    rTics           = double(C{1});
    tRespiration    = rTics*dtRespiration ...
        - log_files.relative_start_acquisition;

    rpulse          = find(rSignals);

    nSamples        = numel(C{1});
    racq_codes       = zeros(nSamples,1);

    if ~isempty(rpulse)
        racq_codes(rpulse) = racq_codes(rpulse) + 4;
        rpulse = tRespiration(rpulse);
    end

    acqpulse          = find(extTriggerSignals);

    if ~isempty(acqpulse)
        racq_codes(acqpulse) = racq_codes(acqpulse) + 8;
    end
end

%% Read cardiac file
if hasCardiacFile

    C = tapas_physio_read_columnar_textfiles(log_files.cardiac);

    nColumns = numel(C);

    % different file formats indicated by diff number of columns

    extTriggerSignals = [];
    switch nColumns
        case 3 % Cologne format
            c           = double(C{2});
            cSignals    = double(C{3});
            cTics           = double(C{1});

        case {4,5}
            [cTics, c, cSignals, extTriggerSignals, stringChannels, verbose] = ...
                tapas_physio_split_data_per_channel_siemens_tics(C, ecgChannel, verbose);
    end

    tCardiac        = cTics*dtCardiac ...
        - log_files.relative_start_acquisition;

    nSamples         = numel(c);
    cacq_codes       = zeros(nSamples,1);

    cpulse          = find(cSignals);

    if ~isempty(cpulse)
        isOxy = any(strfind(upper(log_files.cardiac), 'PULS')); % different codes for PPU

        cacq_codes(cpulse) = cacq_codes(cpulse) + 1 + isOxy; %+1 for ECG, +2 for PULS
        cpulse = tCardiac(cpulse);
    end

    acqpulse          = find(extTriggerSignals);

    if ~isempty(acqpulse)
        cacq_codes(acqpulse) = cacq_codes(acqpulse) + 8;
    end
end

%% plot raw data so far
if DEBUG
    fh = tapas_physio_plot_raw_physdata_siemens_tics(tCardiac, c, tRespiration, r, ...
        hasCardiacFile, hasRespirationFile, cpulse, rpulse, cacq_codes, ...
        racq_codes);
    verbose.fig_handles(end+1) = fh;
end

%% If only one file exists, take t and acq_codes from that accordingly
hasOnlyCardiacFile = hasCardiacFile && ~hasRespirationFile;
if hasOnlyCardiacFile
    t = tCardiac;
    acq_codes = cacq_codes;
end

hasOnlyRespirationFile = ~hasCardiacFile && hasRespirationFile;
if hasOnlyRespirationFile
    t = tRespiration;
    acq_codes = racq_codes;
end

%% Merge acquisition codes, if both files exist, but on same sampling grid
if hasCardiacFile && hasRespirationFile

    haveSameSampling = isequal(tRespiration, tCardiac);
    if haveSameSampling
        acq_codes = cacq_codes + racq_codes;
    else
        %% interpolate to greater precision, if both files exist and
        % 2 different sampling rates are given
        %interpolate acq_codes and trace with lower sampling rate to higher
        %rate

        dtCardiac = tCardiac(2)-tCardiac(1);
        dtRespiration = tRespiration(2) - tRespiration(1);

        isHigherSamplingCardiac = dtCardiac < dtRespiration;
        if isHigherSamplingCardiac
            t = tCardiac;
            rInterp = interp1(tRespiration, r, t);
            racq_codesInterp = interp1(tRespiration, racq_codes, t, 'nearest');
            acq_codes = cacq_codes + racq_codesInterp;

            if DEBUG
                fh = plot_interpolation(tRespiration, r, t, rInterp, ...
                    {'respiratory', 'cardiac'});
                verbose.fig_handles(end+1) = fh;
            end
            r = rInterp;

        else
            t = tRespiration;
            cInterp = interp1(tCardiac, c, t);
            cacq_codesInterp = interp1(tCardiac, cacq_codes, t, 'nearest');
            acq_codes = racq_codes + cacq_codesInterp;

            if DEBUG
                fh = plot_interpolation(tCardiac, c, t, cInterp, ...
                    {'cardiac', 'respiratory'});
                verbose.fig_handles(end+1) = fh;
            end
            c = cInterp;

        end
    end
end
end

%% Local function to plot interpolation result
function fh = plot_interpolation(tOrig, yOrig, tInterp, yInterp, ...
    stringOrigInterp)
fh = tapas_physio_get_default_fig_params();
stringTitle = sprintf('Read-In: Interpolation of %s signal', stringOrigInterp{1});
set(fh, 'Name', stringTitle);
plot(tOrig, yOrig, 'go--');  hold all;
plot(tInterp, yInterp,'r.');
legend({
    sprintf('after interpolation to %s timing', ...
    stringOrigInterp{1}), ...
    sprintf('original %s time series', stringOrigInterp{2}) });
title(stringTitle);
xlabel('time (seconds');
end

author: translationalneuromodeling

repo: tapas


text:
function f = polyterm_value_3d ( n, x )

%*****************************************************************************80
%
%% POLYTERM_VALUE_3D evaluates a single polynomial term in 3D.
%
%  Discussion:
%
%    The polynomial term has the form:
%
%      F(X) = X(1)^E(1) * X(2)^E(2) * X(3)^E(3)
%
%    The exponents E(1:3) are set by calling POLYTERM_EXPONENT_SET.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input, real X(3,N), the points where the polynomial term
%    is to be evaluated.
%
%    Output, real F(N), the value of the polynomial term.
%
  e = [];
  e = polyterm_exponent ( 'GET', e );

  f(1:n,1) = 1.0;

  for i = 1 : 3

    if ( e(i) ~= 0 )
      f(1:n,1) = f(1:n,1) .* ( x(i,1:n).^e(i) )';
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function mObject3 = timeFrequencySTRAIGHTmorphing(mObject1,mObject2,mRate,mixMethod);
%   Morphing based on STRAIGHT parameters
%   mObject3 = timeFrequencySTRAIGHTmorphing(mObject1,mObject2,mRate,mixMethod);

%   Designed and coded by Hideki Kawahara
%   28/Feb./2005
%   Copyright(c) 2005, Hideki Kawahara
%   14/March/2005 bug fix on sampling frequency
%   01/Oct./2005 bug fix on similarity check
%   04/Oct./2005 partial morphing extension
%   18/Oct./2005 direct differential manipulation and API cange
%   29/Jan./2006 bug fix on boundary conditions

switch nargin
    case 0
        mObject3.morphingObject = createMobject;
        mixRate.F0 = 0;
        mixRate.spectrum = 0;
        mixRate.aperiodicity = 0;
        mixRate.coordinate = 0;
        mObject3.mixRate = mixRate;
        mObject3.mixMethods = {'linear','log','differentialLogarithm'};
        return
end;
mObject3 = checkForSimilarity(mObject1,mObject2);
mixRate = checkForMorphingConditions(mRate);
mObject1 = checkForIntegrity(mObject1);
mObject2 = checkForIntegrity(mObject2);
fs = mObject1.samplingFrequency;
if length(mObject3) ==0;return;end;
dtFrame = mObject1.frameUpdateInterval;
endLocation1 = (length(mObject1.F0)-1)*dtFrame; % in ms
endLocation2 = (length(mObject2.F0)-1)*dtFrame; % in ms
timeAnchor1 = [0;mObject1.anchorTimeLocation;endLocation1];
timeAnchor2 = [0;mObject2.anchorTimeLocation;endLocation2];
locationOn1 = (0:length(mObject1.F0)-1)*dtFrame;
locationOn2 = (0:length(mObject2.F0)-1)*dtFrame;
mapFrom1to2 = interp1(timeAnchor1,timeAnchor2,locationOn1);
[nr1,nc1] = size(mObject1.spectrogram);
[nr2,nc2] = size(mObject2.spectrogram);

%---- initialize frequency mapping function
fmapFrom1to2OnTime1 = generateFrequencyMap(mObject1,mObject2);

%---- mixing on mObject1's time axis
nAxis1 = length(locationOn1);
nAxis2 = length(locationOn2);
morphedF0 = zeros(nAxis1,1);
morphedAp = zeros(nr1,nAxis1);
morphedSgram = zeros(nr1,nAxis1);
weightSumF0 = zeros(nAxis1,1);
for ii=1:nAxis1
    mappedIndexOn2 = mapFrom1to2(ii)/dtFrame+1;
    iFloor = floor(mappedIndexOn2);
    iFraction = mappedIndexOn2-iFloor;
    fIndex = floor(fmapFrom1to2OnTime1(:,ii)/fs*2*(nr1-1))+1;
    dAp = iFraction*(mObject2.aperiodicityIndex(:,min(iFloor+1,nAxis2))-mObject2.aperiodicityIndex(:,min(iFloor,nAxis2)));
    ap2on2faxis = mObject2.aperiodicityIndex(:,min(iFloor,nAxis2))+dAp;
    ap2on1faxis = ap2on2faxis(fIndex);
    morphedAp(:,ii) = (1-mixRate.aperiodicity)*mObject1.aperiodicityIndex(:,ii)+mixRate.aperiodicity*ap2on1faxis; %04/Oct/2005 HK
    switch mixMethod
        case 'linear'
            dSgram = iFraction*(mObject2.spectrogram(:,min(iFloor+1,nAxis2))-mObject2.spectrogram(:,min(iFloor,nAxis2)));
            sgram2on2faxis = mObject2.spectrogram(:,min(iFloor,nAxis2))+dSgram;
            sgram2on1faxis = sgram2on2faxis(fIndex);
            morphedSgram(:,ii) = (1-mixRate.spectrum)*mObject1.spectrogram(:,ii)+mixRate.spectrum*sgram2on1faxis;
        case 'log'
            dSgram = iFraction*(log(mObject2.spectrogram(:,min(iFloor+1,nAxis2)))-log(mObject2.spectrogram(:,min(iFloor,nAxis2))));
            sgram2on2faxis = log(mObject2.spectrogram(:,min(iFloor,nAxis2)))+dSgram;
            sgram2on1faxis = sgram2on2faxis(fIndex);
            tmp = (1-mixRate.spectrum)*log(mObject1.spectrogram(:,ii))+mixRate.spectrum*sgram2on1faxis;
            morphedSgram(:,ii) = exp(tmp);
        case 'differentialLogarithm'
            dSgram = iFraction*(mObject2.spectrogram(:,min(iFloor+1,nAxis2))-mObject2.spectrogram(:,min(iFloor,nAxis2)));
            sgram2on2faxis = mObject2.spectrogram(:,min(iFloor,nAxis2))+dSgram;
            sgram2on1faxis = sgram2on2faxis(fIndex);
            tmp = (1-mixRate.spectrum)*log(mObject1.spectrogram(:,ii))+mixRate.spectrum*sgram2on1faxis;
            morphedSgram(:,ii) = exp(tmp);
    end;
    if mObject1.F0(ii)>0
        morphedF0(ii) = (1-mixRate.F0)*log(mObject1.F0(ii));
        weightSumF0(ii) = (1-mixRate.F0);
    end;
    if (mObject2.F0(iFloor)>0) & (mObject2.F0(min(iFloor+1,nAxis2))>0)
        dF0 = iFraction*(log(mObject2.F0(min(iFloor+1,nAxis2)))-log(mObject2.F0(min(iFloor,nAxis2))));
        morphedF0(ii) = mixRate.F0*(log(mObject2.F0(min(iFloor,nAxis2)))+dF0)+morphedF0(ii);
        weightSumF0(ii) = weightSumF0(ii)+mixRate.F0;
    end;
end;
morphedF0(weightSumF0>0) = exp(morphedF0(weightSumF0>0)./weightSumF0(weightSumF0>0));

%----- mapping back onto morphed time axis
timeAnchorMorph = (1-mixRate.coordinate)*timeAnchor1 + mixRate.coordinate*timeAnchor2;
locationOnMorph = (0:(timeAnchorMorph(end)/dtFrame))*dtFrame;
mapFormMorphTo1 = interp1(timeAnchorMorph,timeAnchor1,locationOnMorph);
nAxisMorph = length(locationOnMorph);
morphedApOnMorph = zeros(nr1,nAxisMorph);
morphedSgramOnMorph = zeros(nr1,nAxisMorph);
morphedF0onMorph = zeros(nAxisMorph,1);
%----- set place holders
mObject3.samplingFrequency = fs;
mObject3.F0 = morphedF0onMorph;
mObject3.aperiodicityIndex = morphedApOnMorph;
mObject3.spectrogram = morphedSgramOnMorph;
mObject3.anchorTimeLocation = timeAnchorMorph(2:end-1);
mObject3.anchorFrequency = (1-mixRate.coordinate)*mObject1.anchorFrequency+mixRate.coordinate*mObject2.anchorFrequency;
%------ nitialize frequency mapping function
fmapFromMorphto1OnTimeMorph = generateFrequencyMap(mObject3,mObject1);
for ii=1:nAxisMorph
    mappedIndexOn1 = mapFormMorphTo1(ii)/dtFrame+1;
    iFloor = floor(mappedIndexOn1);
    iFraction = mappedIndexOn1-iFloor;
    fIndex = floor(fmapFromMorphto1OnTimeMorph(:,ii)/fs*2*(nr1-1))+1;
    morphedApOnMorph(:,ii) = morphedAp(fIndex,iFloor) ...
        +iFraction*(morphedAp(fIndex,min(iFloor+1,nAxis1))-morphedAp(fIndex,iFloor));
    morphedSgramOnMorph(:,ii) = morphedSgram(fIndex,iFloor) ...
        +iFraction*(morphedSgram(fIndex,min(iFloor+1,nAxis1))-morphedSgram(fIndex,iFloor));
    if (morphedF0(iFloor)>0) & (morphedF0(min(iFloor+1,nAxis1))>0)
        dF0 = iFraction*(morphedF0(min(iFloor+1,nAxis1))-morphedF0(iFloor));
        morphedF0onMorph(ii) = morphedF0(iFloor)+dF0;
    end;
end;
mObject3.F0 = morphedF0onMorph;
mObject3.aperiodicityIndex = morphedApOnMorph;
mObject3.spectrogram = morphedSgramOnMorph;
mObject3.anchorTimeLocation = timeAnchorMorph(2:end-1);
%mObject3.anchorFrequency = (1-mRate)*mObject1.anchorFrequency+mRate*mObject2.anchorFrequency;
%mObject3 = fmapFromMorphto1OnTimeMorph; % This line is a dummy.
return;

%%% ------ Internal function to check for object's similarity
function mObject3 = checkForSimilarity(mObject1,mObject2)
mObject3 = [];
if mObject1.samplingFrequency ~= mObject2.samplingFrequency;mObject3 = [];return;end;
if mObject1.frameUpdateInterval ~= mObject2.frameUpdateInterval;mObject3 = [];return;end;
if length(mObject1.anchorTimeLocation) ~= length(mObject2.anchorTimeLocation);mObject3 = [];return;end;
nAnchor = length(mObject1.anchorTimeLocation);
for ii=1:nAnchor % check for similarity of anchor structure
    frequencyAnchor1 = mObject1.anchorFrequency(ii,:)';% 01/Oct./2005 by HK
    frequencyAnchor2 = mObject2.anchorFrequency(ii,:)';% 01/Oct./2005 by HK
    if (sum(frequencyAnchor1>0) ~= sum(frequencyAnchor2>0)) | ...
            (sum(frequencyAnchor1<0) ~= sum(frequencyAnchor2<0))
        display('Warning!! Object structures are inconsistent!'); % 01/Oct./2005 by HK
        return;
    end;
end;
mObject3 = createMobject;
m0bject3.samplingFrequency = mObject1.samplingFrequency;
m0bject3.frameUpdateInterval = mObject1.frameUpdateInterval;
return;

%%%--------
function mixRate = checkForMorphingConditions(mRate);
%   04/Oct./2005 added by HK

if ~isstruct(mRate)
    mixRate.F0 = mRate;
    mixRate.spectrum = mRate;
    mixRate.aperiodicity = mRate;
    mixRate.coordinate = mRate;
    return;
end;
mixRate.F0 = mRate.F0;
mixRate.spectrum = mRate.spectrum;
mixRate.aperiodicity = mRate.aperiodicity;
mixRate.coordinate = mRate.coordinate;
return;

%%%--------
function fmapFrom1to2OnTime1 = generateFrequencyMap(mObject1,mObject2);

dtFrame = mObject1.frameUpdateInterval;
endLocation1 = (length(mObject1.F0)-1)*dtFrame; % in ms
timeAnchor1 = [0;mObject1.anchorTimeLocation;endLocation1];
locationOn1 = (0:length(mObject1.F0)-1)*dtFrame;
fs = mObject1.samplingFrequency;
[nr1,nc1] = size(mObject1.spectrogram);
nAnchor = length(mObject1.anchorTimeLocation);
fmapFrom1to2 = zeros(nr1,nAnchor);
frequencyAxis = (0:nr1-1)'/(nr1-1)*fs/2;
numberOfFrequencyAnchors = zeros(nAnchor,1);
for ii=1:nAnchor
    frequencyAnchor1 = mObject1.anchorFrequency(ii,:)';
    frequencyAnchor1 = [0;frequencyAnchor1(frequencyAnchor1>0);fs/2];
    numberOfFrequencyAnchors(ii) = length(frequencyAnchor1(frequencyAnchor1>0));
    frequencyAnchor2 = mObject2.anchorFrequency(ii,:)';
    frequencyAnchor2 = [0;frequencyAnchor2(frequencyAnchor2>0);fs/2];
    fmapFrom1to2(:,ii) = interp1(frequencyAnchor1,frequencyAnchor2,frequencyAxis);
end;
for ii=1:nAnchor
    if numberOfFrequencyAnchors(ii) == 1
        if numberOfFrequencyAnchors(min(ii+1,nAnchor)) > 1
            fmapFrom1to2(:,ii) = fmapFrom1to2(:,min(ii+1,nAnchor));
        elseif numberOfFrequencyAnchors(max(ii-1,1)) > 1
            fmapFrom1to2(:,ii) = fmapFrom1to2(:,max(ii-1,1));
        end;
    end;
end;
fmapFrom1to2 = [fmapFrom1to2(:,1) fmapFrom1to2 fmapFrom1to2(:,nAnchor)];
fmapFrom1to2OnTime1 = interp1(timeAnchor1,fmapFrom1to2',locationOn1)';
return;

%%%-------
function cleanedUpObject = checkForIntegrity(inputObject);

maximumIndex = max([length(inputObject.F0), ...
    size(inputObject.spectrogram,2) ...
    size(inputObject.aperiodicityIndex,2)]);
if length(inputObject.F0) < maximumIndex
    inputObject.F0 = [inputObject.F0(:);inputObject.F0(end)*ones(maximumIndex - length(inputObject.F0),1)];
end;
if size(inputObject.spectrogram,2) < maximumIndex
    numberOfFillIn = maximumIndex-size(inputObject.spectrogram,2);
    inputObject.spectrogram = [inputObject.spectrogram inputObject.spectrogram(:,end)*ones(1,numberOfFillIn)];
end;
if size(inputObject.aperiodicityIndex,2) < maximumIndex
    numberOfFillIn = maximumIndex-size(inputObject.aperiodicityIndex,2);
    inputObject.aperiodicityIndex = [inputObject.aperiodicityIndex inputObject.aperiodicityIndex(:,end)*ones(1,numberOfFillIn)];
end;
cleanedUpObject = inputObject;

author: HidekiKawahara


text:
function compareSampleTraj(rxnNames, samples, models, nBins)
% Compares flux histograms for two or more samples
% for one or more reactions
%
% USAGE:
%
%    compareSampleTraj(rxnNames, samples, models, nBins)
%
% INPUTS:
%    rxnNames:    List of reaction names to compare
%    samples:     Samples to compare
%    models:      Cell array containing COBRA model structures
%
% OPTIONAL INPUTS:
%    nBins:       Number of bins (Default = `nSamples` / 25)

if (nargin < 4)
  [tmp,nSamples] = size(samples{1});
  bins = round(nSamples/25);
end

if (~iscell(rxnNames))
  rxnNameList{1} = rxnNames;
else
  rxnNameList = rxnNames;
end

nRxns = length(rxnNameList);
nX = ceil(sqrt(nRxns));

nY = ceil(nRxns/nX);

for j = 1:nRxns

  clear counts;
  currLB = 1e6;
  currUB = -1e6;
  for i = 1:length(models)
    id = findRxnIDs(models{i}, rxnNameList{j});
    if (isempty(id))
      id = findRxnIDs(models{i}, [rxnNameList{j} '_r']);
      if (isempty(id))
        error('Reaction does not exist');
      end
    end
    currLB = min(currLB, min(samples{i}(id, :)'));
    currUB = max(currUB, max(samples{i}(id, :)'));
  end

  bins = linspace(currLB, currUB, nBins);

  for i = 1:length(models)
    sampleSign = 1;
    id = findRxnIDs(models{i}, rxnNameList{j});
    if (isempty(id))
      id = findRxnIDs(models{i}, [rxnNameList{j} '_r']);
      sampleSign = -1;
    end
    traj(:,i) = sampleSign * samples{i}(id, :)';
  end

  subplot(nY, nX, j);
  plot(traj, '-');
  axis([0 length(traj) currLB currUB]);
  %text((currUB+currLB)/2,max(max(counts))+20,rxnNameList{j});
end

author: opencobra

repo: cobratoolbox


text:
function [cc] = ml2cc(ml)
% Convert volume from milliliters to cubic centimeters.
% Chad Greene 2012
cc = ml;

author: Sable

repo: mcbench-benchmarks


text:
function [fMaskList_matchZ,fMask_matchZ_ZList] = ScreenAllfMasksWithzMasks(hfig)

% fMASKs.MaskDatabase = MaskDatabase;
%     fMASKs.mask_clusID = mask_clusID;
%     fMASKs.i_fish = i_fish;
%     fMASKs.note = 'Auto_defS_woA_Master0.5';

data_masterdir = GetCurrentDataDir();

range_fish =  1:4; % range_fish = GetFishRange();

%% custom params here:
MASKs = getappdata(hfig,'MASKs');
nzMasks = size(MASKs.MaskDatabase,2);
thres_ratio = 0.33;
thres_maxsize = 6*10^5; % ~ size of subpallium

%%
fMaskList_matchZ = cell(1,18);
fMask_matchZ_ZList = cell(1,18);
for i_fishnum = 1:length(range_fish),
    i_fish = range_fish(i_fishnum);
    disp(i_fish);

    %% Load all fMasks from this fish
    name = ['fMASKs_',num2str(i_fish),'.mat'];
    load(fullfile(data_masterdir,'fMASKs',name),'fMASKs');
    nfMasks = size(fMASKs.MaskDatabase,2);
    mask_clusID = fMASKs.mask_clusID;
    %%
    Grid_ratio = zeros(nzMasks,nfMasks);
    for i_zMask = 1:nzMasks,
       zMsk = MASKs.MaskDatabase(:,i_zMask);
       if length(find(zMsk))<thres_maxsize,
           % compare this Z-Brain mask with each fMask
           for i_fMask = 1:nfMasks,
               [tf,ix] = ismember(i_fMask,mask_clusID);
               if tf,
                   fMsk = fMASKs.MaskDatabase(:,ix);
                   if length(find(fMsk))<thres_maxsize,
                       IX_intersect = find(zMsk.*fMsk);
                       r1 = length(IX_intersect)/length(find(zMsk));
                       r2 = length(IX_intersect)/length(find(fMsk));
                       Grid_ratio(i_zMask,i_fMask) = mean([r1,r2]);%max(r1,r2);
                   end
               end

           end
       end
    end

    % Pool from Grid
    [A,IX_zMask] = max(Grid_ratio,[],1);
    IX = find(A>thres_ratio);
    fMaskList_matchZ{i_fish} = IX;
    fMask_matchZ_ZList{i_fish} = IX_zMask(IX);

    %% save
    save(fullfile(data_masterdir,'fMASKs','fMASK_matchZ_Output.mat'),...
        'fMaskList_matchZ','fMask_matchZ_ZList');
    timestamp  = datestr(now,'mmddyy_HHMM');
    backupname = ['fMASK_matchZ_Output_',timestamp,'.mat'];
    save(fullfile(data_masterdir,'fMASKs',backupname),...
        'fMaskList_matchZ','fMask_matchZ_ZList');
end

end

author: xiuyechen

repo: FishExplorer

path: github-repos/MATLAB/xiuyechen-FishExplorer/FishExplorer-c61392cf0835480d64fc03c15f1992935fdc7106/misc scripts/screen clusters with ZBrain masks/ScreenAllfMasksWithzMasks.m


text:
function [ know, x ] = p17_sol ( n )

%*****************************************************************************80
%
%% P17_SOL returns the solution for problem 17.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the problem.  This value
%    is only needed for those problems with variable N.
%
%    Output, integer KNOW.
%    If KNOW is 0, then the solution is not known.
%    If KNOW is positive, then the solution is known, and is returned in X.
%
%    Output, real X(N), the solution, if known.
%
  know = 1;

  x = [ 1.0, 1.0, 1.0, 1.0 ]';

  return
end


text:
function [errY,fracerrY,errR,fracerrR,response] = hmmerror(X,T,hmm,...
    Gamma,test,residuals,actstates)
%
% Computes the error and explained variance for the given data points
%
% INPUT
% X         Observations
% T         Number of time points for each time series
% hmm       hmm data structure
% Gamma     probability of current state cond. on data -
%           inference is run for time points with Gamma=NaN,
% test      sum(T) x 1 vector indicating in which time points
%               the error is to be computed;
%           common use will be to set test(t) = 1 when Gamma is NaN at this point
% residuals     in case we train on residuals, the value of those.
% actstates     Kx1 vector indicating which states were effectively
%               used in the training
%
% OUTPUT
% errY             mean quadratic error
% fracerrY       fractional quadratic error
% errR             mean quadratic error of the residuals
% fracerrR       fractional quadratic error of the residuals
% response      mean of the predictive response distribution for test(t)=1
%
% Author: Diego Vidaurre, OHBA, University of Oxford

N = length(T);
train = hmm.train;

setstateoptions;

if nargin<5
    test = ones(sum(T),1);
end
if nargin<6 || isempty(residuals)
    residuals = ...
        getresiduals(X,T,train.S,train.maxorder,train.order,train.orderoffset,...
        train.timelag,train.exptimelag,train.zeromean);
end
if nargin<7
    actstates = ones(hmm.K,1);
end

Y = [];
te = [];
for j = 1:N
    t0 = sum(T(1:j-1));
    Y = [Y; X(t0+1+train.maxorder:t0+T(j),:)];
    te = [te; test(t0+1+train.maxorder:t0+T(j))];
end

[response,responseR] = hmmpred(X,T,hmm,Gamma,residuals,actstates);
response = response(te==1,:); responseR = responseR(te==1,:);
errY = mean((Y(te==1,regressed) - response(:,regressed)).^2);
fracerrY = errY ./ mean( (Y(te==1,regressed) - ...
    repmat(mean(Y(te==1,regressed)),sum(te),1) ).^2 );
errR = mean((residuals(te==1,regressed) - responseR(:,regressed)).^2);
fracerrR = errR ./ mean( (residuals(te==1,regressed) - ...
    repmat(mean(residuals(te==1,regressed)),sum(te),1) ).^2 );

author: OHBA-analysis

repo: HMM-MAR


text:
function SimilarityMatrix = similarityEuclidean(x)
    %spatial distance measure
    Distances = sqrt(pdist(x)+3e-6).^-1; % 0.05 best so far

    SimilarityMatrix = squareform(Distances) + eye(size(x, 1));

end

author: TadasBaltrusaitis

repo: OpenFace


text:
% std_pacplot() - Commandline function to plot cluster PAC
%                (phase-amplitude coupling).
%
% Usage:
%   >> [STUDY] = std_pacplot(STUDY, ALLEEG, key1, val1, key2, val2);
%   >> [STUDY itcdata itctimes itcfreqs pgroup pcond pinter] = ...
%                std_pacplot(STUDY, ALLEEG ...);
% Inputs:
%   STUDY      - EEGLAB STUDY set comprising some or all of the EEG datasets in ALLEEG.
%   ALLEEG     - global EEGLAB vector of EEG structures for the datasets in the STUDY.
%                Note: ALLEEG for a STUDY set is typically created using load_ALLEEG().
%
% Additional help:
% Inputs and output of this function are strictly identical to the std_erspplot().
% See the help message of this function for more information.
%
% See also: std_erspplot(), pop_clustedit(), pop_preclust()
%
% Authors: Arnaud Delorme, CERCO, July, 2009-

% Copyright (C) Arnaud Delorme, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [STUDY, allpac, alltimes, allfreqs, pgroup, pcond, pinter] = std_pacplot(STUDY, ALLEEG, varargin)

if nargin < 2
    help std_pacplot;
    return;
end;

[STUDY allpac alltimes allfreqs pgroup pcond pinter ] = std_erspplot(STUDY, ALLEEG, 'datatype', 'pac', varargin{:});

author: PatternRecognition

repo: OpenBMI


text:
function f = p08_f ( n, x )

%*****************************************************************************80
%
%% P08_F evaluates the objective function for problem 8.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the argument of the objective function.
%
%    Output, real F, the value of the objective function.
%
  ap = 0.00001;

  t1 = - 0.25 + sum ( x(1:n).^2 );

  t2 = sum ( ( x(1:n) - 1.0 ).^2 );

  f = ap * t2 + t1 * t1;

  return
end


text:
function I = addgrid(I,gsz)

gridpts = 1:gsz(2):size(I,2);
grid_x(:,1) = gridpts;
grid_x(:,3) = gridpts;
grid_x(:,2) = 1;
grid_x(:,4) = size(I,1);
gridpts = 1:gsz(1):size(I,1);
grid_y(:,2) = gridpts;
grid_y(:,4) = gridpts;
grid_y(:,1) = 1;
grid_y(:,3) = size(I,2);

I = insertShape(I,'Line',[grid_x;grid_y]);

author: kristinbranson

repo: JAABA


text:
% The COBRAToolbox: testSearchModel.m
%
% Purpose:
%     - tests the searchModel function.
%
% Authors:
%     - Thomas Pfau - June 2018
%

% initialize the test
fileDir = fileparts(which('testSearchModel'));
% save the current path
currentDir = cd(fileDir);

% load model
model = getDistributedModel('ecoli_core_model.mat');

%Find exactly 'atp[c]'
result = searchModel(model,'atp[c]','similarity',1,'printLevel',0);
assert(isempty(setxor(fieldnames(result),{'mets'}))) % Nothing but mets found.
assert(isempty(setxor({result.mets(:).id},{'atp[c]'}))) % Nothing but atp[c] found.
assert(strcmp(result.mets(1).matches(1).source,'mets')) %Only found in the mets field.

%Search for Glycolysis in the model (will find all subSystems which are
%Glycolysis/Gluconeogenesis
diary('Output.txt');
result = searchModel(model,'Glycolysis','printLevel',1);
diary('off');
glycolysis = findRxnsFromSubSystem(model,'Glycolysis/Gluconeogenesis');
%TODO repair
%assert(isempty(setxor({result.rxns(:).id},glycolysis))); %Its exactly those reactions.
assert(isempty(setxor(fieldnames(result),{'rxns'}))) % There is nothing else.
printOut = fileread('Output.txt');
%All the reactions were mentioned in the printOut
%assert(all(cellfun(@(x) ~isempty(strfind(printOut,['ID: ' x])),glycolysis)));

%Nothing in the model, so the returned struct is empty.
result = searchModel(model,'WeLookForSomethingOdd','printLevel',1);
assert(isempty(fieldnames(result)));
%We add a field that is looking like an annotation field.
model.geneisSomething = repmat({''},size(model.genes));
model.geneisSomething{2} = 'WeLookForSomethingOdderThanThis;WhateverweFindHereIsThere;Maybe';
result = searchModel(model,'WeLookForSomethingOdd','printLevel',1);
%Now, result should have the second gene.
assert(isempty(setxor(fieldnames(result),{'genes'}))) %There are genes, and only genes found now.
assert(isempty(setxor({result.genes(:).id},model.genes(2)))); % and gene2 was found.
assert(strcmp(result.genes(1).matches(1).source,'geneisSomething')); % and gene2 was found.

%Cleanup
delete('Output.txt')

%Return to original directory
cd(currentDir);


text:
function elist=surfpart(f,loopedge)
%
% elist=surfpart(f,loopedge)
%
% partition a triangular surface using a closed loop defined by existing edges
%
% author: Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
% date: 2012/02/09
%
% input:
%      f: input, surface face element list, dimension (be,3)
%      loopedge: a 2-column array, specifying a closed loop in CCW order
%
% output:
%      elist: list of triangles that is enclosed by the loop
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

elist=[];
if(isempty(f) || isempty(loopedge))
    return;
end

if(size(f,2)==3)
    edges=[f(:,[1,2]);
           f(:,[2,3]);
           f(:,[3,1])];             % create all the edges
elseif(size(f,2)==4)
    edges=[f(:,[1,2]);
           f(:,[2,3]);
           f(:,[3,4]);
           f(:,[4,1])];             % create all the edges
else
    error('surfpart only supports triangular and quadrilateral elements');
end

[elist,front]=advancefront(edges,loopedge);
while(~isempty(front))
	[elist0,front0]=advancefront(edges,front);
	elist=unique([elist;elist0]);
	front=front0;
end

repo: REST


text:
function ecef2mapTrafo(obj, mstruct, varargin)
% ECEF2MAPTRAFO Coordinate transformation from ecef to map coordinates.
%
% Example: definition of mstruct for UTM33N
% mstruct       = defaultm('utm');
% mstruct.zone  = '33n';
% mstruct.geoid = referenceEllipsoid('GRS 80');
% mstruct       = defaultm(mstruct);

% Input parsing ----------------------------------------------------------------

p = inputParser;
p.addRequired( 'mstruct');
p.parse(mstruct, varargin{:});
p = p.Results;
% Clear required inputs to avoid confusion
clear mstruct

% Start ------------------------------------------------------------------------

procHierarchy = {'POINTCLOUD' 'ECEF2MAPTRAFO'};
msg('S', procHierarchy);
msg('I', procHierarchy, sprintf('Point cloud label = ''%s''', obj.label));

% Conversion to map coordinates ------------------------------------------------

[lat, lon, hEll] = ecef2geodetic(obj.X(:,1), obj.X(:,2), obj.X(:,3), p.mstruct.geoid); % lat, lon in radian!

[xm, ym] = mfwdtran(p.mstruct, lat*180/pi, lon*180/pi, hEll); % lat, lon in degrees!

% Update coordinates -----------------------------------------------------------

obj.X = [xm ym hEll];
obj.info;

% End --------------------------------------------------------------------------

msg('E', procHierarchy);

end

author: pglira


text:
% DISPLAY_FILTER_BANK_2D Display all the fine-resolution filters of the filter
%    bank and returns the result to display
%
% Usage
%	big_img = DISPLAY_FILTER_BANK_2D(filters, renorm, n)
%
% Input
%    filters (cell of struct): filter banke from a wavelet factory (for
%    instance).
%    r (bool): renormalize the filters for display or not.
%    n (numeric): cropping size, only for 'fourier_multires' type of
%    function. Not required for spatial support filters.
%
% Output
%    big_img (numeric): displayed filter bank.
%
% Description
%    Display the real and imaginary part of a filter bank and returns it.
%
% See also
%    DISPLAY_FILTER_2D

function big_img = display_filter_bank_2d(filters, renorm, n)

if(isfield(filters,'phi'))
    if (~exist('n','var'))
		n = min(32, floor(min(filters.meta.size_in)/2 -1));
    end
    if (~exist('renorm','var'))
		renorm = 1;
	end

	margin = 4;
    margin_outer = 10;
	N = 2*n+1;
	N_margin = N + margin;
	J = filters.meta.J;
	L = filters.meta.L;
	big_img = ones(N_margin*J, (2*L+1)*N_margin + 2*margin_outer);

	% low pass : first on the left
	filt_for_disp = display_filter_2d(filters.phi.filter,'r', n);
	M = max(abs(filt_for_disp(:)));
	if (renorm == 0)
		M = 1;
	end
	big_img((1:N), (1:N)) = real(filt_for_disp)/M;

	% band pass : following filters
	for p = 1:numel(filters.psi.filter)
		filter = filters.psi.filter{p};
		filt_for_disp = display_filter_2d(filter,'r', n);
		M = max(abs(filt_for_disp(:)));
		if (renorm == 0)
			M = 1;
		end

		j = filters.psi.meta.j(p);
		theta = filters.psi.meta.theta(p);

		big_img((1:N)+ j*N_margin, ...
                (1:N) + (theta)*N_margin + margin_outer) = ...
                    real(filt_for_disp)/M;
		big_img((1:N)+ j*N_margin, ...
                   (1:N) + (theta+L)*N_margin + 2*margin_outer) = ...
		imag(filt_for_disp)/M;
	end

	imagesc(big_img);
elseif(isfield(filters,'g'))

    if (~exist('renorm','var'))
		renorm = 1;
    end

    N=filters.meta.size_filter(1);
	margin = 4;
	nFilter=length(filters.g.filter)+length(filters.h.filter);
	N_margin = N + margin;
	L = filters.meta.L;
	big_img = ones(N_margin*nFilter, 2*N_margin);

	% low pass : first on the left
	filt_for_disp = display_filter_2d(filters.h.filter,'r');
	M = max(abs(filt_for_disp(:)));
	if (renorm == 0)
		M = 1;
	end
	big_img((1:N), (1:N)) = real(filt_for_disp)/M;
	big_img((1:N) ,(1:N)+N_margin) = imag(filt_for_disp)/M;

	% band pass : following filters
	for p = 1:numel(filters.g.filter)
		filter = filters.g.filter{p};
		filt_for_disp = display_filter_2d(filter,'r');
		M = max(abs(filt_for_disp(:)));
		if (renorm == 0)
			M = 1;
		end

		%j = filters.g.meta.j(p);
		theta = filters.g.meta.theta(p);

		big_img(1+(1:N)+ p* N_margin, (1:N) + (1-1)*N_margin) = ...
		real(filt_for_disp)/M;
		big_img(1+(1:N)+ p *N_margin, (1:N) + (2-1)*N_margin) = ...
		imag(filt_for_disp)/M;
	end

	imagesc(big_img);

else
    error('The display of this filter bank has not been implemented');
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function yf = ymdhms_to_decimal ( y, m, d, h, n, s )

%*****************************************************************************80
%
%% YMDHMS_TO_DECIMAL converts a Y/M/D/H/Mn/S date to a Decimal Y.F date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, M, D, H, N, S, the YMDHMS date.
%
%    Output, real YF, the Decimal date.
%

%
%  How many days between January 1st and day D?
%
  day_min = 1;
  days = ymd_dif_common ( y, m, day_min, y, m, d );
%
%  How many days in this year total?
%
  day_max = year_length_common ( y );
%
%  The decimal part of the year is ( D + H/24 + N/24*60 + S/24*60*60 ) / DMAX.
%
  f =       s      / 60.0;
  f = ( f + n    ) / 60.0;
  f = ( f + h    ) / 24.0;
  f = ( f + days ) / day_max;

  yf = y + f;

  return
end


text:
%% Copyright (C) 2014-2016 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defop  Method   @@sym not {(@var{x})}
%% @defopx Operator @@sym {~@var{x}} {}
%% @defopx Operator @@sym {!@var{x}} {}
%% Logical "not" of a symbolic array.
%%
%% Example:
%% @example
%% @group
%% syms x y
%% eqn = 2*x == y
%%   @result{} eqn = (sym) 2⋅x = y
%% not(eqn)
%%   @result{} ans = (sym) 2⋅x ≠ y
%% @end group
%% @end example
%%
%% More briefly:
%% @example
%% @group
%% ~(x==y)
%%   @result{} ans = (sym) x ≠ y
%% !(x==y)
%%   @result{} ans = (sym) x ≠ y
%% @end group
%% @end example
%%
%% Applies to each entry of a matrix:
%% @example
%% @group
%% A = [x < y, 2*x + y >= 0, true]
%%   @result{} A = (sym) [x < y  2⋅x + y ≥ 0  True]  (1×3 matrix)
%% ~A
%%   @result{} ans = (sym) [x ≥ y  2⋅x + y < 0  False]  (1×3 matrix)
%% @end group
%% @end example
%%
%% @seealso{@@sym/eq, @@sym/ne, @@sym/logical, @@sym/isAlways}
%% @end defop

function r = not(x)

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

  r = elementwise_op ('Not', x);

end

%!shared t, f
%! t = sym(true);
%! f = sym(false);

%!test
%! % simple
%! assert (isequal( ~t, f))
%! assert (isequal( ~t, f))

%!test
%! % array
%! w = [t t f t];
%! z = [f f t f];
%! assert (isequal( ~w, z))

%!test
%! % number
%! assert (isequal( ~sym(5), f))
%! assert (isequal( ~sym(0), t))

%!test
%! % output is sym
%! syms x
%! e = ~(x == 4);
%! assert (isa (e, 'sym'))
%! assert (strncmp (sympy(e), 'Unequality', 10))

%!test
%! % output is sym even for scalar t/f (should match other bool fcns)
%! assert (isa (~t, 'sym'))

%!test
%! % symbol ineq
%! syms x
%! a = [t  f  x == 1  x ~= 2  x < 3   x <= 4  x > 5   x >= 6];
%! b = [f  t  x ~= 1  x == 2  x >= 3  x > 4   x <= 5  x < 6];
%! assert (isequal( ~a, b))

%!test
%! syms x
%! y = ~x;
%! s = disp(y, 'flat');
%! assert (strcmp (strtrim (s), '~x') || strcmpi (strtrim (s), 'Not(x)'))

%!error not (sym(1), 2)

repo: octsympy


text:
function view = blurMap(view,weightFlag,lambda)
% view = blurMap(view, [weightFlag], [lambda])
% Just like blurPh but operates on the parameter map rather than the ph
% See blurPh
% ARW 061004
if isempty(view)
  myErrorDlg('Volume window must be open to blur map.');
end
if ~strcmp(view.viewType,'Gray')
  myErrorDlg('Blur map only for gray views.');
end
if isempty(view.map)
   myErrorDlg('Maps not loaded.  Load them from the File menu');
end
if ~exist('weightFlag','var')
	weightFlag = 0;
end
if ~exist('lambda','var')
   lambda = 0.1;
end
curScan = getCurScan(view);
map = view.map{curScan}(:);
% Regularization
newMap = regularizeGray(map,view.nodes,view.edges,lambda);
view.map{curScan}(:) = newMap;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%SOM_DEMO4 Data analysis using the SOM.

% Contributed to SOM Toolbox 2.0, February 11th, 2000 by Juha Vesanto
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 1.0beta juuso 071197
% Version 2.0beta juuso 090200 070600

clf reset;
f0 = gcf;
echo on

clc
%    ==========================================================
%    SOM_DEMO4 - DATA ANALYSIS USING THE SOM
%    ==========================================================

%    In this demo, the IRIS data set is analysed using SOM. First, the
%    data is read from ascii file (please make sure they can be found
%    from the current path), normalized, and a map is
%    trained. Since the data also had labels, the map is labelled.

try
  sD = som_read_data('iris.data');
catch
  echo off

  warning('File ''iris.data'' not found. Using simulated data instead.')

  D = randn(50,4);
  D(:,1) = D(:,1)+5;     D(:,2) = D(:,2)+3.5;
  D(:,3) = D(:,3)/2+1.5; D(:,4) = D(:,4)/2+0.3;
  D2 = randn(100,4); D2(:,2) = sort(D2(:,2));
  D2(:,1) = D2(:,1)+6.5; D2(:,2) = D2(:,2)+2.8;
  D2(:,3) = D2(:,3)+5;   D2(:,4) = D2(:,4)/2+1.5;
  sD = som_data_struct([D; D2],'name','iris (simulated)',...
		       'comp_names',{'SepalL','SepalW','PetalL','PetalW'});
  sD = som_label(sD,'add',[1:50]','Setosa');
  sD = som_label(sD,'add',[51:100]','Versicolor');
  sD = som_label(sD,'add',[101:150]','Virginica');

  echo on
end

sD = som_normalize(sD,'var');
sM = som_make(sD);
sM = som_autolabel(sM,sD,'vote');

pause % Strike any key to visualize the map...

clc
%    VISUAL INSPECTION OF THE MAP
%    ============================

%    The first step in the analysis of the map is visual inspection.
%    Here is the U-matrix, component planes and labels (you may
%    need to enlarge the figure in order to make out the labels).

som_show(sM,'umat','all','comp',[1:4],'empty','Labels','norm','d');
som_show_add('label',sM.labels,'textsize',8,'textcolor','r','subplot',6);

%    From this first visualization, one can see that:
%     - there are essentially two clusters
%     - PetalL and PetalW are highly correlated
%     - SepalL is somewhat correlated to PetalL and PetalW
%     - one cluster corresponds to the Setosa species and exhibits
%       small petals and short but wide sepals
%     - the other cluster corresponds to Virginica and Versicolor
%       such that Versicolor has smaller leaves (both sepal and
%       petal) than Virginica
%     - inside both clusters, SepalL and SepalW are highly correlated

pause % Strike any key to continue...

%    Next, the projection of the data set is investigated. A
%    principle component projection is made for the data, and applied
%    to the map. The colormap is done by spreading a colormap on the
%    projection. Distance matrix information is extracted from the
%    U-matrix, and it is modified by knowledge of zero-hits
%    (interpolative) units. Finally, three visualizations are shown:
%    the color code, with clustering information and the number of
%    hits in each unit, the projection and the labels.

echo off

f1=figure;
[Pd,V,me,l] = pcaproj(sD,2); Pm = pcaproj(sM,V,me); % PC-projection
Code = som_colorcode(Pm); % color coding
hits = som_hits(sM,sD);  % hits
U = som_umat(sM); % U-matrix
Dm = U(1:2:size(U,1),1:2:size(U,2)); % distance matrix
Dm = 1-Dm(:)/max(Dm(:)); Dm(find(hits==0)) = 0; % clustering info

subplot(1,3,1)
som_cplane(sM,Code,Dm);
hold on
som_grid(sM,'Label',cellstr(int2str(hits)),...
	 'Line','none','Marker','none','Labelcolor','k');
hold off
title('Color code')

subplot(1,3,2)
som_grid(sM,'Coord',Pm,'MarkerColor',Code,'Linecolor','k');
hold on, plot(Pd(:,1),Pd(:,2),'k+'), hold off, axis tight, axis equal
title('PC projection')

subplot(1,3,3)
som_cplane(sM,'none')
hold on
som_grid(sM,'Label',sM.labels,'Labelsize',8,...
	 'Line','none','Marker','none','Labelcolor','r');
hold off
title('Labels')

echo on

%    From these figures one can see that:
%     - the projection confirms the existence of two different clusters
%     - interpolative units seem to divide the Virginica
%       flowers into two classes, the difference being in the size of
%       sepal leaves

pause % Strike any key to continue...

%    Finally, perhaps the most informative figure of all: simple
%    scatter plots and histograms of all variables. The species
%    information is coded as a fifth variable: 1 for Setosa, 2 for
%    Versicolor and 3 for Virginica. Original data points are in the
%    upper triangle, map prototype values on the lower triangle, and
%    histograms on the diagonal: black for the data set and red for
%    the map prototype values. The color coding of the data samples
%    has been copied from the map (from the BMU of each sample). Note
%    that the variable values have been denormalized.

echo off

% denormalize and add species information
names = sD.comp_names; names{end+1} = 'species';
D = som_denormalize(sD.data,sD); dlen = size(D,1);
s = zeros(dlen,1)+NaN; s(strcmp(sD.labels,'Setosa'))=1;
s(strcmp(sD.labels,'Versicolor'))=2; s(strcmp(sD.labels,'Virginica'))=3;
D = [D, s];
M = som_denormalize(sM.codebook,sM); munits = size(M,1);
s = zeros(munits,1)+NaN; s(strcmp(sM.labels,'Setosa'))=1;
s(strcmp(sM.labels,'Versicolor'))=2; s(strcmp(sM.labels,'Virginica'))=3;
M = [M, s];

f2=figure;

% color coding copied from the map
bmus = som_bmus(sM,sD); Code_data = Code(bmus,:);

k=1;
for i=1:5, for j=1:5,
    if i<j, i1=i; i2=j; else i1=j; i2=i; end
    subplot(5,5,k); cla
    if i<j,
      som_grid('rect',[dlen 1],'coord',D(:,[i1 i2]),...
	       'Line','none','MarkerColor',Code_data,'Markersize',2);
      title(sprintf('%s vs. %s',names{i1},names{i2}))
    elseif i>j,
      som_grid(sM,'coord',M(:,[i1 i2]),...
	       'markersize',2,'MarkerColor',Code);
      title(sprintf('%s vs. %s',names{i1},names{i2}))
    else
      if i1<5, b = 10; else b = 3; end
      [nd,x] = hist(D(:,i1),b); nd=nd/sum(nd);
      nm = hist(M(:,i1),x); nm = nm/sum(nm);
      h=bar(x,nd,0.8); set(h,'EdgeColor','none','FaceColor','k');
      hold on
      h=bar(x,nm,0.3); set(h,'EdgeColor','none','FaceColor','r');
      hold off
      title(names{i1})
    end
    k=k+1;
  end
end

echo on

%    This visualization shows quite a lot of information:
%    distributions of single and pairs of variables both in the data
%    and in the map. If the number of variables was even slightly
%    more, it would require a really big display to be convenient to
%    use.

%    From this visualization we can conform many of the earlier
%    conclusions, for example:
%     - there are two clusters: 'Setosa' (blue, dark green) and
%       'Virginica'/'Versicolor' (light green, yellow, reds)
%       (see almost any of the subplots)
%     - PetalL and PetalW have a high linear correlation (see
%       subplots 4,3 and 3,4)
%     - SepalL is correlated (at least in the bigger cluster) with
%       PetalL and PetalW (in subplots 1,3 1,4 3,1 and 4,1)
%     - SepalL and SepalW have a clear linear correlation, but it
%       is slightly different for the two main clusters (in
%       subplots 2,1 and 1,2)

pause % Strike any key to cluster the map...

close(f1), close(f2), figure(f0), clf

clc
%    CLUSTERING OF THE MAP
%    =====================

%    Visual inspection already hinted that there are at least two
%    clusters in the data, and that the properties of the clusters are
%    different from each other (esp. relation of SepalL and
%    SepalW). For further investigation, the map needs to be
%    partitioned.

%    Here, the KMEANS_CLUSTERS function is used to find an initial
%    partitioning. The plot shows the Davies-Boulding clustering
%    index, which is minimized with best clustering.

subplot(1,3,1)
[c,p,err,ind] = kmeans_clusters(sM, 7); % find at most 7 clusters
plot(1:length(ind),ind,'x-')
[dummy,i] = min(ind)
cl = p{i};

%    The Davies-Boulding index seems to indicate that there are
%    two clusters on the map. Here is the clustering info
%    calculated previously and the partitioning result:

subplot(1,3,2)
som_cplane(sM,Code,Dm)
subplot(1,3,3)
som_cplane(sM,cl)

%    You could use also function SOM_SELECT to manually make or modify
%    the partitioning.

%    After this, the analysis would proceed with summarization of the
%    results, and analysis of each cluster one at a time.
%    Unfortunately, you have to do that yourself. The SOM Toolbox does
%    not, yet, have functions for those purposes.

pause % Strike any key to continue...

clf
clc
%    MODELING
%    ========

%    One can also build models on top of the SOM. Typically, these
%    models are simple local or nearest-neighbor models.

%    Here, SOM is used for probability density estimation. Each map
%    prototype is the center of a gaussian kernel, the parameters
%    of which are estimated from the data. The gaussian mixture
%    model is estimated with function SOM_ESTIMATE_GMM and the
%    probabilities can be calculated with SOM_PROBABILITY_GMM.

[K,P] = som_estimate_gmm(sM,sD);
[pd,Pdm,pmd] = som_probability_gmm(sD,sM,K,P);

%    Here is the probability density function value for the first data
%    sample (x=sD.data(:,1)) in terms of each map unit (m):

som_cplane(sM,Pdm(:,1))
colorbar
title('p(x|m)')

pause % Strike any key to continue...

%    Here, SOM is used for classification. Although the SOM can be
%    used for classification as such, one has to remember that it does
%    not utilize class information at all, and thus its results are
%    inherently suboptimal. However, with small modifications, the
%    network can take the class into account. The function
%    SOM_SUPERVISED does this.

%    Learning vector quantization (LVQ) is an algorithm that is very
%    similar to the SOM in many aspects. However, it is specifically
%    designed for classification. In the SOM Toolbox, there are
%    functions LVQ1 and LVQ3 that implement two versions of this
%    algorithm.

%    Here, the function SOM_SUPERVISED is used to create a classifier
%    for IRIS data set:

sM = som_supervised(sD,'small');

som_show(sM,'umat','all');
som_show_add('label',sM.labels,'TextSize',8,'TextColor','r')

sD2 = som_label(sD,'clear','all');
sD2 = som_autolabel(sD2,sM);       % classification
ok = strcmp(sD2.labels,sD.labels); % errors
100*(1-sum(ok)/length(ok))         % error percentage (%)

echo off

author: ilarinieminen

repo: SOM-Toolbox


text:
function [sol, infos] = gsp_wiener_l2(G,y, A, At, psd, psd_noise, param)
%GSP_WIENER_l2 Solve wiener optimization problem with l2 fidelity term
%   Usage:  sol = gsp_wiener_l2(G, y, ffid, psd, psd_noise)
%           sol = gsp_wiener_l2(G, y, ffid, psd, psd_noise, param)
%           [sol, infos] = gsp_wiener_l2(...)
%
%   Input parameters:
%         G          : Graph (GSP structure)
%         x0         : Measurements (column vector)
%         A          : Operator (anonymous function)
%         At         : Adjoint operator (anonymous function)
%         psd        : PSD filter (anonymous function)
%         psd_noise  : PSD filter of the noise or single number
%         param      : Optional optimization parameters
%   Output parameters:
%         sol        : Solution
%         infos      : Convergence informations
%
%   This function solves the following wiener optimization problem:
%
%     .. argmin_x || A x - y ||_2^2 + || w(L) x ||_2^2
%
%     .. math:: arg\min_x \| A x - y \|_2^2 + \| w(L) x \|_2^2
%
%   Please refer to the reference for more information about this problem.
%   This function requires the UNLocBox to work.
%
%   Please refer to the function gsp_filter_analysis and solvep to know how
%   *param* can be set.
%
%   * *param.nu* : bound on the norm of the operator A (default: 1), i.e.
%
%     .. ` ||A x||^2 <= nu * ||x||^2
%
%     .. math::  \|A x\|^2 \leq \nu  \|x\|^2
%
%   References: perraudin2016stationary

% Author : Nathanael Perraudin
% Date: 6 January 2016

if nargin<7
    param = struct;
end

if ~isfield(param,'nu'), param.nu = 1; end
if ~isfield(param,'verbose'), param.verbose = 1; end

if isnumeric(psd_noise) && sum(abs(psd_noise))==0
    paramproj.A = A;
    paramproj.At = At;
    paramproj.epsilon = 1e-10;
    paramproj.y = y;
    paramproj.tight = 0;
    paramproj.verbose = param.verbose -1;
    ffid.prox = @(x,T) proj_b2(x,T,paramproj);
    ffid.eval = @(x) eps;

else
    % Fidelity term for Wiener optimization
    ffid.grad = @(x) 2*At(A(x)-y);
    ffid.eval = @(x) norm(A(x)-y,'fro')^2;
    ffid.beta = 2*param.nu;
end

[sol, infos] = gsp_wiener_optimization(G, y, ffid, psd, psd_noise, param);

end

repo: gspbox


text:
function w = weekday_check_common ( w )

%*****************************************************************************80
%
%% WEEKDAY_CHECK_COMMON makes sure the Common weekday number is between 1 and 7.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer W, the weekday index.
%
  w = i4_wrap ( w, 1, 7 );

  return
end


text:
%SimuNF_nl.m
%tic
function [] = plotTimecourse(x0, str, i)

dt = 0.005;
tf = 2;
switch str
	case 'NF'
		isField = false;
		offset = 1;
		K = [632.4555         0   68.3772         0    3.1623         0
			0  632.4555         0   68.3772         0    3.1623];
	case 'VF'
		isField = true;
		offset = 2;
		K = [632.4555         0   68.3772         0    3.1623         0
			0  632.4555         0   68.3772         0    3.1623];
	case 'AL'
		isField = true;
		offset = 3;
		K = [547.6615   -6.6753  100.8891  -19.1842    3.1623         0;
			8.1756  447.1638  -19.1842   96.0024         0    3.1623];
	case 'AE'
		offset = nan;
		isField = false;
		K = [547.6615   -6.6753  100.8891  -19.1842    3.1623         0;
			8.1756  447.1638  -19.1842   96.0024         0    3.1623];
end

[t1,y1] = LocalSDESolver(0,tf,[-x0(1),-x0(2),0,0,0,0,0,0]',dt, K, isField);
[t2,y2] = LocalSDESolver(0,tf,[-x0(1),-x0(2),0,0,0,0,0,0]',dt, K, isField);
[t3,y3] = LocalSDESolver(0,tf,[-x0(1),-x0(2),0,0,0,0,0,0]',dt, K, isField);
[t4,y4] = LocalSDESolver(0,tf,[-x0(1),-x0(2),0,0,0,0,0,0]',dt, K, isField);
[t5,y5] = LocalSDESolver(0,tf,[-x0(1),-x0(2),0,0,0,0,0,0]',dt, K, isField);

figure(1)
plot(y1(:,1)+x0(1),y1(:,2)+x0(2), 'b-', ...
	y2(:,1)+x0(1),y2(:,2)+x0(2), 'b-', ...
	y3(:,1)+x0(1),y3(:,2)+x0(2), 'b-', ...
	y4(:,1)+x0(1),y4(:,2)+x0(2), 'b-', ...
	y5(:,1)+x0(1),y5(:,2)+x0(2), 'b-', ...
	'Linewidth', 1.2)

if ~isnan(offset)
	figure(2)
	subplot(8,3,(i-1)*3 + offset)
	plot(t1,sqrt(y1(:,3).^2+y1(:,4).^2), 'b-', ...
		t2,sqrt(y2(:,3).^2+y2(:,4).^2), 'b-', ...
		t3,sqrt(y3(:,3).^2+y3(:,4).^2), 'b-', ...
		t4,sqrt(y4(:,3).^2+y4(:,4).^2), 'b-', ...
		t5,sqrt(y5(:,3).^2+y5(:,4).^2), 'b-', ...
		'Linewidth', 1.2)
	figure(1)
end
end

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [xt,ut,ft,ht,flg] = ptrim(f,x,u,x0,u0,h,opts)
% function [xt,ut,ft,ht,flg] = ptrim(f,x,u,x0,u0,h,opts)
%
% DESCRIPTION
%   This function solves for trim states and inputs for the polynomial
%   dynamical system
%      dx/dt = f(x,u)
%   The trim values (xt,ut) satisfy f(xt,ut)=0.  FSOLVE is used to
%   solve these nonlinear equations.  Initial guesses for the trim
%   state/input can be passed to FSOLVE.  Additional equality
%   constraints on the trim condition can be specified in the form
%   h(x,u)=0 where h is a polynomial vector.
%
% INPUTS
%   f: Vector field of polynomial system  (Nx-by-1 polynomial)
%   x: State  (Nx-by-1 vector of pvars)
%   u: Input  (Nu-by-1 vector of pvars)
%   x0: Initial guess for trim state [Optional, Default: x0=0]
%   u0: Initial guess for trim input [Optional, Default: u0=0]
%   h: Equality constraints (Nh-by-1 polynomial) [Optional]
%   opts: Options for fsolve. See fsolve help [Optional]
%
% OUTPUTS
%   xt: Trim state (Nx-by-1 vector)
%   ut: Trim input (Nu-by-1 vector)
%   ft: f evaluated at (xt,ut) (Nx-by-1 vector)
%   ht: h evaluated at (xt,ut) (Nh-by-1 vector)
%      If ptrim was successful finding a trim point then ft:=f(xt,ut)
%      and ht:=h(xt,ut) will both be equal to zero
%   flg: Exit flag returned by fsolve
%
% SYNTAX
%   [xt,ut,ft,ht,flg] = ptrim(f,x,u)
%   [xt,ut,ft,ht,flg] = ptrim(f,x,u,x0,u0)
%   [xt,ut,ft,ht,flg] = ptrim(f,x,u,x0,u0,h)
%   [xt,ut,ft,ht,flg] = ptrim(f,x,u,x0,u0,h,opts)
%
% EXAMPLE
%   pvar x1 x2 u;
%   x = [x1;x2];
%   f = [-2*x1+x2+x1^2-7; x1-3*x2+u+u^2+3];
%
%   % Find a trim condition
%   [xt,ut,ft] = ptrim(f,x,u)
%
%   % Find a trim condition with x2=4
%   h = x2-4;
%   x0 = []; u0 = [];
%   [xt,ut,ft,ht] = ptrim(f,x,u,x0,u0,h)
%
% See also fsolve, plinearize

% PJS 10/7/2009   Initial Coding

% Error Checking
if nargin==3
    x0 = [];
    u0 = [];
    h = [];
    opts = [];
elseif nargin ==4
    u0 = [];
    h = [];
    opts = [];
elseif nargin==5
    h = [];
    opts = [];
elseif nargin==6
    opts = [];
end

% Set default initial starting guesses
Nx = length(x);
if isempty(x0)
    x0 = zeros(Nx,1);
end
if isempty(u0)
    Nu = length(u);
    u0 = zeros(Nu,1);
end

% Set default fsolve options
if isempty(opts)
    opts = optimset('Display','Off');
end

% Trim equations: f(xt,ut)=0 and h(xt,ut)=0
trimeq = [f;h];

% Solve for trim states/inputs, z:=[x;u]
fh = @(z) double(subs(trimeq,[x; u],z));
z0 = [x0(:); u0(:)];
[zt,fht,flg] = fsolve(fh,z0,opts);

xt = zt(1:Nx);
ut = zt(Nx+1:end);
ft = fht(1:Nx);
ht = fht(Nx+1:end);


text:
function [Knots_n] = refine_linear_grid_2d(Knots, old_spacing, ds, imgsz, varargin)
    k_down = 0.5;
    if nargin >= 5
        k_down = varargin{1};
    end
    if k_down == 0.5
        Knots_sz = ceil(imgsz ./ old_spacing) + 1;
        Knots_n = linear_disp_2d(Knots, Knots_sz, ds+1);
    else
        interp_type = 0;
        ksz2 = ceil(imgsz ./ old_spacing) + 1;
        tmp = cat(3, imresize_my(Knots(:,:,1), ksz2, interp_type, 'for_upsampling_2'), ...
                     imresize_my(Knots(:,:,2), ksz2, interp_type, 'for_upsampling_2'));
        Knots_n = tmp;
    end

%     if all( fl(Knots_n(:,end, 1)) == 0)
%         Knots_n(:,end,1) = Knots_n(:,end-1,1);
%     end
%     if all( fl(Knots_n(:,end, 2)) == 0)
%         Knots_n(:,end,2) = Knots_n(:,end-1,2);
%     end
%     if all( fl(Knots_n(end, :, 1)) == 0)
%         Knots_n(end, :, 1) = Knots_n(end-1, :, 1);
%     end
%     if all( fl(Knots_n(end, :, 2)) == 0)
%         Knots_n(end, :, 2) = Knots_n(end-1, :, 2);
%     end
%     size(Knots)
%     Knots_sz
%     ds+1

%     K1 = imresize(Knots(:,:,1), Knots_sz, 'cubic');
%     K2 = imresize(Knots(:,:,2), Knots_sz, 'cubic');
%     Knots_n = cat(3, K1, K2);

%     figure; imagesc(Knots(:,:,1)); colorbar;
%     figure; imagesc(Knots_n(:,:,1)); colorbar;
%     pause;

%     T = mex_2d_linear_disp_double(Knots, imgsz, old_spacing);
%     ksz = ceil(imgsz ./ new_spacing) + 1;
%     Knots_n =...
%         T([1, round(1:new_spacing(1):imgsz(1))], ...
%           [1, round(1:new_spacing(2):imgsz(2))], ...
%           :);

end

repo: pTVreg


text:
addpath(genpath('D:\GitHub\KiloSort2')) % path to kilosort folder
addpath('D:\GitHub\npy-matlab')

pathToYourConfigFile = 'D:\GitHub\KiloSort2\configFiles'; % take from Github folder and put it somewhere else (together with the master_file)
run(fullfile(pathToYourConfigFile, 'configFileBench384.m'))

% common options for every probe
ops.chanMap     = 'D:\GitHub\KiloSort2\configFiles\neuropixPhase3A_kilosortChanMap.mat';
% ops.trange      = [3300 Inf]; % TIME RANGE IN SECONDS TO PROCESS
% ops.trange      = [3400 Inf]; % TIME RANGE IN SECONDS TO PROCESS
ops.trange      = [0 Inf]; % TIME RANGE IN SECONDS TO PROCESS

ops.Th       = 10;     % threshold on projections (like in Kilosort1)
ops.lam      = 40^2;   % weighting on the amplitude penalty (like in Kilosort1, but it has to be much larger)

ops.ThS      = [8 8];  % lower bound on acceptable single spike quality
ops.momentum = [20 400]; % number of samples to average over
ops.minFR = 1/50;

ops.sigmaMask  = 30;

ops.Nfilt       = 768; % max number of clusters
ops.nfullpasses = 1; % how many forward backward passes to do
ops.nPCs        = 3; % how many PCs to project the spikes into

ops.useRAM = 0;

dset = 'c46';

root = fullfile('H:\DATA\Spikes\', dset);
fname = sprintf('%s_aligned_npx_raw.bin', dset);
ops.fbinary     = fullfile(root, fname);
ops.fproc       = 'H:\DATA\Spikes\temp_wh.dat'; % residual from RAM of preprocessed data
ops.dir_rez     = 'H:\DATA\Spikes\';

% preprocess data
rez = preprocessDataSub(ops);

%% learnAndSolve7;
rez = learnAndSolve8(rez);

fGTname = sprintf('%s_spike_samples_npx.npy', dset);
sGT  = readNPY(fullfile('H:\DATA\Spikes\', fGTname));

testClu = rez.st3(:,2);
testRes = rez.st3(:,1);

gtClu = zeros(numel(sGT), 1, 'int32');
gtRes = sGT;

[allScores, allFPrates, allMissRates, allMerges] = ...
    compareClustering2(gtClu, gtRes, testClu, testRes, []);

savePath = fullfile('H:\DATA\Spikes\', dset);
mkdir(savePath)
rezToPhy(rez, savePath);

%%

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function [ know, x ] = p06_sol ( m, know )

%*****************************************************************************80
%
%% P06_SOL returns the solution for problem 06.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input/output, integer KNOW.
%    On input, KNOW is 0, or the index of the previously returned solution.
%    On output, KNOW is 0 if there are no more solutions, or it is the
%    index of the next solution.
%
%    Output, real X(M), the solution, if known.
%
  if ( know == 0 )
    know = 1;
    x(1:m,1) = 420.9687;
  else
    know = 0;
    x(1:m,1) = 0.0;
  end

  return
end

repo: jburkardt-m


text:
function Plot_continuous_data(xf_detect, sr_sub,handles)
% Function that plot the segment of continuous data.
% It will try to plot a minute, maximun.

detect = handles.par.detection;
stdmin = handles.par.stdmin;
stdmax = handles.par.stdmax;

lx = length(xf_detect);

noise_std_detect = median(abs(xf_detect))/0.6745;

thr = stdmin * noise_std_detect;        %thr for detection is based on detect settings.
throut  = stdmax * noise_std_detect;    %thrmax
max_ylim = 15 * noise_std_detect;     %ylim for plotting. aprox thrmax

cla(handles.cont_data);
hold(handles.cont_data, 'on');

plot(handles.cont_data, (1:lx)/sr_sub, xf_detect)

switch detect
    case 'pos'
        plot(handles.cont_data, [0 lx/sr_sub], [thr thr],'-r')
        ylim(handles.cont_data, [-max_ylim/2 max_ylim])
    case 'neg'
        plot(handles.cont_data, [0 lx/sr_sub], [-thr -thr],'-r')
        ylim(handles.cont_data, [-max_ylim max_ylim/2])
    case 'both'
        plot(handles.cont_data, [0 lx/sr_sub], [thr thr],'-r')
        plot(handles.cont_data, [0 lx/sr_sub], [-thr -thr],'-r')
        ylim(handles.cont_data, [-max_ylim max_ylim])
end

xlim(handles.cont_data, [0 lx/sr_sub])
yl = ylim(handles.cont_data);
plot(handles.cont_data, [0 lx/sr_sub], [throut throut],':r')
plot(handles.cont_data, [0 lx/sr_sub], [-throut -throut],':r')
ylim(handles.cont_data,yl);

%xlabel('Time (sec)')

author: csn-le