% params
totalNum = 7385291;
bValue = 4000;
lamda = 0.5;
qFile = '/home/yinhuan/Data/mapModel/kitti_10/weightVector.txt';
visFilesDir = '/home/yinhuan/Data/mapModel/kitti_10/visMatrix/';
maxQ = 48;
minQ = 1;

splitLength_0 = 50;
saveResultsDir_0 = '/home/yinhuan/Data/mapModel/kitti_10/iter_b_4000/0/';
saveReIndexFile_0 = '/home/yinhuan/Data/mapModel/kitti_10/iter_b_4000/index/0.txt';
saveNewQFile_0 = '/home/yinhuan/Data/mapModel/kitti_10/iter_b_4000/weight/0.txt';

splitLength_1 = 100;
saveResultsDir_1 = '/home/yinhuan/Data/mapModel/kitti_10/iter_b_4000/1/';
saveReIndexFile_1 = '/home/yinhuan/Data/mapModel/kitti_10/iter_b_4000/index/1.txt';
saveNewQFile_1 = '/home/yinhuan/Data/mapModel/kitti_10/iter_b_4000/weight/1.txt';

splitLength_2 = 200;
saveResultsDir_2 = '/home/yinhuan/Data/mapModel/kitti_10/iter_b_4000/2/';
saveReIndexFile_2 = '/home/yinhuan/Data/mapModel/kitti_10/iter_b_4000/index/2.txt';
saveNewQFile_2 = '/home/yinhuan/Data/mapModel/kitti_10/iter_b_4000/weight/2.txt';

splitLength_3 = 400;
saveResultsDir_3 = '/home/yinhuan/Data/mapModel/kitti_10/iter_b_4000/3/';
saveReIndexFile_3 = '/home/yinhuan/Data/mapModel/kitti_10/iter_b_4000/index/3.txt';
saveNewQFile_3 = '/home/yinhuan/Data/mapModel/kitti_10/iter_b_4000/weight/3.txt';

splitLength_4 = 800;
saveResultsDir_4 = '/home/yinhuan/Data/mapModel/kitti_10/iter_b_4000/4/';
saveReIndexFile_4 = '/home/yinhuan/Data/mapModel/kitti_10/iter_b_4000/index/4.txt';
saveNewQFile_4 = '/home/yinhuan/Data/mapModel/kitti_10/iter_b_4000/weight/4.txt';

%% iteration of optimization

visCells_origin = fromVisDirToCells(visFilesDir);

[ epsilon_soft_0, time_sum_0  ] = Uniform_loopCompress_Cells( lamda, qFile, visCells_origin, splitLength_0, totalNum, bValue, saveResultsDir_0 );
min_cost_0 = get_min_cost(minQ, maxQ, saveResultsDir_0, qFile, lamda, epsilon_soft_0);
[totalNum_0, visCells_0] = deleteZeroPoints_Cells( qFile, visCells_origin, saveResultsDir_0, saveReIndexFile_0, saveNewQFile_0 );

[ epsilon_soft_1, time_sum_1  ] = Uniform_loopCompress_Cells( lamda, saveNewQFile_0, visCells_0, splitLength_1, totalNum_0, bValue, saveResultsDir_1 );
min_cost_1 = get_min_cost(minQ, maxQ, saveResultsDir_1, saveNewQFile_0, lamda, epsilon_soft_1);
[totalNum_1, visCells_1] = deleteZeroPoints_Cells( saveNewQFile_0, visCells_0, saveResultsDir_1, saveReIndexFile_1, saveNewQFile_1 );

[ epsilon_soft_2, time_sum_2  ] = Uniform_loopCompress_Cells( lamda, saveNewQFile_1, visCells_1, splitLength_2, totalNum_1, bValue, saveResultsDir_2 );
min_cost_2 = get_min_cost(minQ, maxQ, saveResultsDir_2, saveNewQFile_1, lamda, epsilon_soft_2);
[totalNum_2, visCells_2] = deleteZeroPoints_Cells( saveNewQFile_1, visCells_1, saveResultsDir_2, saveReIndexFile_2, saveNewQFile_2 );

[ epsilon_soft_3, time_sum_3  ] = Uniform_loopCompress_Cells( lamda, saveNewQFile_2, visCells_2, splitLength_3, totalNum_2, bValue, saveResultsDir_3 );
min_cost_3 = get_min_cost(minQ, maxQ, saveResultsDir_3, saveNewQFile_2, lamda, epsilon_soft_3);
[totalNum_3, visCells_3] = deleteZeroPoints_Cells( saveNewQFile_2, visCells_2, saveResultsDir_3, saveReIndexFile_3, saveNewQFile_3 );

[ epsilon_soft_4, time_sum_4  ] = Uniform_loopCompress_Cells( lamda, saveNewQFile_3, visCells_3, splitLength_4, totalNum_3, bValue, saveResultsDir_4 );
min_cost_4 = get_min_cost(minQ, maxQ, saveResultsDir_4, saveNewQFile_3, lamda, epsilon_soft_4);
[totalNum_4, visCells_4] = deleteZeroPoints_Cells( saveNewQFile_3, visCells_3, saveResultsDir_4, saveReIndexFile_4, saveNewQFile_4 );

%% save the results

compressIndex_3 = anal_reindex_last(saveResultsDir_4, saveReIndexFile_3);
compressIndex_2 = anal_reindex_middle(compressIndex_3, saveReIndexFile_2);
compressIndex_1 = anal_reindex_middle(compressIndex_2, saveReIndexFile_1);
compressIndex_0 = anal_reindex_middle(compressIndex_1, saveReIndexFile_0);

dlmwrite('/home/yinhuan/Data/mapModel/kitti_10/iter_b_4000/compressIndex_final.txt', compressIndex_0, 'precision', '%d');


text:
%%

addpath(genpath('..'));
addpath(genpath('../../../../scripts'));

resDir = getResDir;
dataDir = getDataDir;

challenge = 2; % 2D MOT 2015
% challenge = 5; % MOT16

% eval all
challenge = 2; % 2D MOT 2015
challenge = 5; % MOT16

detectorName = 'acf';
detectorName = 'dpm';
detectorName = 'fasterrcnndet-th-0';
% detectorName = '';

sequences=[];
% allSeq={'TUD-Stadtmitte'};
% sequences='PETS09-S2L1';
% sequences={'TUD-Stadtmitte','PETS09-S2L1'};
allSeq = parseSequences([0 challenge],dataDir);
% allSeq = {'MOT16-01'}

    dt0=[];
    gt0=[];
% for each sequence
allFrCnt=0;
for seq=allSeq

    [seqName, seqFolder, imgFolder, imgExt, F, dirImages] ...
        = getSeqInfo(seq, dataDir);

%     detFolder = [dataDir,seqName,filesep,'det',filesep];
%     detFile = [detFolder,'det.txt'];
    [detFolder, detFile]=getDetInfo(seqName,dataDir,detectorName);

%     if ~exist(detFile,'file'), continue; end

    gtFolder= [dataDir,seqName,filesep,'gt',filesep];
    gtFile = [gtFolder,'gt.txt'];

    detRaw=dlmread(detFile);
    gtRaw=dlmread(gtFile);

    for t=1:F
        allFrCnt=allFrCnt+1;
        exdet=find(detRaw(:,1)==t);
        bbox=detRaw(exdet,3:7);
        dt0{allFrCnt}=bbox;

        exgt=find(gtRaw(:,1)==t & gtRaw(:,8)==1);
        gt0{allFrCnt}=[gtRaw(exgt,3:6) zeros(length(exgt),1)];

    end
end

allSeq = parseSequences([1 challenge],dataDir);
% allSeq = {'MOT16-02'}

    dt1=[];
    gt1=[];

% for each sequence
allFrCnt1=0;
for seq=allSeq

    [seqName, seqFolder, imgFolder, imgExt, F, dirImages] ...
        = getSeqInfo(seq, dataDir);

%     detFolder = [dataDir,seqName,filesep,'det',filesep];
%     detFile = [detFolder,'det.txt'];
    [detFolder, detFile]=getDetInfo(seqName,dataDir,detectorName);
%     if ~exist(detFile,'file'), continue; end

    gtFolder= [dataDir,seqName,filesep,'gt',filesep];
    gtFile = [gtFolder,'gt.txt'];

    detRaw=dlmread(detFile);
    gtRaw=dlmread(gtFile);

    for t=1:F
        allFrCnt1=allFrCnt1+1;
        exdet=find(detRaw(:,1)==t);
        bbox=detRaw(exdet,3:7);
        dt1{allFrCnt1}=bbox;

        exgt=find(gtRaw(:,1)==t & gtRaw(:,8)==1);
        gt1{allFrCnt1}=[gtRaw(exgt,3:6) zeros(length(exgt),1)];

    end
end

% lims=[3.1e-3 1e1 .05 1];
lims=[0 1 0 1];
ref=10.^(-2:.25:0);
ref=0:.1:1;

%%
clf; hold on; box on;

% Training Set
[gt,dt]=bbGt('evalRes',gt0,dt0);
[fp,tp,score,miss] = bbGt('compRoc',gt,dt,0,ref);

% Test Set
[gt,dt]=bbGt('evalRes',gt1,dt1);
[fp1,tp1,score1,miss1] = bbGt('compRoc',gt,dt,0,ref);

set(gca,'FontSize',12)
plot(fp,tp,'linewidth',2);
plot(fp1,tp1,'linewidth',2,'color','r');

opX=fp(end); opY=tp(end);
plot(opX,opY,'.','MarkerSize',25);
opX1=fp1(end); opY1=tp1(end);
plot(opX1,opY1,'.','MarkerSize',25,'color','r');

%%
fs=22;
set(gca,'FontSize',fs);
text(opX+.25,opY+.1,sprintf('Rcll: %.1f %%\nPrc: %.1f %%',opX*100,opY*100),'VerticalAlignment','top','HorizontalAlignment','right','FontSize',fs-4);
text(opX1-.05,opY1-.025,sprintf('Rcll: %.1f %%\nPrc: %.1f %%',opX1*100,opY1*100),'VerticalAlignment','top','HorizontalAlignment','right','FontSize',fs-4);
xlim(lims(1:2)); ylim(lims(3:4));

xlabel('Recall'); ylabel('Precision');
title('Detector Performance');
legend('Training Set','Test Set');

% export_fig 'det-dollar.pdf' -png -pdf;
% saveas(gcf,'det-dollar.pdf');
% saveas(gcf,'det-dollar.png');
% export_fig 'det-dollar.pdf' -transparent;
pdfFileName = sprintf('det-%s.pdf',detectorName);
export_fig(pdfFileName, '-transparent');
% export_fig '../../../../../../papers/2015/arxiv-benchmark/figures/det-dollar.pdf' -transparent;

% miss=exp(mean(log(max(1e-10,1-miss)))); roc=[score fp tp];
% figure(show); plotRoc([fp tp],'logx',1,'logy',0,'xLbl','fppi',...
%   'lims',lims,'color','g','smooth',1,'fpTarget',ref);
% title(sprintf('log-average miss rate = %.2f%%',miss*100));

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function t_test_fcns(quiet)
%T_TEST_FCNS  Test T_OK and T_IS and manually check output of failed tests.

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

if nargin < 1
    quiet = 0;
end

ntests = 5;
npass = 29;
nfail = 18;

g = [];
e = [];

tol = 12;
k = 0;
str = {'EXPECT FAIL'; 'EXPECT PASS'};
f = @(s) str{s+1};

t_begin(npass+nfail, quiet);

%% t_ok
k = k + 1; e(k) = 1;
t = sprintf('%s : t_ok(1, ...)', f(e(k)));
g(k) = t_ok(1, t);

k = k + 1; e(k) = 0;
t = sprintf('%s : t_ok(0, ...)', f(e(k)));
g(k) = t_ok(0, t);

%% t_is
k = k + 1; e(k) = 1;
t = sprintf('%s : t_is([], [], ...)', f(e(k)));
expected = [];
got = [];
g(k) = t_is(got, expected, tol, t);

k = k + 1; e(k) = 0;
t = sprintf('%s : t_is(1, [], ...)', f(e(k)));
expected = [];
got = 1;
g(k) = t_is(got, expected, tol, t);

k = k + 1; e(k) = 1;
t = sprintf('%s : t_is(3, 3, ...)', f(e(k)));
expected = 3;
got = 3;
g(k) = t_is(got, expected, tol, t);

k = k + 1; e(k) = 0;
t = sprintf('%s : t_is(1, 3, ...)', f(e(k)));
expected = 3;
got = 1;
g(k) = t_is(got, expected, tol, t);

k = k + 1; e(k) = 1;
t = sprintf('%s : t_is([3 3; 3 3], 3, ...)', f(e(k)));
expected = 3;
got = [3 3; 3 3];
g(k) = t_is(got, expected, tol, t);

k = k + 1; e(k) = 0;
t = sprintf('%s : t_is([3 4; 3 3], 3, ...)', f(e(k)));
expected = 3;
got = [3 4; 3 3];
g(k) = t_is(got, expected, tol, t);

k = k + 1; e(k) = 1;
t = sprintf('%s : t_is([3 3; 4 3], [3 3; 4 3], ...)', f(e(k)));
expected = [3 3; 4 3];
got = [3 3; 4 3];
g(k) = t_is(got, expected, tol, t);

k = k + 1; e(k) = 0;
t = sprintf('%s : t_is([3 4; 3 3], [3 3; 4 3], ...)', f(e(k)));
expected = [3 3; 4 3];
got = [3 4; 3 3];
g(k) = t_is(got, expected, tol, t);

k = k + 1; e(k) = 1;
t = sprintf('%s : t_is(ones(2,4,3), ones(2,4,3), ...)', f(e(k)));
expected = ones(2,4,3);
got = ones(2,4,3);
g(k) = t_is(got, expected, tol, t);

k = k + 1; e(k) = 0;
t = sprintf('%s : t_is(ones(2,4), ones(2,4,3), ...)', f(e(k)));
expected = ones(2,4,3);
got = ones(2,4);
g(k) = t_is(got, expected, tol, t);

k = k + 1; e(k) = 0;
t = sprintf('%s : t_is(ones(2,4,3), <one-elem-different> > tol, ...)', f(e(k)));
expected = ones(2,4,3);
got = ones(2,4,3);
got(2,3,2) = 1.0005;
g(k) = t_is(got, expected, 4, t);

k = k + 1; e(k) = 1;
t = sprintf('%s : t_is(ones(2,4,3), <one-elem-different> < tol, ...)', f(e(k)));
expected = ones(2,4,3);
got = ones(2,4,3);
got(2,3,2) = 1.0005;
g(k) = t_is(got, expected, 3, t);

k = k + 1; e(k) = 1;
t = sprintf('%s : t_is(NaN, NaN, ...)', f(e(k)));
expected = NaN;
got = NaN;
g(k) = t_is(got, expected, tol, t);

k = k + 1; e(k) = 0;
t = sprintf('%s : t_is([], NaN, ...)', f(e(k)));
expected = NaN;
got = [];
g(k) = t_is(got, expected, tol, t);

k = k + 1; e(k) = 0;
t = sprintf('%s : t_is(NaN, [], ...)', f(e(k)));
expected = [];
got = NaN;
g(k) = t_is(got, expected, tol, t);

k = k + 1; e(k) = 0;
t = sprintf('%s : t_is(1, NaN, ...)', f(e(k)));
expected = NaN;
got = 1;
g(k) = t_is(got, expected, tol, t);

k = k + 1; e(k) = 0;
t = sprintf('%s : t_is(NaN, 3, ...)', f(e(k)));
expected = 3;
got = NaN;
g(k) = t_is(got, expected, tol, t);

k = k + 1; e(k) = 1;
t = sprintf('%s : t_is(NaN(3,2), NaN, ...)', f(e(k)));
expected = NaN;
got = NaN(3,2);
g(k) = t_is(got, expected, tol, t);

k = k + 1; e(k) = 0;
t = sprintf('%s : t_is(int32(0.6), 1.0, ...)', f(e(k)));
expected = 0.6;
got = int32(1.0);
g(k) = t_is(got, expected, tol, t);

%% t_str_match
k = k + 1; e(k) = 1;
t = sprintf('%s : t_str_match(got, expected, ...)', f(e(k)));
expected = 'hello t_str_match';
got = 'hello t_str_match';
g(k) = t_str_match(got, expected, t);

k = k + 1; e(k) = 0;
t = sprintf('%s : t_str_match(got, expected, ...)', f(e(k)));
expected = 'hello t_str_match';
got = 'herro t_str_match';
g(k) = t_str_match(got, expected, t);

k = k + 1; e(k) = 1;
t = sprintf('%s : t_str_match(... reps[re])', f(e(k)));
expected = 'hello t_str_match';
got = 'herro t_str_match';
reps = {{'r{2}', 'll'}};
g(k) = t_str_match(got, expected, t, reps);

k = k + 1; e(k) = 0;
t = sprintf('%s : t_str_match(... reps[re])', f(e(k)));
expected = 'hello t_str_match';
got = 'herro t_str_match';
reps = {{'r{2}', 'LL'}};
g(k) = t_str_match(got, expected, t, reps);

k = k + 1; e(k) = 1;
t = sprintf('%s : t_str_match(... reps[!re])', f(e(k)));
expected = 'hello t_str_match';
got = 'herro t_str_match';
reps = {{'rr', 'll', 0}};
g(k) = t_str_match(got, expected, t, reps);

k = k + 1; e(k) = 0;
t = sprintf('%s : t_str_match(... reps[!re])', f(e(k)));
expected = 'hello t_str_match';
got = 'herro t_str_match';
reps = {{'r{2}', 'll', 0}};
g(k) = t_str_match(got, expected, t, reps);

k = k + 1; e(k) = 1;
t = sprintf('%s : t_str_match(... reps[both])', f(e(k)));
expected = 'hel4o t_str_match';
got = 'he56o t_str_match';
reps = {{'\d', 'l', 1, 1}};
g(k) = t_str_match(got, expected, t, reps);

k = k + 1; e(k) = 0;
t = sprintf('%s : t_str_match(... reps[both])', f(e(k)));
expected = 'heL4o t_str_match';
got = 'he56o t_str_match';
reps = {{'\d', 'l', 1, 1}};
g(k) = t_str_match(got, expected, t, reps);

%% t_file_match
[twd, n, ext] = fileparts(which('t_test_fcns'));
exp_fname = fullfile(twd, 't_have_feature', 'rithmaticker.m');
num = fix(1e7*rand);
fname1 = sprintf('test-output-1-%d.txt', num);
fname2 = sprintf('test-output-2-%d.txt', num);
fname3 = sprintf('test-output-3-%d.txt', num);
ex = fileread(exp_fname);
if isempty(strfind(ex, char([13 10])))
    ex2 = strrep(ex, char(10), char([13 10]));
else
    ex2 = strrep(ex, char([13 10]), char(10));
end
ex3 = regexprep(ex2, '\d', 'N');

k = k + 1; e(k) = 1;
t = sprintf('%s : fileread', f(e(k)));
g(k) = t_str_match(ex(1:8), 'function', t);

%% write files
filewrite(fname1, ex);
filewrite(fname2, ex2);
filewrite(fname3, ex3);

k = k + 1; e(k) = 1;
t = sprintf('%s : filewrite 1', f(e(k)));
ok = exist(fname1, 'file') == 2 && strcmp(ex, fileread(fname1));
g(k) = t_ok(ok, t);

k = k + 1; e(k) = 1;
t = sprintf('%s : filewrite 2', f(e(k)));
ok = exist(fname2, 'file') == 2 && ~strcmp(ex, fileread(fname2));
g(k) = t_ok(ok, t);

k = k + 1; e(k) = 1;
t = sprintf('%s : filewrite 3', f(e(k)));
ok = exist(fname3, 'file') == 2 && ~strcmp(ex, fileread(fname3));
g(k) = t_ok(ok, t);

%% t_file_match
k = k + 1; e(k) = 1;
t = sprintf('%s : t_file_match 1', f(e(k)));
g(k) = t_file_match(fname1, exp_fname, t);
k = k + 1; e(k) = 1;
g(k) = t_ok(exist(fname1, 'file') == 2, [t ' : file not deleted']);

k = k + 1; e(k) = 1;
t = sprintf('%s : t_file_match 1 (del)', f(e(k)));
g(k) = t_file_match(fname1, exp_fname, t, {}, 1);
k = k + 1; e(k) = 1;
g(k) = t_ok(exist(fname1, 'file') == 0, [t ' : file deleted']);

k = k + 1; e(k) = 1;
t = sprintf('%s : t_file_match 2', f(e(k)));
g(k) = t_file_match(fname2, exp_fname, t);
k = k + 1; e(k) = 1;
g(k) = t_ok(exist(fname2, 'file') == 2, [t ' : file not deleted']);

k = k + 1; e(k) = 1;
t = sprintf('%s : t_file_match 2 (del)', f(e(k)));
g(k) = t_file_match(fname2, exp_fname, t, {}, 1);
k = k + 1; e(k) = 1;
g(k) = t_ok(exist(fname2, 'file') == 0, [t ' : file deleted']);

k = k + 1; e(k) = 0;
t = sprintf('%s : t_file_match 3 (del)', f(e(k)));
g(k) = t_file_match(fname3, exp_fname, t, {}, 1);
k = k + 1; e(k) = 1;
g(k) = t_ok(exist(fname3, 'file') == 2, [t ' : file not deleted']);

k = k + 1; e(k) = 0;
t = sprintf('%s : t_file_match 3 (reps,del)', f(e(k)));
reps = {{'\d', 'N', 1}};
g(k) = t_file_match(fname3, exp_fname, t, reps);
k = k + 1; e(k) = 1;
g(k) = t_ok(exist(fname3, 'file') == 2, [t ' : file not deleted']);

k = k + 1; e(k) = 1;
t = sprintf('%s : t_file_match 3 (both,del)', f(e(k)));
reps = {{'[\dN]', 'M', 1, 1}};
g(k) = t_file_match(fname3, exp_fname, t, reps, 1);
k = k + 1; e(k) = 1;
g(k) = t_ok(exist(fname3, 'file') == 0, [t ' : file deleted']);

%%-----  reset tests, if this test passes, everything is as expected  -----
t_begin(ntests, quiet);

t = 'number of tests as expected';
t_is(length(g), npass+nfail, tol, t);
t_is(length(e), npass+nfail, tol, t);

t = 'verify test consistency';
t_is(sum( e), npass, tol, t);
t_is(sum(~e), nfail, tol, t);

t = 'passes and fails all as expected';
t_is(g, e, tol, t);

t_end;

function filewrite(fname, content)
[fid, msg] = fopen(fname, 'w');
if fid == -1
    error('%s\nt_test_fcns/filewrite: could not write file ''%s''', msg, fname);
end
fprintf(fid, '%s', content);
fclose(fid);

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function theRank=rankPartition(thePartition)
%%RANKPARTITION Rank a given partition of the integer n when considering
%               the partitions given in reverse lexicographic order. A
%               partition of the integer n is a set of integers that sum to
%               n.
%
%INPUTS: thePartition An sX1 or 1Xs partition. The values start at 1.
%
%OUTPUTS: the Rank. The integer rank of the partition starting from 1.
%
%Let p(n,m) be the number of partitions of n with up to m parts. It is
%known that p(n,m)=p(n,m-1)+p(n-m,m). This relation is given in Chapter
%7.2.1.4 of [1]. The rank generated is derived using that recurrence
%relation.
%
%EXAMPLE:
%The function getNextPartition generates partitions in reverse
%lexicographic order. Here, we demonstrate that the partitions generated
%are ranked correctly. We partition the integer 6.
% n=6;
% [thePartition,theData]=getNextPartition(n);
% while(~isempty(thePartition))
%     theRank=rankPartition(thePartition)
%     [thePartition,theData]=getNextPartition(n,theData);
% end
%One will see that theRank goes from 1 to 11 in order.
%
%REFERENCES:
%[1] D. E. Knuth, The Art of Computer Programming. Vol. 4, Fascicle 3:
%    Generating all Combinations and Partitions, Upper Saddle River, NJ:
%    Addison-Wesley, 2009.
%
%September 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    m=length(thePartition);
    n=sum(thePartition);

    thePartition=sort(thePartition,'descend');

    theRank=0;
    for idx=1:m
        curParts=thePartition(idx)-1;
        if(curParts>0)
            theRank=theRank+numMPartitions(sum(thePartition(idx:m)),curParts,true);
        end
    end
    theRank=numberOfPartitions(n)-theRank;
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function K = correlation_brownian ( s, t )

%*****************************************************************************80
%
%% CORRELATION_BROWNIAN computes the Brownian correlation function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real S(*), T(*), pairs of argument values.
%
%    Output, real K(*), the correlation function values
%
  K = ones ( size ( s ) );
  i = find ( max ( s, t ) ~= 0.0 );
  K(i) = sqrt ( min ( s(i), t(i) ) ./ max ( s(i), t(i) ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef TestPhase
    %TestPhase

    methods (Static)
        function test_radians
            x = randn(5,6);
            y = randn(5,6);
            ang = cv.phase(x, y);
            validateattributes(ang, {class(x)}, {'size',size(x)});

            x = single(x);
            y = single(y);
            ang = cv.phase(x, y);
            validateattributes(ang, {class(x)}, {'size',size(x)});
        end

        function test_degrees
            x = randn(5,6);
            y = randn(5,6);
            ang = cv.phase(x, y, 'Degrees',true);
            validateattributes(ang, {class(x)}, {'size',size(x)});
        end

        function test_error_argnum
            try
                cv.phase();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function y = atanh_pos(x,rnd)
% rigorous atanh for nonnegative double vector x<1 rounded corresponding to rnd
% for internal use in rigorous atanh, acoth

% written  12/30/98     S.M. Rump
% modified 08/31/99     S.M. Rump  major revision
%

  y = x;

  % atanh(x) = log( 1 + 2*x/(1-x) ) / 2

  index = ( x<0.33 );                    % 0 <= x < .33
  if any(index)
    X = x(index);
    setround(-rnd)
    E = 1 - X;
    setround(rnd)
    E = 2*X./E;
    y(index) = log_1( E , rnd ) / 2;     % 0 <= E < 1
  end

  index = ~index;
  if any(index)                          % .33 <= x <= 1
    X = x(index);
    setround(-rnd)
    E = 1 - X;
    setround(rnd)
    E = 1 + 2*X./E;
    y(index) = log_rnd( E , rnd ) / 2;
  end

repo: OpenMAS


text:
%% SQUAREPOISSON Poisson equation in a square domain.
%
%   squarePoisson computes approximations of the Poisson equation in the
%   unit square on a sequence of meshes obtained by uniform refinement. It
%   plots the approximation error (in L2 norm or H1 norm) vs the number of
%   nodes.
%
% See also Poisson, crack, Lshape
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

close all; clear all;
%% Parameters
maxIt = 4; N = zeros(maxIt,1);
erruIuh = zeros(maxIt,1);

%% Generate an initial mesh
[node,elem] = squaremesh([0 1 0 1], 0.25);
%bdFlag = setboundary(node,elem,'Dirichlet','~(x==0)','Neumann','x==0');
bdFlag = setboundary(node,elem,'Dirichlet','all','Neumann','y==1');
%bdFlag = setboundary(node,elem,'Dirichlet');
for k = 1:3
    [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
%     [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
end

%% Get the data of the pde
pde = sincosdata;
% pde = mixBCdata;
option.solver = 'mg';

%% Finite Element Method
for k = 1:maxIt
    % solve the equation
    [u,Du,eqn] = PoissonWG(node,elem,pde,bdFlag,option);
    N(k) = size(elem,1);
    % compute error
    uI = zeros(N(k)+size(eqn.edge,1),1);
    uI(1:N(k)) = pde.exactu((node(elem(:,1),:)+node(elem(:,2),:)+node(elem(:,3),:))/3);
    uI(N(k)+1:end) = pde.exactu((node(eqn.edge(:,1),:)+node(eqn.edge(:,2),:))/2);
    erruIuh(k) = sqrt((u-uI)'*eqn.A*(u-uI));
    % refine mesh
   [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
%     [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
end

%% Plot convergence rates
figure(2);
showrate(N,erruIuh);
%% Error table
% The optimal rate of convergence of the H1-norm (1st order) and L2-norm
% (2nd order) is observed. The 2nd order convergent rate between two
% discrete functions ||DuI-Duh|| is known as superconvergence.

repo: Aero-matlab


function metrics = track_seq(seq_idx, seq_set, tracker, opt)

is_print = 1;   % set is_print to 1 to display the printed info

opt.is_print = is_print;
opt.exit_thre = 0.3;

if strcmp(seq_set, 'train') == 1
    seq_name = opt.mot2d_train_seqs{seq_idx};
    seq_len = opt.mot2d_train_nums(seq_idx);
else
    seq_name = opt.mot2d_test_seqs{seq_idx};
    seq_len = opt.mot2d_test_nums(seq_idx);
end

% mkdir img_traj
if exist('img_traj') ~= 7
    mkdir('img_traj');
end

if exist(['img_traj/' seq_name]) ~= 7
    mkdir(['img_traj/' seq_name]);
end

% read detections
det_file = fullfile(opt.mot, opt.mot2d, seq_set, seq_name, 'det', 'det.txt');
bboxes_det = read_bboxes(det_file);

if strcmp(seq_set, 'train') == 1
    % read ground truth
    gt_file = fullfile(opt.mot, opt.mot2d, seq_set, seq_name, 'gt', 'gt.txt');
    bboxes_gt = read_bboxes(gt_file);
    index_considered = find(bboxes_gt.r == 1);
    bboxes_gt = sub_bboxes(bboxes_gt, index_considered);
end

% initialize frame settings
frame_path = fullfile(opt.mot, opt.mot2d, seq_set, seq_name, 'img1', sprintf('%06d.jpg', 1));
frame_image = imread(frame_path);
frame_size.x = 1;
frame_size.y = 1;
frame_size.w = size(frame_image, 2);
frame_size.h = size(frame_image, 1);

% intialize tracker
tracker.image_width = frame_size.w;
tracker.image_height = frame_size.h;
tracker.max_width = max(bboxes_det.w);
tracker.max_height = max(bboxes_det.h);
tracker.max_score = max(bboxes_det.r);
tracker.min_score = min(bboxes_det.r);
tracker.num_tracked = 0;
tracker.fps = opt.fps(seq_idx);

% for each frame
trackers = [];
id = 0;
for fr = 1:seq_len
    if is_print
        fprintf('%s/%s frame %d\n', seq_set, seq_name, fr);
    else
        fprintf('.');
        if mod(fr, 100) == 0
            fprintf('\n');
        end
    end

    % read frame image
    frame_path = fullfile(opt.mot, opt.mot2d, seq_set, seq_name, 'img1', sprintf('%06d.jpg', fr));
    frame_image = imread(frame_path);

    % extract detection
    sub_idx = find(bboxes_det.fr == fr);
    bboxes = sub_bboxes(bboxes_det, sub_idx);

    % Detection preprocessing based on the scene context
    if strcmp(seq_set, 'test')
        sub_idx = [];
        for idx = 1:numel(bboxes.x)
            FP_flag = 0;
            % filter detections which are too large to be pedestrians
            if ~isempty(opt.h_max{seq_idx})
                if bboxes.h(idx) > opt.h_max{seq_idx}
                    FP_flag = 1;
                end
            end

            % filter detections which are too small to be pedestrians
            if ~isempty(opt.h_min{seq_idx})
                if bboxes.h(idx) < opt.h_min{seq_idx}
                    FP_flag = 1;
                end
            end

            % filter detections in impossible locations
            if ~isempty(opt.backgrounds{seq_idx})
                cx_det = bboxes.x(idx) + bboxes.w(idx) / 2; % center x of the detection
                cy_det = bboxes.y(idx) + bboxes.h(idx) / 2; % center y of the detection
                for bg_idx = 1:size(opt.backgrounds{seq_idx}, 1)
                    bg = opt.backgrounds{seq_idx}(bg_idx,:);
                    if cx_det > bg(1) && cx_det < bg(3) && cy_det > bg(2) && cy_det < bg(4)
                        FP_flag = 1;
                        continue;
                    end
                end
            end
            if FP_flag == 0
                sub_idx = [sub_idx, idx];
            end
        end
        bboxes = sub_bboxes(bboxes, sub_idx);
    end

    % nms
     boxes = [bboxes.x bboxes.y bboxes.x+bboxes.w bboxes.y+bboxes.h bboxes.r];
     sub_idx = nms(boxes, 0.6);
     bboxes = sub_bboxes(bboxes, sub_idx);

    % sort trackers
    [index1, index2] = sort_trackers(trackers, opt);
    index_processed = [];
    for k = 1:2
        % process trackers in the first class or the second class
        if k == 1
            index_track = index1;
        else
            index_track = index2;
        end

        % process trackers
        for i = 1:numel(index_track)
            ind = index_track(i);
            if trackers{ind}.state == opt.STATE_TRACKED || trackers{ind}.state == opt.STATE_ACTIVATED
                % track target
                [bboxes_tmp, index] = find_candidate_detections(trackers(index_processed), bboxes, opt);
                bboxes_eco_track = sub_bboxes(bboxes_tmp, index);
                trackers{ind} = MOT_track(fr, frame_image, frame_size, bboxes_eco_track, trackers{ind}, opt, seq_name);

                if trackers{ind}.state == opt.STATE_TRACKED || trackers{ind}.state == opt.STATE_ACTIVATED
                    index_processed = [index_processed; ind];
                end
            end
        end
        for i = 1:numel(index_track)
            ind = index_track(i);
            if trackers{ind}.state == opt.STATE_LOST
                % associate target
                [bboxes_tmp, index] = find_candidate_detections(trackers(index_processed), bboxes, opt);
                bboxes_associate = sub_bboxes(bboxes_tmp, index);
                trackers{ind} = MOT_associate(fr, frame_image, frame_size, bboxes_associate, trackers{ind}, opt, seq_name);
                index_processed = [index_processed; ind];
            end
        end
    end

    % find detections for initialization
    [bboxes, index] = find_candidate_detections(trackers, bboxes, opt);

    for i = 1:numel(index)
        % extract features
        bbox = sub_bboxes(bboxes, index(i));

        % filter detections using SVM
        det_feat = extract_det_feat(tracker, bbox);
        label = svmpredict(1, det_feat, tracker.w_active, '-q');
        if label < 0
           continue;
        end

        % reset tracker for the new object identity
        tracker.state = opt.STATE_START;
        id = id + 1;
        trackers{end+1} = initialize_tracker(fr, frame_image, id, bboxes, index(i), tracker, opt);
    end

    % resolve tracker conflict
    trackers = handle_conflicting_trackers(trackers, bboxes, opt);

    bboxes_track = [];
    for i = 1:numel(trackers)
        if isempty(bboxes_track)
            bboxes_track = trackers{i}.bboxes;
        else
            bboxes_track = concatenate_bboxes(bboxes_track, trackers{i}.bboxes);
        end
    end

    for i = 1:numel(trackers)
        if trackers{i}.state == opt.STATE_START || trackers{i}.state == opt.STATE_TRACKED || trackers{i}.state == opt.STATE_ACTIVATED
            bbox = sub_bboxes(trackers{i}.bboxes, numel(trackers{i}.bboxes.fr));
            x1 = floor(max(1, bbox.x));
            y1 = floor(max(1, bbox.y));
            x2 = ceil(min(frame_size.w, bbox.x+bbox.w-1));
            y2 = ceil(min(frame_size.h, bbox.y+bbox.h-1));
            img_traj = frame_image(y1:y2, x1:x2, :);
            if exist(['img_traj/' seq_name '/' num2str(trackers{i}.target_id)]) ~= 7
                mkdir(['img_traj/' seq_name '/' num2str(trackers{i}.target_id)]);
            end
            imwrite(img_traj, ['img_traj/' seq_name '/' num2str(trackers{i}.target_id) '/' num2str(fr) '.jpg']);
        end
    end
end

% write tracking results
file_name = sprintf('%s/%s.txt', opt.results_dir, seq_name);
fprintf('write results: %s\n', file_name);
save_results(seq_name, file_name, bboxes_track, tracker.fps*opt.initialization_thre, opt);


text:
function rs = minus(s1,s2)

%tstoolbox/@signal/minus
%   Syntax:
%     * rs=minus(s, offset)
%     * rs=minus(s1,s2)
%
%   Input arguments:
%     * s, s1, s2 - signal object
%     * offset - scalar value
%
%   Calculate difference of signals s1 and s2 or substract a scalar value
%   from s.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(2,2);

if isa(s2, 'signal')
	c = s1.core - s2.core;
%	d = merge(s1.description, s2.description);
	rs = signal(c, s1);
%	rs.description = d;
	rs = addhistory(rs,  ['Subtracted second signal from first signal']);
	rs = addcommandlines(rs, 's = minus(s, s2');
elseif isa(s2, 'double')
	if size(s2)==[1 1]
		c = core(data(s1)-s2);
		rs = signal(c, s1);
		rs = addhistory(rs,  ['Subtracted ' num2str(s2) ' from signal']);
		rs = addcommandlines(rs, 's = minus(s', s2);
	else
	 	error('second argument must be a signal or a scalar');
	end
else
	error('second argument must be a signal or a scalar');
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function value = r8_atanh ( x )

%*****************************************************************************80
%
%% R8_ATANH returns the inverse hyperbolic tangent of a number.
%
%  Discussion:
%
%    Y = R8_ATANH ( X )
%
%    implies that
%
%    X = TANH(Y) = ( EXP(Y) - EXP(-Y) ) / ( EXP(Y) + EXP(-Y) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the number whose inverse hyperbolic
%    tangent is desired.  The absolute value of X should be less than
%    or equal to 1.
%
%    Output, real R8_ATANH, the inverse hyperbolic tangent of X.
%
  if ( 1.0 <= abs ( x ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_ATANH - Fatal error!\n' );
    fprintf ( 1, '  ABS(X) must be < 1.\n' );
    fprintf ( 1, '  Your input is X = %f\n', x );
    error ( 'R8_ATANH - Fatal error!' );
  end

  value = 0.5 * log ( ( 1.0 + x ) / ( 1.0 - x ) );

  return
end


text:
#!/usr/bin/env octave
%% Machine Learning Online Class
%  Exercise 6 | Support Vector Machines
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  exercise. You will need to complete the following functions:
%
%     gaussianKernel.m
%     dataset3Params.m
%     processEmail.m
%     emailFeatures.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% =============== Part 1: Loading and Visualizing Data ================
%  We start the exercise by first loading and visualizing the dataset.
%  The following code will load the dataset into your environment and plot
%  the data.
%

fprintf('Loading and Visualizing Data ...\n')

% Load from ex6data1:
% You will have X, y in your environment
load('ex6data1.mat');

% Plot training data
plotData(X, y);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ==================== Part 2: Training Linear SVM ====================
%  The following code will train a linear SVM on the dataset and plot the
%  decision boundary learned.
%

% Load from ex6data1:
% You will have X, y in your environment
load('ex6data1.mat');

fprintf('\nTraining Linear SVM ...\n')

% You should try to change the C value below and see how the decision
% boundary varies (e.g., try C = 1000)
C = 1;
model = svmTrain(X, y, C, @linearKernel, 1e-3, 20);
visualizeBoundaryLinear(X, y, model);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =============== Part 3: Implementing Gaussian Kernel ===============
%  You will now implement the Gaussian kernel to use
%  with the SVM. You should complete the code in gaussianKernel.m
%
fprintf('\nEvaluating the Gaussian Kernel ...\n')

x1 = [1 2 1]; x2 = [0 4 -1]; sigma = 2;
sim = gaussianKernel(x1, x2, sigma);

fprintf(['Gaussian Kernel between x1 = [1; 2; 1], x2 = [0; 4; -1], sigma = 0.5 :' ...
         '\n\t%f\n(this value should be about 0.324652)\n'], sim);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =============== Part 4: Visualizing Dataset 2 ================
%  The following code will load the next dataset into your environment and
%  plot the data.
%

fprintf('Loading and Visualizing Data ...\n')

% Load from ex6data2:
% You will have X, y in your environment
load('ex6data2.mat');

% Plot training data
plotData(X, y);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ========== Part 5: Training SVM with RBF Kernel (Dataset 2) ==========
%  After you have implemented the kernel, we can now use it to train the
%  SVM classifier.
%
fprintf('\nTraining SVM with RBF Kernel (this may take 1 to 2 minutes) ...\n');

% Load from ex6data2:
% You will have X, y in your environment
load('ex6data2.mat');

% SVM Parameters
C = 1; sigma = 0.1;

% We set the tolerance and max_passes lower here so that the code will run
% faster. However, in practice, you will want to run the training to
% convergence.
model= svmTrain(X, y, C, @(x1, x2) gaussianKernel(x1, x2, sigma));
visualizeBoundary(X, y, model);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =============== Part 6: Visualizing Dataset 3 ================
%  The following code will load the next dataset into your environment and
%  plot the data.
%

fprintf('Loading and Visualizing Data ...\n')

% Load from ex6data3:
% You will have X, y in your environment
load('ex6data3.mat');

% Plot training data
plotData(X, y);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ========== Part 7: Training SVM with RBF Kernel (Dataset 3) ==========

%  This is a different dataset that you can use to experiment with. Try
%  different values of C and sigma here.
%

% Load from ex6data3:
% You will have X, y in your environment
load('ex6data3.mat');

% Try different SVM Parameters here
[C, sigma] = dataset3Params(X, y, Xval, yval);

% Train the SVM
model= svmTrain(X, y, C, @(x1, x2) gaussianKernel(x1, x2, sigma));
visualizeBoundary(X, y, model);

fprintf('Program paused. Press enter to continue.\n');
pause;

author: SaveTheRbtz

repo: ml-class


text:
function [trh, trv, tE] = trfd(h, v)
% Compute "time"-reversed horizontal and vertical finite differences on the
% current estimate x
%
% Use periodic boundaries

%% Application of "time"-reversed horizontal finite difference
tS = tic;

trh = [h(:,end) h];
trh = imfilter(trh, [1 -1], 'circular');
trh = trh(:, 1:end-1);

%% Application of "time"-reversed vertical finite difference
trv = [v(end,:); v];
trv = imfilter(trv, [1 -1]', 'circular');
trv = trv(1:end-1, :);

tE = toc(tS);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function model = xls2model(fileName, biomassRxnEquation, defaultbound)
% Reads a model from Excel spreadsheet.
%
% USAGE:
%
%    model = xls2model(fileName, biomassRxnEquation, defaultbound)
%
% INPUT:
%    fileName:              xls spreadsheet, with one 'Reaction List' and one 'Metabolite List' tab
%
% OPTIONAL INPUTS:
%    biomassRxnEquation:    .xls may have a 255 character limit on each cell,
%                           so pass the biomass reaction separately if it hits this maximum.
%
%    defaultbound:          the deault bound for lower and upper bounds, if
%                           no bounds are specified in the Excel sheet
% OUTPUT:
%    model:                 COBRA Toolbox model
%
% EXAMPLE:
%
%                   'Reaction List' tab headers (case sensitive):
%
%                     * Required:
%
%                       * 'Abbreviation':      HEX1
%                       * 'Reaction':          `1 atp[c] + 1 glc-D[c] --> 1 adp[c] + 1 g6p[c] + 1 h[c]`
%
%                     * Optional:
%
%                       * 'GPR':               (3098.3) or (80201.1) or (2645.3) or ...
%                       * 'Description':       Hexokinase
%                       * 'Subsystem':         Glycolysis
%                       * 'Reversible':        0 (false) or 1 (true)
%                       * 'Lower bound':       0
%                       * 'Upper bound':       1000
%                       * 'Objective':         0/1
%                       * 'Confidence Score':  0,1,2,3,4
%                       * 'EC Number':         2.7.1.1,2.7.1.2
%                       * 'KEGG ID':           R000001
%                       * 'Notes':             Reaction also associated with EC 2.7.1.2
%                       * 'References':        PMID:2043117,PMID:7150652,...
%
%                   'Metabolite List' tab: Required headers (case sensitive): (needs to be complete list of metabolites,
%                   i.e., if a metabolite appears in multiple compartments it has to be represented in multiple rows.
%                   Abbreviations need to overlap with use in Reaction List
%
%                     * Required
%
%                       * 'Abbreviation':      glc-D or glc-D[c]
%                     * Optional:
%
%                       * 'Charged formula' or formula:   C6H12O6
%                       * 'Charge':                       0
%                       * 'Compartment':                  cytosol
%                       * 'Description':                  D-glucose
%                       * 'KEGG ID':                      C00031
%                       * 'PubChem ID':                   5793
%                       * 'ChEBI ID':                     4167
%                       * 'InChI string':                 InChI=1/C6H12O6/c7-1-2-3(8)4(9)5(10)6(11)12-2/h2-11H,1H2/t2-,3-,4+,5-,6?/m1/s1
%                       * 'SMILES':                       OC[C@H]1OC(O)[C@H](O)[C@@H](O)[C@@H]1O
%                       * 'HMDB ID':                      HMDB00122
%
% NOTE:
%
%    Optional inputs may be required for input on unix machines.
%
% NOTE:
%
%    Find an example Excel sheet at `docs/source/examples/ExcelExample.xlsx`
%
% .. Authors:
%    - Ines Thiele, 01/02/09
%    - Richard Que, 04/27/10, Modified reading of PubChemID and ChEBIID so that if met
%      has multiple IDs, all are passed to model. Confidence Scores
%      PubChemIDs, and ChEBIIDs, are properly passed as cell arrays.
%    - Ronan Fleming, 08/17/10, Support for unix
%    - Hulda S.H., 10/11/10, Modified reading of xls document.
%      Identifies columns by their headers. Added reading of HMDB ID.

warning off

% test if Excel is available
excelInstalled = false;
try
    excelObj = actxserver('Excel.Application');
    excelInstalled = true;
    %h.WorkBooks.Item(fileName).Close;
    fprintf(' > Excel is installed.\n\n');
catch ME
    fprintf(' > Excel is not installed; Using xlread.\n\n');
end

if exist(fileName,'file') == 2
    try
        try
            [~, sheets] = xlsfinfo(fileName);
        catch ME
            if strcmp(ME.identifier,'MATLAB:xlsread:FileDoesNotExist')
                [~, sheets] = xlsfinfo(fullfile(pwd, fileName));
            else
                error(ME);
            end
        end
        if ~all(ismember({'Reaction List','Metabolite List'},sheets))
            error(['The provided Excel Sheet (', fileName,') must contain a "Reaction List" and a "Metabolite List sheet as specified here:' sprintf('\n'),...
                   'https://opencobra.github.io/cobratoolbox/docs/ExcelModelFileDefinition.html']);
        end
    catch ME
        fprintf(strrep(sheets, '\', '\\'));
    end
else
    error('File %s not found',fileName);
end

if ~exist('defaultbound','var')
    defaultbound = 1000;
end

%assumes that one has an xls file with two tabs
if isunix || ~excelInstalled
    [~, Strings, rxnInfo] = xlread(fileName,'Reaction List');
    [~, MetStrings, metInfo] = xlread(fileName,'Metabolite List');
else
    [~, Strings, rxnInfo] = xlsread(fileName, 'Reaction List');
    [~, MetStrings, metInfo] = xlsread(fileName, 'Metabolite List');
end

%trim empty row from Numbers and MetNumbers
rxnInfo = rxnInfo(1:size(Strings,1),:);
metInfo = metInfo(1:size(MetStrings,1),:);

if isunix && isempty(MetStrings)
    error('Save .xls file as Windows 95 version using gnumeric not openoffice!');
end

requiredRxnHeaders = {'Abbreviation','Reaction'};
requiredMetHeaders = {'Abbreviation'};

if ~all(ismember(requiredRxnHeaders,Strings(1,:)))
    error(['Required Headers not present in the "Reaction List" sheet of the provided xls file.', sprintf('\n'),...
           'Note, that headers are case sesnitive!', sprintf('\n'),...
           'Another likely source for this issue is a change in the xls format specification.', sprintf('\n'),...
           'Please have a look at the specification at https://opencobra.github.io/cobratoolbox/docs/ExcelModelFileDefinition.html for the current specifications.']);
end

if ~all(ismember(requiredMetHeaders,MetStrings(1,:)))
    error(['Required Headers not present in the "Metabolite List" sheet of the provided xls file.', sprintf('\n'), ...
           'Note, that headers are case sesnitive!', sprintf('\n'),...
           'Another likely source for this issue is a change in the xls format specification.', sprintf('\n'),...
           'Please have a look at the specification at https://opencobra.github.io/cobratoolbox/docs/ExcelModelFileDefinition.html for the current specifications.']);
end

rxnHeaders = rxnInfo(1,:);

for n = 1:length(rxnHeaders)
    if isnan(rxnHeaders{n})
        rxnHeaders{n} = '';
    end
end

% Assuming first row is header row
rxnAbrList = Strings(2:end,strmatch('Abbreviation',rxnHeaders,'exact'));
if ~isempty(strmatch('Description',rxnHeaders,'exact'))
    rxnNameList = Strings(2:end,strmatch('Description',rxnHeaders,'exact'));
else
    rxnNameList = Strings(2:end,strmatch('Abbreviation',rxnHeaders,'exact'));
end
rxnList = Strings(2:end,strmatch('Reaction',rxnHeaders,'exact'));
if ~isempty(strmatch('GPR',rxnHeaders,'exact'))
    grRuleList = Strings(2:end,strmatch('GPR',rxnHeaders,'exact'));
else
    grRuleList = cell(size(rxnList,1),1);
    grRuleList(:) = {''};
end

if ~isempty(strmatch('Proteins',rxnHeaders,'exact'))
    Protein = Strings(2:end,strmatch('Proteins',rxnHeaders,'exact'));
end

if ~isempty(strmatch('Subsystem',rxnHeaders,'exact'))
    subSystemList = Strings(2:end,strmatch('Subsystem',rxnHeaders,'exact'));
    subSystemList = cellfun(@(x) strsplit(x,';'), subSystemList ,'UniformOutput',0);
else
    subSystemList = cell(size(rxnList,1),1);
    subSystemList(:) = {{''}};
end

% initialization with default values
lowerBoundList = -defaultbound*ones(length(rxnAbrList),1);

if ~isempty(strmatch('Lower bound',rxnHeaders,'exact'))
    tmp = rxnInfo(2:end,strmatch('Lower bound',rxnHeaders,'exact'));
    for i = 1:length(tmp)
        if isnumeric(tmp{i})
            lowerBoundList(i) = tmp{i};
        else
            lowerBoundList(i) = str2num(tmp{i});
        end
    end
    lowerBoundList = columnVector(lowerBoundList); %Default -1000
    lowerBoundList(isnan(lowerBoundList)) = -defaultbound;
end

% initialization with default values
upperBoundList = defaultbound*ones(length(rxnAbrList),1);

if ~isempty(strmatch('Upper bound',rxnHeaders,'exact'))
    tmp = rxnInfo(2:end,strmatch('Upper bound',rxnHeaders,'exact'));
    for i = 1:length(tmp)
        if isnumeric(tmp{i})
            upperBoundList(i) = tmp{i};
        else
            upperBoundList(i) = str2num(tmp{i});
        end
    end
    upperBoundList = columnVector(upperBoundList); %Default 1000;
    upperBoundList(isnan(upperBoundList)) = defaultbound;
end

revFlagList = lowerBoundList<0;

% initialization with default values
Objective = zeros(length(rxnAbrList),1);

if ~isempty(strmatch('Objective',rxnHeaders,'exact'))
    tmp = rxnInfo(2:end,strmatch('Objective',rxnHeaders,'exact'));
    for i = 1:length(tmp)
        if isnumeric(tmp{i})
            Objective(i) = tmp{i};
        else
            Objective(i) = str2num(tmp{i});
        end
    end
    Objective = columnVector(Objective);
    Objective(isnan(Objective)) = 0;
end

model = createModel(rxnAbrList,rxnNameList,rxnList,revFlagList,lowerBoundList,upperBoundList,subSystemList,grRuleList);

if ~isempty(strmatch('Confidence Score',rxnHeaders,'exact'))
    confidenceScores = rxnInfo(2:end,strmatch('Confidence Score',rxnHeaders,'exact'));
    if any(~cellfun(@isnumeric, confidenceScores))
        %replace any non numeric elements by their numbers.
        confidenceScores = cellfun(@str2num,confidenceScores,'UniformOutput',0);
    end
    confidenceScores(cellfun(@isempty,confidenceScores)) = {NaN}; %Replace empty by nan.
    model.rxnConfidenceScores = cell2mat(confidenceScores);
    model.rxnConfidenceScores(isnan(model.rxnConfidenceScores)) = 0;
end
if ~isempty(strmatch('EC Number',rxnHeaders,'exact'))
    %This needs to be changed to the new annotation scheme and putting the
    %ECNumbers there.
    model.rxnECNumbers = Strings(2:end,strmatch('EC Number',rxnHeaders,'exact'));
end
if ~isempty(strmatch('Notes',rxnHeaders,'exact'))
    model.rxnNotes = Strings(2:end,strmatch('Notes',rxnHeaders,'exact'));
end
if ~isempty(strmatch('References',rxnHeaders,'exact'))
    model.rxnReferences = Strings(2:end,strmatch('References',rxnHeaders,'exact'));
    numbers = cellfun(@isnumeric ,model.rxnReferences);
    model.rxnReferences(numbers) = cellfun(@convertNumberToID , model.rxnReferences(numbers),'UniformOutput',0);
    model.rxnReferences = cellfun(@(x) regexprep(x,'PMID:',''), model.rxnReferences,'UniformOutput',0);
end

%fill in opt info for metabolites
if ~isempty(Objective) && length(Objective) == length(model.rxns)
    model.c = Objective;
end

metHeaders = metInfo(1,:);

for n = 1:length(metHeaders)
    if isnan(metHeaders{n})
        metHeaders{n} = '';
    end
end

% case 1: all metabolites in List have a compartment assignement

metCol = strmatch('Abbreviation',metHeaders,'exact');
Compartments = {};
mets = MetStrings(:,metCol);
%Now, we could have a problem, if the reactions are presented without
%compartments. In this instance, we would have to first put a "[c]" id
%behind all metabolites.
metCompAbbrev = cellfun(@(x) regexp(x,'.*\[(.*)\]$','tokens'), mets, 'UniformOutput', 0);
%get those which don't have a compartmentID
noncomps = cellfun(@isempty, metCompAbbrev);
mets(noncomps) = strcat(mets(noncomps),'[c]');

[A,B] = ismember(model.mets,mets);
matchingmets = mets(B(A));

if numel(matchingmets) ~= numel(model.mets)
    fprintf('The following metabolites from the reaction formulas did not have a matching metabolite in the metabolite list:\n');
    disp(setdiff(model.mets,mets));
    error('Not all metabolites could be matched');
end

if isempty(strmatch('Compartment',metHeaders,'exact'))
    %we use default compartments
    [compartmentAbbr,compartments] = getDefaultCompartmentSymbols();
    %lets check if all metabolites do have a standard compartment
    metCompAbbrev = cellfun(@(x) regexp(x,['.*\[(' strjoin(compartmentAbbr,'|') '\]$'],'tokens'), matchingmets, 'UniformOutput', 0);
    noncomps = cellfun(@isempty, metCompAbbrev);
    if any(noncomps)
        %So, there are missing compartment ids.
        %lets move all those metabolites to the cytosol, checking, that we
        %don't generate replicates.
        matchingmets(noncomps) = strcat(matchingmets(noncomps),'[c]');
        if numel(unique(matchingmets)) ~= numel(matchingmets)
            [~,ia] = unique(matchingmets);
            non_unique = matchingmets(setdiff(1:numel(matchinmets),ia));
            disp(unique(non_unique))
            error(['The above metabolites are present both without compartment identifier and with id in the cytosol.\n', ...
                  'Metabolites without compartment id are assumed to be located in the cytosol, and these metabolites would lead to duplicate metabolite ids!']);
        end
    end
    %Now, there should be no metabolites without compartment.
    %lets collect the compartments.
    Comps = cellfun(@(x) x{1} ,cellfun(@(x) regexp(x,['.*\[([' strjoin(compartmentAbbr,'') '])\]$'],'tokens'), matchingmets));
    %matchingmets has already the right order.
    model.comps = columnVector(compartmentAbbr(ismember(compartmentAbbr,Comps)));
    model.compNames = columnVector(compartments(ismember(compartmentAbbr,Comps)));
    model.mets = columnVector(matchingmets);
else
    %if Compartments is present, we will create a translation table
    %(ignoring everything that is empty)
    Compartments = MetStrings(B(A),strmatch('Compartment',metHeaders,'exact'));
    Cytosolname = 'cytosol';
    metCompAbbrev = cellfun(@(x) regexp(x,['.*\[(.*)\]$'],'tokens'), matchingmets, 'UniformOutput', 0);
    noncomps = cellfun(@isempty, metCompAbbrev);
    CytoNames = setdiff(unique(Compartments(~cellfun(@isempty, cellfun(@(x) regexp(x,'.*\[(c)\]$','tokens'), matchingmets, 'UniformOutput', 0)))),'');
    if numel(CytoNames) == 1
        Cytosolname = CytoNames{1};
    else
        CytoNames{end+1} = Cytosolname;
        CytoNames = unique(CytoNames);
    end

    if any(noncomps)
        matchingmets(noncomps) = strcat(matchingmets(noncomps),'[c]');
        Compartments(noncomps) = {Cytosolname};
        if numel(unique(matchingmets)) ~= numel(matchingmets)
            [~,ia] = unique(matchingmets);
            non_unique = matchingmets(setdiff(1:numel(matchinmets),ia));
            disp(unique(non_unique))
            error(['The above metabolites are present both without compartment identifier and with id in the cytosol.\n',...
                  'Metabolites without compartment id are assumed to be located in the cytosol, and these metabolites would lead to duplicate metabolite ids!']);
        end
    end
    metCompAbbrev = cellfun(@(x) x{1}, cellfun(@(x) regexp(x,['.*\[(.*)\]$'],'tokens'), matchingmets),'UniformOutput',false);

    %now reorder them and assign names to the ids.
    [ucomps, origpos] = unique(Compartments);
    [model.comps,~,origin] = unique(metCompAbbrev(origpos));
    %Column Vector
    model.comps = columnVector(model.comps);
    if ischar(model.comps)
        model.comps = cellstr(model.comps);
    end
    for i = 1:numel(model.comps)
        %combine all, ignoring empty entries.
        CompNames{i} = strjoin(setdiff(ucomps(origin==i),''),' or ');
    end
    model.compNames = columnVector(CompNames);
    model.mets = columnVector(matchingmets);
end
%%Set metNames
if ~isempty(strmatch('Description',metHeaders,'exact'))
    model.metNames = columnVector(MetStrings(B(A),strmatch('Description',metHeaders,'exact')));
end
%%Set Formulas
if ~isempty(strmatch('Charged formula',metHeaders,'exact'))
    model.metFormulas = columnVector(MetStrings(B(A),strmatch('Charged formula',metHeaders,'exact')));
end
if ~isempty(strmatch('Formula',metHeaders,'exact'))
    model.metFormulas = columnVector(MetStrings(B(A),strmatch('Formula',metHeaders,'exact')));
end
%%Set Charge
if ~isempty(strmatch('Charge',metHeaders,'exact'))
    model.metCharges = cell2mat(columnVector(metInfo(B(A),strmatch('Charge',metHeaders,'exact'))));
end

if ~isempty(strmatch('SMILES',metHeaders,'exact'))
    model.metSmiles= columnVector(MetStrings(B(A),strmatch('SMILES',metHeaders,'exact')));
end

%% Set annotations. (Has to be updated, once annotation structure is defined)
if ~isempty(strmatch('KEGG ID',metHeaders,'exact'))
    model.metKEGGID = columnVector(MetStrings(B(A),strmatch('KEGG ID',metHeaders,'exact')));
end
if ~isempty(strmatch('InChI string',metHeaders,'exact'))
    model.metInChIString = columnVector(MetStrings(B(A),strmatch('InChI string',metHeaders,'exact')));
end
if ~isempty(strmatch('HMDB ID',metHeaders,'exact'))
    model.metHMDBID = columnVector(MetStrings(B(A),strmatch('HMDB ID',metHeaders,'exact')));
end

if ~isempty(strmatch('PubChem ID',metHeaders,'exact'))
    %This is a litte trickier, as PubChemIDs are numbers. So we have to
    %load them differently
    model.metPubChemID = columnVector(metInfo(B(A),strmatch('PubChem ID',metHeaders,'exact')));
    numbers = cellfun(@isnumeric ,model.metPubChemID);
    model.metPubChemID(numbers) = cellfun(@convertNumberToID , model.metPubChemID(numbers),'UniformOutput',0);
end
if ~isempty(strmatch('ChEBI ID',metHeaders,'exact'))
    model.metChEBIID  = columnVector(metInfo(B(A),strmatch('ChEBI ID',metHeaders,'exact')));
    numbers = cellfun(@isnumeric ,model.metChEBIID);
    model.metChEBIID(numbers) = cellfun(@convertNumberToID , model.metChEBIID(numbers),'UniformOutput',0);
end

[~,fileName,extension] = fileparts(fileName);

model.description = [fileName, extension];

warning on
end

function stringNumber = convertNumberToID(number)
if isnan(number)
    stringNumber = '';
else
    stringNumber = num2str(number);
end
end

author: opencobra

repo: cobratoolbox


text:
% DEMO3B - Chaining files with ncdataset and finding attributes using 'value4key'

echo('on')
o2_ = 'Oxygen';
time_ = 'esecs';

%% ---- Access the OpenDAP datasets
m2006 = ncdataset('http://dods.mbari.org/cgi-bin/nph-nc/data/ssdsdata/deployments/m1/200610/m1_aanderaaoxy_20070105.nc');
m2007 = ncdataset('http://dods.mbari.org/cgi-bin/nph-nc/data/ssdsdata/deployments/m1/200711/m1_aanderaaoxy_20071106.nc');
m2008 = ncdataset('http://dods.mbari.org/cgi-bin/nph-nc/data/ssdsdata/deployments/m1/200810/m1_aanderaaoxy_20081008.nc');
t = [m2006.time(time_); m2007.time(time_); m2008.time(time_)];
o2 = [m2006.data(o2_); m2007.data(o2_); m2008.data(o2_)];

%% ---- Find the units of Oxygen and label the display
attr = m2008.attributes(o2_);
units = value4key(attr, 'units'); % Retrieve the units value
name = value4key(attr, 'long_name'); % Retrieve the long_name value

%% ---- Plot the data
figure;
plot(t, o2);...
datetick('x');...
ylabel([name ' [' units ']']);...
title({'Surface Temperature at M1 Mooring in Monterey Bay',m2006.location,m2007.location,m2008.location},'interpreter','none');...
grid;...
shg
echo('off')

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
% OtQueryNext Move to the next record of the query.
%
% Usage:
%    OtQueryNext(hquery)

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% 2006-02-28 Created (Johan Meijdam, VORtech)

function iret = OtQueryNext(hquery)
   iret = Ot('query.next', hquery);
end

author: Sable

repo: mcbench-benchmarks


text:
function [gcX, gcY, gcZ] = get_gc_initial_pos()
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

    % set the initial position in the grid cell network
    global GC_X_DIM;
    global GC_Y_DIM;
    global GC_Z_DIM;

    gcX = floor(GC_X_DIM / 2);  % in 1:36
    gcY = floor(GC_Y_DIM / 2);  % in 1:36
    gcZ = floor(GC_Z_DIM / 2);  % in 1:36
%     gcX = 1;  % in 1:36
%     gcY = 1;  % in 1:36
%     gcZ = 1;  % in 1:36
end

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
function J = computeCost(X, y, theta)
%COMPUTECOST Compute cost for linear regression
%   J = COMPUTECOST(X, y, theta) computes the cost of using theta as the
%   parameter for linear regression to fit the data points in X and y

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta
%               You should set J to the cost.

predictions = X * theta;
sqrErrors = (predictions - y).^2;

J = 1 / (2*m) * sum(sqrErrors);

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

end

author: merwan

save_path: github-repos/MATLAB/merwan-ml-class


function [ epsilon_soft, time_sum, compressIndex  ] = Uniform_loopCompress( lamda, qFile, filesDir, cut_num, totalNum, bValue, saveDir )
%LOOPCOMPRESS Summary of this function goes here
%   Detailed explanation goes here

    fileExt = '*.txt';
    files = dir(fullfile(filesDir,fileExt));

    splitLength =  floor(length(files) / cut_num);

    q_file_t = fopen(qFile);
    observe_value = fscanf(q_file_t, '%d');
    fclose(q_file_t);
    observe_value(find(observe_value > 100)) = 100;  % remove some unusual counts
    q_value = (1 - mapminmax(observe_value', 0, 1))';  % actually same as in CVPR_W 2013

    time_sum = 0;
    saveCnt = 0;
    epsilon_soft = [];

    for i=0:splitLength:length(files)-1
        disp('-----------------------------------------------------------------------');
        start = i;
        finish = i + splitLength - 1;
        % for the last section
        if finish > (length(files) - splitLength) % last section will be longer than others
            finish = length(files) - 1;
        end
        disp(start);
        disp(finish);

        % core
        [compressIndex, time_sec, epsilon_part] = section_compress(lamda, q_value, filesDir, start, finish, totalNum, bValue);
        time_sum = time_sum + time_sec;
        saveName = [saveDir, num2str(saveCnt), '.txt'];
        dlmwrite(saveName, compressIndex, 'precision', '%d');

        % soft part
        epsilon_soft  = [epsilon_soft;epsilon_part];

        disp('Compressed & Saved')
        saveCnt = saveCnt + 1;

        % added fucked!
        if finish == length(files) - 1
            break; % jump out the fucking looping
        end

    end

    disp('All time of programming:');
    disp(time_sum)

end


text:
function mv = mv_exportAmplitudeMap(mv);
%
% mv = mv_exportAmplitudeMap(mv);
%
% Create multiple parameter maps, with the amplitude for
% each selected condition in a multivoxel UI.
%
%
% ras, 09/2006.
if notDefined('mv'), mv = get(gcf, 'UserData'); end

% get amplitudes for all selected conditions
amps = mv_amps(mv);
nConds = size(amps, 2);
selConds = find(tc_selectedConds(mv));

% initialize a view to get parameters about the
% map volume
fn = sprintf('getSelected%s',mv.roi.viewType);
view = eval(fn);
if isempty(view)
    % no selected view of the proper type -- make a hidden one
    mrGlobals; loadSession; saveFlag = 1;
    fn = sprintf('initHidden%s',roi.viewType); view = eval(fn);
end
mapdims = viewGet(view,'dataSize');
nScans = viewGet(view,'numScans');
scan = mv.params.scans(1);

for c = 1:nConds
    % plug in the values to the map volume:
    mapvol = zeros(mapdims);
    ind = roiIndices(view,mv.coords);
    mapvol(ind) = amps(:,c);

    % save map
    mapName = sprintf('Amplitudes_%s_%s_%s', mv.params.ampType, ...
                       mv.trials.condNames{selConds(c)}, mv.roi.name);
    mapPath = fullfile(dataDir(view),mapName);
    if exist(mapPath,'dir')
        load(mapPath,'map','mapName');
        map{scan} = mapvol;
    else
        map = cell(1, nScans);
    end
    map{scan} = mapvol;
    save(mapPath,'map','mapName');
    fprintf('Saved map %s.\n',mapPath);
end

% plug the map volume into the view
fprintf('Assigning %s data to param map for scan %i\n',mapName,scan);
map = cell(1,nScans);
map{scan} = mapvol;
if ~isequal(view.name,'hidden')
    view = setParameterMap(view,map,mapName);
    refreshScreen(view);
end

% evaluate this in the workspace, so the view
% itself is updated
assignin('base','map',map);
evalin('base',sprintf('%s=setParameterMap(%s,map,''%s'');',...
    view.name,view.name,mapName));

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% Gcone_test.m
% test the Gcone object
% Copyright 2008-1-1, Jeff Fessler, University of Michigan

printm 'todo: check user_source_zs vs hct2 with file_source_z'

pn = jf_protected_names;

f.sf1 = {'sf1', ...
	'sf1,p:ts,b:st', 'sf1,p:st,b:st', ...
	'sf1,p:ts,b:ts', 'sf1,p:st,b:ts'};
f.sf2 = {'sf2', ...
	'sf2,p:ts,b:st', 'sf2,p:st,b:st', ...
	'sf2,p:ts,b:ts', 'sf2,p:st,b:ts'};
systypes = { ... % lots of variations on system models!
	'rf1', ... % RR model
	f.sf1{:}, f.sf2{:}, ... % SF variations
	'sf1', 'sf2', 'sf3', 'sf4', 'sf5' ...
};
%systypes = {'rf1', 'sf1', 'dd2'}
if 0 && exist('dd_ge1_mex') == 3 % UM only
	systypes{end+1} = 'dd1';
	systypes{end+1} = 'dd2'; % todo
end
% systypes = {'nn1', 'pd1'};
%systypes = {'dd2'}; % todo! fails adjoint test!?
%systypes = {'sf1', 'sf2'};
%systypes = f.sf2;
nsys = numel(systypes);

%f.class = 'Fatrix';
f.class = 'fatrix2';

% todo: test both ways for fatrix2, not for Fatrix
%is_zxy = false;
is_zxy = true;
if is_zxy
	to_zxy = @(x) permute(x, [3 1 2]);
else
	to_zxy = @(x) x;
end

% small systems for basic tests
if 1 && (0 || ~isvar('A1')), printm 'setup small'
	f.down = 16;

	dfs_list = [inf inf 0]; % parallel flat arc
	dsd_list = [inf 949.075 949.075]; % parallel flat arc

	if 0 % arc only
		dfs_list = [0];
		dsd_list = [949.075];
	end

	for kk = 1:numel(dfs_list)

		cgs = ct_geom('ge2', 'nt', 320, ...
			'source_z0', -20, 'pitch', 0.5, ... % test helix
			'dfs', dfs_list(kk), ... % arc or flat
			'dsd', dsd_list(kk), ... % fan or parallel beam
			'down', f.down);

		for ii=1:nsys
			systype = systypes{ii};
			printm('testing type %s dfs=%g dsd=%g', ...
				systype, cgs.dfs, cgs.dsd)

			if streq(systype, 'dd1') || streq(systype, 'dd2')
				f.dy = 1; % DD requires square pixels
				if isinf(cgs.dsd) % no parallel for DD
					continue
				end
			else
				f.dy = 0.7; % stress test SF with non-square
			%	f.dy = 'dx';
			end

			igs = image_geom('nx', 512, 'ny', 480, 'nz', 416, ...
				'dx', 1, 'dy', f.dy, 'dz', 0.5, ...
				'offset_x', 2.9, 'offset_y', 3.5, ...
				'offset_z', -3.4, ...
				'mask', 'all-but-edge', ... % trick: for hct2
				'down', f.down);

			if im, cgs.plot(igs); end

			args = {cgs, igs, 'type', systype, 'class', f.class};
			if streq(systype, 'sf4') || streq(systype, 'sf5')
				args = {args{:}, 'mexarg', {single(0)}};
			end

			args = {args{:}, 'zxy', is_zxy};

			A1 = Gcone(args{:}, 'nthread', 1);
			Ac = Gcone(args{:});

			if 1 % adjoint
				tester_tomo2(A1, to_zxy(igs.mask), 'G2', Ac)
				test_adjoint(A1, 'big', 1, 'tol', 5e-5)
				test_adjoint(Ac, 'big', 1, 'tol', 5e-5)
			end

			if streq(systype, 'dd1') || streq(systype, 'dd2')
				continue % don't bother hct2 test for DD
			end

			% hereafter is hct2 test
			Ah = Gcone(args{:}, 'use_hct2', 1);

			if pn.has_hct2
				if streq(systype, 'nn1') || streq(systype, 'pd1')
					thresh = 3e-2; % big due to rounding in nn1
				else
					thresh = 4e-6;
				end

				if 1 % mex vs hct2
					xs = single(igs.mask);
					xs(round(end/3), round(2*end/3), round(end/4)) = 10;
					xs = to_zxy(xs);
					t1 = A1 * xs;
					t2 = Ah * xs;
					equivs(t1, t2, 'thresh', thresh)

					b1 = A1' * t1;
					b2 = Ah' * t1;
%					im plc 1 3, im(1, b1), im(2, b2), im(3, b2-b1)
					equivs(b1, b2, 'thresh', thresh)
				end
				if 0 % block
					B1 = Gblock(A1, 2);
					Bh = Gblock(Ah, 2);
					t1 = B1{2} * xs;
					t2 = Bh{2} * xs;
					equivs(t1, t2, 'thresh', thresh)
				end

				tester_tomo2(Ah, to_zxy(igs.mask), ...
					'equiv_thresh', thresh) % paces
				test_adjoint(Ah, 'big', 1, 'tol', 5e-5)
			end
		end % systype
	end % dfs
end % small

if 0, printm 'rf1' % test rf1 vs dd2 small
	A0 = Gcone(cgs, igs, 'type', 'rf1');
	Ad = Gcone(cgs, igs, 'type', 'dd2');
%	x0 = single(igs.mask);
	x0 = igs.circ;
	im(x0)
	cpu etic
	yd = Ad * x0;
	cpu etoc dd
	if 1
		cpu etic
		y0 = A0 * x0;
		cpu etoc rf1
	end
	im plc 2 2
	im(yd), im(y0), im(y0-yd)
return
end

if ~isvar('x0'), printm 'x0 big'
	f.down = 4;
	igb = image_geom('nx', 512, 'ny', 480, 'nz', 416, ...
		'dx', 1, 'dz', 0.5, ...
		'dy', 0.7, ... % stress test
		'offset_x', 12.9, 'offset_y', 3.5, 'offset_z', -3.4, ...
		'down', f.down);
	ell = [3*igb.dx 5*igb.dx -2*igb.dz ...
		igb.dx*igb.nx/3 igb.dy*igb.ny/4 igb.zfov/4 ...
		0 0 10];
	x0 = ellipsoid_im(igb, ell, 'oversample', 2);
%	clf, im(x0), return
%prompt
end

% big systems for accuracy tests
if 0 || ~isvar('Ab'), printm 'setup big'
	cgb = ct_geom('ge1', 'nt', 320, ...
		'source_z0', -40, 'pitch', 0.5, ... % test helix
	...%	'dfs', inf, ... % flat detector
	...%	'dsd', inf, 'dfs', inf, ... % parallel beam
		'down', f.down);

	clear Ab Ah
	for ii=1:nsys
		systype = systypes{ii};
		if streq(systype, 'dd', 2) && isinf(cgb.dsd)
			Ab{ii} = [];
			continue
		end
		Ab{ii} = Gcone(cgb, igb, 'type', systype, ...
				'zxy', is_zxy, 'class', f.class);
	end
end

if ~isvar('ya'), printm 'analytical projections'
	ya = ellipsoid_proj(cgb, ell, 'oversample', 2);
%	im clf, im(ya)
%prompt
end

if 0 || ~isvar('yb'), printm 'discrete projections'
	nrmse = @(x,y) norm(y(:)-x(:)) / norm(x(:)) * 100;
	for ii=1:nsys
		if ~isempty(Ab{ii})
			cpu etic
			yb{ii} = Ab{ii} * to_zxy(x0);
			f.time(ii) = cpu('etoc');
			printm('%s: time %g nrmse %g %%', ...
				systypes{ii}, f.time(ii), nrmse(ya, yb{ii}))
		else
			f.time(ii) = 0;
			yb{ii} = cgb.zeros;
		end
	end
end

if 0 % toggle analytical vs big discrete
	ii = 1;
	im_toggle(ya(:,end/2,:), yb{ii}(:,end/2,:), ...
		yb{ii}(:,end/2,:) - ya(:,end/2,:))
	im(ya(:,end/2,:) - yb{ii}(:,end/2,:))
return
end

if 0, printm 'look at error in worst views'
	ilist = 1:nsys;
	ilist = [1 2 7];
	im('plc', 2, numel(ilist))
	for jj=1:numel(ilist)
		ii = ilist(jj);
		err = yb{ii} - ya;
		tmp = reshape(err, [], cgb.na);
		tmp = sum(tmp.^2); % error in each view
		ia = imax(tmp); % worst view
		im(jj, err(:,:,ia)), cbar h
		titlef('%s ia=%d', systypes{ii}, ia)
		im('subplot', jj+numel(ilist))
		plot(tmp), axis tight
	end
return
end

if 1, printm 'projection profiles'
	it = cgb.nt;
	it = round(cgb.nt/2); % middle
	it = it + [-2 0 2];
	ia = imin(abs(cgb.ad - 45)); % closest to 45
%	ia = ceil(ia/2);
	pro = @(y) col(y(:,it,ia));
	arg = [];
	for ii=1:numel(systypes)
		arg = [arg pro(yb{ii})];
	end
	if im
		clf, plot([arg pro(ya)])
		text(10, 200, sprintf('ia=%d', ia))
		text(10, 400, sprintf('ang=%g', cgb.ad(ia)))
		ir_legend({systypes{:}, 'true'})
		axisy(0, 1.2 * max(col(pro(ya))))
		grid
	end
%prompt
return
end

if 0 % dd1 vs dd2 - they match well
	i_dd1 = strmatch('dd1', systypes);
	i_dd2 = strmatch('dd2', systypes);
	im clf, im(yb{i_dd1} - yb{i_dd2}), cbar
	equivs(yb{i_dd1}, yb{i_dd2}, 'thresh', 2e-5)
return
end

if 1, printm 'show projections and differences'
	im('plc', 2, 1+nsys)
	im(1, x0)
	ia = round([1 cgb.na/4 cgb.na/2 cgb.na]);
	im(nsys+2, ya(:,:,ia)), axis normal

	for ii=1:nsys
		tmp = yb{ii};
		im(ii+1, tmp(:,:,ia)), axis normal
		xlabel(systypes{ii})
		im(ii+2+nsys, tmp(:,:,ia)-ya(:,:,ia)), axis normal
	end
prompt
end

if 1, printm 'show back-projections and pairwise differences'
	im('plc', nsys, nsys)
	iz = round([1 igb.nz/4 igb.nz/2 igb.nz]);
	iz = 1:2:igb.nz;
	clear yt
	for ii=1:nsys
%		tmp = cgb.ones;
		tmp = cgb.zeros; tmp(:,:,20) = 1;
		yt{ii} = Ab{ii}' * tmp;
	end
	for ii=1:nsys
		for jj=1:nsys
			if ii == jj
				im('subplot', (ii-1)*nsys + jj)
				im(yt{ii}(:,:,iz))
				xlabel(systypes{ii})
			elseif ii > jj
				im('subplot', (ii-1)*nsys + jj)
				im(yt{ii}(:,:,iz) - yt{jj}(:,:,iz))
				xlabelf([systypes{ii} ' - ' systypes{jj}])
			end
		end
	end
end


text:
function [figH,m,cbH] = showMontage(imVol, slices, cmap, crop, numCols, figNum, flip)
% Show a montage of the volume data in imVol
%
% [figH,m,cbH] = showMontage(imVol, [slices=[]], [cmap=gray(256)], [crop=[]], [numCols=[]], figNum=figure, [flip='none'])
%
% Example:
%  t1    = niftiRead('t1.nii.gz');
%  imVol = t1.data; % that is --> : [226x271x226 int16]
%
% flip options are 'none', 'axial'
%
% HISTORY:
% 2007.02.22 RFD wrote it.
%
% Bob (c) VISTASOFT Team
%

%% Handle the Inputs

% If a file was passed in then we try to read it with niftiRead
if ~isnumeric(imVol) && exist(imVol,'file')
    try
        ni = niftiRead(imVol);
        imVol = ni.data;
    catch
        warning('A file was passed in, but that file does not appear to be a a file we can read.');
        figH = [];
        m    = [];
        cbH  = [];
        return
    end

end

if(~exist('slices','var') || isempty(slices))
  slices = [];
end

if(~exist('cmap','var') || isempty(cmap))
  cmap = gray(256);
end

if(exist('crop','var') && ~isempty(crop))
    % newSz = [diff(crop')+1 size(imVol,3)];
    imVol = imVol(crop(1,1):crop(1,2),:,:);
    imVol = imVol(:,crop(2,1):crop(2,2),:);
end

if(~exist('numCols','var')), numCols = []; end

if(~exist('figNum','var'))
    if exist('mrvNewGraphWin','file')
        figH = mrvNewGraphWin;
    else
        figH = figure;
    end
else
    figH = figure(figNum);
end

if(~exist('flip','var') || isempty(flip))
    flip = 'none';
end

if(strcmpi(flip(1),'a'))
    imVol = flip(permute(imVol,[2 1 3 4]),1);
end

%% Make the montage

colormap(cmap);
m = makeMontage(imVol,slices,[],numCols);
imagesc(m);
axis image;
cbH = colorbar;

if(nargout<2),   clear m;    end
if(nargout<1),   clear figH; end

return;


text:
function value = digamma ( x )

%*****************************************************************************80
%
%% DIGAMMA calculates the digamma or Psi function = d ( LOG ( GAMMA ( X ) ) ) / dX
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2004
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    J Bernardo,
%    Psi ( Digamma ) Function,
%    Algorithm AS 103,
%    Applied Statistics,
%    Volume 25, Number 3, pages 315-317, 1976.
%
%  Parameters:
%
%    Input, real X, the argument of the digamma function.
%    0 < X.
%
%    Output, real VALUE, the value of the digamma function at X.
%
  c = 8.5;
  d1 = -0.5772156649;
  s = 0.00001;
  s3 = 0.08333333333;
  s4 = 0.0083333333333;
  s5 = 0.003968253968;
%
%  The argument must be positive.
%
  if ( x <= 0.0 )

    value = 0.0;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DIGAMMA - Fatal error!\n' );
    fprintf ( 1, '  X <= 0.\n' );
    error ( 'DIGAMMA - Fatal error!' );
%
%  Use approximation if argument <= S.
%
  elseif ( x <= s )

    value = d1 - 1.0 / x;
%
%  Reduce the argument to DIGAMMA(X + N) where C <= (X + N).
%
  else

    value = 0.0;
    y = x;

    while ( y < c )
      value = value - 1.0 / y;
      y = y + 1.0;
    end
%
%  Use Stirling's (actually de Moivre's) expansion if C < argument.
%
    r = 1.0 / ( y * y );
    value = value + log ( y ) - 0.5 / y - r * ( s3 - r * ( s4 - r * s5 ) );

  end

  return
end


text:
function [ Qx,qx,theta ] = surface_U_v_direct( W,Hd,Hr,Theta,G,N,K,grt,beta )
    ebs=beta;
    B=Hd*W;
    B=B.';
     %%
    A=zeros(N,1,K,K);
    for i0=1:K
        for k0=1:K
         atp=diag(Hr(k0,:))*G*W(:,i0);
         A(:,:,i0,k0)=atp;
        end
    end
 %%
    theta=diag(Theta');
 %%
    Qx=zeros(N,N);
    for k0=1:K
        for i0=1:K
            Qx=Qx-abs(ebs(k0))^2.*A(:,:,i0,k0)*A(:,:,i0,k0)';
        end
    end
    qx=zeros(N,1);
    for k0=1:K
        tmp=zeros(N,1);
        for i0=1:K
            tmp=tmp+abs(ebs(k0))^2.*(B(i0,k0)')*A(:,:,i0,k0);
        end
        qx=qx+sqrt(grt(k0))*ebs(k0)'*A(:,:,k0,k0)-tmp;
    end
end

author: guohuayan


text:
% @author: Maziar Raissi

function plot_surface_griddata(X_star, u_star, xlab, ylab, tit)

x_l = min(X_star(:,1));
x_r = max(X_star(:,1));

y_l = min(X_star(:,2));
y_r = max(X_star(:,2));

nn = 100;
x = linspace(x_l, x_r, nn)';
y = linspace(y_l, y_r, nn)';
[Xplot, Yplot] = meshgrid(x,y);

Zplot = griddata(X_star(:,1),X_star(:,2),u_star,Xplot,Yplot,'cubic');

surface(Xplot, Yplot, Zplot);
xlabel(xlab);
ylabel(ylab);
title(tit);

axis tight
colormap jet
shading interp
colorbar
set(gca,'FontSize',14);
set(gcf, 'Color', 'w');

author: maziarraissi

repo: HPM

save_path: github-repos/MATLAB/maziarraissi-HPM


text:
%  INTERNAL FUNCTION: creates a grid of unique permutations of derivatives,
%  with non-increasing elements from left to right
%
%  ::
%
%    R1=derivatives_grid(R0)
%
%  Args:
%
%     - **R0** [vector|matrix]: initial permutation of derivatives
%
%  Returns:
%     :
%
%     - **R1** [matrix]: permutation indexes for derivatives of next order
%
%  See also:
%      mygrid
%
%

author: jmaih


text:
clear
close all

plotx2east

load martensite_single_grain

[grains,ebsd('indexed').grainId] =...
    calcGrains(ebsd('indexed'),'angle',3*degree);

cond_grains = ismember(grains.id,unique(ebsd(ids_of_interest).grainId));

%%
job = parentGrainReconstructor(ebsd,grains) % define job
job.useBoundaryOrientations = true;
clearvars grains ebsd

job.p2c = orientation.KurdjumovSachs(job.p2c.CS,job.p2c.SS);
job.calcParent2Child('quantile',0.6,'c2c'); % get representative OR

%%
job.calcHyperGraph3('threshold',5*degree,'c2c','mergeSimilar','mergethreshold',8*degree,'keepGraph')

numIters = [3 10 50];
grainId = 1844;

%%

for k = 1:length(numIters);
    if k == 1
    job.clusterHyperGraph3('numIter',numIters(k),...
        'inflationPower',1,'merged','mergethreshold',8*degree,'keepGraph')
        graphplotfunc(job,cond_grains,grainId,k,numIters)
    else
    job.clusterHyperGraph3('numIter',numIters(k) - numIters(k-1),...
        'inflationPower',1,'merged','mergethreshold',8*degree,'keepGraph')
        graphplotfunc(job,cond_grains,grainId,k,numIters)
    end
end

for k = 1:length(numIters);
    if k == 1
    job.clusterHyperGraph3('numIter',numIters(k),...
        'inflationPower',1.05,'merged','mergethreshold',8*degree,'keepGraph')
        graphplotfunc(job,cond_grains,grainId,k,numIters)
    else
    job.clusterHyperGraph3('numIter',numIters(k) - numIters(k-1),...
        'inflationPower',1.05,'merged','mergethreshold',8*degree,'keepGraph')
        graphplotfunc(job,cond_grains,grainId,k,numIters)
    end
end

%%
function graphplotfunc(job,cond_grains,grainId,k,numIters)
p2c = orientation.NishiyamaWassermann(job.p2c.CS,job.p2c.SS);
p2c = p2c.project2FundamentalRegion(job.p2c);

pOri = variants(p2c, job.grains.meanOrientation, job.votes.parentId(:,1));

nextAxis
plot(job.grains(cond_grains),...
    pOri(cond_grains))
lims = job.ebsd.extent;
xlim([lims(1) lims(2)])
ylim([lims(3) lims(4)])

text(lims(1)+5,lims(4)-10, [num2str(numIters(k)) ' iterations'],...
    'HorizontalAlignment','left',...
    'FontSize',20, 'FontWeight','bold','BackgroundColor', 'white')

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
clear all, close all, clc

load ../DATA/allFaces.mat

% We use the first 36 people for training data
trainingFaces = faces(:,1:sum(nfaces(1:36)));
avgFace = mean(trainingFaces,2);  % size n*m by 1;

% compute eigenfaces on mean-subtracted training data
X = trainingFaces-avgFace*ones(1,size(trainingFaces,2));
[U,S,V] = svd(X,'econ');

figure, axes('position',[0  0  1  1]), axis off
imagesc(reshape(avgFace,n,m)), colormap gray

for i=1:50  % plot the first 50 eigenfaces
    pause(0.1);  % wait for 0.1 seconds
    imagesc(reshape(U(:,i),n,m)); colormap gray;
end

%% Now show eigenface reconstruction of image that was omitted from test set

testFace = faces(:,1+sum(nfaces(1:36))); % first face of person 37
axes('position',[0  0  1  1]), axis off
imagesc(reshape(testFace,n,m)), colormap gray

testFaceMS = testFace - avgFace;
for r=25:25:2275
    reconFace = avgFace + (U(:,1:r)*(U(:,1:r)'*testFaceMS));
    imagesc(reshape(reconFace,n,m)), colormap gray
    title(['r=',num2str(r,'%d')]);
    pause(0.1)
end

%% Project person 2 and 7 onto PC5 and PC6

P1num = 2;  % person number 2
P2num = 7;  % person number 7
P1 = faces(:,1+sum(nfaces(1:P1num-1)):sum(nfaces(1:P1num)));
P2 = faces(:,1+sum(nfaces(1:P2num-1)):sum(nfaces(1:P2num)));
P1 = P1 - avgFace*ones(1,size(P1,2));
P2 = P2 - avgFace*ones(1,size(P2,2));

figure
subplot(1,2,1), imagesc(reshape(P1(:,1),n,m)); colormap gray, axis off
subplot(1,2,2), imagesc(reshape(P2(:,1),n,m)); colormap gray, axis off

% project onto PCA modes 5 and 6
PCAmodes = [5 6];
PCACoordsP1 = U(:,PCAmodes)'*P1;
PCACoordsP2 = U(:,PCAmodes)'*P2;

figure
plot(PCACoordsP1(1,:),PCACoordsP1(2,:),'kd','MarkerFaceColor','k')
axis([-4000 4000 -4000 4000]), hold on, grid on
plot(PCACoordsP2(1,:),PCACoordsP2(2,:),'r^','MarkerFaceColor','r')
set(gca,'XTick',[0],'YTick',[0]);

author: dynamicslab


text:
function [beta_gibbs,sigma_gibbs,beta_mean,sigma_mean,lambda_posterior]=panel4gibbs(N,n,h,T,k,q,Yi,Xi,s0,omegab,v0,It,Bu,pick,pickf)

% compute first  preliminary elements
% compute sbar
sbar=h+s0;
% compute the inverse of omegab
invomegab=diag(1./diag(omegab));
% initiate the Gibbs sampler
% initiate the counting of iterations
count=1;
pickcount=1;
% initiate the record matrices
beta_gibbs=zeros(q,It-Bu,N);
sigma_gibbs=zeros(n^2,It-Bu,N);
beta_mean=zeros(q,It-Bu);
sigma_mean=zeros(n^2,It-Bu);
lambda_posterior=zeros((It-Bu),1);

% step 1: compute initial values
% initial value for beta (use OLS values)
for ii=1:N
beta(:,ii)=bear.vec((Xi(:,:,ii)'*Xi(:,:,ii))\(Xi(:,:,ii)'*Yi(:,:,ii)));
end
% initial value for b
b=(1/N)*sum(beta,2);
% initial value for lambda1
lambda1=0.01;
sigmab=lambda1*omegab;
% initial value for sigma (use OLS values)
for ii=1:N
eps=Yi(:,:,ii)-Xi(:,:,ii)*reshape(beta(:,ii),k,n);
sigma(:,:,ii)=(1/(T-k-1))*eps'*eps;
end

hbar = bear.parfor_progressbar(It-Bu,'Progress of Panel BVAR Gibbs Sampler');  %create the progress bar

% run the Gibbs sampler
while count<=It

% step 2: obtain b
% first compute betam, the mean value of the betas over all units
betam=(1/N)*sum(beta,2);
% draw b from a multivariate normal N(betam,(1/N)*sigmab))
b=betam+chol(bear.nspd((1/N)*sigmab),'lower')*mvnrnd(zeros(q,1),eye(q))';

% step 3: obtain sigmab
% compute first vbar
for ii=1:N
temp(1,ii)=(beta(:,ii)-b)'*invomegab*(beta(:,ii)-b);
end
vbar=v0+sum(temp,2);
% compute lambda1
lambda1=bear.igrandn(sbar/2,vbar/2);
% recover sigmab
sigmab=lambda1*omegab;

% step 4: draw the series of betas
% first obtain the inverse of sigmab
invsigmab=diag(1./diag(sigmab));
% then loop over units
for ii=1:N
% take the choleski factor of sigma of unit ii, inverse it, and obtain from it the inverse of the original sigma
C=bear.trns(chol(bear.nspd(sigma(:,:,ii)),'Lower'));
invC=C\speye(n);
invsigma=invC*invC';
% obtain omegabar
invomegabar=kron(invsigma,Xi(:,:,ii)'*Xi(:,:,ii))+invsigmab;
% invert
C=bear.trns(chol(bear.nspd(invomegabar),'Lower'));
invC=C\speye(q);
omegabar=invC*invC';
% obtain betabar
betabar=omegabar*(kron(invsigma,Xi(:,:,ii)')*bear.vec(Yi(:,:,ii))+invsigmab*b);
% draw beta
beta(:,ii)=betabar+chol(bear.nspd(omegabar),'lower')*mvnrnd(zeros(q,1),eye(q))';
end

% step 5: draw the series of sigmas
% loop over units
for ii=1:N
% compute Stilde
Stilde=(Yi(:,:,ii)-Xi(:,:,ii)*reshape(beta(:,ii),k,n))'*(Yi(:,:,ii)-Xi(:,:,ii)*reshape(beta(:,ii),k,n));
% draw sigma
sigma(:,:,ii)=bear.iwdraw(Stilde,T);
end

   % record phase
   % if the burn-in sample phase is not yet over
   if count<=Bu
   % simply add 1 to the iteration count
   count=count+1;
   % on the other hand, if the burn-in sample phase is over
   elseif count>Bu
   % adding one iteration to the count will depend on wether post-burn selection applies
      % if there is no post burn selection
      if pick==0
      % record the draw
         % loop over units
         for ii=1:N
         beta_gibbs(:,count-Bu,ii)=beta(:,ii);
         sigma_gibbs(:,count-Bu,ii)=bear.vec(sigma(:,:,ii));
         end
      % and add one to the count
      count=count+1;
      % if there is post burn selection, only one draw over 'fpick' draws will be retained
      elseif pick==1
         % if the iteration does not correspond to fpick, don't record the results, don't increase the regular count, but do increase pickcount by 1
         if pickcount~=pickf
         pickcount=pickcount+1;
         % on the other hand, if the iteration does correspond to fpick
         elseif pickcount==pickf
         % do record the results
            % loop over units
            beta_mean(:,count-Bu)=b;
            sigma_mean(:,count-Bu)=bear.vec(mean(sigma,3));
            lambda_posterior(count-Bu)=lambda1;
            for ii=1:N
            beta_gibbs(:,count-Bu,ii)=beta(:,ii);
            sigma_gibbs(:,count-Bu,ii)=bear.vec(sigma(:,:,ii));
            end
         % then increase the regular count by 1 and re-initialise pickcount
         count=count+1;
         pickcount=1;
         end
      end
   end

   hbar.iterate(1);   % update progress by one iteration

end

close(hbar);   %close progress bar

author: european-central-bank

repo: BEAR-toolbox


text:
maxmag = ceil(10*max(newt2.Magnitude))/10;
mima = min(newt2.Magnitude);
if mima > 0 ; mima = 0 ; end

[bval,xt2] = hist(newt2.Magnitude,(mima:0.1:maxmag));
% normalise to annula rates
bval = bval/(max(newt2.Date)-min(newt2.Date));
bvalsum = cumsum(bval); % N for M <=
bval2 = bval(length(bval):-1:1);
bvalsum3 = cumsum(bval(length(bval):-1:1));    % N for M >= (counted backwards)
xt3 = (maxmag:-0.1:mima);

backg_ab = log10(bvalsum3);

figure_w_normalized_uicontrolunits(bfig);delete(gca);delete(gca); delete(gca); delete(gca)
rect = [0.22,  0.3, 0.65, 0.6];           % plot Freq-Mag curves
axes('position',rect);

%%
% plot the cum. sum in each bin  %%
%%

%pl =semilogy(xt3,bvalsum3,'sb');
%set(pl,'LineWidth',1.0,'MarkerSize',6,...
%    'MarkerFaceColor','w','MarkerEdgeColor','k');
%hold on
pl1 =semilogy(xt3,bval2,'^b');
set(pl1,'LineWidth',1.0,'MarkerSize',4,...
    'MarkerFaceColor',[0.7 0.7 .7],'MarkerEdgeColor','k');

bv2 = [];bv3 = [] ; me = [];BV = [];
ni2 = 50;

for i = 1:ni2/1:length(newt2)-ni2
    [bv magco stan ] =  bvalca2(newt2(i:i+ni2,:));
    l = newt2(i:i+ni2,:) >= magco;
    nn2 = newt2(i:i+ni2,l);

    [bval,xt2] = hist(nn2(:,6),(mima:0.1:maxmag));
    % normalise to annula rates
    bval = bval/(max(nn2(:,3))-min(nn2(:,3)));
    bvalsum = cumsum(bval); % N for M <=
    bval2 = bval(length(bval):-1:1);
    bvalsum3 = cumsum(bval(length(bval):-1:1));    % N for M >= (counted backwards)

    hold on
    pl1 =semilogy(xt3,bval2,'^b');
    set(pl1,'LineWidth',1.0,'MarkerSize',4,...
        'MarkerFaceColor',[0.7 0.7 .7],'MarkerEdgeColor','k');
    pause

end

bv2 = [bv2 ; magco newt2(i+ni2/2,3)];
BV = [BV ; magco newt2(i,3) ; magco newt2(i+ni2,3) ; inf inf];

author: CelsoReyes


text:
function h = p28_h ( n, x )

%*****************************************************************************80
%
%% P28_H evaluates the Hessian for problem 28.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real H(N,N), the N by N Hessian matrix.
%
  h = zeros ( n, n );

  r = sqrt ( x(1)^2 + x(2)^2 );

  rx1 = x(1) / r;
  rx2 = x(2) / r;

  rx1x1 = x(2)^2 / r^3;
  rx1x2 = - x(1) * x(2) / r^3;
  rx2x1 = - x(1) * x(2) / r^3;
  rx2x2 = x(1)^2 / r^3;
%
%  F = A * B
%  dFdX1 = ( Ar * B + A * Br ) * Rx1
%  d2FdX1dX1 = ( Arr * B + Ar * Br ) * Rx1^2 + ( Ar * B + A * Br ) * Rx1x1
%  etc
%
  a = sqrt ( r );
  ar = 0.5 / sqrt ( r );
  arr = - 0.25 / sqrt ( r^3 );

  b = 1.0 + ( sin ( 50.0 * r^0.2 ) )^2;
  br = 10.0 * sin ( 100.0 * r^0.2 ) * r^(-0.8);
  brr = 200.0 * cos ( 100.0 * r^0.2 ) * r^(-1.6) ...
    - 10.0 * sin ( 100.0 * r^0.2 ) * 0.8 * r^(-1.8);

  h(1,1) = ( arr * b + 2.0 * ar * br + a * brr ) * rx1 * rx1 ...
    + ( ar * b + a * br ) * rx1x1;

  h(1,2) = ( arr * b + 2.0 * ar * br + a * brr ) * rx1 * rx2 ...
    + ( ar * b + a * br ) * rx1x2;

  h(2,1) = ( arr * b + 2.0 * ar * br + a * brr ) * rx2 * rx1 ...
    + ( ar * b + a * br ) * rx2x1;

  h(2,2) = ( arr * b + 2.0 * ar * br + a * brr ) * rx2 * rx2 ...
    + ( ar * b + a * br ) * rx2x2;

  return
end

author: johannesgerer


text:
function test2
%TEST2 test cs_sparse, cs_permute, cs_pvec, cs_ipvec, cs_symperm
%
% Example:
%   test2
% See also: testall

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

rand ('state', 0)
% clf

for trial = 1:100
    m = fix (10 * rand (1)) ;
    n = fix (10 * rand (1)) ;
    nz = fix (100 * rand (1)) ;

    i = 1 + fix (m * rand (nz,1)) ;
    j = 1 + fix (n * rand (nz,1)) ;
    x = rand (nz,1) ;

    if (~ispc)
        if (mod (trial, 2) == 1)
            x = x + 1i * (2*rand(nz,1)-1) ;
        end
    end

    A = sparse (i,j,x) ;
    B = cs_sparse (i,j,x) ;
    D = cs_sparse2 (i,j,x) ;
    fprintf ('%3d %3d %6d : %6d %6d : %d\n', ...
        m, n, nz, nnz (A), nnz(B), nz-nnz(A)) ;

    err = norm (A-B,1) / max (1, norm (A,1)) ;
    if (err > 0)
        disp ('err = ') ;
        disp (err) ;
    end
    if (err > 1e-14)
        error ('!') ;
    end

    if (nnz (B-D) > 0)
        error ('!') ;
    end

    if (nnz (A) ~= nnz (B))
        error ('nz!') ;
    end

    if (max (1,nnz (B)) ~= max (1,nzmax (B)))
        nnz (B)
        nzmax (B)
        error ('nzmax!') ;
    end
    % pack

    [m n] = size (A) ;
    p = randperm (m) ;
    q = randperm (n) ;
    C1 = A (p,q) ;
    C2 = cs_permute (A,p,q) ;
    err = norm (C1-C2,1) ;
    if (err > 0)
        error ('!') ;
    end

%    subplot (1,2,1) ; spy (A)
%    subplot (1,2,2) ; spy (C2)
%    drawnow

    x = rand (m,1) ;

    if (~ispc)
        if (mod (trial, 2) == 1)
            x = x + 1i * (2*rand(m,1)-1) ;
        end
    end

    x1 = x (p) ;
    x2 = cs_pvec (x, p) ;

    err = norm (x1-x2,1) ;
    if (err > 0)
        error ('!') ;
    end

    x1 = zeros (m,1) ;
    x1 (p) = x ;                                                            %#ok
    x2 = cs_ipvec (x, p) ;                                                  %#ok

    n = min (m,n) ;
    B = A (1:n, 1:n) ;
    p = randperm (n) ;
    B = B+B' ;

    C1 = triu (B (p,p)) ;
    C2 = cs_symperm (B,p) ;

    try
        pp = amd (C2) ;                                                     %#ok
    catch
        pp = symamd (C2) ;                                                  %#ok
    end

    err = norm (C1-C2,1) ;
    if (err > 0)
        error ('!') ;
    end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function GameOfLife
% This is a simple simulation of Conway Game of life GoL
% it is good for understanding Cellular Automata (CA) concept
% GoL Rules:
%   1. Survival: an alive cell live if it has 2 or 3 alive neighbors
%   2. Birth: a dead cell will be alive if it has 3 alive neighbors
%   3. Deaths:
%        Lonless: alive cell dies if it has 0 or 1 alive neighbors
%        Overcrowding: alive cell dies if it has 4 or more alive neighbors
% Any questions related to CA are welcome
% By: Ibraheem Al-Dhamari

clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%      Initialization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% size= 500x500
% different random initial values
%A= rand(500,500);
%A= ones(500,500);

% periodic configuration
%A= zeros(500,500);
% A(100,200)=1;
% A(100,201)=1;
% A(100,202)=1;

% initial from image
 A=imread('CA02.JPG');

% convert to binary--> % states={0,1}
A=im2bw(A);
% visualize the initial states
disp('the binary image')
imshow(A);
pause
% boundary type: 0= reflection
%                1= doublication
%                2= null, zeros

% this step enlarge A with 4 virtual vectors
A=Bnd(A,0);
[d1,d2]=size(A);
disp('the extended binary image')
whos A
imshow(A);
pause
B=A;
t=0;
stp=false; % to stop when if no new configrations
%B is the CA in time t
%A is the CA in time t+1
%t is the number of generations

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Play ^_^
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while ~stp & (t<10) % repeat for 10 generations
    % for each cell in the CA
    for i=2:d1-1
        for j=2:d2-1
            % apply Game of life rule
            A(i,j)=GOL(A,B,i,j)            ;
        end
    end
    % visualize what happened
    disp('the CA image')
    imshow(A);
    drawnow;
%    pause
    % save B
    if A==B
       stp=true; % no more new states
    end
    B=A;
    t=t+1
end

%==========================================
%   Game of Life Rules
%==========================================
function s=GOL (A,B,i,j)
% game of life rule
sm=0;
% count number of alive neighbors
sm=sm+ B(i-1,j-1)+B(i-1,j)+B(i-1,j+1);
sm=sm+ B(i,j-1)+           B(i,j+1);
sm=sm+ B(i+1,j-1)+B(i+1,j)+B(i+1,j+1);

% compute the new state of the current cell
s=B(i,j);
if B(i,j)==1
    if (sm>1)&&(sm<4)
        s=1;
    else
        s=0 ;
    end
else
    if sm==3
       s=1;
    end
end

%==========================================
%   Boundary Type
%==========================================
function bA= Bnd(A,k)
% add new four vectors based on boundary type
[d1, d2]=size(A);
d1=d1+2; d2=d2+2;
X=ones(d1,d2);
X=im2bw(X);
X(2:d1-1,2:d2-1)=A;
imshow(X);
whos A X
if k==0 % Reflection
   X(  1  , 2:d2-1)=A(end , :);
   X(  d1 , 2:d2-1)=A( 1  , :);
   X( 2:d1-1 , 1  )=A(: , end);
   X( 2:d1-1 , d2 )=A(: ,  1 );

   X(1,1)    =A(end,end);
   X(1,end)  =A(end,1);
   X(end,1)  =A(1,end);
   X(end,end)=A(1,1);
elseif k==1 % Double
   X(  1  , 2:d2-1)=A( 1  , :);
   X(  d1 , 2:d2-1)=A(end , :);
   X( 2:d1-1 , 1  )=A(: ,  1 );
   X( 2:d1-1 , d2 )=A(: , end);

   X(1,1)    =A(end,1);
   X(1,end)  =A(end,end);
   X(end,1)  =A(1,1);
   X(end,end)=A(1,end);

else % k==2 % zeros
   X(  1  ,:)=0;
   X( end ,:)=0;
   X(: ,  1  )=0;
   X(: , end )=0;
end
bA=X;


text:
clc;
clear;
close all;

T0 = 10000000 ; % initial temperature
r = 0.999 ; % temperature damping rate
Ts = 0.001 ; % stop temperature

model = initModel();

% initialization
while(1)
route = randomSol(model);
if(isFeasible(route,model))
    break;
end
end

cost = calculateCost(route,model);
T = T0;
min = cost;

cnt = 1;

% SA
while(T > Ts)
    flag = '#';
    mode = randi([1 3]);
    newRoute = createNeibor(route,model,mode);
    newCost = calculateCost(newRoute,model);
    delta = newCost - cost;

    if(delta < 0)
        cost = newCost;
        route = newRoute;
        flag = '*';
    else
        p=exp(-delta/T);
        if rand <= p
             cost = newCost;
             route = newRoute;
             flag = '^';
        end
    end

     if cost < min
         min = cost;
     end

     costArr(cnt) = cost;

    T = T*r; %  annealing
    disp([flag 'Iteration ' num2str(cnt) ': Best Cost = ' num2str(cost) ' T = ' num2str(T)]);
    cnt = cnt+1;

end
disp(min);
plot(costArr);

author: lzane

repo: VRP-using-SA-with-Matlab


text:
% -*- INTERNAL UNDOCUMENTED FUNCTION -*-
function [geometry, msh, space, u] = ...
              solve_laplace_3d_iso (problem_data, method_data)

warning ('geopdes:obsolete', 'Function SOLVE_LAPLACE_3D_ISO is obsolete. Using SOLVE_LAPLACE_ISO instead')

[geometry, msh, space, u] = solve_laplace_iso (problem_data, method_data);

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function msm_to_mm_coordinate_integer_skewsymmetric ( output_filename, a )

%*****************************************************************************80
%
%% MSM_TO_MM_COORDINATE_INTEGER_SKEWSYMMETRIC writes a "matrix coordinate integer skew-symmetric" Matrix Market file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the name of the file to which the information
%    should be written.
%
%    Input, sparse matrix A, the NROW by NCOL matrix, stored in MATLAB sparse
%    matrix format, which is to be written to the file.
%
  [ nrow, ncol ] = size ( a );
  nnzeros = nnz ( a );

  fid = fopen ( output_filename, 'wt+' );

  if ( fid < 0 );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MSM_TO_MM_COORDINATE_INTEGER_SKEWSYMMETRIC - Fatal error!\n' );
    fprintf ( 1, '  Cannot open the output file.\n' );
    error ( 'MSM_TO_MM_COORDINATE_INTEGER_SKEWSYMMETRIC - Fatal error!');
  end;

  fprintf ( fid, '%%%%MatrixMarket matrix coordinate integer skew-symmetric\n');
  fprintf ( fid, '%%%%  Created by MSM_TO_MM_COORDINATE_INTEGER_SKEWSYMMETRIC.M\n' );
  fprintf ( fid, '  %d  %d  %d\n', nrow, ncol, nnzeros );

  for j = 1 : ncol
    [ rows, temp, vals ] = find ( a(:,j) );
    sz = size ( rows );
    for k2 = 1 : sz
      if ( j < rows(k2) )
        fprintf ( fid, '  %d  %d  %d\n', rows(k2), j, vals(k2) );
      end
    end

  end

  fclose ( fid );

  return
end

repo: jburkardt-m


text:
function A = makehatch_plus(hatch,n,m)
%MAKEHATCH_PLUS Predefined hatch patterns
%
% Modification of MAKEHATCH to allow for selection of matrix size. Useful whe using
%   APPLYHATCH_PLUS with higher resolution output.
%
% input (optional) N    size of hatch matrix (default = 6)
% input (optional) M    width of lines and dots in hatching (default = 1)
%
%  MAKEHATCH_PLUS(HATCH,N,M) returns a matrix with the hatch pattern for HATCH
%   according to the following table:
%      HATCH        pattern
%     -------      ---------
%        /          right-slanted lines
%        \          left-slanted lines
%        |          vertical lines
%        -          horizontal lines
%        +          crossing vertical and horizontal lines
%        x          criss-crossing lines
%        .          square dots
%        c          circular dots
%        w          Just a blank white pattern
%        k          Just a totally black pattern
%
%  See also: APPLYHATCH_PLUS

%  By Ben Hinkle, bhinkle@mathworks.com
%  This code is in the public domain.

% Modified Brian FG Katz    8-aout-03
% Modified David M Kaplan    19-fevrier-08

if ~exist('n','var'), n = 6; end
if ~exist('m','var'), m = 1; end
n=round(n);

switch (hatch)
  case '\'
    [B,C] = meshgrid( 0:n-1 );
    B = B-C;
    clear C
    A = abs(B) <= m/2;
    A = A | abs(B-n) <= m/2;
    A = A | abs(B+n) <= m/2;
  case '/'
    A = fliplr(makehatch_plus('\',n,m));
  case '|'
    A=zeros(n);
    A(:,1:m) = 1;
  case '-'
    A = makehatch_plus('|',n,m);
    A = A';
  case '+'
    A = makehatch_plus('|',n,m);
    A = A | A';
  case 'x'
    A = makehatch_plus('\',n,m);
    A = A | fliplr(A);
  case '.'
    A=zeros(n);
    A(1:2*m,1:2*m)=1;
  case 'c'
    [B,C] = meshgrid( 0:n-1 );
    A = sqrt(B.^2+C.^2) <= m;
    A = A | fliplr(A) | flipud(A) | flipud(fliplr(A));
  case 'w'
    A = zeros(n);
  case 'k'
    A = ones(n);
  otherwise
    error(['Undefined hatch pattern "' hatch '".']);
end