text:
function varargout = classify_naive_bayes(meth, varargin)
% Naive Bayes classifier
%
% :Usage:
% ::
%
%     % set up model structure
%     bayes_model = classify_naive_bayes('setup', Y, Xi, [activation_cutoff, selectivity_cutoff]);
%
% :Inputs:
%
%   **Y:**
%        is full (not sparse); empty features will be eliminated
%
%   **Xi:**
%        is obs x classes, an indicator matrix of 1's and 0's
%
% :TEST: test classifier; make a prediction about classes from data
% ::
%
%    [class_est log_joint best_log map p_obs_act_given_class] = classify_naive_bayes('test', Y, bayes_model);
%
% :EVAL: evaluate classification accuracy
% ::
%
%    [prop_correct, confusion_mtx, misclass, prop_correct_by_class, chance, chance_95_ci] = classify_naive_bayes('eval', true_class, class_est, wh_obs);
%
% :APPARENT: apparent classification; with full dataset
% ::
%
%     bayes_model = classify_naive_bayes('apparent', Y, bayes_model);
%
% :XVAL: crossvalidate
% ::
%
%    xval = classify_naive_bayes('xval', Y, Xi);
%
%    bayes_model = classify_naive_bayes('write', bayes_model, Y, volInfo, conditionnames)
%    bayes_model = classify_naive_bayes('write', bayes_model, Y, MC_Setup.volInfo, MC_Setup.Xinms);
%
% To add feature abstraction step within xval:
% ::
%
%    xval = classify_naive_bayes('xval', Y, Xi, 0, .9, .05, 1, volInfo );
%
% :PLOT:
% ::
%
%    classify_naive_bayes('plot', bayes_model, ['pa|t', 'lr', 'lr surface',  'map plot', or 'class plot']);
%
% :Optional Inputs: (any order)
%
% Threshold (abs. value), and colors in cell array
% ::
%
%    classify_naive_bayes('plot', bayes_model, 'lr', .10, {[1 .7 0] [0 0 1]});
%
%
% :Examples:
% ::
%
%    [bayes_model, Y] = classify_naive_bayes('setup', Y, Xi);
%
%    % Test obs. 2
%    tic, [class_est log_joint best_log map] = classify_naive_bayes('test',Y(2,:), bayes_model); toc
%
%    % Get apparent classification rate and look at confusion matrix
%    tic, bayes_model = classify_naive_bayes('apparent', Y, bayes_model); toc
%    bayes_model.apparent.confusion_mtx
%
%    % Cross-validate
%    bayes_model.xval = classify_naive_bayes('xval', Y, Xi);
%    bayes_model.xval.confusion_mtx
%
% :Example 2:
% ::
%
%     Select features, and do apparent and cross-validated classification
%    Y = MC_Setup.unweighted_study_data';
%    wh = sum(Y) > 5;
%    Y = Y(:, wh);
%    whos Y
%    [bayes_model, Y] = classify_naive_bayes('setup', Y, Xi);
%    bayes_model = classify_naive_bayes('apparent', Y, bayes_model);
%    bayes_model.apparent.prop_correct, bayes_model.apparent.confusion_mtx
%    bayes_model.xval = classify_naive_bayes('xval', Y, Xi);
%    bayes_model.xval.prop_correct, bayes_model.xval.confusion_mtx
%
% :Example 3:
%    create_figure('hist'); hist(bayes_model.pa1_given_t, 100);
%    xval = classify_naive_bayes('xval', Y, Xi, .05, .3);
%
% Get results from key regions and run classifier only on those regions:
% ::
%
%     cl = classify_naive_bayes('plot', bayes_model, 'lr', .10, {[1 .7 0] [0 0 1]});
%    [studybyroi,studybyset] = Meta_cluster_tools('getdata',cl{1},dat',volInfo);
%    bayes_model_regions = classify_naive_bayes('setup', studybyroi, Xi);
%    bayes_model_regions = classify_naive_bayes('apparent', studybyroi,bayes_model_regions);
%    disp('Apparent confusion')
%    disp(bayes_model_regions.apparent.confusion_mtx)
%
%    bayes_model_regions.xval = classify_naive_bayes('xval', studybyroi, Xi);
%    disp('Cross-validated confusion')
%    disp(bayes_model_regions.xval.confusion_mtx)
%
%    fprintf('Proportion correct: Apparent: %3.0f%%  Xval: %3.0f%%\n', 100*bayes_model_regions.apparent.prop_correct, 100*bayes_model_regions.xval.prop_correct);
%    fprintf('Proportion correct by class: \t'); fprintf('%3.0f%%\t', 100*bayes_model_regions.xval.prop_correct_by_class);
%    fprintf('\n');

%    sz = cat(1,cl{1}(:).numVox);
%    cl{1}(sz < 10) = [];
%
%    subcl = subclusters_from_local_max(cl{1}, 10);
%
% ..
%    Tor Wager, Sept 2007
% ..

switch meth

    case 'setup'
        Y = full(varargin{1});
        Xi = varargin{2};

        activation_cutoff = [];
        selectivity_cutoff = [];
        return_full_ps = 1;         % for faster processing if not all fields are needed (i.e., in xval)

        if length(varargin) > 2

            activation_cutoff = varargin{3};
            selectivity_cutoff = varargin{4};

        end

        if length(varargin) > 4
            return_full_ps = varargin{5};
        end

        k = 1; % default k regularization param
        if length(varargin) > 5
            k = varargin{6};
        end

        % setup
        [Y, whkeep] = eliminate_empty_Y(Y);

        [bayes_model, Y] = classify_naive_bayes_setup(Y, Xi, activation_cutoff, selectivity_cutoff, whkeep, return_full_ps, k);

        gam = 1; % default gam prior weighting param
        if length(varargin) > 6
            gam = varargin{7};
        end
        bayes_model.params.gam = gam;

        varargout{1} = bayes_model;
        varargout{2} = Y;

    case 'apparent'

        Y = full(varargin{1});
        bayes_model = varargin{2};

        bayes_model.true_class = indic2condf(bayes_model.Xi);

        % test apparent classification rate
        [bayes_model.apparent.class_est bayes_model.apparent.log_joint bayes_model.apparent.best_log bayes_model.apparent.map bayes_model.apparent.pact_given_task] = ...
            classify_naive_bayes_test(Y, bayes_model);

        [bayes_model.apparent.prop_correct, bayes_model.apparent.confusion_mtx, bayes_model.apparent.misclass, ...
            bayes_model.apparent.prop_correct_by_class, bayes_model.apparent.chance, bayes_model.apparent.chance_95_ci] = ...
            classify_naive_bayes_eval(bayes_model.true_class, bayes_model.apparent.class_est);

        varargout{1} = bayes_model;

    case 'test'
        Y = full(varargin{1});
        bayes_model = varargin{2};

        [class_est log_joint best_log map p_obs_act_given_class] = classify_naive_bayes_test(Y, bayes_model);

        varargout{1} = class_est;
        varargout{2} = log_joint;
        varargout{3} = best_log;
        varargout{4} = map;
        varargout{5} = p_obs_act_given_class;

    case 'eval'
        bayes_model = varargin{1};
        class_est = varargin{2};
        wh_obs = [];
        if length(varargin) > 2, wh_obs = varargin{3}; end

        [prop_correct, confusion_mtx, misclass, prop_correct_by_class, chance, chance_95_ci] = classify_naive_bayes_eval(bayes_model, class_est, wh_obs);

        varargout{1} = prop_correct;
        varargout{2} = confusion_mtx;
        varargout{3} = misclass;
        varargout{4} = prop_correct_by_class;
        varargout{5} = chance;
        varargout{6} = chance_95_ci;

        case 'abstract'
          %bayes_model = classify_naive_bayes('abstract', Y, thresh, bayes_model, volInfo);

        Y = full(varargin{1});
        thresh = varargin{2};

        bayes_model = varargin{3};
        volInfo = varargin{4};

        reducedY = bayes_meta_feature_abstract(Y, thresh, bayes_model, volInfo);

        [bayes_model, Y] = classify_naive_bayes('setup', reducedY, bayes_model.Xi, bayes_model.params.activation_cutoff, bayes_model.params.selectivity_cutoff, 1, bayes_model.params.k);

        varargout{1} = bayes_model;
        varargout{2} = Y;

    case 'xval'
        Y = full(varargin{1});
        Xi = varargin{2};

        activation_cutoff = [];
        selectivity_cutoff = [];
        if length(varargin) > 2

            activation_cutoff = varargin{3};
            selectivity_cutoff = varargin{4};

        end

        k = .5;
        gam = 1;
        if length(varargin) > 4

            k = varargin{5};
            gam = varargin{6};

        end

        if length(varargin) > 6
            volInfo = varargin{7};
            do_feature_abstract = 1;
        end

        % initialize
        if do_feature_abstract, Yorig = Y; end

        [Y, whkeep] = eliminate_empty_Y(Y);

        Y = sparse(Y);

        nobs = size(Y, 1);
        fprintf('\nCross-validating %3.0f observations: 00000', nobs);

        log_joint = zeros(nobs, size(Xi, 2));
        map = zeros(nobs, size(Xi, 2));
        class_est = zeros(nobs, 1);
        best_log = zeros(nobs, 1);

        true_class = indic2condf(Xi);

        % Update volinfo for feature abstraction
% %         volInfo.wh_inmask = double(whkeep');
% %         volInfo.n_inmask = sum(whkeep);
% %         volInfo.image_indx(volInfo.image_indx) = whkeep';
% %         volInfo.xyzlist(~whkeep',:) = [];

        % run
        for i = 1:nobs

            if mod(i,10) == 0, fprintf(1,'\b\b\b\b\b%05d',i); end

            include_in_training = true(nobs, 1);
            include_in_training(i) = 0;

            % setup: does feature-selection as well, if requested; no
            % need to return Y, because ...
            bayes_model = classify_naive_bayes_setup(Y(include_in_training, :), Xi(include_in_training, :), activation_cutoff, selectivity_cutoff, whkeep, 0, k);

            bayes_model.params.k = k;
            bayes_model.params.gam = gam;

            testdata = Y(i, bayes_model.whkeep_from_notempty);

            % feature abstraction step
            if do_feature_abstract
                [reducedY, cl, cluster_indx] = bayes_meta_feature_abstract(Yorig(include_in_training, :), .1, .1, bayes_model, volInfo, 0);
                cluster_indx = cluster_indx(bayes_model.whkeep);

                bayes_model = classify_naive_bayes('setup', reducedY, Xi(include_in_training, :), 0, 0, gam, k);

                wh_cl = cluster_indx(testdata);
                wh_cl = unique(wh_cl); wh_cl(wh_cl == 0) = [];
                testdata = false(1, bayes_model.nfeatures);
                testdata(wh_cl) = 1;

            end

            [class_est(i) log_joint(i,:) best_log(i) map(i,:)] = classify_naive_bayes_test(testdata, bayes_model);
        end

        fprintf(1, '\n');

        [prop_correct, confusion_mtx, misclass, prop_correct_by_class, chance, chance_95_ci] = classify_naive_bayes_eval(true_class, class_est);

        varargout{1} = struct('class_est', class_est, 'log_joint', log_joint, 'best_log', best_log, 'map', map, ...
            'prop_correct', prop_correct, 'prop_correct_by_class', prop_correct_by_class, 'confusion_mtx', confusion_mtx, ...
            'misclass', misclass, 'chance', chance, 'chance_95_ci', chance_95_ci, 'params', bayes_model.params);

    case {'write', 'images', 'write_images'}
        bayes_model = varargin{1};
        Y = varargin{2};
        volInfo = varargin{3};
        conditionnames = varargin{4};

        bayes_model = write_images(bayes_model, Y, volInfo, conditionnames);
        varargout{1} = bayes_model;

    case {'plot', 'output'}
        bayes_model = varargin{1};
        disptype = varargin{2};

        cl = display_output(bayes_model, disptype, varargin{3:end});

        varargout{1} = cl;

    otherwise
        error('Unknown method (invalid first argument).');
end

end   % end main function

function [bayes_model, Y] = classify_naive_bayes_setup(Y, Xi, activation_cutoff, selectivity_cutoff, whkeep, return_full_ps, k)

% Y is data matrix, nobs observations x nfeatures variables
% in a meta-analysis, e.g., this is studies x voxels
%
% Each row of Y is an observation; if Y has one row, we're classifying an
% observation (e.g., one brain map)
%
% Xi are indicators for tasks (classes), ntasks columns of nobs elements
% (i.e., an nobs x ntasks matrix)

% eliminate features with no 'on' values (activations)

% k is regularization param, biases towards 0.5
if nargin < 7, k = 1; end

whkeep_from_notempty = whkeep;  % indicator of features after feature-selection, in index of not-empties

% feature selection, if requested
if ~isempty(selectivity_cutoff) && ~isempty(activation_cutoff)

    [Y, whkeep, whkeep_from_notempty] = select_features(Y, Xi, activation_cutoff, selectivity_cutoff, whkeep);

end

[nobs, nfeatures] = size(Y);
nclasses = size(Xi, 2);

if return_full_ps
    [priors, pa1_given_t, pa0_given_t, pt_given_act1, pt_given_act0, pa1_given_not_t] = ...
        bayes_get_probabilities(Y, Xi, k);

else
    % key results (not all) only; saves computation time in
    % crossvalidation.
    [priors, pa1_given_t, pa0_given_t] = bayes_get_probabilities(Y, Xi, k);

end

true_class = indic2condf(Xi);

bayes_model = struct('nobs', nobs, 'nfeatures', nfeatures, 'nclasses', nclasses, ...
    'whkeep', whkeep, 'whkeep_from_notempty', whkeep_from_notempty, ...
    'params', [], ...
    'priors', priors, 'pa1_given_t', pa1_given_t, 'pa0_given_t', pa0_given_t, ...
    'Xi', Xi, 'true_class', true_class);

bayes_model.params.gam = 1;  % gamma, weighting on prior vs. joint evidence
bayes_model.params.k = k;      % k, regularization of pa|t, bias towards 0.5.
bayes_model.params.activation_cutoff = activation_cutoff;
bayes_model.params.selectivity_cutoff = selectivity_cutoff;

if return_full_ps
    bayes_model.pt_given_act1 = pt_given_act1;
    bayes_model.pt_given_act0 = pt_given_act0;
    bayes_model.pa1_given_not_t = pa1_given_not_t;

end

end

function [class_est log_joint best_log map p_obs_act_given_class] = classify_naive_bayes_test(Y, bayes_model)

% Estimate classes, log joint p, and MAP estimates for a test data set

% P(task | act)
% The var. below is all that is needed to choose the best class.
% argmax log_joint
% do this for each observation

[nobs, nfeatures] = size(Y);

if nfeatures ~= bayes_model.nfeatures
    %  try eliminating zero voxels
    Y = Y(:, bayes_model.whkeep);
    [nobs, nfeatures] = size(Y);
end

if nfeatures ~= bayes_model.nfeatures
    error('Number of features (i.e., voxels/regions) in Y must equal that in bayes_model structure or that of original training data.');
end

log_joint = zeros(nobs, bayes_model.nclasses);
map = zeros(nobs, bayes_model.nclasses);

class_est = zeros(nobs, 1);
best_log = zeros(nobs, 1);

gam = bayes_model.params.gam;    % prior weighting factor; from 0 to nfeatures.

logpriors = log(bayes_model.priors);
%logpriors = logpriors(ones(nfeatures, bayes_model.nclasses));

for i = 1:nobs

    wh_active = Y(i,:)' > 0;

    % sum of log prob. of activation at each active voxel, plus sum of
    % log prob. of NOT activation at each inactive voxel (1 - p(activation))
    % Same as ln of product of all p(act state = observed state |
    % class) for all features

    % prob. that activity state is the observed one in the test vector,
    % given each class.
    % p(Fi = fi | C = c)
    % p(ACTi = acti | T = t)

    % % % %     p_obs_act_given_class(wh_active, :) = bayes_model.pt_given_act1(wh_active, :);         % p(active | t, obs active)
    % % % %     p_obs_act_given_class(~wh_active, :) = bayes_model.pt_given_act0(~wh_active, :);

    p_obs_act_given_class = zeros(nfeatures, bayes_model.nclasses);
    p_obs_act_given_class(wh_active, :) = bayes_model.pa1_given_t(wh_active, :);         % p(active | t, obs active)
    p_obs_act_given_class(~wh_active, :) = bayes_model.pa0_given_t(~wh_active, :);   % p(not active | t, obs not active)

    loglikelihood = sum(log(p_obs_act_given_class));
    log_joint(i,:) = gam * log(bayes_model.priors) + loglikelihood;      % log of joint prob. p(T, Y), p(task, act1, act2, act3, etc.)

    %log_joint(i,:) = sum(log(p_obs_act_given_class)); % proportional to the posterior.

    [best_log(i) class_est(i)] = max(log_joint(i,:));

    % Divide joint by evidence
    % get scaled estimate of p(task | activation)
    % These are Max. a posteriori estimates, and values may be > 1 because
    % we're assuming independence when that assumption is not likely to be true
    % * something still wrong with this
    % % % %
    % % % %         % evidence
    % % % %         log_evidence = sum(bayes_model.log_evidence_act(wh_active)) + sum(bayes_model.log_evidence_notact(~wh_active));
    % % % %         % because something is wrong, exp() gives Inf values a lot, so isn't
    % very sensible.
    % % % %         map(i,:) = log_joint(i,:) - log_evidence;  %exp(log_joint(i,:) - log_evidence);

    map(i,:) = exp(log_joint(i,:) ./ (gam + nfeatures) );  %exp(log_joint(i, :)); % will be very small
    %sum(log_joint) + nfeatures * logpriors;

    %map = map ./ sum(map);

end

end

function [prop_correct, m, misclass, prop_correct_by_class, chance, chance_95_ci] = classify_naive_bayes_eval(true_class, class_est, wh_obs)

if nargin < 3 || isempty(wh_obs)
    % assume we have all observations.
    wh_obs = 1:length(true_class);
end

[m,dp,corr,far,misclass] = confusion_matrix(true_class(wh_obs),class_est);

totaln = sum(m(:));
prop_correct = trace(m) ./ totaln;

nclasses = length(unique(true_class(true_class ~= 0)));

prop_correct_by_class = diag(m) ./ sum(m,2);

chance = 1 ./ nclasses;
chance_95_ci = [chance - sqrt(chance * (1 - chance) ./ totaln)   chance + sqrt(chance * (1 - chance) ./ totaln)];

end

function [Y, whkeep] = eliminate_empty_Y(Y)

fprintf('Eliminating any empty features. \n');

% using whkeep helps avoid out of memory errors.
whzero = (sum(Y) < eps);
whkeep = ~whzero;
Y = Y(:, whkeep);

end

function [Y, whkeep, whomsave] = select_features(Y, Xi, activation_cutoff, selectivity_cutoff, whkeep)

% select features
% -----------------------------------------------------------------

[nobs, nfeatures] = size(Y);
ntasks = size(Xi,2);

sumtasks = sum(Y); sumtasks(sumtasks==0) = 1;  % avoid warning for empty tasks

% p(task | activation) estimates
ptask = (Xi' * Y) ./ repmat(sumtasks,ntasks,1);

maxp = max(ptask);  % max prob of task across vars

pt = (Xi' * Y) ./ repmat(sum(Xi)',1,nfeatures);  % proportion of each task that activated in an area
ptmax = max(pt);
whomsave = maxp > (selectivity_cutoff) & ptmax > activation_cutoff;

if sum(whomsave) == 0
    disp('Warning: No variables meet feature selection criteria.');
    whomsave = (maxp == max(maxp));
end

whkeep(whkeep) = whomsave;  % indices of saved voxels in whole-brain (original) index vector
%whkeep = whkeep(whomsave);  % overall indices of saved voxels
Y = Y(:,whomsave);

%fprintf(1,'Selected: %3.0f ',nvars);
% get weights on voxels based on n activations.... ******

end

function bayes_model = write_images(bayes_model, Y, volInfo, conditionnames)
% bayes_model = write_images(bayes_model, MC_Setup.volInfo, MC_Setup.Xinms)

% Xi = bayes_model.Xi;

%[priors, pa1_given_t, pa0_given_t, varargout] = bayes_get_probabilities(Y, Xi, bayes_model.params.k);

nfeatures_orig = length(bayes_model.whkeep);

pa1_given_t = NaN * zeros(nfeatures_orig, bayes_model.nclasses);
pa1_given_not_t = NaN * zeros(nfeatures_orig, bayes_model.nclasses);

pt_given_act1 = NaN * zeros(nfeatures_orig, bayes_model.nclasses);
pt_given_act0 = NaN * zeros(nfeatures_orig, bayes_model.nclasses);

pa1_given_t(bayes_model.whkeep,:) = bayes_model.pa1_given_t;
pa1_given_not_t(bayes_model.whkeep,:) = bayes_model.pa1_given_not_t;

pt_given_act1(bayes_model.whkeep,:) = bayes_model.pt_given_act1;
pt_given_act0(bayes_model.whkeep,:) = bayes_model.pt_given_act0;

% now done in bayes_get_probabilities
% % pa1_given_not_t = (double(~(Xi)') * Y)' ./ repmat(sum(~Xi), size(Y, 2), 1);
% % pa1_given_not_t(pa1_given_not_t < eps) = 1 ./ bayes_model.nobs;
% % pa1_given_not_t(pa1_given_not_t == 1) = 1 - (1 ./ bayes_model.nobs);

%%%%pa1_given_not_t(~bayes_model.whkeep, :) = NaN;

disp('Writing images to disk in the current directory.');

for i = 1:bayes_model.nclasses
    pa1_given_tname{i} = ['pa|t_' conditionnames{i} '.img'];
    iimg_reconstruct_3dvol(pa1_given_t(:,i), volInfo, 'outname', pa1_given_tname{i});

    disp(pa1_given_tname{i})

    pa_given_nottname{i} = ['pa|nott_' conditionnames{i} '.img'];
    iimg_reconstruct_3dvol(pa1_given_not_t(:,i), volInfo, 'outname', pa_given_nottname{i});

    disp(pa_given_nottname{i})

    pt_given_aname{i} = ['pt|a_' conditionnames{i} '.img'];
    iimg_reconstruct_3dvol(pt_given_act1(:,i), volInfo, 'outname', pt_given_aname{i});

    disp(pt_given_aname{i})

    pt_given_notaname{i} = ['pt|nota_' conditionnames{i} '.img'];
    iimg_reconstruct_3dvol(pt_given_act0(:,i), volInfo, 'outname', pt_given_notaname{i});

    disp(pt_given_notaname{i})
end

% likelihood ratio at each voxel
% should penalize this in some way for low activations...
%pa1_given_t(pa1_given_t == 0) = 1;
%pa1_given_not_t(pa1_given_not_t == 0) = 1;

% higher numbers mean more association between activation and class.
% lr for 2 classes can both be positive IF not all obs. fall within class 1
% or class 2!!!
evidence_yes = (sum(Y) ./ bayes_model.nobs)';
evidence_yes = evidence_yes(:, ones(1, bayes_model.nclasses));

% % priors = NaN * zeros(nfeatures_orig, bayes_model.nclasses);
% % priors(bayes_model.whkeep,:) = bayes_model.priors(ones(bayes_model.nfeatures, bayes_model.nclasses));

% add .1 to regularize...this could be k param
lr = ( (pa1_given_t + .1) ./ (evidence_yes + .1) ) - 1;
%pt_given_act1 ./ pt_given_act0;
lr(lr < .05 & lr > -.05) = NaN;
%lr = log(lr);

for i = 1:bayes_model.nclasses
    lrname{i} = ['likeratio_' conditionnames{i} '.img'];
    iimg_reconstruct_3dvol(lr(:,i), volInfo, 'outname', lrname{i});

    disp(lrname{i})
end

bayes_model.image_output.volInfo = volInfo;
bayes_model.image_output.pa1_given_tname = pa1_given_tname;
bayes_model.image_output.pa_given_nottname = pa_given_nottname;

bayes_model.image_output.pt_given_aname = pt_given_aname;
bayes_model.image_output.pt_given_notaname = pt_given_notaname;

bayes_model.image_output.likerationame = lrname;

bayes_model.image_output.pa1_given_t = pa1_given_t;
bayes_model.image_output.pa1_given_not_t = pa1_given_not_t;
bayes_model.image_output.likeratio = lr;

end

function cl = display_output(bayes_model, disptype, varargin)
% display_output(bayes_model, disptype, varargin)

switch spm('Ver')
    case 'SPM2'
        % spm_defaults is a script
        disp('WARNING: spm defaults not set for spm2. Make sure your defaults are set correctly');

    case 'SPM5'
        % spm_defaults is a function
        spm_defaults()

    case 'SPM8'
        % spm_defaults is a function
        spm_defaults()

    otherwise
        % unknown SPM
        disp('Unknown version of SPM!');
        spm_defaults()
end

cl = [];
colors1 = [1 .7 0];
colors2 = [0 0 1];

for i = 1:length(varargin)
    if iscell(varargin{i})
        colors1 = varargin{i}{1};
        colors2 = varargin{i}{2};
    elseif length(varargin{i}) == 1
        thresh = varargin{i};
    end
end

switch disptype

    case {'pa|t', 'lr', 'pt|a'}

        overlay = which('spm2_single_subj_T1_scalped.img');
        spm_check_registration(repmat(overlay, bayes_model.nclasses, 1));

        for i = 1:bayes_model.nclasses

            switch disptype
                case 'pa|t'
                    cl{i} = mask2clusters(bayes_model.image_output.pa1_given_tname{i});
                    fprintf('Displayed: p(a | t) for each task.\n');

                    spm_orthviews_name_axis(bayes_model.image_output.pa1_given_tname{i}, i);

                case 'lr'
                    cl{i} = mask2clusters(bayes_model.image_output.likerationame{i});
                    fprintf('Displayed: regions with positive likelihood ratio for each task.\n');

                    spm_orthviews_name_axis(bayes_model.image_output.likerationame{i}, i);

                case 'pt|a'
                    cl{i} = mask2clusters(bayes_model.image_output.pt_given_aname{i});
                    fprintf('Displayed: p(t | a) for each task.\n');

                    spm_orthviews_name_axis(bayes_model.image_output.pt_given_aname{i}, i);

            end

            % threshold
            if exist('thresh', 'var')
                for j = 1:length(cl{i})
                    wh_omit = abs( cat(2, cl{i}(j).Z) ) < thresh;
                    cl{i}(j).XYZ(:, wh_omit) = [];
                    cl{i}(j).XYZmm(:, wh_omit) = [];
                    cl{i}(j).Z(:, wh_omit) = [];
                end

                clu = clusters2CLU(cl{i});
                cl{i} = tor_extract_rois([], clu, clu);

            end

            cluster_orthviews(cl{i}, 'handle', i, 'add');

        end

        switch disptype
            case 'pa|t'
                %spm_orthviews_change_colormap([.2 .2 .4], [1 1 0], [.9 .6 .1]);  %slate to orange to yellow
                spm_orthviews_change_colormap([.5 0 0], [1 1 0], [1 .5 0]);  %blue to yellow

            case  'lr'
                spm_orthviews_change_colormap(colors2, colors1);  %blue to white to yellow

            case 'pt|a'
                spm_orthviews_change_colormap(colors2, colors1);

        end

    case 'lr surface'

        poscm = colormap_tor([.2 .2 .4], [1 1 0], [.9 .6 .1]);  %slate to orange to yellow
        negcm = colormap_tor([0 0 1], [0 .3 1]);  % light blue to dark blue (shouldn't be needed!)

        for i = 1:bayes_model.nclasses

            cl{i} = mask2clusters(bayes_model.image_output.likerationame{i});

            fprintf('Displayed: regions with positive likelihood ratio for each task.\n');

            % Left medial and lateral surface
            % ----------------------------------------------------------
            create_figure('Brain Surface'); scn_export_papersetup(800);

            cluster_surf(cl{i}, 4, 'heatmap', 'colormaps', poscm, negcm, 'hires left');
            sh = cluster_surf(cl{i}, 4, 'heatmap', 'colormaps', poscm, negcm, 'brainstem');
            saveas(gcf, [bayes_model.image_output.likerationame{i}(1:end-4) '_left_medial'], 'png');

            view(270,0); lightRestoreSingle;
            saveas(gcf, [bayes_model.image_output.likerationame{i}(1:end-4) '_left_lateral'], 'png');

            % Right medial and lateral surface
            % ----------------------------------------------------------
            create_figure('Brain Surface'); cluster_surf(cl{i}, 4, 'heatmap', 'colormaps', poscm, negcm, 'hires right');
            sh = cluster_surf(cl{i}, 4, 'heatmap', 'colormaps', poscm, negcm, 'brainstem');
            scn_export_papersetup(800); saveas(gcf, [bayes_model.image_output.likerationame{i}(1:end-4) '_right_lateral'], 'png');

            view(270,0); lightRestoreSingle;
            saveas(gcf, [bayes_model.image_output.likerationame{i}(1:end-4) '_right_medial'], 'png');

            % Limbic
            % ---------------------------------------------------------
            create_figure('Brain Surface'); cluster_surf(cl{i}, 4, 'heatmap', 'colormaps', poscm, negcm, 'limbic');
            saveas(gcf, [bayes_model.image_output.likerationame{i}(1:end-4) '_limbic'], 'png');
            saveas(gcf, [bayes_model.image_output.likerationame{i}(1:end-4) '_limbic'], 'fig');
        end

    case 'map plot'

        % MAP plot

        colors = {'b' 'r' 'k' 'g' 'm'};
        create_figure('MAP Plot across observations');
        for i = 1:bayes_model.nclasses

            trueclassi = find(bayes_model.Xi(:,i) == 1);
            myclassmap = bayes_model.apparent.map(trueclassi, i);

            plot(bayes_model.apparent.map(:, i), 'color', colors{i});
            plot(trueclassi, bayes_model.apparent.map(find(bayes_model.Xi(:,i) == 1), i), 'o', 'Color', colors{i},'MarkerFaceColor',colors{i});

            incorrect_obs = ~(myclassmap == max(bayes_model.apparent.map(trueclassi, :), [], 2));
            incorrect_obs_vec = zeros(bayes_model.nobs, 1);
            incorrect_obs_vec(trueclassi(incorrect_obs)) = 1;

            plot(find(incorrect_obs_vec), myclassmap(incorrect_obs), 'kx', 'MarkerSize', 12, 'LineWidth', 2);
        end
        xlabel('Observations'); ylabel('log(MAP) (up to scaling constant)');

    case 'class plot'

        pr = bayes_model.priors;
        fr = [linspace(1, 0, 10)' linspace(0, 1, 10)'];   % frequencies of guesses, from always A to always B

        xfr = linspace(1, 0, 10);

        %hr_byclass = fr .* repmat(pr, size(fr, 1), 1)
        colors = {'b' 'r' 'k' 'g' 'm'};

        hr = fr * pr';

        create_figure('Summary Line Plot');
        %plot(xfr, hr_byclass);
        plot(xfr, hr, 'k', 'LineWidth', 2);

        phighestfreq = sum(bayes_model.xval.class_est == wh) ./ length(bayes_model.xval.class_est);  % p chosen for highest-frequency class
        plot(phighestfreq, bayes_model.xval.prop_correct, 'ko', 'LineWidth', 2, 'MarkerFaceColor', [.5 .5 .5], 'MarkerSize', 12);

        xx = repmat(phighestfreq, 2, bayes_model.nclasses);
        yy = [bayes_model.xval.prop_correct_by_class'; bayes_model.xval.prop_correct_by_class'];

        for i = 1:size(yy, 2)
            plot(xx(:,i), yy(:,i), 'o', 'Color', colors{i}, 'LineWidth', 2, 'MarkerFaceColor', colors{i}, 'MarkerSize', 12);
        end

        % Standard error of hit rate
        % var(hits for each category) = p*(1-p)*N for that category
        % where p is the prior p(correct) for that category, and N is the number of
        % choices in that category given the frequency of sampling that category.

        se_hr = sqrt( hr .* (1 - hr) ./ bayes_model.nobs );
        ci95_hr = 1.96 * se_hr;

        plot(xfr, hr + ci95_hr, 'k--', 'LineWidth', 1);
        plot(xfr, hr - ci95_hr, 'k--', 'LineWidth', 1);
        fill_around_line(hr, ci95_hr, 'k', xfr);

        ylabel('Correct classification rate');
        xlabel('Prop. choices from highest class');

        % N = fr .* bayes_model.nobs;     % for each of k categories, in a row

        % pq = repmat(pr .* (1 - pr), size(fr,1), 1);
        % varsum = sum(pq .* N, 2)
        %
        %
        % se_hr = sqrt( sum(repmat(pr .* (1 - pr), size(fr,1), 1) ./ (bayes_model.nobs .* fr), 2) )
        % se_hr(isinf(se_hr)) = 0;

    otherwise
        error('Unknown display option.')
end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% timewarp()   - Given two event marker vectors, computes a matrix
%                that can be used to warp a time series so that its
%                evlatencies match newlatencies. Values of the warped
%                timeserie that falls between two frames in the original
%                timeserie will be linear interpolated.
% Usage:
%   >> warpmat = timewarp(evlatency, newlatency)
%
% Necessary inputs:
%   evlatency  - [vector] event markers in the original time series, in frames
%                Markers must be ordered by increasing frame latency.
%                If you want to warp the entire time-series, make sure
%                the first (1) and last frames are in the vector.
%   newlatency - [vector] desired warped event time latencies. The original
%                time series will be warped so that the frame numbers of its
%                events (see evlatency above) match the frame numbers in
%                newlatency. newlatency frames must be sorted by ascending
%                latency. Both vectors must be the same length.
%
% Optional outputs:
%      warpmat - [matrix] Multiplying this matrix with the original
%                time series (column) vectors performs the warping.
%
% Example:
%      % In 10-frame vectors, warp frames 3 and 5 to frames 4 and 8,
%      % respectively. Generate a matrix to warp data epochs in
%      % variable 'data' of size (10,k)
%      >> warpmat = timewarp([1 3 5  10], [1 4 8 10])
%      >> warped_data = warpmat*data;
%
% Authors: Jean Hausser, SCCN/INC/UCSD, 2006
%
% See also: angtimewarp(), phasecoher(), erpimage()
%
function M=timewarp(evLatency, newLatency)
  M = [0];
  if min(sort(evLatency) == evLatency) == 0
    error('evLatency should be in ascending order');
    return;
  end
  if min(sort(newLatency) == newLatency) == 0
    error('newLatency should be in ascending order');
    return;
  end
  if length(evLatency) ~= length(newLatency)
    error('evLatency and newLatency must have the same length.');
    return;
  end
  if length(evLatency) < 2 | length(newLatency) < 2
    error(['There should be at least two events in evlatency and ' ...
          'newlatency (e.g., "begin" and "end")'] );
    return;
  end
  if evLatency(1) ~= 1
    disp(['Assuming old and new time series beginnings are synchronized.']);
    disp(['Make sure you have defined an ending event in both the old and new time series!']);
    evLatency(end+1)=1;
    newLatency(end+1)=1;
    evLatency = sort(evLatency);
    newLatency = sort(newLatency);
  end

  t = 1:max(evLatency);

  for k=1:length(evLatency)-1
    for i=evLatency(k):evLatency(k+1)-1
      tp(i) = (t(i)-evLatency(k)) * ...
              (newLatency(k+1) - newLatency(k))/...
              (evLatency(k+1) - evLatency(k)) + ...
              newLatency(k);
    end
  end

  % Check what's going on at tp(max(newLatency)), should equal t(max(evLatency))
  tp(max(evLatency)) = max(newLatency);
  ts = tp-min(newLatency)+1;

% $$$   M = sparse(max(newLatency)-min(newLatency)+1, max(evLatency));
  M = zeros(max(newLatency)-min(newLatency)+1, max(evLatency));

  k = 0;
  for i=1:size(M,1)
    while i > ts(k+1)
      k = k+1;
    end
% $$$     k = k-1;

    if k == 0
      % Check wether i == ts(1) and i == 1
      % In that case, M(1,1) = 1
      M(1,1) = 1;
    else
      M(i,k) = 1 - (i-ts(k))/(ts(k+1)-ts(k));
      M(i,k+1) = 1 - (ts(k+1)-i)/(ts(k+1)-ts(k));
    end
  end

author: PatternRecognition

repo: OpenBMI


text:
% Build a spectrum-like harmonic model using a given sinusoids set
%
% Input
%  sins     : Nx3 (or Nx4), frequency[Hz], amplitude, phase, (harmnb)
%             Only positive frequencies, with DC.
%  opt      : Options (see below)
%
% Output
%  M        : Spectrum built from sins.
%             (With same format as a DFT result)
%
% Copyright (c) 2012 University of Crete - Computer Science Department
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file 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 Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%
% TODO
%  Manage holes in the harmonic structure
%

function [M, opt] = sin2shm(sins, opt)

    if nargin<2
        opt.nbh             = Inf; % Number of harmonics to be use in the model
                                   % (excluding DC).
        opt.max_freq        = 8000;% Consider harmonics up to this limit.

        opt.prevpow2        = false;% if true, reduce the harmonic model to the
                                  % biggest possible power of 2 length
        opt.dc_extrap       = 1;  % 0: Use the current value
                                  % 1; Copy the amplitude of the first harm
                                  % 2: Use linear extrapolation on log scale
                                  % both 1 and 2 assume angle(DC)=0
    end
    if nargin==0; M=opt; return; end

    % TODO manage holes in the harmonic structure

    if ~isempty(opt.max_freq) && opt.max_freq>0
        inds = sins(1,:)<=opt.max_freq;
        sins = sins(:,inds);
    end
    sins = sins(:,1:min(size(sins,2), 1+opt.nbh));
    if opt.prevpow2
        % in order to speed up the computation of the min phase
        harmsize = 2^(nextpow2(size(sins,2))-1) + 1;
        sins = sins(:,1:min(size(sins,2), harmsize));
    end

    M = [sins(2,:).'.*exp(1i*sins(3,:).')];

    M(1) = sign(real(M(1)))*abs(M(1));        % set DC phase to 0 or pi

    M(end) = sign(real(M(end)))*abs(M(end));  % set the Nyquist phase to 0 or pi

    M = hspec2spec(M);

    % Extrapolate the model DC
    % Doesn't keep the DC phase !
    % Angle(DC) never makes sens in acoustic measurement
    % For technical reasons, we can therefore assume any constant (here zero).
    if opt.dc_extrap==1;     M(1) = abs(M(2));
    elseif opt.dc_extrap==2; M(1) = db2lin(lin2db(M(2))+(lin2db(M(2))-lin2db(M(3))));   end

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function net = dsift_net(binSize, varargin)
% Define a CNN equivalent to dense SIFT

opts.numOrientations = 4 ;
opts = vl_argparse(opts, varargin) ;
opts.binSize = binSize ;

% Spatial derivatives along NO directions
NO = opts.numOrientations ;
dx = [0 0 0 ; -1 0 1 ; 0  0 0]/2 ;
dy = dx' ;
for i=0:2*NO-1
  t = (2*pi)/(2*NO)*i ;
  spatialDer{i+1} = cos(t)*dx+sin(t)*dy;
end
spatialDer = single(cat(4, spatialDer{:})) ;

% Spatial bilienar binning
a = 1 - abs(((1:2*opts.binSize) - (2*opts.binSize+1)/2)/opts.binSize);
bilinearFilter = repmat(single(a'*a), [1, 1, 1, 2*NO]) ;

% Stacking of SIFT cells into 4x4 blocks.

sigma = 1.5 ;
mask = {} ;
t = 0 ;
for i=1:4
  for j=1:4
    for o=1:2*NO
      t=t+1 ;
      mask{t} = zeros(4,4,2*NO) ;
      mask{t}(i,j,o) = exp(-0.5*((i-2.5).^2 + (j-2.5).^2) / sigma^2) ;
    end
  end
end
mask = single(cat(4, mask{:})) ;

net.layers = {} ;
net.layers{end+1} = struct('type','conv', ...
  'filters', spatialDer, ...
  'biases', zeros(size(spatialDer,4),1,'single'), ...
  'stride', 1, 'pad', 0) ;
if 0
  net.layers{end+1} = struct('type','noffset', 'param', [.5*cos(2*pi/8), .5]) ;
  net.layers{end+1} = struct('type','relu') ;
else
  net.layers{end+1} = get_hog_binning_layer(opts) ;
end
net.layers{end+1} = struct('type','conv', ...
  'filters', bilinearFilter, ...
  'biases', [], ...
  'stride', binSize, 'pad', 0) ;
net.layers{end+1} = struct('type','conv', ...
  'filters', mask, ...
  'biases', [], ...
  'stride', 3, 'pad', 0) ;
net.layers{end+1} = struct('type', 'normalize', 'param', [128*2, 0.000001, 1, .5]) ;
net.layers{end+1} = get_hog_clamp_layer(opts) ;

% -------------------------------------------------------------------------
function l = get_hog_clamp_layer(opts)
% -------------------------------------------------------------------------
l.type = 'custom' ;
l.forward = @hog_clamp_forward ;
l.backward = @hog_clamp_backward ;

function res_ = hog_clamp_forward(ly,res,res_)
res_.x = min(res.x,single(0.2)) ;
%res_.x = res.x ;

function res = hog_clamp_backward(ly,res,res_)
res.dzdx = res_.dzdx .* (res.x <= 0.2) ;

% -------------------------------------------------------------------------
function l = get_hog_binning_layer(opts)
% -------------------------------------------------------------------------
l.type = 'custom' ;
l.NO = opts.numOrientations ;
l.forward = @hog_binning_forward ;
l.backward = @hog_binning_backward ;

function res_ = hog_binning_forward(ly,res,res_)
x = res.x ;
n2 = sum(x.^2,3)/ly.NO ;
n = sqrt(n2) ;
cs = bsxfun(@rdivide, x, max(n, 1e-10)) ;
cs = max(min(cs,1),-1) ;
delta = 1 - (2*ly.NO)/(2*pi)*acos(cs) ;
w = max(0, delta) ;
res_.x = bsxfun(@times, n, w) ;

function res = hog_binning_backward(ly,res,res_)
% forward computations
x = res.x ;
n2 = sum(x.^2,3)/ly.NO ;
n = sqrt(n2) ;
cs = bsxfun(@rdivide, x, max(n, 1e-10)) ;
cs = max(min(cs,1),-1) ;
delta = 1 - (2*ly.NO)/(2*pi)*acos(cs) ;
w = max(0, delta) ;

dn = bsxfun(@rdivide, x/ly.NO, max(n, 1e-10)) ;
dwdcs = ((2*ly.NO)/(2*pi) ./ sqrt(1.0000001 - cs.^2)) .* (delta > 0) .* res_.dzdx ;
dw = ...
  + dwdcs .* repmat(1./n, [1 1 2*ly.NO]) ...
  - bsxfun(@times, sum(bsxfun(@rdivide, dwdcs.*x/ly.NO, n.*n2),3), x) ;
res.dzdx = bsxfun(@times, sum(res_.dzdx .* w, 3), dn) +  bsxfun(@times, n, dw) ;

author: aravindhm

repo: deep-goggle


text:
function [A,C,newcenters] = manually_refine_components(Y,A,C,centers,img,sx,options)
% function that allows to add or remove components.
% Usage:
% after runnning use left click to remove and right click to add components
% at specific locations. Hit enter when done
% Inputs
% ------
% centers: coordinates of the centers of the components
% img: image used to identify components (for instance NCA
%
% Returns:
% --------
% newcenters: updated component centers

defoptions = CNMFSetParms;
if nargin < 7 || isempty(options); options = defoptions; end
if ~isfield(options,'d1') || isempty(options.d1); options.d1 = input('What is the total number of rows? \n'); end          % # of rows
if ~isfield(options,'d2') || isempty(options.d2); options.d2 = input('What is the total number of columns? \n'); end       % # of columns
if ~isfield(options,'cont_threshold') || isempty(options.cont_threshold); cont_threshold = defoptions.cont_threshold; else cont_threshold = options.cont_threshold; end          % # of rows
if nargin < 6 || isempty(sx)
    sx = 5;
end
if nargin < 5 || isempty(img)
    img = std(Y,[],3);
end
if nargin < 4 || isempty(centers)
    centers = com(A,options.d1,options.d2);
end

min_distance_point_selection=2;
x=1;
T = size(C,2);
newcenters=[centers];
% ident_point=[0,0];
fig = figure;
imagesc(img,[min(img(:)),max(img(:))]);
    axis equal; axis tight; hold all;
    scatter(newcenters(:,2),newcenters(:,1),'mo'); hold on;
    title('Center of ROIs found from initialization algorithm');
    xlabel({'Press left click to add new component, right click to remove existing component'; 'Press enter to exit'},'fontweight','bold');
    drawnow;
    cmap = colormap;
for i = 1:size(A,2)
    a_srt = sort(A(:,i),'descend');
    ff = find(cumsum(a_srt.^2) >= cont_threshold*sum(a_srt.^2),1,'first');
    contour(reshape(A(:,i),options.d1,options.d2),[0,0]+a_srt(ff),'Linecolor',[1,0,1]/2);
    hold on;
end

while ~isempty(x)
%     scatter(ident_point(1),ident_point(2),'go');
    [x,y,button]=ginput(1);
    disp(button); hold on;
    if ~isempty(x)
        pixel=round([x y]);
        if button==1
            disp(['Adding pixel at:' num2str(fliplr(pixel))])
            newcenters=[newcenters; fliplr(pixel)];
            int_x = round(newcenters(end,1)) + (-sx:sx);
            if int_x(1)<1
                int_x = int_x + 1 - int_x(1);
            end
            if int_x(end)>options.d1
                int_x = int_x - (int_x(end)-options.d1);
            end
            int_y = round(newcenters(end,2)) + (-sx:sx);
            if int_y(1)<1
                int_y = int_y + 1 - int_y(1);
            end
            if int_y(end)>options.d2
                int_y = int_y - (int_y(end)-options.d2);
            end
            Ypatch = reshape(Y(int_x,int_y,:),(2*sx+1)^2,size(C,2));
            Ypatch = bsxfun(@minus, Ypatch, median(Ypatch,2));
            [INT_x,INT_y] = meshgrid(int_x,int_y);
            coor = sub2ind([options.d1,options.d2],INT_x(:),INT_y(:));
            Y_res = Ypatch - A(coor,:)*C;
            [atemp, ctemp, ~, ~, newcenter, ~] = greedyROI(reshape(Y_res,2*sx+1,2*sx+1,T), 1, options);
            %[atemp, ctemp] = initialize_components(reshape(Y_res,2*sx+1,2*sx+1,T), 1,sx,options);  % initialize
            % find contour
            a_srt = sort(atemp,'descend');
            ff = find(cumsum(a_srt.^2) >= cont_threshold*sum(a_srt.^2),1,'first');
            A(coor,end+1) = atemp/norm(atemp);
            C(end+1,:) = ctemp*norm(atemp);
            new_center = com(A(:,end),options.d1,options.d2);
            newcenters(end,:) = new_center;
            scatter(new_center(2),new_center(1),'mo'); hold on;
            contour(reshape(A(:,end),options.d1,options.d2),[0,0]+a_srt(ff),'Linecolor',[1,0,1]/2);
            hold on;
            colormap(fig,cmap);

            drawnow;
%             atemp = max(mean(Y_res,2),0);
%             for i = 1:10
%                ctemp = max(atemp'*Y_res,0)/norm(atemp)^2;
%                atemp = max(Y_res*ctemp',0)/norm(ctemp)^2;
%             end

        elseif button==3
            [m,id]=min(sum(bsxfun(@minus,newcenters,[y,x]).^2,2));
            ident_point=[newcenters(id,2),newcenters(id,1)];
            if m<=min_distance_point_selection
                disp(['Removing point:' num2str(ident_point)])
                newcenters(id,:)=[];
                A(:,id) = [];
                C(id,:) = [];
                % replot after removing
                clf;
                imagesc(img,[min(img(:)),max(img(:))]);
                    axis equal; axis tight; hold all;
                    scatter(newcenters(:,2),newcenters(:,1),'mo'); hold on;
                    title('Center of ROIs found from initialization algorithm');
                    xlabel({'Press left click to add new component, right click to remove existing component'; 'Press enter to exit'},'fontweight','bold');
                    drawnow;
                    cmap = colormap;
                for i = 1:size(A,2)
                    a_srt = sort(A(:,i),'descend');
                    ff = find(cumsum(a_srt.^2) >= cont_threshold*sum(a_srt.^2),1,'first');
                    contour(reshape(A(:,i),options.d1,options.d2),[0,0]+a_srt(ff),'Linecolor',[1,0,1]/2);
                    hold on;
                end
            else
                disp('Selection too far from any point')
            end
        end
    end
end

end

author: zhoupc


text:
function [ ndp, xdp, ydp ] = dif_deriv ( nd, xd, yd )

%*****************************************************************************80
%
%% DIF_DERIV computes the derivative of a polynomial in divided difference form.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 June 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carl deBoor,
%    A Practical Guide to Splines,
%    Springer, 2001,
%    ISBN: 0387953663,
%    LC: QA1.A647.v27.
%
%  Parameters:
%
%    Input, integer ND, the size of the input table.
%
%    Input, real XD(ND), the abscissas for the divided
%    difference table.
%
%    Input, real YD(ND), the divided difference table.
%
%    Output, integer NDP, the size of the output table, which is ND-1.
%
%    Input, real XDP(NDP), the abscissas for the divided
%    difference table for the derivative.
%
%    Output, real YDP(NDP), the divided difference
%    table for the derivative.
%

%  Using a temporary copy of the difference table, shift the
%  abscissas to zero.
%
  xd_temp(1:nd) = xd(1:nd);
  yd_temp(1:nd) = yd(1:nd);

  [ xd_temp, yd_temp ] = dif_shift_zero ( nd, xd_temp, yd_temp );
%
%  Construct the derivative.
%
  ndp = nd - 1;

  xdp(1:ndp) = 0.0;

  for i = 1 : ndp
    ydp(i) = i * yd_temp(i+1);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [huff entropy avglength redundancy]=huffman(alpha,prob)
s=sum(prob(:));
s=roundn(s,-4);
% Calculate length of source and probability
la=length(alpha);
lp=length(prob);
if (la==lp & s==1)
       %Calculate the Entropy
         entropy=prob.*log2(prob);
         entropy=-sum(entropy(:));
         pos=1:lp;
         [prs idx]=sort(prob,'descend');
         npos=pos(idx);
         %Minimum Variance
         idx=find(prs==min(prs(:)));
         tp=npos(idx);
         tp=sort(tp,'descend');
         npos(idx)=tp;
         %
         codebook(1:lp)={''};
         ps=npos;
         np=lp;
         cb=zeros([lp-1 3]);
         cnt=lp+1;
         prb=prs;
         for i=1:lp-1
              fst=ps(np-1);
              sec=ps(np);

                if fst<=lp
                   codebook(fst)=strcat('0',codebook(fst));
                else
                  codebook=encod(fst,cb,codebook,'0',lp);

                end

                if sec<=lp
                    codebook(sec)=strcat('1',codebook(sec));
                else
                   codebook=encod(sec,cb,codebook,'1',lp);
                end
              cb(i,1)=cnt;
              cb(i,2)=fst;
              cb(i,3)=sec;
              if np>2
                  ps=ps(1:np-2);
                  ps(np-1)=cnt;
                  cnt=cnt+1;
                  prbt=prb(1:np-2);
                  prbt(np-1)=prb(np-1)+prb(np);
                  prb=prbt;
                  [prb idx]=sort(prb,'descend');
                  ps=ps(idx);
                  %Minimum Variance
                  idx=find(prb==prbt(np-1));
                  tp=ps(idx);
                  tp=sort(tp,'descend');
                  ps(idx)=tp;
                  %
                  np=np-1;
              end

         end
         for i=1:lp
            huff(i).sym=alpha(i);
            huff(i).prob=prob(i);
            huff(i).code=codebook(i);
         end
           avglength=0;
           for i=1:lp
               avglength=avglength+huff(i).prob*length(cell2mat(huff(i).code));
           end
           redundancy=((avglength-entropy)/entropy)*100;
else
    display('Error in input.....');
    huff=[];
end

function codebook=encod(fs,cb,codebook,str,lp)
idx=find(cb(:,1)==fs);
x=cb(idx,2);
y=cb(idx,3);
 if x<=lp & y<=lp
     codebook(x)=strcat(str,codebook(x));
     codebook(y)=strcat(str,codebook(y));
     return
 elseif x<=lp
     codebook(x)=strcat(str,codebook(x));
     codebook=encod(y,cb,codebook,str,lp);
 elseif y<=lp
     codebook(y)=strcat(str,codebook(y));
     codebook=encod(x,cb,codebook,str,lp);
 else
     codebook=encod(x,cb,codebook,str,lp);
     codebook=encod(y,cb,codebook,str,lp);
 end

author: Sable

repo: mcbench-benchmarks


text:
Network Resnet50 {
	Layer CONV1 {
		Type: CONV
		Stride { X: 2, Y: 2 }
		Dimensions { K: 64, C: 3, R: 7, S: 7, Y:224, X:224 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV2_1_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 64, C: 64, R: 1, S: 1, Y: 56, X: 56 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV2_1_2 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 64, C: 64, R: 3, S: 3, Y: 56, X: 56 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV2_1_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 256, C: 64, R: 1, S: 1, Y: 56, X: 56 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV2_1_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 256, R: 1, S: 1, Y: 56, X: 56 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV2_2_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 64, C: 256, R: 1, S: 1, Y: 56, X: 56 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV2_2_2 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 64, C: 64, R: 3, S: 3, Y: 56, X: 56 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV2_2_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 256, C: 64, R: 1, S: 1, Y: 56, X: 56 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV2_2_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 256, R: 1, S: 1, Y: 56, X: 56 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV2_3_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 64, C: 256, R: 1, S: 1, Y: 56, X: 56 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV2_3_2 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 64, C: 64, R: 3, S: 3, Y: 56, X: 56 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV2_3_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 256, C: 64, R: 1, S: 1, Y: 56, X: 56 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV2_3_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 256, R: 1, S: 1, Y: 56, X: 56 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV3_1_1 {
		Type: CONV
		Stride { X: 2, Y: 2 }
		Dimensions { K: 128, C: 256, R: 1, S: 1, Y: 56, X: 56 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV3_1_2 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 128, C: 128, R: 3, S: 3, Y: 28, X: 28 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV3_1_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 512, C: 128, R: 1, S: 1, Y: 28, X: 28 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV3_1_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 512, R: 1, S: 1, Y: 28, X: 28 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV3_2_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 128, C: 256, R: 1, S: 1, Y: 28, X: 28 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV3_2_2 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 128, C: 128, R: 3, S: 3, Y: 28, X: 28 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV3_2_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 512, C: 128, R: 1, S: 1, Y: 28, X: 28 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV3_1_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 512, R: 1, S: 1, Y: 28, X: 28 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV3_3_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 128, C: 256, R: 1, S: 1, Y: 28, X: 28 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV3_3_2 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 128, C: 128, R: 3, S: 3, Y: 28, X: 28 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV3_3_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 512, C: 128, R: 1, S: 1, Y: 28, X: 28 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV3_3_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 512, R: 1, S: 1, Y: 28, X: 28 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV3_4_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 128, C: 256, R: 1, S: 1, Y: 28, X: 28 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV3_4_2 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 128, C: 128, R: 3, S: 3, Y: 28, X: 28 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV3_4_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 512, C: 128, R: 1, S: 1, Y: 28, X: 28 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV3_4_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 512, R: 1, S: 1, Y: 28, X: 28 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_1_1 {
		Type: CONV
		Stride { X: 2, Y: 2 }
		Dimensions { K: 256, C: 512, R: 1, S: 1, Y: 28, X: 28 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_1_2 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 256, C: 256, R: 3, S: 3, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_1_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1024, C: 256, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_1_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 1024, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_2_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 256, C: 1024, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_2_2 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 256, C: 256, R: 3, S: 3, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_2_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1024, C: 256, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_2_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 1024, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_3_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 256, C: 1024, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_3_2 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 256, C: 256, R: 3, S: 3, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_3_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1024, C: 256, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_3_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 1024, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer CONV4_4_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 256, C: 1024, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_4_2 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 256, C: 256, R: 3, S: 3, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_4_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1024, C: 256, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_4_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 1024, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_5_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 256, C: 1024, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_5_2 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 256, C: 256, R: 3, S: 3, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_5_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1024, C: 256, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_5_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 1024, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_6_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 256, C: 1024, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_6_2 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 256, C: 256, R: 3, S: 3, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_6_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1024, C: 256, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV4_6_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 1024, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV5_1_1 {
		Type: CONV
		Stride { X: 2, Y: 2 }
		Dimensions { K: 512, C: 1024, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV5_1_2 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 512, C: 512, R: 3, S: 3, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV5_1_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 2048, C: 512, R: 1, S: 1, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV5_1_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 2048, R: 1, S: 1, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV5_2_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 512, C: 2048, R: 1, S: 1, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV5_2_2 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 512, C: 512, R: 3, S: 3, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV5_2_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 2048, C: 512, R: 1, S: 1, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV5_2_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 2048, R: 1, S: 1, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV5_3_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 512, C: 2048, R: 1, S: 1, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV5_3_2 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 512, C: 512, R: 3, S: 3, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV5_3_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 2048, C: 512, R: 1, S: 1, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer CONV5_3_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 2048, R: 1, S: 1, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

	Layer FC1000 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1000, C: 2048, R: 7, S: 7, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}

	}

}

author: maestro-project

repo: maestro


text:
function hypervolume = hypervolume3D(F,ub,lb)
% Copyright (c) 2011, Johannes
% 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
%
% 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.
%
% Method for ND objective function values as described in:
%
% 'M. Fleischer. The measure of Pareto Optima Applications to
%  Multi-objective Metaheuristics. EMO 2003, LNCSS 2632
%  519-533, 2003.'
%
% Author: Johannes W. Kruisselbrink
% Last modified: March 17, 2011

	hypervolume = 0;
    F  = -F' + ones(size(F'));
	[M, l] = size(F); ub = ub + ones(1,M); ub = ub';

	% Remove the duplicates from F and compute Lebesque measure of L
	L = unique(F','rows')';

	while l >= 1
		if (l > 1)
			b = zeros(M,1);
			spawn_vector = repmat(L(:,1), 1, M);
			for i = 1:M
				% Bound b(i) is either the least upper bound of the i-th value of the
				% other points or it is the value of the absolute upper bound ub(i)
				difL = (L(i,2:end) - L(i,1));
				lub = find(difL > 0);
				if (length(lub) > 0)
					b(i) = min(L(i,lub+1));
				else
					b(i) = ub(i);
				end

				b(i) = min((difL > 0) .* L(i,2:end) + (difL <= 0) * ub(i));
				% Update i-th spawn vector
				spawn_vector(i,i) = b(i);
			end

			% Compute lop-Off volume and update lebesgue measure
			lov = prod(b - L(:,1));
			hypervolume = hypervolume + lov;

			% Remove L(:,1) from L
			L = L(:,2:end);

			% Add the spawn_vector to L, but first filter dominated
			% solutions and the solutions that touch the upper bounds
			% from the spawn_vector
			L = nd_filter(L, spawn_vector, ub);

		else
			lov = prod(ub - L(:,1));
			hypervolume = hypervolume + lov;
			L = [];
		end
		% Update l
		[M, l] = size(L);
    end
end

function L = nd_filter(L, spawn_vector, ub)
% Implementation of the filter routine as described in:
%
% 'M. Fleischer. The measure of Pareto Optima Applications
%  to Multi-objective Metaheuristics. EMO 2003, LNCSS 2632
%  519-533, 2003.'
%
% Author: Johannes W. Kruisselbrink
% Last modified: March 17, 2011

	[M, l_L] = size(L);
	[M, l_sp] = size(spawn_vector);
	do_assign = zeros(1, l_sp);

	for i = 1 : l_sp

		% Find if the spawnvector hits the upper bound
		at_ub = false;
		for j = 1:M
			if (spawn_vector(j,i) == ub(j))
				at_ub = true;
				break;
			end
		end

		% For this if statement, the following would be more elegant
		% (replacing the loop above), but less efficient:
		% if (all(spawn_vector(:,i) ~= ub))
		if (at_ub == false)
			do_assign(i) = 1;
			for j = 1 : l_L
				if (weakly_dominates(L(:,j), spawn_vector(:,i)))
					do_assign(i) = 0;
					break;
				end
			end
		end
	end
	L = [spawn_vector(:,find(do_assign == 1)), L];
end

function d = weakly_dominates(fA, fB)
% [d] = weakly_dominates(fA, fB)
%
% Compares two solutions A and B given their objective function
% values fA and fB. Returns whether A weakly dominates B.
%
% Input:
% - fA					- The objective function values of solution A
% - fB					- The objective function values of solution B
%
% Output:
% - d					- d is 1 if fA dominates fB, otherwise d is 0
%
% Author: Johannes W. Kruisselbrink
% Last modified: March 17, 2011

	% Elegant, but not very efficient
	%d = (all(fA <= fB) && any(fA < fB));

	% Not so elegant, but more efficient
	d = true;
	for i = 1:length(fA)
		if (fA(i) > fB(i))
			d = false;
			return
		end
	end
end

author: Eric-Bradford

repo: TS-EMO


text:
function in=tril(in)

precision=in(1).precision;
for ii=1:size(in,1)
 if ii+1<=size(in,2)
  in(ii,ii+1:end)=mp(0,precision);
 end
end

author: opencobra

repo: cobratoolbox


text:
function r8ge_print ( m, n, a, title )

%*****************************************************************************80
%
%% R8GE_PRINT prints a R8GE matrix.
%
%  Discussion:
%
%    The R8GE storage format is used for a general M by N matrix.  A storage
%    space is made for each logical entry.  The two dimensional logical
%    array is mapped to a vector, in which storage is by columns.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows of the matrix.
%    M must be positive.
%
%    Input, integer N, the number of columns of the matrix.
%    N must be positive.
%
%    Input, real A(M,N), the R8GE matrix.
%
%    Input, string TITLE, a title to be printed.
%
  r8ge_print_some ( m, n, a, 1, 1, m, n, title );

  return
end


text:
function [Gal] = nmps22Gal(nmps2)
% Convert acceleration from nanometers per square-second to galileos
% Chad A. Greene 2012
Gal = nmps2*1e-7;


text:
function [varargout]=rhombicDodecahedronMesh(varargin)

% function [E,V,C,F,CF,CFF]=rhombicDodecahedronMesh(r,nCopies)
% ------------------------------------------------------------------------
% Creates a rhombic dodecahedron mesh where r sets the radias and nCopies
% (a 1x3 vector) sets the number of copies in the x, y, and z direction.
% The output consists of:
%
% Fc_Q, Fc_T: the quadrilateral and triangular face cell arrays (1 cell
% entry per element).
%
% Ft_Q, Ft,T: the quadrilateral and triangular face arrays
%
% Ct_Q, Ct,T: color/label data for the face arrays
%
% Vt: the vertex array
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2019/02/08 Created
% 2019/10/13 Changed orientation for gridding to create more regular grid
% ------------------------------------------------------------------------

%% Parse input

switch nargin
    case 1
        r=varargin{1};
        nCopies=2;
    case 2
        r=varargin{1};
        nCopies=varargin{2};
end

if isempty(nCopies)
    nCopies=2;
end

if numel(nCopies)==1
    nCopies=nCopies*ones(1,3);
end

%%

% Get single rhombic dodecahedron
[~,Vs]=rhombicDodecahedron(r);
V_offsets=2*r*eye(3,3); %Offset vectors

% Copy over rhombic dodecahedron in grid like fashion 4 times and shift
% some of the copies to nest them together.
nCopies1=nCopies;
[Vg1]=gridClone(Vs,nCopies1,V_offsets);

nCopies2=nCopies;
nCopies2([1,2])=nCopies2([1,2])-1;
[Vg2]=gridClone(Vs,nCopies2,V_offsets);

nCopies3=nCopies;
nCopies3([1,3])=nCopies3([1,3])-1;
[Vg3]=gridClone(Vs,nCopies3,V_offsets);

nCopies4=nCopies;
nCopies4([2,3])=nCopies4([2,3])-1;
[Vg4]=gridClone(Vs,nCopies4,V_offsets);

%%

Vg2(:,1)=Vg2(:,1)+r;
Vg2(:,2)=Vg2(:,2)+r;

Vg3(:,1)=Vg3(:,1)+r;
Vg3(:,3)=Vg3(:,3)+r;

Vg4(:,2)=Vg4(:,2)+r;
Vg4(:,3)=Vg4(:,3)+r;

V=[Vg1;Vg2;Vg3;Vg4];

%% Create element and face arrays
E=reshape((1:1:size(V,1)),14,size(V,1)/14)';
C=(1:1:size(E,1))'; %Element colors
[F,CF,CFF]=element2patch(E,C);

%% Merging point and fix element and face indices
[F,V,~,indFix]=mergeVertices(F,V);
E=indFix(E);

%% Collect output
varargout{1}=E;
varargout{2}=V;
varargout{3}=C;
varargout{4}=F;
varargout{5}=CF;
varargout{6}=CFF;

end

%%

function [Vg]=gridClone(Vs,nCopies,V_offsets)

nTotal=prod(nCopies); %Total number of copies
indC=ones(size(Vs,1),1)*(1:1:nTotal);
indC=indC(:);
[I,J,K] = ind2sub(nCopies,indC);
I=I-1; J=J-1; K=K-1;

%Defining offsets
D1=I*V_offsets(1,:);
D2=J*V_offsets(2,:);
D3=K*V_offsets(3,:);

%Defining vertices matrix
Vg=repmat(Vs,nTotal,1)+(D1+D2+D3);

end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function Y = inv_chol(L)
% Matrix Inversion using Cholesky Decomposition
%
% Finds the inverse of the matrix X, given its (lower triangular) Cholesky
% Decomposition; i.e. X = LL', according to the paper 'Matrix Inversion
% Using Cholesky Decomposition', Aravindh Krishnamoorthy, Deepak Menon,
% arXiv:1111.4144.
%

% Version 0.1, 2013-05-25, Aravindh Krishnamoorthy
% e-mail: aravindh.k@ieee.org

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initializations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = size(L, 1) ;
Y = zeros(N, N) ;
% Work with the upper triangular matrix
R = L' ;
% Construct the auxillary diagonal matrix S = 1/rii
S = inv(diag(diag(R))) ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j=N:-1:1
    for i=j:-1:1
        Y(i,j) = S(i,j) - R(i,i+1:end)*Y(i+1:end,j) ;
        Y(i,j) = Y(i,j)/R(i,i) ;
        % Write out the symmetric element
        Y(j,i) = conj(Y(i,j)) ;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Test
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% N = 4 ;
% X = complex(randn(N, N), randn(N, N)) ;
% X = X'*X ;
% L = chol(X, 'lower') ;
% I = inv_chol(L) ;
% norm(I-inv(X))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


text:
function [y] = symextend(x,Nnum)

y(:,1:Nnum) = fliplr(x(:,1:Nnum));
y = [fliplr(x(:,1:Nnum)) x x(:,end:-1:end-Nnum+1)];
y = [flipud(y(1:Nnum,:)); y ;y(end:-1:end-Nnum+1,:)];

author: thomaskuestner


text:
%IM_MINF Fixed mapping for minimum filter (DIP_Image) (DIP_Image)
%
%	B = IM_MINF(A,SIZE,SHAPE)
%	B = A*IM_MINF([],SIZE,SHAPE)
%	B = A*IM_MINF(SIZE,SHAPE)
%
% INPUT
%   A        Dataset with object images dataset (possibly multi-band)
%   SIZE     Filter width in pixels, default SIZE = 7
%   SHAPE    String with shape:'rectangular', 'elliptic', 'diamond'
%            Default: elliptic
%
% OUTPUT
%   B        Dataset with filtered images
%
% SEE ALSO (PRTools Guide)
% DATASETS, DATAFILES, DIP_IMAGE, MINF

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function b = im_minf(varargin)

	argin = shiftargin(varargin,'scalar');
  argin = setdefaults(argin,[],7,'elliptic');
  if mapping_task(argin,'definition')
    b = define_mapping(argin,'fixed');
    b = setname(b,'Minimum filter');
  else
    [a,size,shape] = deal(argin{:});
    if isa(a,'prdataset') % allows datafiles too
      isobjim(a);
      b = filtim(a,mfilename,{size,shape});
    elseif isa(a,'double') || isa(a,'dip_image') % here we have a single image
      if checktoolbox('dipimage')
        a = 1.0*dip_image(a);
        b = minf(a,size,shape);
      else
        diplibwarn
        %prwarning(1,'Rectangular shape only')
        shape = ones(size,size);
        b = ordfilt2(a,1,shape);
      end
    end
  end

return

author: marianux

repo: ecg-kit

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


text:
classdef Filter_PDE_Density < Filter

    properties (Access = private)
        epsilon
        Acomp
        Anodal2Gauss
        M
        x_reg
        LHS
        bc
        quadrature
    end

    methods (Access = public)

        function obj = Filter_PDE_Density(cParams)
            obj.init(cParams);
            obj.createQuadrature();
            obj.computeBoundaryConditions();
            obj.createMassMatrix();
            obj.epsilon = cParams.mesh.computeMeanCellSize();
            obj.Anodal2Gauss = obj.computeA();
            lhs = obj.createProblemLHS();
            obj.LHS = decomposition(lhs);
        end

        function x0 = getP0fromP1(obj,x)
            obj.x_reg =  obj.getP1fromP1(x);
            x0 = zeros(obj.mesh.nelem,obj.quadrature.ngaus);
            for igaus = 1:obj.quadrature.ngaus
                x0(:,igaus) = obj.Anodal2Gauss{igaus}*obj.x_reg;
            end
        end

        function RHS = integrate_L2_function_with_shape_function(obj,x)
            RHS = obj.M*x;
        end

        function obj = updateEpsilon(obj,epsilon)
            if obj.hasEpsilonChanged(epsilon)
                obj.epsilon = epsilon;
                lhs = obj.createProblemLHS();
                obj.LHS = decomposition(lhs);
            end
        end

        function x_reg = getP1fromP1(obj,x)
            RHS = obj.integrate_L2_function_with_shape_function(x);
            x_reg = obj.solveFilter(RHS);
        end

        function x_reg = getP1fromP0(obj,x0)
            quad = Quadrature.set(obj.mesh.type);
            quad.computeQuadrature('LINEAR');
            s.dV = obj.mesh.computeDvolume(quad)';
            s.x        = x0;
            RHS        = obj.Acomp.integrateP1FunctionWithShapeFunction(s);
            x_reg      = obj.solveFilter(RHS);
        end

    end

    methods (Access = private)

        function createQuadrature(obj)
            q = Quadrature.set(obj.mesh.type);
            q.computeQuadrature('LINEAR');
            obj.quadrature = q;
        end

        function createMassMatrix(obj)
            s.type   = 'MassMatrix';
            s.mesh   = obj.mesh;
            s.fun    = P1Function.create(obj.mesh, 1);
            s.quadratureOrder = 'QUADRATICMASS';
            lhs = LHSintegrator.create(s);
            obj.M = lhs.compute();
        end

        function A_nodal_2_gauss = computeA(obj)
            p1f = P1Function.create(obj.mesh,1);
            s.nnode   = obj.mesh.nnodeElem;
            s.nelem   = obj.mesh.nelem;
            s.npnod   = obj.mesh.nnodes;
            s.ngaus   = obj.quadrature.ngaus;
            s.connec  = obj.mesh.connec;
            s.shape   = p1f.computeShapeFunctions(obj.quadrature);
            obj.Acomp = Anodal2gausComputer(s);
            obj.Acomp.compute();
            A_nodal_2_gauss = obj.Acomp.A_nodal_2_gauss;
        end

        function itHas = hasEpsilonChanged(obj,eps)
            if isempty(obj.epsilon)
                obj.epsilon = 0;
            end
            var = abs(eps - obj.epsilon)/eps;
            itHas = var > 1e-15;
        end

        function x_reg = solveFilter(obj,RHS)
            RHS = obj.bc.fullToReducedVector(RHS);
            s.type = 'DIRECT';
            Solv = Solver.create(s);
            x = Solv.solve(obj.LHS,RHS);
            x_reg = obj.bc.reducedToFullVector(x);
        end

        function computeBoundaryConditions(obj)
            s.scale        = obj.femSettings.scale;
            s.mesh         = obj.mesh;
            s.bc{1}.dirichlet = [];
            s.bc{1}.pointload = [];
            s.bc{1}.ndimf     = [];
            s.bc{1}.ndofs     = [];
            s.ndofs        = obj.mesh.nnodes;
            obj.bc         = BoundaryConditions(s);
            obj.bc.compute();
        end

        function lhs = createProblemLHS(obj)
            s          = obj.femSettings;
            s.mesh     = obj.mesh;
            s.type     = obj.LHStype;
            problemLHS = LHSintegrator.create(s);
            lhs        = problemLHS.compute(obj.epsilon);
            lhs        = obj.bc.fullToReducedMatrix(lhs);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function ms_regular_training(video, method, frames)

opts.expDir = ['net/' method '/' num2str(frames) '/' video] ;

opts.train.batchSize = 5 ;
opts.train.numEpochs = 20;
opts.train.continue = false ;
opts.train.useGpu = true ;
opts.train.learningRate    = 1e-3;
opts.train.expDir = opts.expDir ;

scales = [1, 0.75, 0.5];

%opts = vl_argparse(opts, varargin) ;

% --------------------------------------------------------------------
%                                                         Prepare data
% --------------------------------------------------------------------

imgDir = ['../CDNetDataset/' method '/' num2str(frames) 'frames/' video '/input'];
labelDir = ['../CDNetDataset/' method '/' num2str(frames) 'frames/' video '/GT'];
imdb = getImdb(imgDir,labelDir);

mask = imread(['../CDNetDataset/' method '/' num2str(frames) 'frames/' video '/ROI.bmp']);
mask = mask(:,:,1);
A = max(max(mask));
mask(mask == A) = 1;
if size(mask,1) > 400 || size(mask,2) >400
    mask = imresize(mask, 0.5, 'nearest');
end
imdb.mask = single(double(mask));

imdb.half_size = 15;

%%%%%%Yi%%%%%% redefined the net
load('net');

net.layers{end-1} = struct('type', 'conv', ...
    'filters', 0.1*randn(1,1,64,1, 'single'), ...
    'biases', zeros(1, 1, 'single'), ...
    'stride', 1, ...
    'pad', 0) ;
net.layers{end} = struct('type', 'sigmoidcrossentropyloss');

load('meanPixel.mat');
imdb.meanPixel = meanPixel;
imdb.scales = scales;

[net,info] = cnn_train_adagrad_ms(net, imdb, @getBatch,...
    opts.train,'errorType','euclideanloss',...
    'conserveMemory', true);

end

function [im, labels] = getBatch(imdb, batch)
% --------------------------------------------------------------------

half_size = imdb.half_size;
meanPixel = imdb.meanPixel;

for ii = 1:numel(batch)

    imagename = imdb.images.name{batch(ii)};
    im_ii = single(imread(imagename));
    labelname = imdb.images.labels{batch(ii)};
    roi = imread(labelname);
    labels_ii = zeros(size(roi, 1), size(roi, 2));
    labels_ii( roi == 50 )  = 0.25;       %shade
    labels_ii( roi == 170 ) = 0.75;       %object boundary
    labels_ii( roi == 255 ) = 1;          %foreground

    % resize the image to half size
    if size(im_ii,1) > 400 || size(im_ii,2) >400
        im_ii = imresize(im_ii, 0.5, 'nearest');
        labels_ii = imresize(labels_ii, 0.5, 'nearest');
    end

    im(:,:,:,ii) = im_ii;
    labels(:,:,1,ii) = labels_ii;
end
end

function imdb = getImdb(imgDir, labelDir)

files = dir([imgDir '/*.jpg']);
label_files = dir([labelDir '/*.png']);
names = {};labels = {};

for ii = 1:numel(files)
    names{end+1} = [imgDir '/' files(ii).name];
    labels{end+1} = [labelDir '/' label_files(ii).name];
end

imdb.images.set = ones(1,numel(names));
imdb.images.name = names ;
imdb.images.labels = labels;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function f = or(f, g)
%|   CLASSICFUN logical OR.
%   F | G performs a logical OR of the CLASSICFUN objects F and G and returns a CLASSICFUN
%   containing elements set to either logical 1 (TRUE) or logical 0 (FALSE).
%   An element of the output CLASSICFUN is set to 1 if either input CLASSICFUN contains a
%   non-zero element at that point, otherwise it is set to 0.  F and G must
%   either be identically zero or have roots in their domains.  If this is not
%   the case, garbage is returned with no warning.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

f.onefun = f.onefun | g.onefun;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%
% Implementation of Bertsekas' auction algorithm [1] for a very fast
% solution of the linear assignment problem, i.e. it is able
% to solve a 1000-by-1000 problem in ~1.5s.
%
% In practice, with respect to running time, the auction algorithm
% outperforms the Kuhn-Munkres (or Hungarian) algorithm significantly.
% The auction algorithm has average run-time complexity O(N^2*log(N)) and
% worst-case complexity O(N^3), whereas the Kuhn-Munkres algorithm has
% average- and worst-case complexity  O(N^3).
%
% Note that only global optima are found for integer-valued benefit matrices.
% For real-valued benefit matrices a scaling of the values needs to be applied
% by multiplication with a large number. This scaling factor depends on the
% desired accuracy, as the global solution is found for the integral part
% of the benefit matrix, whilst there is no guarantee that the fractional part
% of the benefits are properly taken into account. However, in practical cases
% it seems advantageous to not round the resulting benefit matrix and retain
% the fractional parts in the benefit matrix. Also note that larger scalings
% of the benefit matrix increase the run-time, so a problem-specific trade-off
% between runtime and accuracy must be chosen.
%
% Runtime can be slightly decreased by installing the FastSet library,
% available on
% http://www.levmuchnik.net/Content/ProgrammingTips/MatLab/FastSet/FastSet.html
%
% Input:
%               A                       N-by-N benefit matrix (higher
%                                       values indicate a better match)
%               [epsilon]               Initial value of epsilon (optional)
%               [epsilonDecreaseFactor] Decreas factor of epsilon
%                                       (optional)
%
% Output:
%               assignments             Resulting assignments
%               [prices]                Prices used during auctions
%                                       (optional)
%
% Example:      See the function test() below for a usage example.
%               Typically only the benefit matrix A is given as input and the
%               first output argument is relevant. epsilon and
%               epsilonDecreaseFactor can be used to heuristically adapt
%               runtime. The example below can also be used for testing
%               whether a sufficient amount of scaling has been performed.
%               This is done by solving the assignment problem using CVX
%               (available on http://cvxr.com/cvx/download/ ) and comparing
%               its result to the solution with the auction algorithm.
%
%
% [1]	Bertsekas, D.P. 1998. Network Optimization Continuous and Discrete Models.
%
% Implementation by Florian Bernard ( f.bernardpi [at] gmail [dot] com )
%
% Created on 25/07/2014

function [assignments, prices] = ...
    assignmentProblemAuctionAlgorithm(A, epsilon, epsilonDecreaseFactor)

N = size(A,1);
assignments = nan(N,1);
prices = ones(N,1);

% check which unique function we use
if ( exist('fast_unique', 'file') )
    % requires FastSet, available on
    % http://www.levmuchnik.net/Content/ProgrammingTips/MatLab/FastSet/FastSet.html
    uniqueFcn = @fast_unique;
else
    uniqueFcn = @unique;
end

% heuristic for setting epsilon
A = A*(N+1);
if ( ~exist('epsilon', 'var') || isempty(epsilon) )
    maxAbsA = max(abs(A(:)));
    epsilon = maxAbsA/25;
%     if ( sum(abs(A(:)) > 0.0001) > 0.7*N*N ) % non-sparse
%         ...
%     else
%         epsilon = 0.5*((N*maxVal)/5 + N*maxVal); % see page 260 in [1]
%     end
end

if ( ~exist('epsilonDecreaseFactor', 'var') )
    epsilonDecreaseFactor = 0.2;
end

while (epsilon >= 1)
    % the outer loop performs epsilon-scaling in order to speed up execution
    % time. In particular, an updated prices array is computed in each
    % round, which speeds up further rounds with lower values of epsilon.
    assignments = nan(N,1);
    while (any(isnan(assignments)))
    %% forward-auction algorithm
        %% bidding phase
        % find unassigned indices
        unassignedIdx = find(isnan(assignments));
        nUnassigned = numel(unassignedIdx);

        % find best and second best objects
        AijMinusPj = bsxfun(@minus, A(unassignedIdx,:), prices');

        [~,viIdx] = max(AijMinusPj, [], 2);

        for i=1:nUnassigned
            AijMinusPj(i,viIdx(i)) = -inf;
        end
        wi = max(AijMinusPj, [], 2);

        % compute bids
        bids = nan(nUnassigned,1);
        for i=1:nUnassigned
            bids(i) = A(unassignedIdx(i),viIdx(i)) - wi(i) + epsilon;
        end

        %% assignment phase
        objectsThatHaveBeenBiddedFor = uniqueFcn(viIdx);
        for uniqueObjIdx=1:numel(objectsThatHaveBeenBiddedFor)
            currObject = objectsThatHaveBeenBiddedFor(uniqueObjIdx);
            personssWhoGaveBidsForJ = find(viIdx==currObject);

            [prices(currObject), idx] = max(bids(personssWhoGaveBidsForJ));
            personWithHighestBid = unassignedIdx(personssWhoGaveBidsForJ(idx));

            % remove previous assignment and store new assignment (the person
            % with highest bid)
            assignments(assignments==currObject) = nan;
            assignments(personWithHighestBid) = currObject;

        end
    end
    epsilon = epsilon * epsilonDecreaseFactor; % refine epsilon
end
end

function test()
N = 1000;

X = rand(N,3);
reference = rand(N,3);

A = X*reference';
%% Apply Auction algorithm
tic
    % scale A such that the integer version of A has sufficient accuracy
    % if S1 and P are not equal, the scaling factor needs to be increased
    scalingFactor = 10^6;
    Ascaled = A*scalingFactor;

    [assignments] = assignmentProblemAuctionAlgorithm(Ascaled);

    % create permutation matrix from assignments
    linIdx = sub2ind(size(A),assignments',1:size(A,1));
    P = sparse(size(A,1),size(A,2));
    P(linIdx) = 1;
toc

%% Use CVX for checking correctness of result
%  (Useful in particular for non-integer matrices A)
%  Note that this is really slow as CVX needs to create the canonical form
%  of the LP before it can call Gurobi. Creating the canonical form by hand
%  and directly calling Gurobi is much faster.
E = ones(N,1);

tic
    cvx_begin % quiet
        cvx_solver gurobi % uncomment this line if Gurobi is not available
        variable S1(N,N) nonnegative;
        maximize (trace(S1*A));
        S1*E==E;
        E'*S1==E';
    cvx_end
toc

% check if results are equal
if ( sum(sum(full(S1)~=full(P))) == 0 )
    disp('Correct results using Auction algorithm');
else
    error('Incorrect results using Auction algorithm. You probably need to scale A by a larger factor!');
end
end

author: OshriHalimi


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

% Script chap::2::script
% Density Bachelier Model

t = 5;                                     % maturity
x0 = 100;                                  % spot value
r= 0;                                       % riskless rate
q=0;                                        % dividend yield

figname = 'Risk Neutral Density - Bachelier';

x = 95:.1:105;                            % range

sigma_base = 0.2;                           % volatility base scenario

y_base = pg(x, x0, t, r, q, sigma_base);
legendname_base = 'Base';
%% Changing sigma_base
sigma_low = 0.15;
sigma_high = 0.25;

y_low = pg(x, x0, t, r, q, sigma_low);
y_high = pg(x, x0, t, r, q, sigma_high);

legendname_low = 'Changing \sigma low value';
legendname_high = 'Changing \sigma high value';

createfigure_density(x,y_base,y_low,y_high,...
    figname, legendname_base,legendname_low,legendname_high);

x = 95:.1:105;                            % range

sigma_base = 0.2;                           % volatility base scenario

y_base = pg(x, x0, t, r, q, sigma_base);
legendname_base = 'Base';


text:
function [ x, y, z, w ] = ld0014 ( )

%*****************************************************************************80
%
%% LD0014 computes the 14 point Lebedev angular grid.
%
%  Modified:
%
%    14 September 2010
%
%  Author:
%
%    Dmitri Laikov
%
%  Reference:
%
%    Vyacheslav Lebedev, Dmitri Laikov,
%    A quadrature formula for the sphere of the 131st
%    algebraic order of accuracy,
%    Russian Academy of Sciences Doklady Mathematics,
%    Volume 59, Number 3, 1999, pages 477-481.
%
%  Parameters:
%
%    Output, real X(N), Y(N), Z(N), W(N), the coordinates
%    and weights of the points.
%
  n = 0;
  x = zeros(14,1);
  y = zeros(14,1);
  z = zeros(14,1);
  w = zeros(14,1);
  a = 0.0;
  b = 0.0;
  v = 0.6666666666666667E-01;
  [ n, x, y, z, w ] = gen_oh ( 1, n, a, b, v, x, y, z, w );
  v = 0.7500000000000000E-01;
  [ n, x, y, z, w ] = gen_oh ( 3, n, a, b, v, x, y, z, w );

  return
end


text:
function [ a, b ] = p20_lim ( )

%*****************************************************************************80
%
%% P20_LIM returns the integration limits for problem 20.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the limits of integration.
%
  a = -1.0;
  b = 1.0;

  return
end


text:
function [ t_stop, y_stop ] = p26_stop ( neqn )

%*****************************************************************************80
%
%% P26_STOP returns the stopping point for problem p26.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_STOP, Y_STOP(NEQN), the final data.
%
  y_stop = zeros ( neqn, 1 );

  t_stop = 20.0;
  y_stop = [ ...
   -1.294460621213470e+01; ...
   -2.208575158908672e-15 ];

  return
end