text:
function n_errors = gsp_test_graph_learning

tol = 1e-5;

n_errors = test_graph_l2(tol) + test_graph_log(tol); %+ test_graph_LX(tol*100) ;
% n_errors = n_errors + test_graph_LX(tol*100);

end

function n_errors = test_graph_LX(tol)

n_errors = 0;

n = 20;
b = 1;

X = rand(20, 5);

params.verbosity = 0;
params.maxit = 10000;

W = gsp_learn_graph_LX_fro(X, b, params);

% zero diagonal
if not(all(diag(W) == 0))
    n_errors = n_errors + 1;
    warning('W has non zero diagonal');
end

% non-negative weights
if not(all(W(:) >= -tol))
    n_errors = n_errors + 1;
    warning('W has non-negative values');
end

% sum equal to n
if (norm(sum(W) - b) / n) > tol
    n_errors = n_errors + 1;
    warning('W has degrees smaller than the limit requested');
end

% W symmetric
if (norm(W - W', 'fro') / n) > tol
    n_errors = n_errors + 1;
    warning('W not symmetric');
end

end

function n_errors = test_graph_l2(tol)

n_errors = 0;

n = 20;

Z = gsp_distanz(randn(n)).^2;

params.maxit = 100000;
params.tol = 0.1 * tol;
W = gsp_learn_graph_l2_degrees(1*Z, 1, params);

% zero diagonal
if not(all(diag(W) == 0))
    n_errors = n_errors + 1;
end

% non-negative weights
if not(all(W(:) >= -tol))
    n_errors = n_errors + 1;
end

% sum equal to n
if abs(sum(W(:) - n) / n > tol)
    n_errors = n_errors + 1;
end

% % step-size correct
% params.step_size = .9;
% W2 = gsp_learn_graph_l2_degrees(10*Z, 10, params);
% if not(norm(W - W2, 'fro')/norm(W, 'fro') < 10*tol)
%     n_errors = n_errors + 1;
% end

end

function n_errors = test_graph_log(tol)

n_errors = 0;

n = 20;

Z = gsp_distanz(rand(n)).^2;

% We'll get a very sparse graph and check the degrees quality later
W = gsp_learn_graph_log_degrees(Z, 1, 1);

params.tol = tol * 1e-3;
params.maxit = 50000;
params.step_size = .9;
params.verbosity = 1;
W = gsp_learn_graph_log_degrees(Z, .031, .003, params);

% zero diagonal
if not(all(diag(W) == 0))
    n_errors = n_errors + 1;
end

% non-negative weights
if not(all(W(:) >= -tol))
    n_errors = n_errors + 1;
end

% positive degrees
if not(all(sum(W) >= tol * max(sum(W))))
    n_errors = n_errors + 1;
end

end

repo: gspbox


text:
function s = get_limited_categories(categorical_obj, lenlimit)
    % get_limited_categories gets the categories for a categorical, but designed for display
    % which means you limit the character length.
    %
    % s = get_limited_categories(categorical_obj) will get a string contianning the first few and
    % last few categories with ellipses in between.
    % default length is 80

    c = categories(categorical_obj);
    if nargin == 1
        lenlimit = 80;
    end

    s = strjoin(c,', ');
    if numel(s) > lenlimit
        commas = find(s==',');
        break1=max(commas(commas<25));
        break2=min(commas(commas>(length(s)-25)));
        s=[s(1:break1),'...',s(break2:end)];
    end
end

author: CelsoReyes


text:
function T = isotopes(varargin)
%ISOTOPES Returns NMR properties of the isotopes
% T = isotopes(i1, i2, ...)
%
% EXAMPLE
%     H = isotopes('1H'), T = isotopes('3H'), or if the atomic mass is not
%     specified, i.e. I = isotopes('H'), the function returns array of
%     structures with data for all isotopes of the specified element.
%     isotopes('n') returns NMR properties for neutron. Isotopes called
%     without arguments returns the whole table.
%
% The output structure has the following fields:
%     .name - element name;
%     .isotope - element symbol and mass number, i.e. 19F;
%     .Z - atomic number;
%     .mass - mass number;
%     .I - nuclear spin;
%     .gamma - magnetogyric ratio in rad/(T s) calculated from the resonant
%              frequency: 2*pi*F(@1T);
%     .ResFreq - resonant frequency in Hertz for an applied field B0 of 1 tesla
%                (in cgs units, 10 kilogauss);
%     .RelSensB - sensitivity relative to 1H (=1) assuming an equal number
%                 of nuclei, constant temperature, and constant B0:
%                 0.0076508(mu/mu_N)^3(I + 1)/I^2;
%     .RelSensF - sensitivity relative to 1H (=1) assuming an equal number
%                 of nuclei, constant temperature, and constant resonance
%                 frequency: O.23871(mu/mu_N)(I + 1);
%     .u - nuclear magnetic moment in units of the nuclear magneton, (mu/mu_N);
%     .Q - nuclear quadrupole moment in units of femtometers squared
%          (1 fm2 = 10^-2 barn)
%
%
% REFERENCES
% 1. "Nuclear spins, moments, and other data related to NMR spectroscopy" (9-93)
%    CRC Handbook of Chemistry and Physics, 83th Ed., CRC Press, Boca Raton, FL,
%    2002.
% 2. Stone, N. J., <www.nndc.bnl.gov/nndc/stone_moments/>
% 3. http://www.hbcpnetbase.com/

% Copyright (c) 2002 by Peter L. Volegov (volegov@unm.edu)
% All Rights Reserved.

%% Table of the Isotopes
% Element, Name, Z, AtomicMass, Spin, F@1T/MHz, RelSens(H), RelSens(F), mu/mu_N, Q/fm2
TheTable = {
'n', 'neutron', 1, 1, 1/2, 29.1647, 0.32139, 0.6850, -1.913042720, 0.0000; ...
'H', 'Hydrogen', 1, 1, 1/2, 42.5775, 1.00000, 1.0000, 2.792847337, 0.0000; ...
'H', 'Deuterium', 1, 2, 1, 6.5359, 0.00965, 0.4093, 0.857438228, 0.2860; ...
'H', 'Tritium', 1, 3, 1/2, 45.4148, 1.21354, 1.0666, 2.978962500, 0.0000; ...
'He', 'Helium', 2, 3, 1/2, 32.4360, 0.44212, 0.7618, -2.127624800, 0.0000; ...
'Li', 'Lithium', 3, 6, 1, 6.2661, 0.00850, 0.3925, 0.822046700, -0.0808; ...
'Li', 'Lithium', 3, 7, 3/2, 16.5483, 0.29356, 1.9433, 3.256440000, -4.0100; ...
'Be', 'Beryllium', 4, 9, 3/2, 5.9842, 0.01388, 0.7027, -1.177600000, 5.2880; ...
'B', 'Boron', 5, 10, 3, 4.5752, 0.01985, 1.7193, 1.800645000, 8.4590; ...
'B', 'Boron', 5, 11, 3/2, 13.6630, 0.16522, 1.6045, 2.688649000, 4.0590; ...
'C', 'Carbon', 6, 13, 1/2, 10.7084, 0.01591, 0.2515, 0.702411800, 0.0000; ...
'N', 'Nitrogen', 7, 14, 1, 3.0777, 0.00101, 0.1928, 0.403761000, 2.0440; ...
'N', 'Nitrogen', 7, 15, 1/2, 4.3173, 0.00104, 0.1014, -0.283188800, 0.0000; ...
'O', 'Oxygen', 8, 17, 5/2, 5.7742, 0.02910, 1.5822, -1.893790000, -2.5580; ...
'F', 'Fluorine', 9, 19, 1/2, 40.0776, 0.83400, 0.9413, 2.628868000, 0.0000; ...
'Ne', 'Neon', 10, 21, 3/2, 3.3631, 0.00246, 0.3949, -0.661797000, 10.1550; ...
'Na', 'Sodium', 11, 23, 3/2, 11.2688, 0.09270, 1.3233, 2.217522000, 10.4000; ...
'Mg', 'Magnesium', 12, 25, 5/2, 2.6083, 0.00268, 0.7147, -0.855450000, 19.9400; ...
'Al', 'Aluminum', 13, 27, 5/2, 11.1031, 0.20689, 3.0424, 3.641507000, 14.6600; ...
'Si', 'Silicon', 14, 29, 1/2, 8.4655, 0.00786, 0.1988, -0.555290000, 0.0000; ...
'P', 'Phosphorus', 15, 31, 1/2, 17.2515, 0.06652, 0.4052, 1.131600000, 0.0000; ...
'S', 'Sulfur', 16, 33, 3/2, 3.2717, 0.00227, 0.3842, 0.643821200, -6.7800; ...
'Cl', 'Chlorine', 17, 35, 3/2, 4.1765, 0.00472, 0.4905, 0.821874300, -8.1650; ...
'Cl', 'Chlorine', 17, 37, 3/2, 3.4765, 0.00272, 0.4083, 0.684123600, -6.4350; ...
'Ar', 'Argone', 18, 37, 3/2, 5.8190, 0.01276, 0.6833, 1.145000000, 0.0000; ...
'Ar', 'Argone', 18, 39, 7/2, 3.4600, 0.01130, 1.7079, -1.590000000, 0.0000; ...
'K', 'Potassium', 19, 39, 3/2, 1.9893, 0.00051, 0.2336, 0.391466200, 5.8500; ...
'K', 'Potassium', 19, 40, 4, 2.4737, 0.00523, 1.5493, -1.298100000, -7.3000; ...
'K', 'Potassium', 19, 41, 3/2, 1.0919, 8.00000, 0.1282, 0.214870100, 7.1100; ...
'Ca', 'Calcium', 20, 43, 7/2, 2.8688, 0.00642, 1.4150, -1.317260000, -4.0800; ...
'Sc', 'Scandium', 21, 45, 7/2, 10.3591, 0.30244, 5.1093, 4.756487000, -22.0000; ...
'Ti', 'Titanium', 22, 47, 5/2, 2.4041, 0.00210, 0.6587, -0.788480000, 30.2000; ...
'Ti', 'Titanium', 22, 49, 7/2, 2.4048, 0.00378, 1.1861, -1.104170000, 24.7000; ...
'V', 'Vanadium', 23, 50, 6, 4.2505, 0.05571, 5.5904, 3.345689000, 21.0000; ...
'V', 'Vanadium', 23, 51, 7/2, 11.2133, 0.38360, 5.5306, 5.148705700, -5.2000; ...
'Cr', 'Chromium', 24, 53, 3/2, 2.4115, 0.00091, 0.2832, -0.474540000, -15.0000; ...
'Mn', 'Manganese', 25, 55, 5/2, 10.5763, 0.17881, 2.8980, 3.468720000, 33.0000; ...
'Fe', 'Iron', 26, 57, 1/2, 1.3816, 0.00003, 0.0324, 0.090623000, 0.0000; ...
'Co', 'Cobalt', 27, 59, 7/2, 10.0770, 0.27841, 4.9702, 4.627000000, 42.0000; ...
'Ni', 'Nickel', 28, 61, 3/2, 3.8114, 0.00359, 0.4476, -0.750020000, 16.2000; ...
'Cu', 'Copper', 29, 63, 3/2, 11.2982, 0.09342, 1.3268, 2.223290000, -22.0000; ...
'Cu', 'Copper', 29, 65, 3/2, 12.1030, 0.11484, 1.4213, 2.381670000, -20.4000; ...
'Zn', 'Zinc', 30, 67, 5/2, 2.6694, 0.00287, 0.7314, 0.875479000, 15.0000; ...
'Ga', 'Gallium', 31, 69, 3/2, 10.2478, 0.06971, 1.2034, 2.016590000, 17.1000; ...
'Ga', 'Gallium', 31, 71, 3/2, 13.0208, 0.14300, 1.5291, 2.562270000, 10.7000; ...
'Ge', 'Germanium', 32, 73, 9/2, 1.4897, 0.00141, 1.1546, -0.879467700, -19.6000; ...
'As', 'Arsenic', 33, 75, 3/2, 7.3150, 0.02536, 0.8590, 1.439475000, 31.4000; ...
'Se', 'Selenium', 34, 77, 1/2, 8.1571, 0.00703, 0.1916, 0.535060000, 0.0000; ...
'Br', 'Bromine', 35, 79, 3/2, 10.7042, 0.07945, 1.2570, 2.106400000, 31.3000; ...
'Br', 'Bromine', 35, 81, 3/2, 11.5384, 0.09951, 1.3550, 2.270562000, 26.2000; ...
'Kr', 'Krypton', 36, 83, 9/2, 1.6442, 0.00190, 1.2744, -0.970669000, 25.9000; ...
'Rb', 'Rubidium', 37, 85, 5/2, 4.1254, 0.01061, 1.1304, 1.353030000, 27.6000; ...
'Rb', 'Rubidium', 37, 87, 3/2, 13.9811, 0.17703, 1.6418, 2.751240000, 13.3500; ...
'Sr', 'Strontium', 38, 87, 9/2, 1.8525, 0.00272, 1.4358, -1.093603000, 33.5000; ...
'Y', 'Yttrium', 39, 89, 1/2, 2.0949, 0.00012, 0.0492, -0.137415400, 0.0000; ...
'Zr', 'Zirconium', 40, 91, 5/2, 3.9748, 0.00949, 1.0891, -1.303620000, -17.6000; ...
'Nb', 'Niobium', 41, 93, 9/2, 10.4523, 0.48821, 8.1011, 6.170500000, -32.0000; ...
'Mo', 'Molybdenum', 42, 95, 5/2, 2.7874, 0.00327, 0.7638, -0.914200000, -2.2000; ...
'Mo', 'Molybdenum', 42, 97, 5/2, 2.8463, 0.00349, 0.7799, -0.933500000, 25.5000; ...
'Tc', 'Technetium', 43, 99, 9/2, 9.6294, 0.38174, 7.4633, 5.684700000, -12.9000; ...
'Ru', 'Ruthenium', 44, 99, 5/2, 1.9553, 0.00113, 0.5358, -0.641300000, 7.9000; ...
'Ru', 'Ruthenium', 44, 101, 5/2, 2.1916, 0.00159, 0.6005, -0.718800000, 45.7000; ...
'Rh', 'Rhodium', 45, 103, 1/2, 1.3477, 0.00003, 0.0317, -0.088400000, 0.0000; ...
'Pd', 'Palladium', 46, 105, 5/2, 1.9570, 0.00113, 0.5364, -0.642000000, 66.0000; ...
'Ag', 'Silver', 47, 107, 1/2, 1.7331, 0.00007, 0.0407, -0.113679600, 0.0000; ...
'Ag', 'Silver', 47, 109, 1/2, 1.9924, 0.00010, 0.0468, -0.130690600, 0.0000; ...
'Cd', 'Cadmium', 48, 111, 1/2, 9.0692, 0.00966, 0.2130, -0.594886100, 0.0000; ...
'Cd', 'Cadmium', 48, 113, 1/2, 9.4871, 0.01106, 0.2228, -0.622300900, 0.0000; ...
'In', 'Indium', 49, 113, 9/2, 9.3655, 0.35121, 7.2588, 5.528900000, 79.9000; ...
'In', 'Indium', 49, 115, 9/2, 9.3856, 0.35348, 7.2744, 5.540800000, 81.0000; ...
'Sn', 'Tin', 50, 115, 1/2, 14.0077, 0.03561, 0.3290, -0.918830000, 0.0000; ...
'Sn', 'Tin', 50, 117, 1/2, 15.2610, 0.04605, 0.3584, -1.001040000, 0.0000; ...
'Sn', 'Tin', 50, 119, 1/2, 15.9660, 0.05273, 0.3750, -1.047280000, 0.0000; ...
'Sb', 'Antimony', 51, 121, 5/2, 10.2551, 0.16302, 2.8100, 3.363400000, -36.0000; ...
'Sb', 'Antimony', 51, 123, 7/2, 5.5532, 0.04659, 2.7389, 2.549800000, -49.0000; ...
'Te', 'Tellurium', 52, 123, 1/2, 11.2349, 0.01837, 0.2639, -0.736947800, 0.0000; ...
'Te', 'Tellurium', 52, 125, 1/2, 13.5454, 0.03220, 0.3181, -0.888505100, 0.0000; ...
'I', 'Iodine', 53, 127, 5/2, 8.5778, 0.09540, 2.3504, 2.813273000, -71.0000; ...
'Xe', 'Xenon', 54, 129, 1/2, 11.8604, 0.02162, 0.2786, -0.777976300, 0.0000; ...
'Xe', 'Xenon', 54, 131, 3/2, 3.5159, 0.00282, 0.4129, 0.691861900, -11.4000; ...
'Cs', 'Cesium', 55, 133, 7/2, 5.6234, 0.04838, 2.7735, 2.582025000, -0.3430; ...
'Ba', 'Barium', 56, 135, 3/2, 4.2582, 0.00500, 0.5001, 0.837943000, 16.0000; ...
'Ba', 'Barium', 56, 137, 3/2, 4.7634, 0.00700, 0.5594, 0.937365000, 24.5000; ...
'La', 'Lanthanum', 57, 138, 5, 5.6615, 0.09404, 5.3188, 3.713646000, 45.0000; ...
'La', 'Lanthanum', 57, 139, 7/2, 6.0612, 0.06058, 2.9895, 2.783045500, 20.0000; ...
'Ce', 'Cerium', 58, 137, 3/2, 4.8800, 0.00752, 0.5729, 0.960000000, 0.0000; ...
'Ce', 'Cerium', 58, 139, 3/2, 5.3900, 0.01012, 0.6326, 1.060000000, 0.0000; ...
'Ce', 'Cerium', 58, 141, 7/2, 2.3700, 0.00364, 1.1708, 1.090000000, 0.0000; ...
'Pr', 'Praseodymium', 59, 141, 5/2, 13.0359, 0.33483, 3.5720, 4.275400000, -5.9000; ...
'Nd', 'Neodymium', 60, 143, 7/2, 2.3190, 0.00339, 1.1440, -1.065000000, -63.0000; ...
'Nd', 'Neodymium', 60, 145, 7/2, 1.4290, 0.00079, 0.7047, -0.656000000, -33.0000; ...
'Pm', 'Promethium', 61, 143, 5/2, 11.5900, 0.23510, 3.1748, 3.800000000, 0.0000; ...
'Pm', 'Promethium', 61, 147, 7/2, 5.6200, 0.04827, 2.7714, 2.580000000, 74.0000; ...
'Sm', 'Samarium', 62, 147, 7/2, 1.7748, 0.00152, 0.8753, -0.814900000, -26.0000; ...
'Sm', 'Samarium', 62, 149, 7/2, 14631.0000, 0.00085, 0.7216, -0.671800000, 7.4000; ...
'Eu', 'Europium', 63, 151, 5/2, 10.5856, 0.17929, 2.9006, 3.471800000, 90.3000; ...
'Eu', 'Europium', 63, 153, 5/2, 4.6745, 0.01544, 1.2809, 1.533100000, 241.0000; ...
'Gd', 'Gadolinium', 64, 155, 3/2, 1.3120, 0.00015, 0.1541, -0.258200000, 127.0000; ...
'Gd', 'Gadolinium', 64, 157, 3/2, 1.7200, 0.00033, 0.2020, -0.338500000, 135.0000; ...
'Tb', 'Terbium', 65, 159, 3/2, 10.2300, 0.06945, 1.2019, 2.014000000, 143.2000; ...
'Dy', 'Dysprosium', 66, 161, 5/2, 1.4654, 0.00048, 0.4015, -0.480600000, 250.7000; ...
'Dy', 'Dysprosium', 66, 163, 5/2, 2.0508, 0.00130, 0.5619, 0.672600000, 265.0000; ...
'Ho', 'Holmium', 67, 165, 7/2, 9.0883, 0.20423, 4.4825, 4.173000000, 358.0000; ...
'Er', 'Erbium', 68, 167, 7/2, 1.2281, 0.00050, 0.6057, -0.563900000, 356.5000; ...
'Tm', 'Thulium', 69, 169, 1/2, 3.5310, 0.00057, 0.0829, -0.231600000, 0.0000; ...
'Yb', 'Ytterbium', 70, 171, 1/2, 7.5261, 0.00552, 0.1768, 0.493670000, 0.0000; ...
'Yb', 'Ytterbium', 70, 173, 5/2, 2.0730, 0.00135, 0.5680, -0.679890000, 280.0000; ...
'Lu', 'Lutetium', 71, 175, 7/2, 4.8626, 0.03128, 2.3983, 2.232700000, 349.0000; ...
'Lu', 'Lutetium', 71, 176, 7, 3.4510, 0.03975, 6.0516, 3.169000000, 497.0000; ...
'Hf', 'Hafnium', 72, 177, 7/2, 1.7282, 0.00140, 0.8524, 0.793500000, 336.5000; ...
'Hf', 'Hafnium', 72, 179, 9/2, 1.0856, 0.00055, 0.8414, -0.640900000, 379.3000; ...
'Ta', 'Tantalum', 73, 180, 9, 4.0870, 0.10610, 11.5175, 4.825000000, 0.0000; ...
'Ta', 'Tantalum', 73, 181, 7/2, 5.1627, 0.03744, 2.5463, 2.370500000, 317.0000; ...
'W', 'Tungsten', 74, 183, 1/2, 1.7957, 0.00008, 0.0422, 0.117784800, 0.0000; ...
'Re', 'Rhenium', 75, 185, 5/2, 9.7176, 0.13870, 2.6627, 3.187100000, 218.0000; ...
'Re', 'Rhenium', 75, 187, 5/2, 9.8170, 0.14300, 2.6900, 3.219700000, 207.0000; ...
'Os', 'Osmium', 76, 187, 1/2, 0.9856, 0.00001, 0.0231, 0.064651890, 0.0000; ...
'Os', 'Osmium', 76, 189, 3/2, 3.3536, 0.00244, 0.3938, 0.659933000, 85.6000; ...
'Ir', 'Iridium', 77, 191, 3/2, 0.7658, 0.00003, 0.0899, 0.150700000, 81.6000; ...
'Ir', 'Iridium', 77, 191, 3/2, 0.8319, 0.00004, 0.0977, 0.163700000, 75.1000; ...
'Pt', 'Platinum', 78, 195, 1/2, 9.2922, 0.01039, 0.2182, 0.609520000, 0.0000; ...
'Au', 'Gold', 79, 197, 3/2, 0.7406, 0.00003, 0.0870, 0.145746000, 54.7000; ...
'Hg', 'Mercury', 80, 199, 1/2, 7.7123, 0.00594, 0.1811, 0.505885500, 0.0000; ...
'Hg', 'Mercury', 80, 201, 3/2, 2.8469, 0.00149, 0.3343, -0.560225700, 38.6000; ...
'Tl', 'Thallium', 81, 203, 1/2, 24.7316, 0.19598, 0.5809, 1.622257900, 0.0000; ...
'Tl', 'Thallium', 81, 205, 1/2, 24.9749, 0.20182, 0.5866, 1.638214600, 0.0000; ...
'Pb', 'Lead', 82, 207, 1/2, 9.0340, 0.00955, 0.2122, 0.592580000, 0.0000; ...
'Bi', 'Bismuth', 83, 209, 9/2, 6.9630, 0.14433, 5.3967, 4.110600000, 51.6000; ...
'Po', 'Polonium', 84, 209, 1/2, 11.7000, 0.02096, 0.2757, 0.770000000, 0.0000; ...
'Rn', 'Radon', 86, 211, 1/2, 9.1600, 0.00997, 0.2152, 0.601000000, 0.0000; ...
'Fr', 'Francium', 87, 223, 3/2, 5.9500, 0.01362, 0.6982, 1.170000000, 117.0000; ...
'Ra', 'Radium', 88, 223, 3/2, 1.3746, 0.00017, 0.1614, 0.270500000, 125.0000; ...
'Ra', 'Radium', 88, 225, 1/2, 11.1870, 0.01814, 0.2627, -0.733800000, 0.0000; ...
'Ac', 'Actinium', 89, 227, 3/2, 5.6000, 0.01131, 0.6564, 1.100000000, 170.0000; ...
'Th', 'Thorium', 90, 229, 5/2, 1.4000, 0.00042, 0.3843, 0.460000000, 430.0000; ...
'Pa', 'Protactinium', 91, 231, 3/2, 10.2000, 0.06903, 1.1995, 2.010000000, -172.0000; ...
'U', 'Uranium', 92, 235, 7/2, 0.8300, 0.00015, 0.4082, -0.380000000, 493.6000; ...
'Np', 'Neptunium', 93, 237, 5/2, 9.5700, 0.13264, 2.6234, 3.140000000, 388.6000; ...
'Pu', 'Plutonium', 94, 239, 1/2, 3.0900, 0.00038, 0.0727, 0.203000000, 0.0000; ...
'Am', 'Americium', 95, 243, 5/2, 4.6000, 0.01446, 1.2532, 1.500000000, 421.0000; ...
};

%% Find isotopes
if nargin == 0
    % Return the whole table
    T = GetDataFromTheTable(TheTable);

 else
    T = [];
    for j = 1:nargin
        curIn = varargin{j};
        if ischar(curIn)
            nn = isstrprop(curIn, 'alpha');
            curElem = curIn(nn);

            % Find the element symbol
            ie = strmatch(curElem, {TheTable{:, 1}}, 'exact');

            % Isotope mass number
            nn = isstrprop(curIn,  'digit');
            curMass = curIn;
            curMass(~nn) = ' ';
            curMass = strtok(curMass);
            if ~isempty(curMass)
                curMass = str2double(curMass);
                if ~isnan(curMass)
                    im = find(curMass == [TheTable{:, 4}]);
                    ie = intersect(ie, im);
                end
            end

            % Get the isotop(s) properties
            if isempty(ie)
                warning('MATLAB:isotopes:CanNotFindTheElement', ...
                        ['Can not find the element: ' curIn]);
                curT = GetDataFromTheTable;
                curT.name = '<< not found >>';
                curT.isotop = curIn;
                curT.mass = curMass;
                T = [T; curT(:)];
            else
                curT = GetDataFromTheTable(TheTable, ie);
                T = [T; curT(:)];
            end

        else
            error('Input argument(s) should be char, i.e. ''19F''');
        end
    end
end

%
function T = GetDataFromTheTable(TheTable, TblInd)
%Returns data in a structure

% Factor to Hz
fctHz = 1e6;

% Create the output structure template
StructTmpl = struct('name', '', 'isotop', '', 'Z', 0, 'mass', 0, 'I', 0, ...
                    'gamma', 0, 'ResFreq', 0, 'RelSensB', 0, ...
                    'RelSensF', 0, 'u', 0, 'Q', 0);
if nargin < 1
    T = StructTmpl;
    return;

elseif (nargin < 2)
    TblInd = (1:size(TheTable, 1));

elseif isempty(TblInd)
    T = StructTmpl;
    return;

end

nIndNum = length(TblInd);
T = repmat(StructTmpl, nIndNum, 1);
for i = 1:nIndNum
    j = TblInd(i);
    T(i).name = TheTable{j, 2};
    T(i).isotop = sprintf('%.0d%s', TheTable{j, 4}, TheTable{j, 1});
    T(i).Z = TheTable{j, 3};
    T(i).mass = TheTable{j, 4};
    T(i).I = TheTable{j, 5};
    T(i).gamma = 2*pi*fctHz*TheTable{j, 6};
    T(i).ResFreq = TheTable{j, 6}*fctHz;
    T(i).RelSensB = TheTable{j, 7};
    T(i).RelSensF = TheTable{j, 8};
    T(i).u = TheTable{j, 9};
    T(i).Q = TheTable{j, 10};
end

author: Sable

repo: mcbench-benchmarks


text:
function pass = test_flipshiftrotate( )
% Test diskfun flip, rotate, flipxy and circshift functions.

tol = 1000*chebfunpref().cheb2Prefs.chebfun2eps;

% See the random number generator.

%gaussian shifted
f1 = @(x, y) 10*exp(-20*(x-.5).^2-20*(y+.5).^2) ;
h1 = @(x, y) 10*exp(-20*(-x-.5).^2-20*(y+.5).^2) ;
g1 = @(x, y) 10*exp(-20*(x-.5).^2-20*(-y+.5).^2) ;
f = diskfun(f1);
h = diskfun(h1);
g = diskfun(g1);
w = diskfun(@(x,y)x).*(g +h +f)+diskfun(@(x,y) y);
r1 = @(x,y) x.*(f1(x,y) + h1(x,y) + g1(x,y))+y;
ud = diskfun(@(x,y) r1(x,-y));
lr = diskfun(@(x,y) r1(-x, y));
tp = diskfun(@(x,y) r1(y,x));
pass(1) =  ( norm( fliplr(f) - h )  < tol);
pass(2) = ( norm( flipud(f) - g)  < tol) ;
pass(3) = ( norm( flipud(w) - ud) < tol);
pass(4) = ( norm( fliplr(w) -lr) < tol);
pass(5) = ( norm( flipxy(w)-tp) < tol);
pass(6) = ( norm( flipdim(w,1)-ud) < tol);
pass(7) = ( norm( flipdim(w,2)-lr) < tol);

%test rotate/circshift
th = [0, pi/2, 5*pi/7,  -pi/3, -pi/4, -5*pi ];
[theta, rad] = meshgrid(trigpts(30, [-pi, pi]), chebpts(30, [0, 1]));
w = @(t, r) exp(-20*(r.*cos(t)-.5).^2-20.*(r.*sin(t)).^2 );
f = diskfun(w, 'polar');
for j = 1:length(th)
    g = rotate(f, th(j));
    m = circshift(f, th(j));
    pass(j+7) = ( norm(w(theta-th(j), rad)-g(theta, rad, 'polar'))  < tol ...
        || norm(g-m) ==0);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% (Internal) Design the wavelet decomposition filters for wavedet algorithm
%
% Prototype:
% ----------
% q_filters = qs_filter_design(scales, fs, N);
%
% Description:
% ------------
% Mimics the transfer function of the filters used for ECG delineation in
% [Martinez et al. 2004] for an arbitrary sampling frequency and filter order N.
%
% Martinez et al. "A Wavelet-Based ECG Delineator: Evaluation on Standard
% Databases" IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 51, NO. 4,
% APRIL 2004.
%
% WARNING :
% ---------
% As this routines iterates through several configurations in order to
% converge, the user should check the transfer functions of the filters. My
% suggestion is once you obtain a desired filter bank for a given Fs or
% configuration, save or cache it in a .mat file in order to use it during
% operation. For example, if you usually work with signals sampled at 360
% Hz, a good choice is to have a cached version of the filters for this Fs
% in a .mat file called "wt_filters_6 scales_360 Hz.mat". You can use this
% function on-line with your algorithm at your own risk.
%
% Examples :
% ----------
%
% q_filters = qs_filter_design(4, 250);
%
% q_filters = qs_filter_design(5, 360);
%
% q_filters = qs_filter_design(6, 1000);
%
% you can check filter characteristics using:
%
% fss = [250 360 500 1000]; %Hz
% for fs  = fss
%     fvtool(q_filters, 'fs', fs )
% end
%
% Author: Mariano Llamedo Soria (llamedom at frba.utn.edu.ar)
% Version: 0.1 beta
% Birthdate: 17/2/11
% Last update: 22/02/13
% Copyright 2008-2015
%
function q_filters = qs_filter_design(scales, fs, N)

design_iter_attemps = 10;

if( nargin < 2 || isempty(fs) )
    fs = 250; %Hz
end

if( nargin < 3 || isempty(N) )
    %valor emp�rico obtenido de varios dise�os.
    N = max(10, round(fs*4/30 + 16+2/3));
end

%Pruebo que N no sea demasiado distinto a lo recomendable.
recommended_N = max(10, round(fs*4/30 + 16+2/3));
if( abs(N - recommended_N) > 0.1*N )
    warning(['Check the transfer functions of the differentiator filters designed. Recommended order N = ' num2str(recommended_N) ]);
end

%frecuencia a la que est� dise�ado el delineador, y que se toma para
%referencia para que las escalas signifiquen lo mismo a cualquier Fs.
f_ref = 250; %Hz
f_ratio = f_ref/fs;

% Funciones de transferencia correctas para el dise�o de los filtros utilizadas
% por el wavedet a 250 Hz.
empirical_tf = { ...
    [2;-2] ...
    [0.250000000000000;0.750000000000000;0.500000000000000;-0.500000000000000;-0.750000000000000;-0.250000000000000;] ...
    [0.0312500000000000;0.0937500000000000;0.187500000000000;0.312500000000000;0.343750000000000;0.281250000000000;0.125000000000000;-0.125000000000000;-0.281250000000000;-0.343750000000000;-0.312500000000000;-0.187500000000000;-0.0937500000000000;-0.0312500000000000;] ...
    [0.00390625000000000;0.0117187500000000;0.0234375000000000;0.0390625000000000;0.0585937500000000;0.0820312500000000;0.109375000000000;0.140625000000000;0.160156250000000;0.167968750000000;0.164062500000000;0.148437500000000;0.121093750000000;0.0820312500000000;0.0312500000000000;-0.0312500000000000;-0.0820312500000000;-0.121093750000000;-0.148437500000000;-0.164062500000000;-0.167968750000000;-0.160156250000000;-0.140625000000000;-0.109375000000000;-0.0820312500000000;-0.0585937500000000;-0.0390625000000000;-0.0234375000000000;-0.0117187500000000;-0.00390625000000000;] ...
    [0.000488281250000000;0.00146484375000000;0.00292968750000000;0.00488281250000000;0.00732421875000000;0.0102539062500000;0.0136718750000000;0.0175781250000000;0.0219726562500000;0.0268554687500000;0.0322265625000000;0.0380859375000000;0.0444335937500000;0.0512695312500000;0.0585937500000000;0.0664062500000000;0.0727539062500000;0.0776367187500000;0.0810546875000000;0.0830078125000000;0.0834960937500000;0.0825195312500000;0.0800781250000000;0.0761718750000000;0.0708007812500000;0.0639648437500000;0.0556640625000000;0.0458984375000000;0.0346679687500000;0.0219726562500000;0.00781250000000000;-0.00781250000000000;-0.0219726562500000;-0.0346679687500000;-0.0458984375000000;-0.0556640625000000;-0.0639648437500000;-0.0708007812500000;-0.0761718750000000;-0.0800781250000000;-0.0825195312500000;-0.0834960937500000;-0.0830078125000000;-0.0810546875000000;-0.0776367187500000;-0.0727539062500000;-0.0664062500000000;-0.0585937500000000;-0.0512695312500000;-0.0444335937500000;-0.0380859375000000;-0.0322265625000000;-0.0268554687500000;-0.0219726562500000;-0.0175781250000000;-0.0136718750000000;-0.0102539062500000;-0.00732421875000000;-0.00488281250000000;-0.00292968750000000;-0.00146484375000000;-0.000488281250000000;] ...
    [6.10351562500000e-05;0.000183105468750000;0.000366210937500000;0.000610351562500000;0.000915527343750000;0.00128173828125000;0.00170898437500000;0.00219726562500000;0.00274658203125000;0.00335693359375000;0.00402832031250000;0.00476074218750000;0.00555419921875000;0.00640869140625000;0.00732421875000000;0.00830078125000000;0.00933837890625000;0.0104370117187500;0.0115966796875000;0.0128173828125000;0.0140991210937500;0.0154418945312500;0.0168457031250000;0.0183105468750000;0.0198364257812500;0.0214233398437500;0.0230712890625000;0.0247802734375000;0.0265502929687500;0.0283813476562500;0.0302734375000000;0.0322265625000000;0.0339965820312500;0.0355834960937500;0.0369873046875000;0.0382080078125000;0.0392456054687500;0.0401000976562500;0.0407714843750000;0.0412597656250000;0.0415649414062500;0.0416870117187500;0.0416259765625000;0.0413818359375000;0.0409545898437500;0.0403442382812500;0.0395507812500000;0.0385742187500000;0.0374145507812500;0.0360717773437500;0.0345458984375000;0.0328369140625000;0.0309448242187500;0.0288696289062500;0.0266113281250000;0.0241699218750000;0.0215454101562500;0.0187377929687500;0.0157470703125000;0.0125732421875000;0.00921630859375000;0.00567626953125000;0.00195312500000000;-0.00195312500000000;-0.00567626953125000;-0.00921630859375000;-0.0125732421875000;-0.0157470703125000;-0.0187377929687500;-0.0215454101562500;-0.0241699218750000;-0.0266113281250000;-0.0288696289062500;-0.0309448242187500;-0.0328369140625000;-0.0345458984375000;-0.0360717773437500;-0.0374145507812500;-0.0385742187500000;-0.0395507812500000;-0.0403442382812500;-0.0409545898437500;-0.0413818359375000;-0.0416259765625000;-0.0416870117187500;-0.0415649414062500;-0.0412597656250000;-0.0407714843750000;-0.0401000976562500;-0.0392456054687500;-0.0382080078125000;-0.0369873046875000;-0.0355834960937500;-0.0339965820312500;-0.0322265625000000;-0.0302734375000000;-0.0283813476562500;-0.0265502929687500;-0.0247802734375000;-0.0230712890625000;-0.0214233398437500;-0.0198364257812500;-0.0183105468750000;-0.0168457031250000;-0.0154418945312500;-0.0140991210937500;-0.0128173828125000;-0.0115966796875000;-0.0104370117187500;-0.00933837890625000;-0.00830078125000000;-0.00732421875000000;-0.00640869140625000;-0.00555419921875000;-0.00476074218750000;-0.00402832031250000;-0.00335693359375000;-0.00274658203125000;-0.00219726562500000;-0.00170898437500000;-0.00128173828125000;-0.000915527343750000;-0.000610351562500000;-0.000366210937500000;-0.000183105468750000;-6.10351562500000e-05;] ...
                };

Grid_size = 1024;

% Creo una grilla de muestreo en frecuencia logaritmica para que la
% zona en que derivan los filtros sea una recta de pendiente constante.
F_log = logspace(-3, min(0, log10(fs/f_ref)), Grid_size) * pi;

filter_count = 1;

for ii = rowvec(scales)

    [g_amp F] = freqz(empirical_tf{ii},1, F_log);
    F = F * f_ratio / pi;
    g_amp = abs(g_amp);

    %averiguo hasta qu� muestra se comporta como un derivador, ya que ser�
    %un par�metro de dise�o.
    slope_aux = diff( log10(g_amp) );
    end_diff_idx = find(slope_aux < 0.95*slope_aux(1), 1, 'first');

    %dise�o un derivador hasta dicha frecuencia, con una banda de
    %transici�n dada por la expresion , cuando sea posible.
    ftrans = min(70, 251 * exp(-0.63*ii));
    diff_order = round(N*1.8.^((ii*0.4 -0.6)));
    %Los N tienen que ser necesariamente pares para el dise�o del derivador.
    if( rem(diff_order,2) ~= 0 )
        diff_order = diff_order + 1;
    end

    [msgstr, msgid] = lastwarn;
    %itero hasta que se dise�a correctamente.
    jj = 0;
    effective_order = diff_order;
    aux_seq = linspace(0.5,1.1,design_iter_attemps);
    while( jj < design_iter_attemps )
        d = fdesign.differentiator('n,fp,fst', effective_order, F(end_diff_idx), min( 0.95, F(end_diff_idx)+(ftrans*2/fs) ) );
        try
            Hd = design(d,'equiripple');
            [~, msgid] = lastwarn('','');
        catch MException
            %on error, force iteration with different config.
            msgid = 'signal:firpm:DidNotConverge';
        end
        if(strcmpi(msgid,'signal:firpm:DidNotConverge'))
            %fallo en la convergencia, iteramos de nuevo.
            %recorro linealmente por el rango +10:-50%
            effective_order = round(diff_order * aux_seq(jj+1));
            %Los effective_order tienen que ser necesariamente pares para el dise�o del derivador.
            if( rem(effective_order,2) ~= 0 )
                effective_order = effective_order + 1;
            end
            jj = jj + 1;
        else
            jj = design_iter_attemps;
        end
    end

    if(strcmpi(msgid,'signal:firpm:DidNotConverge'))
    %fallo en la convergencia, reportamos el error.
        error('qs_filter_design:Impossible2Design', 'Impossible to design the differentiator filter. Please try another N value, or check the filters transfer functions manually.')
    end

    %Vuelvo a ver la transferencia real del derivador dise�ado y del filtro
    %a emular para que tengan una respuesta similar en la zona en que se
    %comporta como derivador. Averiguo el factor de escala entre ambas
    %transferencias.
    g_amp_emp = freqz(empirical_tf{ii},1, F_log);
    g_amp_emp = abs(g_amp_emp);

    g_amp_real = freqz(Hd, F * pi);
    g_amp_real = abs(g_amp_real);

    aux_idx = round(end_diff_idx/2);
    aux_scale = g_amp_emp(aux_idx)/g_amp_real(aux_idx);

    %Escalo la zona del derivador.
    Hd.Numerator = Hd.Numerator * aux_scale;
    g_amp_real = g_amp_real * aux_scale;

    %ahora dise�o un filtro de compensaci�n para que la zona en que no se
    %deriva sea similar al filtro de referencia. Esta ir� desde el m�ximo
    %de la transferencia pasobanda hasta donde haya una amplitud mayor a
    %0.5 veces.
    [~, max_idx] = max(g_amp);

    aux_3db_unique = find( g_amp_real > 0.5 & g_amp_emp > 0.5, 1, 'last');

    if( max_idx >= aux_3db_unique)
        max_idx = end_diff_idx;
    end
    aux_idx = max_idx:aux_3db_unique;

    %Creo una grilla de frecuencia - magnitud arbitraria para el dise�o del
    %filtro. Esta transferencia no deber� afectar la zona derivadora (H(w)
    %= 1) y compensar la zona indicada por aux_idx, emulando la
    %transferencia de referencia y teniendo un valor de atenuaci�n
    %considerable en Nyquist (H(w) = 1e-3)
    F_comp = [0 max(0.9*F(end_diff_idx), (F(max_idx)-F(end_diff_idx))/2) F(aux_idx) 1];
    g_amp_comp = [1 1 g_amp(aux_idx)./g_amp_real(aux_idx) 1e-3];
    W = ones(1, length(F_comp));
    W(3:end-1) = 10;

    %Se dise�a el filtro
    d = fdesign.arbmag('N,F,A', diff_order, F_comp, g_amp_comp);
    Hd_comp = design(d,'firls', 'weights', W, 'FilterStructure', 'dfsymfir');

    % y se cascadea al original.
    Hd = dfilt.cascade(Hd, Hd_comp);

    q_filters(filter_count) = Hd;

    filter_count = filter_count + 1;

end

author: marianux

repo: ecg-kit

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


text:
function calpak_test615 ( )

%*****************************************************************************80
%
%% CALPAK_TEST615 tests YEAR_LENGTH_REPUBLICAN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST615\n' );
  fprintf ( 1, '  For the Republican calendar:\n' );
  fprintf ( 1, '  YEAR_LENGTH_REPUBLICAN determines the length of a year.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Year  Length\n' );
  fprintf ( 1, '\n' );

  for y = 1 : 6
    sy = y_to_s_republican ( y );
    fprintf ( 1, '  %10s  %3d\n', sy, year_length_republican ( y ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function obj = Jacob_pTop_func_new(in1,in2,in3,in4,in5)
coef_f0_q_sym1 = in1(:,1);
coef_f0_q_sym2 = in1(:,2);
coef_f0_q_sym3 = in1(:,3);
coef_f0_q_sym4 = in1(:,4);
coef_f0_q_sym5 = in1(:,5);
coef_f0_q_sym6 = in1(:,6);
coef_f0_q_sym7 = in1(:,7);
coef_f0_q_sym8 = in1(:,8);
coef_f0_q_sym9 = in1(:,9);
coef_f0_q_sym10 = in1(:,10);
coef_f0_q_sym11 = in1(:,11);
coef_f1_q_sym1 = in2(:,1);
coef_f0_q_sym12 = in1(:,12);
coef_f1_q_sym2 = in2(:,2);
coef_f0_q_sym13 = in1(:,13);
coef_f1_q_sym3 = in2(:,3);
coef_f0_q_sym14 = in1(:,14);
coef_f1_q_sym4 = in2(:,4);
coef_f0_q_sym15 = in1(:,15);
coef_f1_q_sym5 = in2(:,5);
coef_f0_q_sym16 = in1(:,16);
coef_f1_q_sym6 = in2(:,6);
coef_f0_q_sym17 = in1(:,17);
coef_f1_q_sym7 = in2(:,7);
coef_f0_q_sym18 = in1(:,18);
coef_f1_q_sym8 = in2(:,8);
coef_f0_q_sym19 = in1(:,19);
coef_f1_q_sym9 = in2(:,9);
coef_f0_q_sym20 = in1(:,20);
coef_f0_q_sym21 = in1(:,21);
coef_f2_q_sym1 = in3(:,1);
coef_f0_q_sym22 = in1(:,22);
coef_f2_q_sym2 = in3(:,2);
coef_f0_q_sym23 = in1(:,23);
coef_f2_q_sym3 = in3(:,3);
coef_f0_q_sym24 = in1(:,24);
coef_f2_q_sym4 = in3(:,4);
coef_f2_q_sym5 = in3(:,5);
coef_f2_q_sym6 = in3(:,6);
coef_f2_q_sym7 = in3(:,7);
coef_f2_q_sym8 = in3(:,8);
coef_f2_q_sym9 = in3(:,9);
coef_f3_q_sym1 = in4(:,1);
coef_f3_q_sym2 = in4(:,2);
coef_f3_q_sym3 = in4(:,3);
coef_f3_q_sym4 = in4(:,4);
coef_f3_q_sym5 = in4(:,5);
coef_f3_q_sym6 = in4(:,6);
coef_f3_q_sym7 = in4(:,7);
coef_f3_q_sym8 = in4(:,8);
coef_f3_q_sym9 = in4(:,9);
coef_f1_q_sym10 = in2(:,10);
coef_f1_q_sym11 = in2(:,11);
coef_f1_q_sym12 = in2(:,12);
coef_f1_q_sym13 = in2(:,13);
coef_f1_q_sym14 = in2(:,14);
coef_f1_q_sym15 = in2(:,15);
coef_f1_q_sym16 = in2(:,16);
coef_f1_q_sym17 = in2(:,17);
coef_f1_q_sym18 = in2(:,18);
coef_f1_q_sym19 = in2(:,19);
coef_f1_q_sym20 = in2(:,20);
coef_f1_q_sym21 = in2(:,21);
coef_f1_q_sym22 = in2(:,22);
coef_f1_q_sym23 = in2(:,23);
coef_f1_q_sym24 = in2(:,24);
coef_f2_q_sym10 = in3(:,10);
coef_f2_q_sym11 = in3(:,11);
coef_f2_q_sym12 = in3(:,12);
coef_f2_q_sym13 = in3(:,13);
coef_f2_q_sym14 = in3(:,14);
coef_f2_q_sym15 = in3(:,15);
coef_f2_q_sym16 = in3(:,16);
coef_f2_q_sym17 = in3(:,17);
coef_f2_q_sym18 = in3(:,18);
coef_f2_q_sym19 = in3(:,19);
coef_f2_q_sym20 = in3(:,20);
coef_f2_q_sym21 = in3(:,21);
coef_f2_q_sym22 = in3(:,22);
coef_f2_q_sym23 = in3(:,23);
coef_f2_q_sym24 = in3(:,24);
coef_f3_q_sym10 = in4(:,10);
coef_f3_q_sym11 = in4(:,11);
coef_f3_q_sym12 = in4(:,12);
coef_f3_q_sym13 = in4(:,13);
coef_f3_q_sym14 = in4(:,14);
coef_f3_q_sym15 = in4(:,15);
coef_f3_q_sym16 = in4(:,16);
coef_f3_q_sym17 = in4(:,17);
coef_f3_q_sym18 = in4(:,18);
coef_f3_q_sym19 = in4(:,19);
coef_f3_q_sym20 = in4(:,20);
coef_f3_q_sym21 = in4(:,21);
coef_f3_q_sym22 = in4(:,22);
coef_f3_q_sym23 = in4(:,23);
coef_f3_q_sym24 = in4(:,24);
q0 = in5(1,:);
q1 = in5(2,:);
q2 = in5(3,:);
q3 = in5(4,:);
t2 = coef_f0_q_sym11.*q0;
t3 = coef_f0_q_sym18.*q1;
t4 = coef_f0_q_sym22.*q2;
t5 = coef_f0_q_sym24.*q3;
t6 = q0.^2;
t7 = q0.^3;
t8 = q1.^2;
t9 = q1.^3;
t10 = q2.^2;
t11 = q2.^3;
t12 = q3.^2;
t13 = q3.^3;
t24 = coef_f0_q_sym6.*q0.*q1.*q2.*2.0;
t25 = coef_f0_q_sym7.*q0.*q1.*q3.*2.0;
t26 = coef_f0_q_sym9.*q0.*q2.*q3.*2.0;
t27 = coef_f0_q_sym16.*q1.*q2.*q3.*2.0;
t14 = coef_f0_q_sym1.*t7;
t15 = coef_f0_q_sym12.*t9;
t16 = coef_f0_q_sym19.*t11;
t17 = coef_f0_q_sym23.*t13;
t18 = coef_f0_q_sym2.*q1.*t6;
t19 = coef_f0_q_sym3.*q2.*t6;
t20 = coef_f0_q_sym5.*q0.*t8;
t21 = coef_f0_q_sym4.*q3.*t6;
t22 = coef_f0_q_sym8.*q0.*t10;
t23 = coef_f0_q_sym10.*q0.*t12;
obj = reshape([coef_f0_q_sym11.*q1-coef_f1_q_sym11.*q0.*2.0-coef_f1_q_sym18.*q1-coef_f1_q_sym22.*q2-coef_f1_q_sym24.*q3+coef_f0_q_sym5.*t9-coef_f1_q_sym1.*t7.*4.0-coef_f1_q_sym12.*t9-coef_f1_q_sym19.*t11-coef_f1_q_sym23.*t13+coef_f0_q_sym1.*q1.*t6.*3.0+coef_f0_q_sym2.*q0.*t8.*2.0+coef_f0_q_sym6.*q2.*t8+coef_f0_q_sym7.*q3.*t8+coef_f0_q_sym8.*q1.*t10-coef_f1_q_sym2.*q1.*t6.*3.0-coef_f1_q_sym3.*q2.*t6.*3.0+coef_f0_q_sym10.*q1.*t12-coef_f1_q_sym4.*q3.*t6.*3.0-coef_f1_q_sym5.*q0.*t8.*2.0-coef_f1_q_sym8.*q0.*t10.*2.0-coef_f1_q_sym10.*q0.*t12.*2.0-coef_f1_q_sym13.*q2.*t8-coef_f1_q_sym14.*q3.*t8-coef_f1_q_sym15.*q1.*t10-coef_f1_q_sym17.*q1.*t12-coef_f1_q_sym20.*q3.*t10-coef_f1_q_sym21.*q2.*t12+coef_f0_q_sym3.*q0.*q1.*q2.*2.0+coef_f0_q_sym4.*q0.*q1.*q3.*2.0+coef_f0_q_sym9.*q1.*q2.*q3-coef_f1_q_sym6.*q0.*q1.*q2.*2.0-coef_f1_q_sym7.*q0.*q1.*q3.*2.0-coef_f1_q_sym9.*q0.*q2.*q3.*2.0-coef_f1_q_sym16.*q1.*q2.*q3,coef_f0_q_sym11.*q2-coef_f2_q_sym11.*q0.*2.0-coef_f2_q_sym18.*q1-coef_f2_q_sym22.*q2-coef_f2_q_sym24.*q3+coef_f0_q_sym8.*t11-coef_f2_q_sym1.*t7.*4.0-coef_f2_q_sym12.*t9-coef_f2_q_sym19.*t11-coef_f2_q_sym23.*t13+coef_f0_q_sym1.*q2.*t6.*3.0+coef_f0_q_sym3.*q0.*t10.*2.0+coef_f0_q_sym5.*q2.*t8+coef_f0_q_sym6.*q1.*t10+coef_f0_q_sym9.*q3.*t10+coef_f0_q_sym10.*q2.*t12-coef_f2_q_sym2.*q1.*t6.*3.0-coef_f2_q_sym3.*q2.*t6.*3.0-coef_f2_q_sym4.*q3.*t6.*3.0-coef_f2_q_sym5.*q0.*t8.*2.0-coef_f2_q_sym8.*q0.*t10.*2.0-coef_f2_q_sym10.*q0.*t12.*2.0-coef_f2_q_sym13.*q2.*t8-coef_f2_q_sym14.*q3.*t8-coef_f2_q_sym15.*q1.*t10-coef_f2_q_sym17.*q1.*t12-coef_f2_q_sym20.*q3.*t10-coef_f2_q_sym21.*q2.*t12+coef_f0_q_sym2.*q0.*q1.*q2.*2.0+coef_f0_q_sym4.*q0.*q2.*q3.*2.0+coef_f0_q_sym7.*q1.*q2.*q3-coef_f2_q_sym6.*q0.*q1.*q2.*2.0-coef_f2_q_sym7.*q0.*q1.*q3.*2.0-coef_f2_q_sym9.*q0.*q2.*q3.*2.0-coef_f2_q_sym16.*q1.*q2.*q3,coef_f0_q_sym11.*q3-coef_f3_q_sym11.*q0.*2.0-coef_f3_q_sym18.*q1-coef_f3_q_sym22.*q2-coef_f3_q_sym24.*q3+coef_f0_q_sym10.*t13-coef_f3_q_sym1.*t7.*4.0-coef_f3_q_sym12.*t9-coef_f3_q_sym19.*t11-coef_f3_q_sym23.*t13+coef_f0_q_sym1.*q3.*t6.*3.0+coef_f0_q_sym4.*q0.*t12.*2.0+coef_f0_q_sym5.*q3.*t8+coef_f0_q_sym7.*q1.*t12+coef_f0_q_sym8.*q3.*t10+coef_f0_q_sym9.*q2.*t12-coef_f3_q_sym2.*q1.*t6.*3.0-coef_f3_q_sym3.*q2.*t6.*3.0-coef_f3_q_sym4.*q3.*t6.*3.0-coef_f3_q_sym5.*q0.*t8.*2.0-coef_f3_q_sym8.*q0.*t10.*2.0-coef_f3_q_sym10.*q0.*t12.*2.0-coef_f3_q_sym13.*q2.*t8-coef_f3_q_sym14.*q3.*t8-coef_f3_q_sym15.*q1.*t10-coef_f3_q_sym17.*q1.*t12-coef_f3_q_sym20.*q3.*t10-coef_f3_q_sym21.*q2.*t12+coef_f0_q_sym2.*q0.*q1.*q3.*2.0+coef_f0_q_sym3.*q0.*q2.*q3.*2.0+coef_f0_q_sym6.*q1.*q2.*q3-coef_f3_q_sym6.*q0.*q1.*q2.*2.0-coef_f3_q_sym7.*q0.*q1.*q3.*2.0-coef_f3_q_sym9.*q0.*q2.*q3.*2.0-coef_f3_q_sym16.*q1.*q2.*q3,q0.*2.0,t2+t3.*2.0+t4+t5+t14+t15.*4.0+t16+t17+t18.*2.0+t19+t20.*3.0+t21+t22+t23+t24+t25+t27-coef_f1_q_sym18.*q0-coef_f1_q_sym2.*t7-coef_f1_q_sym5.*q1.*t6.*2.0+coef_f0_q_sym13.*q2.*t8.*3.0-coef_f1_q_sym6.*q2.*t6+coef_f0_q_sym14.*q3.*t8.*3.0+coef_f0_q_sym15.*q1.*t10.*2.0-coef_f1_q_sym7.*q3.*t6+coef_f0_q_sym17.*q1.*t12.*2.0+coef_f0_q_sym20.*q3.*t10+coef_f0_q_sym21.*q2.*t12-coef_f1_q_sym12.*q0.*t8.*3.0-coef_f1_q_sym15.*q0.*t10-coef_f1_q_sym17.*q0.*t12+coef_f0_q_sym9.*q0.*q2.*q3-coef_f1_q_sym13.*q0.*q1.*q2.*2.0-coef_f1_q_sym14.*q0.*q1.*q3.*2.0-coef_f1_q_sym16.*q0.*q2.*q3,coef_f0_q_sym18.*q2-coef_f2_q_sym18.*q0+coef_f0_q_sym15.*t11-coef_f2_q_sym2.*t7+coef_f0_q_sym2.*q2.*t6+coef_f0_q_sym6.*q0.*t10+coef_f0_q_sym12.*q2.*t8.*3.0+coef_f0_q_sym13.*q1.*t10.*2.0+coef_f0_q_sym16.*q3.*t10+coef_f0_q_sym17.*q2.*t12-coef_f2_q_sym5.*q1.*t6.*2.0-coef_f2_q_sym6.*q2.*t6-coef_f2_q_sym7.*q3.*t6-coef_f2_q_sym12.*q0.*t8.*3.0-coef_f2_q_sym15.*q0.*t10-coef_f2_q_sym17.*q0.*t12+coef_f0_q_sym5.*q0.*q1.*q2.*2.0+coef_f0_q_sym7.*q0.*q2.*q3+coef_f0_q_sym14.*q1.*q2.*q3.*2.0-coef_f2_q_sym13.*q0.*q1.*q2.*2.0-coef_f2_q_sym14.*q0.*q1.*q3.*2.0-coef_f2_q_sym16.*q0.*q2.*q3,coef_f0_q_sym18.*q3-coef_f3_q_sym18.*q0+coef_f0_q_sym17.*t13-coef_f3_q_sym2.*t7+coef_f0_q_sym2.*q3.*t6+coef_f0_q_sym7.*q0.*t12+coef_f0_q_sym12.*q3.*t8.*3.0+coef_f0_q_sym14.*q1.*t12.*2.0+coef_f0_q_sym15.*q3.*t10+coef_f0_q_sym16.*q2.*t12-coef_f3_q_sym5.*q1.*t6.*2.0-coef_f3_q_sym6.*q2.*t6-coef_f3_q_sym7.*q3.*t6-coef_f3_q_sym12.*q0.*t8.*3.0-coef_f3_q_sym15.*q0.*t10-coef_f3_q_sym17.*q0.*t12+coef_f0_q_sym5.*q0.*q1.*q3.*2.0+coef_f0_q_sym6.*q0.*q2.*q3+coef_f0_q_sym13.*q1.*q2.*q3.*2.0-coef_f3_q_sym13.*q0.*q1.*q2.*2.0-coef_f3_q_sym14.*q0.*q1.*q3.*2.0-coef_f3_q_sym16.*q0.*q2.*q3,q1.*2.0,coef_f0_q_sym22.*q1-coef_f1_q_sym22.*q0-coef_f1_q_sym3.*t7+coef_f0_q_sym13.*t9+coef_f0_q_sym3.*q1.*t6+coef_f0_q_sym6.*q0.*t8-coef_f1_q_sym6.*q1.*t6+coef_f0_q_sym15.*q2.*t8.*2.0-coef_f1_q_sym8.*q2.*t6.*2.0+coef_f0_q_sym16.*q3.*t8-coef_f1_q_sym9.*q3.*t6+coef_f0_q_sym19.*q1.*t10.*3.0+coef_f0_q_sym21.*q1.*t12-coef_f1_q_sym13.*q0.*t8-coef_f1_q_sym19.*q0.*t10.*3.0-coef_f1_q_sym21.*q0.*t12+coef_f0_q_sym8.*q0.*q1.*q2.*2.0+coef_f0_q_sym9.*q0.*q1.*q3+coef_f0_q_sym20.*q1.*q2.*q3.*2.0-coef_f1_q_sym15.*q0.*q1.*q2.*2.0-coef_f1_q_sym16.*q0.*q1.*q3-coef_f1_q_sym20.*q0.*q2.*q3.*2.0,t2+t3+t4.*2.0+t5+t14+t15+t16.*4.0+t17+t18+t19.*2.0+t20+t21+t22.*3.0+t23+t24+t26+t27-coef_f2_q_sym22.*q0-coef_f2_q_sym3.*t7+coef_f0_q_sym13.*q2.*t8.*2.0+coef_f0_q_sym14.*q3.*t8+coef_f0_q_sym15.*q1.*t10.*3.0+coef_f0_q_sym17.*q1.*t12+coef_f0_q_sym20.*q3.*t10.*3.0-coef_f2_q_sym6.*q1.*t6+coef_f0_q_sym21.*q2.*t12.*2.0-coef_f2_q_sym8.*q2.*t6.*2.0-coef_f2_q_sym9.*q3.*t6-coef_f2_q_sym13.*q0.*t8-coef_f2_q_sym19.*q0.*t10.*3.0-coef_f2_q_sym21.*q0.*t12+coef_f0_q_sym7.*q0.*q1.*q3-coef_f2_q_sym15.*q0.*q1.*q2.*2.0-coef_f2_q_sym16.*q0.*q1.*q3-coef_f2_q_sym20.*q0.*q2.*q3.*2.0,coef_f0_q_sym22.*q3-coef_f3_q_sym22.*q0+coef_f0_q_sym21.*t13-coef_f3_q_sym3.*t7+coef_f0_q_sym3.*q3.*t6+coef_f0_q_sym9.*q0.*t12+coef_f0_q_sym13.*q3.*t8+coef_f0_q_sym16.*q1.*t12+coef_f0_q_sym19.*q3.*t10.*3.0+coef_f0_q_sym20.*q2.*t12.*2.0-coef_f3_q_sym6.*q1.*t6-coef_f3_q_sym8.*q2.*t6.*2.0-coef_f3_q_sym9.*q3.*t6-coef_f3_q_sym13.*q0.*t8-coef_f3_q_sym19.*q0.*t10.*3.0-coef_f3_q_sym21.*q0.*t12+coef_f0_q_sym6.*q0.*q1.*q3+coef_f0_q_sym8.*q0.*q2.*q3.*2.0+coef_f0_q_sym15.*q1.*q2.*q3.*2.0-coef_f3_q_sym15.*q0.*q1.*q2.*2.0-coef_f3_q_sym16.*q0.*q1.*q3-coef_f3_q_sym20.*q0.*q2.*q3.*2.0,q2.*2.0,coef_f0_q_sym24.*q1-coef_f1_q_sym24.*q0-coef_f1_q_sym4.*t7+coef_f0_q_sym14.*t9+coef_f0_q_sym4.*q1.*t6+coef_f0_q_sym7.*q0.*t8-coef_f1_q_sym7.*q1.*t6+coef_f0_q_sym16.*q2.*t8-coef_f1_q_sym9.*q2.*t6+coef_f0_q_sym17.*q3.*t8.*2.0+coef_f0_q_sym20.*q1.*t10+coef_f0_q_sym23.*q1.*t12.*3.0-coef_f1_q_sym10.*q3.*t6.*2.0-coef_f1_q_sym14.*q0.*t8-coef_f1_q_sym20.*q0.*t10-coef_f1_q_sym23.*q0.*t12.*3.0+coef_f0_q_sym9.*q0.*q1.*q2+coef_f0_q_sym10.*q0.*q1.*q3.*2.0+coef_f0_q_sym21.*q1.*q2.*q3.*2.0-coef_f1_q_sym16.*q0.*q1.*q2-coef_f1_q_sym17.*q0.*q1.*q3.*2.0-coef_f1_q_sym21.*q0.*q2.*q3.*2.0,coef_f0_q_sym24.*q2-coef_f2_q_sym24.*q0+coef_f0_q_sym20.*t11-coef_f2_q_sym4.*t7+coef_f0_q_sym4.*q2.*t6+coef_f0_q_sym9.*q0.*t10+coef_f0_q_sym14.*q2.*t8+coef_f0_q_sym16.*q1.*t10+coef_f0_q_sym21.*q3.*t10.*2.0-coef_f2_q_sym7.*q1.*t6+coef_f0_q_sym23.*q2.*t12.*3.0-coef_f2_q_sym9.*q2.*t6-coef_f2_q_sym10.*q3.*t6.*2.0-coef_f2_q_sym14.*q0.*t8-coef_f2_q_sym20.*q0.*t10-coef_f2_q_sym23.*q0.*t12.*3.0+coef_f0_q_sym7.*q0.*q1.*q2+coef_f0_q_sym10.*q0.*q2.*q3.*2.0+coef_f0_q_sym17.*q1.*q2.*q3.*2.0-coef_f2_q_sym16.*q0.*q1.*q2-coef_f2_q_sym17.*q0.*q1.*q3.*2.0-coef_f2_q_sym21.*q0.*q2.*q3.*2.0,t2+t3+t4+t5.*2.0+t14+t15+t16+t17.*4.0+t18+t19+t20+t21.*2.0+t22+t23.*3.0+t25+t26+t27-coef_f3_q_sym24.*q0-coef_f3_q_sym4.*t7+coef_f0_q_sym13.*q2.*t8+coef_f0_q_sym14.*q3.*t8.*2.0+coef_f0_q_sym15.*q1.*t10+coef_f0_q_sym17.*q1.*t12.*3.0+coef_f0_q_sym20.*q3.*t10.*2.0+coef_f0_q_sym21.*q2.*t12.*3.0-coef_f3_q_sym7.*q1.*t6-coef_f3_q_sym9.*q2.*t6-coef_f3_q_sym10.*q3.*t6.*2.0-coef_f3_q_sym14.*q0.*t8-coef_f3_q_sym20.*q0.*t10-coef_f3_q_sym23.*q0.*t12.*3.0+coef_f0_q_sym6.*q0.*q1.*q2-coef_f3_q_sym16.*q0.*q1.*q2-coef_f3_q_sym17.*q0.*q1.*q3.*2.0-coef_f3_q_sym21.*q0.*q2.*q3.*2.0,q3.*2.0],[4,4]);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [flatView, allDist] = measureCorticalDistanceLineROI(flatView, coords,plotFlag)
%
% USAGE: measureCorticalDistanceLineROI(flatView, coords)
%
% AUTHOR:  Dougherty
% DATE:    2002.04.10
% PURPOSE:
%   Compute the shortest cortical manifold distance between
%   points. The coords is passed in are assume to be a line ROI
%   where the first coord is at one end of the line and the last
%   coord is at the other end.
%
% HISTORY
%
% 7/16/02 djh, replaced mrSESSION.vAnatomyPath with global vANATOMYPATH

sampleDist = 5;

global vANATOMYPATH;
mmPerPix = readVolAnatHeader(vANATOMYPATH);

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

% Get a gray structure because we need the gray nodes.
grayView = getSelectedGray;
if isempty(grayView)
    grayView = initHiddenGray;
end

% the third coordinate is the 'slice', which, for flat views, means left or right hemisphere.
slice = coords(3,1);
if (slice==1)
    nodes = grayView.allLeftNodes;
    edges = grayView.allLeftEdges;
else
    nodes = grayView.allRightNodes;
    edges = grayView.allRightEdges;
end

nCoords = size(coords,2);
disp(['Finding nodes for ',num2str(nCoords),' coords...']);
allNodeIndices = zeros(1,nCoords);
% Find the nearest gray node for each coordinate.
%
for(ii=1:nCoords)
    % get nearest flat coordinate (not all points on the flat correspond to flat coordinates)
    flatDistances = (flatView.coords{slice}(1,:) - coords(1,ii)).^2 + ...
        (flatView.coords{slice}(2,:) - coords(2,ii)).^2;
    % There is a one-to-many mapping of flatCoords to grayCoords, but we ignore that
    % here by using 'min', which will always reuturn one value, even if there are several
    % identical minima.
    % FIX THIS- we should always grab layer 1, or something more consistent
    % than relying on min's arbitrary sort.
    [val,coordIndex] = min(flatDistances);

    grayNode = find(nodes(2,:) == flatView.grayCoords{slice}(1,coordIndex) & ...
            nodes(1,:) == flatView.grayCoords{slice}(2,coordIndex) & ...
            nodes(3,:) == flatView.grayCoords{slice}(3,coordIndex));
    % This should produce exactly one index.

    % Catch errors.
    if(isempty(grayNode))
        myErrorDlg('No gray nodes were found!');
    end
    if(length(grayNode)>1)
        disp([mfilename,': WARNING- coord ',num2str(ii),'- more than one grayNode found!']);
        grayNode = grayNode(1);
    end
    allNodeIndices(ii) = grayNode;
end

sampleNodes = allNodeIndices(1);
nodeIndices = allNodeIndices(2:end);
done = 0;
while(~done)
    % Now, compute the manifold distance to all other points from the given 'start' point.
    allDist = mrManDist(nodes, edges, sampleNodes(end), mmPerPix, -1, 0);
    drop = intersect(find(allDist<=sampleDist), nodeIndices);
    if(~isempty(drop))
        for(ii=1:length(drop))
            nodeIndices = nodeIndices(nodeIndices~=drop(ii));
        end
    end
    if(~isempty(nodeIndices))
        % *** WE SHOULD GET THE NEXT NEAREST
        nearest = find(allDist(nodeIndices)==min(allDist(nodeIndices)));
        nearest = nearest(1);
        sampleNodes = [sampleNodes, nodeIndices(nearest)];
        nodeIndices = nodeIndices(nodeIndices~=nodeIndices(nearest));
    end
    if(isempty(nodeIndices))
        done = 1;
    end
end

% This would be much simpler if we kept track of which coords we were sampling in the loop above.
for(ii=1:length(sampleNodes))
    thisOne = find(allNodeIndices==sampleNodes(ii));
    sampleCoords(:,ii) = coords(:, thisOne(1));
end

% We loop for the number of line segements, which is the number of coords - 1.
for(ii=1:size(sampleCoords,2)-1)
    % Draw a line for each measured segment
    h(ii) = line(sampleCoords(2,ii:ii+1), sampleCoords(1,ii:ii+1), 'Color', 'r', 'LineWidth', 2);

    % Now, compute the manifold distance between these points.
    % mrManDist returns the distance to all other points from the given 'start' point.
    allDist = mrManDist(nodes, edges, sampleNodes(ii), mmPerPix, -1, 0);
    % We just want the distance from the start point to the end point, so we
    % pull that out by providing the index of the end point.
    dist(ii) = allDist(sampleNodes(ii+1));
    disp(['Cortical distance of segment ',num2str(ii),': ',num2str(dist(ii)),' mm.']);
end
disp(['Total cortical distance: ',num2str(sum(dist)),' mm.']);
if (plotFlag)
uiwait(msgbox(['Total cortical distance: ',num2str(sum(dist)),' mm.'], ...
        'Cortical Distance', 'modal'));
end

    for ii=1:length(h)
    delete(h(ii))
end
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function quenchWriteStateFromFiberGroup(fg,t1File,outName)
%
%   quenchWriteStateFromFiberGroup([fg],[t1File],[outName])
%
% OVERVIEW: Takes a mrDiffusion fiber group and associated t1File and
%           creates a Quench '.qst' state file using the fg.subgroup field
%           to identify different group assignments in Quench.
%
% NOTE:
%           This function saves an accompanying pathway database file
%           (.pdb) that will be loaded in quench when the state file is
%           loaded. The T1 file can be stored anywhere as the full path is
%           written to the sate file - the pdb file and the qst files are
%           saved, and should be stored, in the same directory.
%
% INPUTS:
%       fg      -  a mrDiffusion fiber structure: can be a filename (.mat)
%                  with a subgroup field wherein each distinct fiber is
%                  given an idex unique to the fibergroup it belongs to.
%                  The field fg.subgroupNames should have the names for
%                  each of the unique fiber groups indexed.
%
%       t1File  -  full path to t1 file. Important: Do not use symbolic
%                  links. If you use a symbolic link in linux then we can
%                  resovle the link, but in windows or MAC you're on your
%                  own.
%
%       outName -  the stem for the qst and pdb filenames. If not provided
%                  the name will be taken from fg.name.
%
% WEB RESOURCES:
%       mrvBrowseSVN('quenchWriteStateFromFiberGroup');
%
%
% EXAMPLE USAGE:
%
%   Default usage: Provide path to fiber group and t1 file - take name
%                  from fiber group.
%       >> fg = 'MoriGroups.mat';
%       >> t1File = '/some/Path/t1File.nii.gz';
%       >> quenchWriteStateFromFiberGroup(fg,t1File,[]);
%
% (C) Stanford University, VISTA Lab 2012
%

%% Check inputs and open file for writing

% Fiber group
if notDefined('fg')
    fg = mrvSelectFile('r','*.mat','QST: Select fiber group',pwd);
    if isempty(fg); return; end
end

% Check that fg is a struct or a file that we can read in
if ~isstruct(fg) && exist(fg,'file')
    fprintf('Reading fiber group structure from file...');
    fg = fgRead(fg);
    fprintf('done.\n')
elseif isfield(fg,'fg.fibers')
    disp('Fiber group structure already loaded.');
end

% Check for subroups and subgroup names fields. If they don't exist however
% the fiber group structure contains an array of fiber groups then convert
% it into the desired format. Otherwise exit the function
if ~isfield(fg,'subgroup') && ~isfield(fg,'subgroupNames') && length(fg) > 1
    fg = dtiFgArrayToFiberGroup(fg);
elseif ~isfield(fg,'subgroup') && length(fg) <=1
    error('No subgroup field found in fibergroup structure!');
elseif ~isfield(fg,'subgroupNames') && length(fg) <=1
    error('No subgroupNames field found in fiber structure!');
end

% T1 File
if notDefined('t1File')
    t1File = mrvSelectFile('r','*.nii*','QST: Select t1File',pwd);
    if isempty(t1File); return; end
else
    [p f e] = fileparts(t1File);
    if isempty(p)
        disp('WARNING: t1File must be a full path to the t1File - NOT A SYMBOLIC LINK!!!!');
        t1File = mrvSelectFile(r,'.nii.gz','Choose t1File',pwd);
        if isempty(t1File); return; end
    end
end

% Check that t1File is not a soft link:
if isunix
    [status t1File] = system((sprintf('readlink -f %s',t1File)));
    if status ~= 0
        disp('Could not resolve real path to t1File');
    end
else
    disp('NON-UNIX system detected: Could not verify the path of the t1File - make sure that t1File is not a symbolic Link.')
end

% If outname is not supplied then we use the fg.name field and make sure
% that the extension is not yet attached
if notDefined('outName')
    outName = fg.name;
end

[path file ext] = fileparts(outName);

if ~isempty(ext)
    outName = fullfile(path,file);
end

% Remove any spaces in the file name and replace with '_'
outName = regexprep(outName,'\W','_');

%% Set color and name assignments

% Get the information we need to do the assignments
fgInds    = fg.subgroup;
numGroups = numel(unique(fgInds));

disp([num2str(numGroups) ' fiber groups will be written to the state file:']);

% Number of unique colors used by default in Quench = 16
num       = 1:16;
cinds     = repmat(num,1,(round(numGroups/16)+1));

% Build cell array of colors for color assignments
c{1}='0.0784314 0.352941 0.784314 1';
c{2}='0.596078 0.305882 0.639216 1';
c{3}='1 1 0.2 1';
c{4}='0.301961 0.686275 0.290196 1';
c{5}='0.745098 0.156863 0.156863 1';
c{6}='1 0.498039 0 1';
c{7}='0.470588 0.705882 0.705882 1';
c{8}='0.470588 0.392157 0.196078 1';
c{9}='0.564706 0.705882 1 1';
c{10}='0.831373 0.682353 0.854902 1';
c{11}='0.964706 0.964706 0.690196 1';
c{12}='0.678431 0.87451 0.670588 1';
c{13}='0.913725 0.607843 0.607843 1';
c{14}='0.976471 0.784314 0.611765 1';
c{15}='0.756863 0.882353 0.882353 1';
c{16}='0.780392 0.737255 0.635294 1';

% Loop over the groups and assign the color
for ii = 1:numGroups
    color{ii} = c{cinds(ii)};
end

% Name assignments: Check group names field and build group name cell array
% for displaying each of the names
for ii = 1:numGroups
    if numGroups == length(fg.subgroupNames)
        gname{ii} = fg.subgroupNames(ii).subgroupName;
        % Remove any spaces in the file names and replace with '_'
        gname{ii} = regexprep(gname{ii},'\W','_');
        fprintf(' %s - %s\n',num2str(ii),gname{ii});
    else
        error('Length of subgroupNames is not == to the number of subgroups');
    end
end

%% Start writing the huge, annoying text file
%  Note that quench is INCREDIBLY picky about line endings and spaces.
%  Things have to be just as they are below or Quench will CRASH - hard.

fprintf('Writing Quench state file [%s.qst] ... ', outName);

% Open the state file for writing
qst = fopen([outName '.qst'],'wt');

        fprintf(qst,'Version 2\n');

        fprintf(qst,'--- Main Window ---\n');
        fprintf(qst,'Position 1865 52\n');
        fprintf(qst,'Size 1088 915\n');
        fprintf(qst,'Is Maximized 0\n');
        fprintf(qst,'Current background image index 0\n');

        fprintf(qst,'\n--- Volumes ---\n');
        fprintf(qst,'Num Volumes :1\n');
        fprintf(qst,'%s\n', t1File); % This has to be the full path to the t1.

        fprintf(qst,'--- PDB Info ---\n');
        fprintf(qst,'Num PDB''s: 1\n');
        fprintf(qst,'%s\n',[outName '.pdb']); % This could be the full path to the fiber group

        fprintf(qst,'--- Pathway Assignment ---\n');
        fprintf(qst,'Locked: 0\n');
        fprintf(qst,'Selected Group: 1\n');
        fprintf(qst,'Num Assigned: %s\n',num2str(numel(fgInds)));

        % Convert the fgInds matrix to a string and remove the brackets
        indstr      = mat2str(fgInds);
        indstr(1)   = [];
        indstr(end) = [];
        fprintf(qst,'%s \n',indstr); % the trailing space is important!!!

        fprintf(qst,'\n--- Pathway Groups ---\n');
        fprintf(qst,'Num Groups: %s\n',num2str(numGroups+1)); % This should be +1
        fprintf(qst,'Name: Trash\n');
        fprintf(qst,'Color: 0.501961 0.501961 0.501961 1\n');
        fprintf(qst,'Visible: 1\n');
        fprintf(qst,'Active: 1\n\n');

        % Loop over the groups and write out the unique names and colors
        for ii = 1:numGroups
            fprintf(qst,'Name: %s\n', gname{ii});
            fprintf(qst,'Color: %s\n', color{ii});
            fprintf(qst,'Visible: 1\n');
            fprintf(qst,'Active: 1\n\n');
        end

        % Not sure that these positions will work for every screen
        fprintf(qst,'Camera Position: -503.393,-16.364,26.7371\n');
        fprintf(qst,'Camera View Up: 0.0488575,-0.000488639,0.998806\n');
        fprintf(qst,'Camera Focal Point: 0.062871,-13.8796,2.11129\n');

        fprintf(qst,'--- Move tool ---\n');
        fprintf(qst,'Position: (0,0,0)\n');
        fprintf(qst,'Scale :(19.8361,19.8361,19.8361)\n');
        fprintf(qst,'Visible : 0\n');

        fprintf(qst,'--- Scale tool ---\n');
        fprintf(qst,'Position: (0,0,0)\n');
        fprintf(qst,'Scale :(19.8361,19.8361,19.8361)\n');
        fprintf(qst,'Visible : 0\n');

        fprintf(qst,'\n--- Volume State ---\n');
        fprintf(qst,'Volume Section"s Position: 0,-18,18\n');
        fprintf(qst,'Volume Section"s Visibility: 1,1,1\n');
        fprintf(qst,'Active Image :1\n');

        fprintf(qst,'\n--- Overlays ---\n');
        fprintf(qst,'Num overlay items :1\n');
        fprintf(qst,'Opacity :1\n');
        fprintf(qst,'Range :0 32767\n');
        fprintf(qst,'ColorMapIndex :0\n');
        fprintf(qst,'Visible :0\n');

        fprintf(qst,'\n--- Gesture Panel ---\n');
        fprintf(qst,'Visible :1\n');
        fprintf(qst,'Position :(755,4)\n');
        fprintf(qst,'SelectMode :0\n');

        fprintf(qst,'\n--- Selection Panel ---\n');
        fprintf(qst,'Visible :1\n');
        fprintf(qst,'Position :(857,4)\n');
        fprintf(qst,'SelectMode :0\n');

        fprintf(qst,'\n--- Refine Selection panel ---\n');
        fprintf(qst,'Position -1 -1\n');
        fprintf(qst,'Size 450 260\n');
        fprintf(qst,'Visible 0\n');
        fprintf(qst,'Num ROIs 0\n');

        fprintf(qst,'\n--- Background and overlay panel ---\n');
        fprintf(qst,'Position -1 -1\n');
        fprintf(qst,'Size 458 250\n');
        fprintf(qst,'Visible 0\n');

        fprintf(qst,'\n--- ROIs ---\n');
        fprintf(qst,'Num ROIs: 0\n\n');

fclose(qst);

fprintf('done.\n');

%% Save the fiber group struct as a pdb

fprintf('Writing pathway database [%s.pdb] ... ', outName);

fgWrite(fg,outName,'pdb')

fprintf('done! \n');

return

%%
%#ok<*ASGLU>
%#ok<*AGROW>
%#ok<*NASGU>


text:
% demo.m

close all;
clear;
clc;

%% set parameters
k = 5;

%% load data
load('./dataset/ORL_Face_img_cov.mat');

%% perform RCM k-NN classifier with
% GRCM2 with eigenvalue-based distance
grcm_accuracy = rcm_knn_classifier(TrainSet, TestSet,'GRCM', '2', 'EV', k);
% RCM4 with eigenvalue-based distance
rcm_accuracy = rcm_knn_classifier(TrainSet, TestSet, 'RCM', '4', 'EV', k);

fprintf('\n');
fprintf('# GRCM2 Accuracy = %5.2f\n', grcm_accuracy);
fprintf('# RCM4 Accuracy = %5.2f\n', rcm_accuracy);

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function value = r8_beta_pdf ( alpha, beta, rval )

%*****************************************************************************80
%
%% R8_BETA_PDF evaluates the PDF of a beta distribution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2013
%
%  Author:
%
%    Original FORTRAN90 version by Guannan Zhang.
%    MATLAB version by John Burkardt.
%
%  Parameters:
%
%    Input, real ALPHA, BETA, shape parameters.
%    0.0 < ALPHA, BETA.
%
%    Input, real RVAL, the point where the PDF is evaluated.
%
%    Output, real VALUE, the value of the PDF at RVAL.
%
  if ( alpha <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_BETA_PDF - Fatal error!\n' );
    fprintf ( 1, '  Parameter ALPHA is not positive.\n' );
    error ( 'R8_BETA_PDF - Fatal error!' );
  end

  if ( beta <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_BETA_PDF - Fatal error!\n' );
    fprintf ( 1, '  Parameter BETA is not positive.\n' );
    error ( 'R8_BETA_PDF - Fatal error!' );
  end

  if ( rval < 0.0 || 1.0 < rval )

    value = 0.0;

  else

    temp = r8_gamma_log ( alpha + beta ) - r8_gamma_log ( alpha ) ...
      - r8_gamma_log ( beta );

    value = exp ( temp ) * rval ^ ( alpha - 1.0 ) ...
      * ( 1.0 - rval ) ^ ( beta - 1.0 );

  end

  return
end


text:
function handle = imageVisualiseVargplvm(imageVals, imageSize, transpose, negative, ...
				 scale, thresh)

% IMAGEVISUALISE Helper code for showing an image during 2-D visualisation.
% FORMAT
% DESC is a helper function for plotting image data using latent
% variable models.
% ARG imageValues : the values to set the image data to.
% ARG imageSize : the size of the image.
% ARG transpose : whether the resized image needs to be transposed
% (default 1, which is yes).
% ARG negative : whether to display the negative of the image
% (default 0, which is no).
% ARG scale : whether or not to use the imagesc function (defaults
% to 1, which is yes).
% ARG thresh : An array, thresh = [minVal, threshDown, threshUp, maxVal]
% which gives value minVal to all image elements smaller than threshDown
% and value maxVal to all image elements larger than threshUp.
% RETURN handle : a the handle to the image data.
%
% COPYRIGHT : Neil D. Lawrence, 2003, 2004, 2006
%
% SEEALSO : imageModify, lvmResultsDynamic

% SHEFFIELDML

if nargin < 3 || isempty(transpose)
  transpose = 1;
end
if nargin< 4 || isempty(negative)
  negative = 0;
end
if nargin < 5 || isempty(scale)
  scale = 1;
end
if nargin < 6
    thresh = [];
end
if negative
  imageVals = -imageVals;
end
imageData = reshape(imageVals, imageSize(1), imageSize(2));
if transpose
  imageData = imageData';
end
if ~isempty(thresh)
    imageData(imageData < thresh(2)) = thresh(1);
    imageData(imageData >= thresh(3)) = thresh(4);
end
if scale
  handle = imagesc(imageData);
else
  handle = image(imageData);
end
colormap gray

author: SheffieldML

repo: vargplvm


text:
%% This file is the PDE simulation file of SINDy-PI discovered Belousov-Zhabotinsky reaction.
% Coded By: K
% Last Updated: 2019/09/19
%%
function [rhs,x_t,z_t,s_t,u_t,x,z,s,u,x_x,z_x,s_x,u_x,x_y,z_y,s_y,u_y,x_xx,z_xx,s_xx,u_xx,x_yy,z_yy,s_yy,u_yy,x_lap,z_lap,s_lap,u_lap]=...
    BZ_Reaction_SINDy_PI_PDE(t,xzsut,Kx,Kxx,Ky,Kyy,K22,n,N,NeedDev)

% Calculate u and v terms
xt=reshape((xzsut(1:N)),n,n);
zt=reshape((xzsut((N+1):(2*N))),n,n);
st=reshape((xzsut(2*N+1:3*N)),n,n);
ut=reshape((xzsut((3*N+1):(4*N))),n,n);

x=real(ifft2(xt));
z=real(ifft2(zt));
s=real(ifft2(st));
u=real(ifft2(ut));

% Reaction Terms (You need to manully code up the equation discovered by DL-SINDy)
xtrhs=reshape((fft2(  (0.24667*x + 0.33333*s + 0.5*z + 3.3333*x.*s - 5.0*x.*z + 2.1333*x.^2 - 3.3333*x.^3)./(x + 0.1)  )),N,1);
ztrhs=reshape((fft2(  x + 0.4*u - 1.3*z  )),N,1);
strhs=reshape((fft2(  0.17333*x-0.66667*s  )),N,1);
utrhs=reshape((fft2(  100*z-133.33*u  )),N,1);

rhs=[-K22.*xzsut(1:N)+xtrhs
     -0.1*K22.*xzsut(N+1:2*N)+ztrhs
     -K22.*xzsut(2*N+1:3*N)+strhs
     -K22.*xzsut(3*N+1:4*N)+utrhs
     ];

 % If you don't need to extract the value, don't run the following code to
 % speed up
if NeedDev==1
    % Get the derivative you want
    x_x=real(ifft2(reshape((1j*Kx).*xzsut(1:N),n,n)));
    z_x=real(ifft2(reshape((1j*Kx).*xzsut(N+1:2*N),n,n)));
    s_x=real(ifft2(reshape((1j*Kx).*xzsut(2*N+1:3*N),n,n)));
    u_x=real(ifft2(reshape((1j*Kx).*xzsut(3*N+1:4*N),n,n)));
    %
    x_y=real(ifft2(reshape((1j*Ky).*xzsut(1:N),n,n)));
    z_y=real(ifft2(reshape((1j*Ky).*xzsut(N+1:2*N),n,n)));
    s_y=real(ifft2(reshape((1j*Ky).*xzsut(2*N+1:3*N),n,n)));
    u_y=real(ifft2(reshape((1j*Ky).*xzsut(3*N+1:4*N),n,n)));
    %
    x_xx=real(ifft2(reshape((1j*Kxx).^2.*xzsut(1:N),n,n)));
    z_xx=real(ifft2(reshape((1j*Kxx).^2.*xzsut(N+1:2*N),n,n)));
    s_xx=real(ifft2(reshape((1j*Kxx).^2.*xzsut(2*N+1:3*N),n,n)));
    u_xx=real(ifft2(reshape((1j*Kxx).^2.*xzsut(3*N+1:4*N),n,n)));
    %
    x_yy=real(ifft2(reshape((1j*Kyy).^2.*xzsut(1:N),n,n)));
    z_yy=real(ifft2(reshape((1j*Kyy).^2.*xzsut(N+1:2*N),n,n)));
    s_yy=real(ifft2(reshape((1j*Kyy).^2.*xzsut(2*N+1:3*N),n,n)));
    u_yy=real(ifft2(reshape((1j*Kyy).^2.*xzsut(3*N+1:4*N),n,n)));
    %
    x_t=real(ifft2(reshape(rhs(1:N),n,n)));
    z_t=real(ifft2(reshape(rhs(N+1:2*N),n,n)));
    s_t=real(ifft2(reshape(rhs(2*N+1:3*N),n,n)));
    u_t=real(ifft2(reshape(rhs(3*N+1:4*N),n,n)));
    %
    x_lap=real(ifft2(reshape(-K22.*xzsut(1:N),n,n)));
    z_lap=real(ifft2(reshape(-K22.*xzsut(N+1:2*N),n,n)));
    s_lap=real(ifft2(reshape(-K22.*xzsut(2*N+1:3*N),n,n)));
    u_lap=real(ifft2(reshape(-K22.*xzsut(3*N+1:4*N),n,n)));

end

author: dynamicslab

repo: SINDy-PI


text:
function [formula, protons] = getFormulaFromInChI(InChI)
% Extracts the chemical formula of a given compound from
% the InChI string provided
%
% USAGE:
%
%    [formula, protons] = getFormulaFromInChI(InChI)
%
% INPUT:
%    InChI:      The Inchi String of the chemical formula (e.g. InChI=
%                extract formula from `InChI = 1S/C3H4O3/c1-2(4)3(5)6/h1H3, (H,5,6)/p-1` for pyruvate
%
% OUTPUTS:
%    formula:    The chemical formula (including the protonation state
%    protons:    The total number of protons

[token,rem] = strtok(InChI, '/');
formula=strtok(rem, '/');

%This could be a composite formula, so combine it.
tokens = strsplit(formula,'.');

%The protonation state can also modify the formula! To get it, we remove
%any reconstruction fields, as they do not influence it.
InChI = regexprep(InChI,'/r.*','');
p_layer = regexp(InChI,'/p(.*?)/|/p(.*?)$','tokens');
protonationProtons = 0;
if ~isempty(p_layer)
    individualProtons = cellfun(@(x) {strsplit(x{1},';')},p_layer);
    protonationProtons = cellfun(@(x) sum(cellfun(@(y) eval(y) , x)), individualProtons);
end

%Calc the coefs for all formulas
if (numel(tokens) > 1) || (~isempty(regexp(formula,'(^[0-9]+)'))) || (~isempty(p_layer))
    CoefLists = cellfun(@(x) calcFormula(x), tokens,'UniformOutput',0);
    if ~isempty(p_layer)
        %CoefLists = [CoefLists;{{'H';protonationProtons}}];%was crashing
        %with formula C62H90N13O14P.Co.H2O
        CoefLists{end+1} = {'H';protonationProtons};
    end
    %and now, combine them.
    Elements = {};
    Coefficients = [];
    for i = 1:numel(CoefLists)
        if isempty(CoefLists{i})
            %This should only happen, if there was no actual formula.
            continue
        end
        currentForm = CoefLists{i};
        Elements = [Elements,setdiff(currentForm(1,:),Elements)];
        current_coefs = cell2mat(currentForm(2,:));
        [A,B] = ismember(Elements,currentForm(1,:));
        %Extend the coefficients if necessary
        Coefficients(end+1:numel(Elements)) = 0;
        Coefficients(A) = Coefficients(A)+current_coefs;
    end

    Coefs = num2cell(Coefficients);
    Coefs(cellfun(@(x) x == 1, Coefs)) = {[]};
    Coefs = cellfun(@(x) num2str(x) , Coefs,'UniformOutput',0);
    if nargout > 1
        protons = Coefficients(ismember(Elements,'H'));
    end
    formula = strjoin([Elements , {''}],Coefs);
else
    %had to add this for some inchi, e.g.
    %InChI=1/C21H30O4/c1-19-8-5-14(23)11-13(19)3-4-15-16(19)6-9-20(2)17(15)7-10-21(20,25)18(24)12-22/h11,15-17,22,25H,3-10,12H2,1-2H3/t15-,16+,17+,19+,20+,21+/m1/s1
    protons = numAtomsOfElementInFormula(formula, 'H',0);
end

zero='0';
indZero=strfind(formula,zero);
if ~isempty(indZero)
    if isletter(formula(indZero-1))
        if 0
            warning('Formula contains a zero with a letter preceeding it, replacing with letter.')
            fprintf('%s%s\n','Formula before:, ', formula)
        end
        formula = strrep(formula, formula(indZero-1:indZero), formula(indZero-1));
        if 0
            fprintf('%s%s\n','Formula after:, ', formula)
        end
    end
end

end

function [CoefList] = calcFormula(Formula)
multiplier = 1;
isReplicated = regexp(Formula,'(^[0-9]+)','tokens');
ElementTokens = regexp(Formula,'([A-Z][a-z]?)([0-9]*)','tokens');
Elements = cellfun(@(x) x{1}, ElementTokens,'UniformOutput',0);
Coefs = cellfun(@(x) str2num(x{2}), ElementTokens,'UniformOutput',0);
Coefs(cellfun(@isempty, Coefs)) = {1};

if ~isempty(isReplicated)
    multiplier = str2num(isReplicated{1}{1});
    Coefs = cellfun(@(x) x*multiplier, Coefs,'UniformOutput',0);
end

CoefList = [Elements;Coefs];
end

author: opencobra

repo: cobratoolbox


text:
function p = predict(Theta1, Theta2, X)
%PREDICT Predict the label of an input given a trained neural network
%   p = PREDICT(Theta1, Theta2, X) outputs the predicted label of X given the
%   trained weights of a neural network (Theta1, Theta2)

% Useful values
m = size(X, 1);
num_labels = size(Theta2, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned neural network. You should set p to a
%               vector containing labels between 1 to num_labels.
%
% Hint: The max function might come in useful. In particular, the max
%       function can also return the index of the max element, for more
%       information see 'help max'. If your examples are in rows, then, you
%       can use max(A, [], 2) to obtain the max for each row.
%

a1 = [ones(m, 1) X];

z2 = a1 * Theta1';
a2 = sigmoid(z2);

a2 = [ones(m, 1) a2];
z3 = a2 * Theta2';
hypo = sigmoid(z3);
[probability indices] = max(hypo');
p = indices';

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

end

author: AvaisP


text:
function x = bernoulli_cdf_inv ( cdf, a )

%*****************************************************************************80
%
%% BERNOULLI_CDF_INV inverts the Bernoulli CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real CDF, the value of the CDF.
%    0.0D+00 <= CDF <= 1.0.
%
%    Input, real A, the parameter of the PDF.
%    0.0D+00 <= A <= 1.0.
%
%    Output, integer X, the corresponding argument.
%
  if ( cdf < 0.0 | 1.0 < cdf )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BERNOULLI_CDF_INV - Fatal error!\n' );
    fprintf ( 1, '  CDF < 0 or 1 < CDF.\n' );
    error ( 'BERNOULLI_CDF_INV - Fatal error!' );
  end

  if ( cdf <= 1.0 - a )
    x = 0;
  else
    x = 1;
  end

  return
end


text:
function [H,WtW,WtX] = nnls_fpgm(X,W,options)

% Computes an approximate solution of the following nonnegative least
% squares problem (NNLS)
%
%           min_{H >= 0} ||X-WH||_F^2
%
% using a fast gradient method;
% See Nesterov, Introductory Lectures on Convex Optimization: A Basic
% Course, Kluwer Academic Publisher, 2004.
%
% Input / Output; see nnls_input_output.m
%
% + options.proj allows to use a contraints on the columns or rows of H so
%   that the entries in each column/row sum to at most one
%   options.proj = 0: no projection (default).
%   options.proj = 1: projection of the columns on {x|x>=0, sum(x) <= 1}
%   options.proj = 2: projection of the rows {x|x>=0, sum(x) = 1}
%
% + options.alpha0 is the FPGM  extrapolation parameter (default=0.05)
%
% Code modified from https://sites.google.com/site/nicolasgillis/code
%
% This file has been ported from
%       nnls_FPGM.m at https://gitlab.com/ngillis/nmfbook/-/tree/master/algorithms
%       by Nicolas Gillis (nicolas.gillis@umons.ac.be)

    if nargin <= 2
        options = [];
    end
    if ~isfield(options,'delta')
        options.delta = 1e-6; % Stopping condition depending on evolution of the iterate V:
        % Stop if ||V^{k}-V^{k+1}||_F <= delta * ||V^{0}-V^{1}||_F
        % where V^{k} is the kth iterate.
    end
    if ~isfield(options,'inner_max_epoch')
        options.inner_max_epoch = 500;
    end
    if ~isfield(options,'proj') % Projection on the unit simplex and the origin
        options.proj = 0;
    end
    if ~isfield(options,'alpha0') % Parameter for FPGM ~ extrapolation parameter
        options.alpha0 = 0.05;
    end

    W = full(W);
    [m, n] = size(X);
    [m, r] = size(W);
    WtW = W'*W;
    WtX = W'*X;

    % If no initial matrices are provided, H is initialized as follows:
    if ~isfield(options,'init') || isempty(options.init)
        H = nnls_init(X,W,WtW,WtX);
    else
        H = options.init;
    end

    % Hessian and Lipschitz constant
    L = norm(WtW,2);
    % Linear term
    WtX = W'*X;
    alpha0 = options.alpha0; % Parameter of FPGM, can be tuned.
                             % If options.alpha0 = 0 --> no acceleration, PGM
    alpha(1) = alpha0;
    if options.proj == 1
        H = SimplexProj( H ); % Project columns of H onto the simplex and origin
    elseif options.proj == 0
        H = max(H,0);
    elseif options.proj == 2
        H = SimplexColProj(H'); % Project rows of H onto the simplex
        H = H';
    end
    Y = H; % second sequence
    i = 1;
    % Stop if ||V^{k}-V^{k+1}||_F <= delta * ||V^{0}-V^{1}||_F
    eps0 = 0; eps = 1;
    while i <= options.inner_max_epoch && eps >= options.delta*eps0
        % Previous iterate
        Hp = H;
        % FGM Coefficients; see Nesterov's book
        alpha(i+1) = ( sqrt(alpha(i)^4 + 4*alpha(i)^2 ) - alpha(i)^2) / (2);
        beta(i) = alpha(i)*(1-alpha(i))/(alpha(i)^2+alpha(i+1));
        % Projection step
        H = Y - (WtW*Y-WtX) / L;
        if options.proj == 1
            H = SimplexProj( H ); % Project columns of H onto the set {x|x>=0, sum(x) <= 1}
        elseif options.proj == 0
            H = max(H,0);
        elseif options.proj == 2
            H = SimplexColProj(H'); % Project rows of H onto the simplex
            H = H';
        end
        % `Optimal' linear combination of iterates
        Y = H + beta(i)*(H-Hp);
        if i == 1
            eps0 = norm(H-Hp,'fro');
        end
        eps = norm(H-Hp,'fro');
        i = i + 1;
    end
end

repo: NMFLibrary


function [ icpPoseNew ] = poseSparsification_yq( savePoseName, saveIndexName, icpPose, expDis )
%POSESPARCIFICATION Summary of this function goes here
%   Detailed explanation goes here

    % YQ dataset, data preparation
    % same for park dataset

    keepIndex = [1];
    icpPoseNew(1,:) = icpPose(1,:); % for plotting

    cnt = 1;
    for i = 2:length(icpPose)
        pose1 = [icpPose(i,4), icpPose(i,8)];
        pose2 = [icpPose(keepIndex(cnt),4), icpPose(keepIndex(cnt),8)];
        dis = norm(pose1 - pose2);

        if dis > expDis
            cnt = cnt + 1;
            keepIndex(cnt,:) = i;
            icpPoseNew(cnt,:) = icpPose(i,:);
        end
    end

    keepIndex = keepIndex - 1; % from zero

    length(keepIndex)

    % save the index to the file
    dlmwrite(saveIndexName, keepIndex, 'delimiter', '\t');

    dlmwrite(savePoseName, icpPoseNew, 'delimiter', '\t');

end


text:
function cyclic_reduction_test01 ( )

%*****************************************************************************80
%
%% CYCLIC_REDUCTION_TEST01 tests C83_CR_FA, C83_CR_SL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2010
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CYCLIC_REDUCTION_TEST01\n' );
  fprintf ( 1, '  C83_CR_FA factors a complex tridiagonal matrix;\n' );
  fprintf ( 1, '  C83_CR_SL solves a factored system.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
%
%  Set the matrix values.
%
  a(1,1) = 0.0;
  for j = 2 : n
    a(1,j) = -1.0  - ( j - 1 ) * i;
  end

  for j = 1 : n
    a(2,j) = 2.0 + 2.0 * j * i;
  end

  for j = 1 : n - 1
    a(3,j) = -1.0 - ( j + 1 ) * i;
  end
  a(3,n) = 0.0;
%
%  Set the desired solution.
%
  for j = 1 : n
    x(j) = j + 10 * j * i;
  end
%
%  Compute the corresponding right hand side.
%
  b = c83_mxv ( n, a, x );
%
%  Factor the matrix.
%
  a_cr = c83_cr_fa ( n, a );
%
%  Solve the linear system.
%
  x = c83_cr_sl ( n, a_cr, b );

  c8vec_print ( n, x, '  Solution:' );

  return
end


text:
function scalar = array2scalar(F_array, grid, ge, ft, axis, osc, physQ, intercept)

chkarg(istypesizeof(grid, 'Grid2d') || istypesizeof(grid, 'Grid3d'), ...
	'"grid" should be instance of Grid2d or Grid3d.');
chkarg(istypesizeof(ge, 'GT'), '"ge" should be instance of GT.');
chkarg(istypesizeof(ft, 'FT'), '"ft" should be instance of FT.');
chkarg(istypesizeof(axis, 'Axis'), '"axis" should be instance of Axis.');
chkarg(istypesizeof(osc, 'Oscillation'), '"osc" should be instance of Oscillation.');
chkarg(istypesizeof(physQ, 'PhysQ'), '"physQ" should be instance of PhysQ.');

is3D = true;
if istypesizeof(grid, 'Grid2d')
	is3D = false;
end

if ft == FT.e
    gt = ge;  % grid type for E-field
else  % ft == FT.h
    gt = alter(ge);  % grid type for H-field
end

gt_array = gt(ones(1, Axis.count));
gt_array(axis) = alter(gt);

if is3D
	grid3d = grid;
	chkarg(istypesizeof(F_array, 'arbitrary', grid3d.N), '"F_array" should be %d-by-%d-by-%d array.', grid3d.Ncell{:});

	V = F_array;
	for w = Axis.elems
		V = attach_extra_F(V, ft, gt_array(w), axis, grid3d.bc(w), w, grid3d.kBloch(w)*grid3d.L(w));
	end

	scalar = Scalar3d(V, grid3d, gt_array, osc, physQ, [physQ.symbol, '_', char(axis)]);
else  % grid is Grid2d
	if nargin < 7  % no intercept
		intercept = NaN;
	end
	chkarg(istypesizeof(intercept, 'real'), '"intercept" should be real.');  % NaN is real

	grid2d = grid;
	chkarg(istypesizeof(F_array, 'arbitrary', grid2d.N), '"F_array" should be %d-by-%d array.', grid2d.Ncell{:});

	gt_array = gt_array(grid2d.axis);

	if axis == grid2d.axis(Dir.h)
		axis_temp = Axis.x;
	elseif axis == grid2d.axis(Dir.v)
		axis_temp = Axis.y;
	else  % axis == grid2d.normal_axis
		axis_temp = Axis.z;
	end

	C = F_array;
	C = attach_extra_F(C, ft, gt_array(Dir.h), axis_temp, grid2d.bc(Dir.h), Axis.x, grid2d.kBloch(Dir.h)*grid2d.L(Dir.h));  % treat C as 3D array with one element in 3rd dimension
	C = attach_extra_F(C, ft, gt_array(Dir.v), axis_temp, grid2d.bc(Dir.v), Axis.y, grid2d.kBloch(Dir.v)*grid2d.L(Dir.v));  % treat C as 3D array with one element in 3rd dimension

	scalar = Scalar2d(C, grid2d, gt_array, osc, physQ, [physQ.symbol, '_', char(axis)], intercept);
end

function Fw_array = attach_extra_F(Fw_array, ft, gt, w, bc_v, v, kvLv)
% Augment the Fw array with ghost boundary elements in the v-axis.
% Surprisingly, this function combines attach_extra_E() and attach_extra_H().

% chkarg(istypesizeof(Fw_array, 'complex', [0 0 0]), '"Fw_array" should be 3D array with complex elements.');
chkarg(istypesizeof(ft, 'FT'), '"ft" should be instance of FT.');
chkarg(istypesizeof(gt, 'GT'), '"gt" should be instance of GT.');
chkarg(istypesizeof(w, 'Axis'), '"w" should be instance of Axis.');
chkarg(istypesizeof(bc_v, 'BC'), '"bc_v" should be instance of BC.');
chkarg(istypesizeof(v, 'Axis'), '"v" should be instance of Axis.');
chkarg(istypesizeof(kvLv, 'real'), '"kvLv" should be real.');

ind_n = {':',':',':'};
ind_p = {':',':',':'};
Nv = size(Fw_array, int(v));
if gt == GT.prim
	if bc_v == BC.p
		ind_p{v} = 1;
		Fw_array = cat(int(v), Fw_array, Fw_array(ind_p{:})*exp(-1i*kvLv));
	else  % bc_v == BC.e or BC.m
		% At the ghost point tangential primary fields are always zero, so
		% normal dual fields are also zero.
		size_extra = size(Fw_array);
		size_extra(v) = 1;
		Fw_array = cat(int(v), Fw_array, zeros(size_extra));
	end
else  % gt == GT.dual
	if bc_v == BC.p
		ind_n{v} = Nv;
		ind_p{v} = 1;
		Fw_array = cat(int(v), Fw_array(ind_n{:})*exp(1i*kvLv), Fw_array, Fw_array(ind_p{:})*exp(-1i*kvLv));
	else  % bc_v == BC.e or BC.m
		ind_n{v} = 1;
		ind_p{v} = Nv;
		ft_match_bc = (ft == FT.e && bc_v == BC.e) || (ft == FT.h && bc_v == BC.m);
		if (v ~= w && ft_match_bc) || (v == w && ~ft_match_bc)
			Fw_array = cat(int(v), -Fw_array(ind_n{:}), Fw_array, Fw_array(ind_p{:}));
		else  % (v == w && ft_match_bc) || (v ~= w && ~ft_match_bc)
			Fw_array = cat(int(v), Fw_array(ind_n{:}), Fw_array, Fw_array(ind_p{:}));
		end
	end
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
%% inpainting RobustPCA example: moon picture corrupted with some text
addpath('../');

% read image and add the mask
img = double(imread('moon.tif'))/255;
img = img(141:140+256, 51:50+256);
msk = zeros(size(img));
msk(65:192,65:192) = imresize(imread('text.png'), 0.5);
img_corrupted = img;
img_corrupted(msk > 0) = nan;

fprintf(1, '%d corrupted entries\n', nnz(isnan(img_corrupted)));

% create a matrix X from overlapping patches
ws = 16; % window size
no_patches = size(img, 1) / ws;
X = zeros(no_patches^2, ws^2);
k = 1;
for i = (1:no_patches*2-1)
    for j = (1:no_patches*2-1)
        r1 = 1+(i-1)*ws/2:(i+1)*ws/2;
        r2 = 1+(j-1)*ws/2:(j+1)*ws/2;
        patch = img_corrupted(r1, r2);
        X(k,:) = patch(:);
        k = k + 1;
    end
end

% apply Robust PCA
lambda = 0.02; % close to the default one, but works better
tic
[L, S] = RobustPCA(X, lambda, 1.0, 1e-5);
toc

% reconstruct the image from the overlapping patches in matrix L
img_reconstructed = zeros(size(img));
img_noise = zeros(size(img));
k = 1;
for i = (1:no_patches*2-1)
    for j = (1:no_patches*2-1)
        % average patches to get the image back from L and S
        % todo: in the borders less than 4 patches are averaged
        patch = reshape(L(k,:), ws, ws);
        r1 = 1+(i-1)*ws/2:(i+1)*ws/2;
        r2 = 1+(j-1)*ws/2:(j+1)*ws/2;
        img_reconstructed(r1, r2) = img_reconstructed(r1, r2) + 0.25*patch;
        patch = reshape(S(k,:), ws, ws);
        img_noise(r1, r2) = img_noise(r1, r2) + 0.25*patch;
        k = k + 1;
    end
end
img_final = img_reconstructed;
img_final(~isnan(img_corrupted)) = img_corrupted(~isnan(img_corrupted));

% show the results
figure;
subplot(2,2,1), imshow(img_corrupted), title('Corrupted image')
subplot(2,2,2), imshow(img_final), title('Recovered image')
subplot(2,2,3), imshow(img_reconstructed), title('Recovered low-rank')
subplot(2,2,4), imshow(img_noise), title('Recovered sparse')

fprintf(1, 'ws=%d\tlambda=%f\trank(L)=%d\tcard(S)=%d\terr=%f\n', ...
       ws, lambda, rank(L), nnz(S), norm(img - img_final, 'fro'));

author: dlaptev

repo: RobustPCA

save_path: github-repos/MATLAB/dlaptev-RobustPCA


text:
function [newimage neworigin] = ...
    imagepatch(Ibig,Ismall,origin,offset,padval)
% IMAGEPATCH writes one image over another at location specified by offset.
% I = imagepatch(image1,image2,[originx originy],[offsetx offsety]) returns
% a new image, I, containing image1 overwritten by image 2 at the location
% specified by the offset and origin coordinates.
%
% [I origin] = imagepatch(im1,im2,[origin_x origin_y],[offset_x offset_y])
% returns a new image, I, and a vector, origin, containing the pixel
% coordinates corresponding to the origin in the new image.
%
% I = imagepatch(im1,im2,[origin_x origin_y],[offset_x offset_y],padval)
% returns a new image, where padval specifies the padding value used by
% PADARRAY if it is necessary to expand im1 to accomodate im2 at the offset
% coordinates.
%
% Example:
%             moon = imread('moon.tif');
%             pout = imread('pout.tif');
%             %centered image
%             origin = round(flipdim(size(moon),2)./2);
%             offset = [200 250];
%             imshow(imagepatch(moon,pout,origin,offset));
%
% See also IMAGESC, IMAGE, IMREAD, IMWRITE

% written by David Sedarsky
% July, 2007

if nargin == 0
    %if true, show example of use
    padval = 0;
    help imagepatch
    Ibig = ones(10);        %initialize example bigframe
    Ibig(1:2:end,:)=0.7;    %draw some stripes for fun
    Ismall = 0.5 .* ones(6);%initialize example subframe
    origin = [5 5];         %set example origin coordinates
    offset = [5 -4];        %set example offset for subframe
elseif nargin < 4
    error(['4 inputs, (imagematrix,imagematrix,xycoords,xyoffsets)'...
        ', required.']);
elseif nargin == 4
    padval = 0;  %set default pad value for use of PADARRAY below
elseif nargin > 5
    %crude padval check
    if ischar(padval)
        if ~strcmpi(padval,'circular')||strcmp(padval,'c')
            if ~strcmpi(padval,'replicate')||strcmp(padval,'r')
                if ~strcmpi(padval,'symmetric')||strcmp(padval,'s')
                    error(['Invalid pad value. See HELP PADARRAY for'...
                        ' valid pad value input format.']);
                end
            end
        end
    elseif isscalar(padval)
        %scalar value is fine
    else
        error(['Invalid pad value. See HELP PADARRAY for'...
                        ' valid pad value input format.']);
    end
end

bigframe = Ibig;

%get image dimensions
[ydim xdim] = size(Ismall);
[ybig xbig] = size(Ibig);

%get origin
xorig = origin(1);
yorig = origin(2);

%get offsets and offset directions
xoff = offset(1);
if xoff > 0
    xdir = 'positive';
elseif xoff < 0
    xdir = 'negative';
else
    xdir = 'zero';
end
yoff = offset(2);
if yoff > 0
    ydir = 'positive';
elseif yoff < 0
    ydir = 'negative';
else
    ydir = 'zero';
end

%spacer is the distance between origin and edge of subframe
xspacer = abs(xoff)-floor(xdim./2);
yspacer = abs(yoff)-floor(ydim./2);

%distance: Origin to side of new frame in pixels
if strcmp(xdir,'positive')
    right_side = xspacer + xdim;
    left_side = xspacer + 1;
elseif strcmp(xdir,'negative')
    right_side = -xspacer;
    left_side = -xspacer - (xdim-1)
else
    %do zero case
%     right_side = floor(xdim./2);
%     left_side = -floor(xdim./2);
    right_side = xspacer + xdim;
    left_side = xspacer + 1;
end
if strcmp(ydir,'positive')
    top_side = yspacer + 1;
    bottom_side = yspacer + ydim;
elseif strcmp(ydir,'negative')
    top_side = -yspacer - (ydim-1);
    bottom_side = -yspacer;
else
    %do zero case
%     top_side = -floor(ydim./2);
%     bottom_side = floor(ydim./2);
    top_side = yspacer + 1;
    bottom_side = yspacer + ydim;
end

%make a copy of the origin to modify later when we change bigframe
neworigin = origin;

%Check that subframe lies within bigframe, pad bigframe if it's too small
%checking RIGHT
if xorig + right_side > xbig
    %if true, need to expand bigframe
    right_padsize = [0 xorig+right_side-xbig];
    bigframe = padarray(bigframe,right_padsize,padval,'post');
end
%checking LEFT
if xorig + left_side < 1
    %if true, need to expand bigframe
    left_padsize = [0 abs(xorig+left_side)+1];
    bigframe = padarray(bigframe,left_padsize,padval,'pre');
    %adjust origin -- leftside padding moves the origin wrt first pixel
    neworigin = neworigin + left_padsize(end:-1:1)
end
%checking BOTTOM
if yorig + bottom_side > ybig
    %if true, need to expand bigframe
    bottom_padsize = [yorig+bottom_side-ybig 0];
    bigframe = padarray(bigframe,bottom_padsize,padval,'post');
end
%checking TOP
if yorig + top_side < 1
    %if true, need to expand bigframe
    top_padsize = [abs(yorig+top_side)+1 0];
    bigframe = padarray(bigframe,top_padsize,padval,'pre');
    %adjust origin -- topside padding moves the origin wrt first pixel
    neworigin = neworigin + top_padsize(end:-1:1);
end

nxorig = neworigin(1);
nyorig = neworigin(2);

%distance: new Origin to side of new frame in pixels
right = nxorig + right_side;
left = nxorig + left_side;
top = nyorig + top_side;
bottom = nyorig + bottom_side;

%place new frame inside expanded bigframe
bigframe(top:bottom,left:right) = Ismall;

if nargin == 0
    figure('Name','''imagepatch'' demonstration image.');
    imshow(bigframe,'InitialMagnification', 'fit');
    title(['offset:  x = ' num2str(xoff) ' ,  y = ' num2str(yoff)])
else
    newimage = bigframe;
end


text:
function [data_sets] = fiff_find_evoked(fname)
%
% [data_sets] = fiff_find_evoked(fname)
%
% Find all evoked data sets in a fif file and create a list of descriptors
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%

global FIFF;
if isempty(FIFF)
    FIFF = fiff_define_constants();
end

me = 'MNE:fiff_find_evoked';

%
%   Open the file
%
[ fid, tree ] = fiff_open(fname);
data_sets = struct('comment',{},'aspect_kind',{},'aspect_name',{});
%
%   Find all evoked data sets
%
evoked = fiff_dir_tree_find(tree, FIFF.FIFFB_EVOKED);
if length(evoked) == 0
    fclose(fid);
    return
end
%
%   Identify the aspects
%
naspect = 0;
for k = 1:length(evoked)
    sets(k).aspects  = fiff_dir_tree_find(evoked(k), FIFF.FIFFB_ASPECT);
    sets(k).naspect = length(sets(k).aspects);
    naspect = naspect + sets(k).naspect;
end
%
%   Collect the desired information
%
count = 1;
for k = 1:length(evoked)
    evoked_comment = find_tag(evoked(k), FIFF.FIFF_COMMENT);
    for a = 1:sets(k).naspect
        aspect_comment = find_tag(sets(k).aspects(a), FIFF.FIFF_COMMENT);
        aspect_kind    = find_tag(sets(k).aspects(a), FIFF.FIFF_ASPECT_KIND);
        if ~isempty(aspect_comment)
            data_sets(count).comment = aspect_comment.data;
        elseif ~isempty(evoked_comment)
            data_sets(count).comment = evoked_comment.data;
        else
            data_sets(count).comment = 'No comment';
        end
        if ~isempty(aspect_kind)
            data_sets(count).aspect_kind = aspect_kind.data;
        else
            data_sets(count).aspect_kind = -1;
        end
        switch data_sets(count).aspect_kind
        case FIFF.FIFFV_ASPECT_AVERAGE
            data_sets(count).aspect_name = 'Average';
        case FIFF.FIFFV_ASPECT_STD_ERR
            data_sets(count).aspect_name = 'Standard error';
        case FIFF.FIFFV_ASPECT_SINGLE
            data_sets(count).aspect_name = 'Single';
        case FIFF.FIFFV_ASPECT_SUBAVERAGE
            data_sets(count).aspect_name = 'Subaverage';
        case FIFF.FIFFV_ASPECT_ALTAVERAGE
            data_sets(count).aspect_name = 'Alt. subaverage';
        case FIFF.FIFFV_ASPECT_SAMPLE
            data_sets(count).aspect_name = 'Sample';
        case FIFF.FIFFV_ASPECT_POWER_DENSITY
            data_sets(count).aspect_name = 'Power density spectrum';
        case FIFF.FIFFV_ASPECT_DIPOLE_WAVE
            data_sets(count).aspect_name = 'Dipole source waveform';
            otherwise
            data_sets(count).aspect_name = 'Unknown';
        end
        count = count + 1;
    end
end

fclose(fid);

return

function [tag] = find_tag(node, findkind)

for p = 1:node.nent
    kind = node.dir(p).kind;
    pos  = node.dir(p).pos;
    if kind == findkind
        tag = fiff_read_tag(fid, pos);
        return
    end
end
tag = [];
return
end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function Calls = Automerge_Callback(Calls1,Calls2,AudioFile)
%% Merges two detection files into one

Calls=[Calls1; Calls2];

% Audio info
audio_info = audioinfo(AudioFile);

%% Merge overlapping boxes
Calls = merge_boxes(Calls.Box, Calls.Score, Calls.Type, audio_info, 1, 0, 0);

author: DrCoffey

repo: DeepSqueak

save_path: github-repos/MATLAB/DrCoffey-DeepSqueak


text:
function y = tmat_mxp2 ( a, n, x )

%*****************************************************************************80
%
%% TMAT_MXP2 multiplies a geometric transformation matrix times N points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Foley, van Dam, Feiner, Hughes,
%    Computer Graphics, Principles and Practice,
%    Addison Wesley, Second Edition, 1990.
%
%  Parameters:
%
%    Input, real A(4,4), the geometric transformation matrix.
%
%    Input, integer N, the number of points to be multiplied.
%
%    Input, real X(3,N), the points to be multiplied.
%
%    Output, real Y(3,N), the transformed points.  Each product is
%    accumulated in a temporary vector, and then assigned to the
%    result.  Therefore, it is legal for X and Y to share memory.
%
  for k = 1 : n
    y(1:3,k) = a(1:3,4) + a(1:3,1:3) * x(1:3,k);
  end

  return
end