text:
function [results_obj, stats, indx] = searchlight(dat, varargin)
% Run searchlight multivariate prediction/classification on an image_vector
% or fmri_data object OR two objects, for cross-prediction.
%
% :Usage:
% ::
%
%    [list outputs here] = function_name(list inputs here, [optional inputs])
%    [results_obj, stats, indx] = searchlight(dat, [optional inputs])
%
%
% :Features:
%  - Runs searchlight with standard, pre-defined algorithms
%  - Custom-entry definition of holdout sets
%  - Can re-use searchlight spheres after initial definition
%  - Custom-entry definition of any spheres/regions of interest
%  - Uses Matlab's parallel processing toolbox (parfor)
%
% Type help image_vector.searchlight to display this help information
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2014  Tor Wager and...
%
%     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/>.
% ..
%
%
% :Inputs:
%
%   **dat:**
%        image_vector or fmri_data object with data
%
%   **dat.Y:**
%        required: true outcomes for each observation (image) in dat
%
% :Optional Inputs:* Keyword followed by input variable:
%
%   **r:**
%        searchlight radius, voxels
%   **dat2:**
%        second dataset, for cross-prediction
%   **indx:**
%        sparse logical matrix. each COLUMN is index of inclusion sets for each region/sphere in searchlight
%        This takes a long time to calculate, but can be saved and
%        re-used for a given mask
%   **holdout_set:**
%        Followed by integer vector of which observations belong to which
%        holdout set, for cross-validation. This is passed into fmri_data.predict.m.  Default is
%        empty.
%
% :Outputs:
%
%   **results_obj:**
%        fmri_data object with results maps
%
%   **stats:**
%        selected statistics for each sphere in searchlight
%
%   **indx:**
%        sparse logical matrix. each COLUMN is index of inclusion sets for each region/sphere in searchlight
%        * this can be re-used for all data with the same mask/structure. *
%
%
% :Examples:
% ::
%
%    % Define a sensible gray-matter mask:
%    dat = fmri_data(which('scalped_avg152T1_graymatter.img'));
%    dat = threshold(dat, [.8 Inf], 'raw-between');
%    dat = trim_mask(dat);
%
%    % Create fake data and holdout indicator index vector
%    dat.dat = randn(dat.volInfo.n_inmask, 30);
%    dat.Y = dat.dat(111111, :)' + .3 * randn(30, 1);
%    holdout_set = ones(6, 1); for i = 2:5, holdout_set = [holdout_set; i*ones(6, 1)]; end
%
%    % Run, and run again with existing indx
%    pool = parpool(12);  % initialize parallel processing (12 cores)
%    [results_obj, stats, indx] = searchlight(dat, 'holdout_set', holdout_set);
%    results_obj = searchlight(dat, 'holdout_set', holdout_set, 'indx', indx);
%
% :See also:
% region.m, fmri_data.predict.m
%
% ..
%    Programmers' notes:
%    List dates and changes here, and author of changes
% ..

% ..
% DEFAULTS AND INPUTS
% ..

r = 5; % For defining regions
indx = [];      % cell; with logical index of inclusion sets for each region/sphere in searchlight

% For prediction/classification
algorithm_name = 'cv_lassopcr';
holdout_set = [];
results_obj = [];

% process variable arguments
for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}
            % functional commands
            case {'r', 'indx', 'algorithm_name', 'dat2', 'holdout_set'}
                str = [varargin{i} ' = varargin{i + 1};'];
                eval(str)
                varargin{i} = [];
                varargin{i + 1} = [];

            case {'cross_predict'}
                algorithm_name = varargin{i};

            otherwise
                % use all other inputs as optional arguments to specific
                % functions, to be interpreted by them
                predfun_inputs{end + 1} = varargin{i};

                if (i+1) <= length(varargin) && ~ischar(varargin{i + 1})
                    predfun_inputs{end + 1} = varargin{i + 1};
                end

        end
    end
end

n = dat.volInfo.n_inmask;

%% Set up indices for spherical searchlight, if not entered previously

if isempty(indx)

    indx = searchlight_sphere_prep(dat, r);

else

    fprintf('Using input indx to define regions/spheres...\n');

end

%% Check for data and return if empty

if ~isa(dat, 'fmri_data') || isempty(dat.Y)

    fprintf('Returning indx only: No data in dat.Y to predict or data is not an fmri_data object');

    return

end

%% Run cross-predict (or other) function

% Get rough time estimate first
predict_time_estimate(dat, indx, algorithm_name, holdout_set);

t = tic;
fprintf('Running prediction in each region...');

% parfor i = 1:n
%
%     indxcell{i} = indx(:, i);
%
% end

output_val = cell(1, n);

parfor i = 1:n

    output_val{i} = predict_wrapper(dat, indx(:, i), algorithm_name, holdout_set);

end

e = toc(t);

[hour, minute, second] = sec2hms(e);
fprintf(1,'Done in %3.0f hours %3.0f min %2.0f sec\n',hour, minute, second);

%% Save results map(s) in object

results_obj = mean(dat);
results_obj.dat = cat(1, output_val.cverr{:});

end % function

% -------------------------------------------------------------------------
% -------------------------------------------------------------------------
%
% Sub-functions
%
% -------------------------------------------------------------------------
% -------------------------------------------------------------------------

% SPHERE PREP
% -------------------------------------------------------------------------

function indx = searchlight_sphere_prep(dat, r)

n = dat.volInfo.n_inmask;
indx = cell(1, n);

t = tic;
fprintf('Setting up seeds...');

parfor i = 1:n
    seed{i} = dat.volInfo.xyzlist(i, :);
end

e = toc(t);
fprintf('Done in %3.2f sec\n', e);

% Set up indices for spherical searchlight
% -------------------------------------------------------------------------
% These could be indices for ROIs, user input, previously saved indices...

% First, a rough time estimate:
% -------------------------------------------------------------------------
fprintf('Searchlight sphere construction can take 20 mins or more! (est: 20 mins with 8 processors/gray matter mask)\n');
fprintf('It can be re-used once created for multiple analyses with the same region definitions\n');
fprintf('Getting a rough time estimate for how long this will take...\n');

n_to_run = min(500, n);
t = tic;
parfor i = 1:n_to_run

    mydist = sum([dat.volInfo.xyzlist(:, 1) - seed{i}(1) dat.volInfo.xyzlist(:, 2) - seed{i}(2) dat.volInfo.xyzlist(:, 3) - seed{i}(3)] .^ 2, 2);
    indx{i} = mydist <= r.^2;

end
e = toc(t);
estim = e * n / n_to_run;

[hour, minute, second] = sec2hms(estim);
fprintf(1,'\nEstimate for whole brain = %3.0f hours %3.0f min %2.0f sec\n',hour, minute, second);

% Second, do it for all voxels/spheres:
% -------------------------------------------------------------------------

t = tic;

fprintf('Constructing spheres for each seed...');

%infdist = inf * ones(size(dat.volInfo.xyzlist(:, 1), 1), 1);

parfor i = 1:n

    mydist = sum([dat.volInfo.xyzlist(:, 1) - seed{i}(1) dat.volInfo.xyzlist(:, 2) - seed{i}(2) dat.volInfo.xyzlist(:, 3) - seed{i}(3)] .^ 2, 2);
    indx{i} = mydist <= r.^2;

end

e = toc(t);
fprintf('Done in %3.2f sec\n', e);

% Make sparse matrix
t = tic;
fprintf('Sparsifying region indices...');

indxcat = cat(2, indx{:});
indxcat = sparse(indxcat);

e = toc(t);
fprintf('Done in %3.2f sec\n', e);

% NOTE: IdenTICAL but twice as slow
% indx1 = cell(1, n);
% tic
% parfor i = 1:1000
%     indx1{i} = logical(iimg_xyz2spheres(seed{i}, dat.volInfo.xyzlist, r));
% end
% toc

end % searchlight_sphere_prep

% PREDICTION
% -------------------------------------------------------------------------

function output_val = predict_wrapper(dat, indx, algorithm_name, holdout_set)
% This function wraps predict.m to use its functionality.
% It could be extended to handle optional inputs, etc.
% Right now it is basic.
% Note: pass in indx{i}, which becomes indx here, for ONE
% region/searchlight sphere

dat.dat = dat.dat(indx, :);
dat.removed_voxels(~indx) = true;

[cverr, stats] = predict(dat, 'algorithm_name', algorithm_name, 'nfolds', holdout_set, 'useparallel', 0, 'verbose', 0);

% Parse output
% Only storing important info and full data weight map - we will probably
% want xval weight maps too at some point. We will also probably want a
% p-value and standard errors for accuracy too.

output_val = struct;
switch algorithm_name
    case 'cv_svm'
        output_val.cverr = stats.cverr;
        output_val.dist_from_hyperplane_xval = stats.dist_from_hyperplane_xval;
        output_val.weight_obj = stats.weight_obj;
        output_val.intercept = stats.other_output{3};

    case 'cv_lassopcr'
        output_val.yfit = stats.yfit;
        output_val.pred_outcome_r = stats.pred_outcome_r;
        output_val.weight_obj = stats.weight_obj;
        output_val.intercept = stats.other_output{3};
end

end % subfunction

% PREDICTION time estimate
% -------------------------------------------------------------------------

function predict_time_estimate(dat, indx, algorithm_name, holdout_set)

t = tic;
fprintf('Getting rough time estimate: Running prediction for up to 500 voxels...');

n = dat.volInfo.n_inmask;
output_val = cell(1, n);

n_to_run = min(500, n);

parfor i = 1:n_to_run

    output_val{i} = predict_wrapper(dat, indx(:, i), algorithm_name, holdout_set);

end

e = toc(t);
fprintf('Done in %3.2f sec\n', e);

estim = e * n / n_to_run;

[hour, minute, second] = sec2hms(estim);
fprintf(1,'Estimate for whole brain = %3.0f hours %3.0f min %2.0f sec\n',hour, minute, second);

end

% Convert time
% -------------------------------------------------------------------------

function [hour, minute, second] = sec2hms(sec)
%SEC2HMS  Convert seconds to hours, minutes and seconds.
%
%   [HOUR, MINUTE, SECOND] = SEC2HMS(SEC) converts the number of seconds in
%   SEC into hours, minutes and seconds.

hour   = fix(sec/3600);      % get number of hours
sec    = sec - 3600*hour;    % remove the hours
minute = fix(sec/60);        % get number of minutes
sec    = sec - 60*minute;    % remove the minutes
second = sec;
end

% % WANI'S CROSS-CLASSIFICATION FUNCTION - to edit and integrate
% % -------------------------------------------------------------------------
%
% function [acc, p, se, stats1 stats2] = cv_svm_cross_subf(data1, data2, cv_assign, dobalanced, balanced_ridge, doscale)
%
% % run cv_svm on train data and test data
%
% if doscale
%     data1 = rescale(data1, 'zscoreimages');
%     data2 = rescale(data2, 'zscoreimages');
% end
%
% % preallocate some variables
% predicted11 = NaN(size(data1.Y));
% predicted12 = NaN(size(data2.Y));
% predicted21 = NaN(size(data1.Y));
% predicted22 = NaN(size(data2.Y));
%
% dist_from_hyper11 = NaN(size(data1.Y));
% dist_from_hyper12 = NaN(size(data2.Y));
% dist_from_hyper21 = NaN(size(data1.Y));
% dist_from_hyper22 = NaN(size(data2.Y));
%
% % get cv assignment
% u = unique(cv_assign);
% nfolds = length(u);
% [trIdx, teIdx] = deal(cell(1, nfolds));
%
% for i = 1:length(u)
%     teIdx{i} = cv_assign == u(i);
%     trIdx{i} = ~teIdx{i};
% end
%
% stats1.Y = data1.Y;
% stats2.Y = data2.Y;
%
% % Use all
% trainobj1 = data(double(data1.dat'),data1.Y);
% trainobj2 = data(double(data2.dat'),data2.Y);
%
% svmobj1 = svm({'optimizer="andre"','C=1','child=kernel'});
% if dobalanced, svmobj1.balanced_ridge = balanced_ridge; end
% [~, svmobj1] = train(svmobj1, trainobj1);
%
% svmobj2 = svm({'optimizer="andre"','C=1','child=kernel'});
% if dobalanced, svmobj2.balanced_ridge = balanced_ridge; end
% [~, svmobj2] = train(svmobj2, trainobj2);
%
% % 1) weights
% stats1.all{1} = get_w(svmobj1)';
% stats2.all{1} = get_w(svmobj2)';
% % 2) intercepts
% stats1.all{2} = svmobj1.b0;
% stats2.all{2} = svmobj2.b0;
%
% stats1.all_descrip = {'1:weights 2:intercept 3:distance from hyperplane for test data1 4: dist for data2'};
% stats2.all_descrip = {'1:weights 2:intercept 3:distance from hyperplane for test data2 4: dist for data1'};
%
% % Cross-val loop starts here
%
% for i = 1:numel(teIdx)
%
%     % Prepare the hot-warm and Rej-friend data in Spider format
%     trainobj1 = data(double(data1.dat'),data1.Y);
%     trainobj2 = data(double(data2.dat'),data2.Y);
%
%     % Select training and test data - trainobj, testobj
%     testobj1 = trainobj1;
%     testobj2 = trainobj2;
%
%     trainobj1.X = trainobj1.X(trIdx{i}, :);
%     trainobj1.Y = trainobj1.Y(trIdx{i}, :);
%
%     trainobj2.X = trainobj2.X(trIdx{i}, :);
%     trainobj2.Y = trainobj2.Y(trIdx{i}, :);
%
%     testobj1.X = testobj1.X(teIdx{i}, :);
%     testobj1.Y = testobj1.Y(teIdx{i}, :);
%
%     testobj2.X = testobj2.X(teIdx{i}, :);
%     testobj2.Y = testobj2.Y(teIdx{i}, :);
%
%
%     %-------------------------------------------------------------------
%     % Train on data1 and test on data1 and data2
%
%     % Set up an SVM object
%     svmobj1 = svm({'optimizer="andre"','C=1','child=kernel'});
%
%     if dobalanced
%         svmobj1.balanced_ridge = balanced_ridge;
%     end
%
%     % train on hw data set
%     [~, svmobj1] = train(svmobj1, trainobj1);
%
%     % Test on both HW and RF test set
%
%     testobj11 = test(svmobj1, testobj1);
%     testobj12 = test(svmobj1, testobj2);
%
%     w = get_w(svmobj1)';
%     b0 = svmobj1.b0;
%
%     dist11 = testobj1.X * w + b0;
%     dist12 = testobj2.X * w + b0;
%
%     if ~isequal((dist11>0)*2-1,testobj11.X)
%         error('something is wrong');
%     elseif ~isequal((dist12>0)*2-1,testobj12.X)
%         error('something is wrong');
%     end
%
%     % Save predictions (cross-val)
%     predicted11(teIdx{i}) = testobj11.X;
%     predicted12(teIdx{i}) = testobj12.X;
%
%
%     %-------------------------------------------------------------------
%     % Train on RF and test on HW and RF
%
%     % Set up an SVM object
%     svmobj2 = svm({'optimizer="andre"','C=1','child=kernel'});
%
%     if dobalanced
%         svmobj2.balanced_ridge = balanced_ridge;
%     end
%     % train on rf data set
%     [~, svmobj2] = train(svmobj2, trainobj2);
%
%     % Test on both HW and RF test set
%
%     testobj21 = test(svmobj2, testobj1);
%     testobj22 = test(svmobj2, testobj2);
%
%     w = get_w(svmobj2)';
%     b0 = svmobj2.b0;
%
%     dist21 = testobj1.X * w + b0;
%     dist22 = testobj2.X * w + b0;
%
%     if ~isequal((dist21>0)*2-1,testobj21.X)
%         error('something is wrong');
%     elseif ~isequal((dist22>0)*2-1,testobj22.X)
%         error('something is wrong');
%     end
%
%     % Save predictions (cross-val)
%     predicted21(teIdx{i}) = testobj21.X;
%     predicted22(teIdx{i}) = testobj22.X;
%
%     % get stats from svmobj
%     % 1) weights
%     stats1.cvoutput{i,1} = get_w(svmobj1)';
%     stats2.cvoutput{i,1} = get_w(svmobj2)';
%     % 2) intercepts
%     stats1.cvoutput{i,2} = svmobj1.b0;
%     stats2.cvoutput{i,2} = svmobj2.b0;
%     % 3) distance from hyperplane for test data 1 and 2
%     stats1.cvoutput{i,3} = testobj1.X * stats1.cvoutput{i,1} + stats1.cvoutput{i,2};
%     stats1.cvoutput{i,4} = testobj2.X * stats1.cvoutput{i,1} + stats1.cvoutput{i,2};
%
%     stats2.cvoutput{i,3} = testobj1.X * stats2.cvoutput{i,1} + stats2.cvoutput{i,2};
%     stats2.cvoutput{i,4} = testobj2.X * stats2.cvoutput{i,1} + stats2.cvoutput{i,2};
%
%     stats1.cvoutput_descrip = {'1:weights 2:intercept 3:distance from hyperplane for test data1 4: dist for data2'};
%     stats2.cvoutput_descrip = {'1:weights 2:intercept 3:distance from hyperplane for test data2 4: dist for data1'};
%
%     dist_from_hyper11(teIdx{i}) = stats1.cvoutput{i,3};
%     dist_from_hyper12(teIdx{i}) = stats1.cvoutput{i,4};
%     dist_from_hyper21(teIdx{i}) = stats2.cvoutput{i,3};
%     dist_from_hyper22(teIdx{i}) = stats2.cvoutput{i,4};
% end
%
% stats1.yfit = predicted11;
% stats2.yfit = predicted22;
% stats1.testfit = predicted12;
% stats2.testfit = predicted21;
%
% stats1.all{3} = dist_from_hyper11;
% stats1.all{4} = dist_from_hyper12;
% stats2.all{3} = dist_from_hyper22;
% stats2.all{4} = dist_from_hyper21;
%
% % RO on H vs. W
% % res_temp = binotest(predicted11(1:(length(predicted11)/2)) ==
% %data1.Y(1:(length(predicted11)/2)), 0.5);
% % test_results.accuracy(1,1) = res_temp.prop;
% % test_results.accuracy_se(1,1) = res_temp.SE;
% % test_results.accuracy_p(1,1) = res_temp.p_val;
% %
% % res_temp = binotest(predicted12(1:(length(predicted12)/2)) ==
% %data2.Y(1:(length(predicted12)/2)), 0.5);
% % test_results.accuracy(1,2) = res_temp.prop;
% % test_results.accuracy_se(1,2) = res_temp.SE;
% % test_results.accuracy_p(1,2) = res_temp.p_val;
% %
% % res_temp = binotest(predicted21(1:(length(predicted21)/2)) ==
% %data1.Y(1:(length(predicted21)/2)), 0.5);
% % test_results.accuracy(1,3) = res_temp.prop;
% % test_results.accuracy_se(1,3) = res_temp.SE;
% % test_results.accuracy_p(1,3) = res_temp.p_val;
% %
% % res_temp = binotest(predicted22(1:(length(predicted22)/2)) ==
% %data2.Y(1:(length(predicted22)/2)), 0.5);
% % test_results.accuracy(1,4) = res_temp.prop;
% % test_results.accuracy_se(1,4) = res_temp.SE;
% % test_results.accuracy_p(1,4) = res_temp.p_val;
%
% subjn = numel(teIdx);
% outcome = [true(subjn,1); false(subjn,1)];
%
% hvsw_dist11 = stats1.all{3}(1:(subjn*2));
% rvsf_dist12 = stats1.all{4}(1:(subjn*2));
% hvsw_dist21 = stats2.all{4}(1:(subjn*2));
% rvsf_dist22 = stats2.all{3}(1:(subjn*2));
%
%
% % ROC = roc_plot(hvsw_dist11, outcome, 'threshold', 0);
% % test_results.accuracy_thresh0(1) = ROC.accuracy;
% % test_results.accuracy_thresh0_se(1) = ROC.accuracy_se;
% % test_results.accuracy_thresh0_p(1) = ROC.accuracy_p;
%
% try
%     ROC = roc_plot(hvsw_dist11, outcome, 'twochoice');
%     acc(1) = ROC.accuracy;
%     se(1) = ROC.accuracy_se;
%     p(1) = ROC.accuracy_p;
% catch err
%     % ROC = roc_plot(hvsw_dist11, outcome, 'twochoice');
%     acc(1) = NaN;
%     se(1) = NaN;
%     p(1) = NaN;
% end
%
% %
% % ROC = roc_plot_wani(rvsf_dist12, outcome, 'threshold', 0);
% % test_results.accuracy_thresh0(2) = ROC.accuracy;
% % test_results.accuracy_thresh0_se(2) = ROC.accuracy_se;
% % test_results.accuracy_thresh0_p(2) = ROC.accuracy_p;
% %
% try
%     ROC = roc_plot(rvsf_dist12, outcome, 'twochoice');
%     acc(2) = ROC.accuracy;
%     se(2) = ROC.accuracy_se;
%     p(2) = ROC.accuracy_p;
% catch err
%     acc(2) = NaN;
%     se(2) = NaN;
%     p(2) = NaN;
% end
% %
% %
% % ROC = roc_plot(hvsw_dist21, outcome, 'threshold', 0);
% % test_results.accuracy_thresh0(3) = ROC.accuracy;
% % test_results.accuracy_thresh0_se(3) = ROC.accuracy_se;
% % test_results.accuracy_thresh0_p(3) = ROC.accuracy_p;
% %
% try
%     ROC = roc_plot(hvsw_dist21, outcome, 'twochoice');
%     acc(3) = ROC.accuracy;
%     se(3) = ROC.accuracy_se;
%     p(3) = ROC.accuracy_p;
% catch err
%     acc(3) = NaN;
%     se(3) = NaN;
%     p(3) = NaN;
% end
%
% % ROC = roc_plot(rvsf_dist22, outcome, 'threshold', 0);
% % test_results.accuracy_thresh0(4) = ROC.accuracy;
% % test_results.accuracy_thresh0_se(4) = ROC.accuracy_se;
% % test_results.accuracy_thresh0_p(4) = ROC.accuracy_p;
% try
%     ROC = roc_plot(rvsf_dist22, outcome, 'twochoice');
%     acc(4) = ROC.accuracy;
%     se(4) = ROC.accuracy_se;
%     p(4) = ROC.accuracy_p;
% catch err
%     acc(4) = NaN;
%     se(4) = NaN;
%     p(4) = NaN;
% end
%
% close;

% end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function pinv = perm_inverse ( n, p )

%*****************************************************************************80
%
%% PERM_INVERSE computes the inverse of a permutation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer N, the number of values being permuted.
%    N must be positive.
%
%    Input, integer P(N), describes the permutation.
%    P(I) is the item which is permuted into the I-th place
%    by the permutation.
%
%    Output, integer PINV(N), the inverse permutation.
%

%
%  Check.
%
  perm_check ( n, p );

  pinv(p(1:n)) = 1:n;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function f=fal(e,a,d)
    if abs(e)<d
        f=e*d^(a-1);
    else f=(abs(e))^a*sign(e);
    end

author: TianfaYao

repo: ADRC

save_path: github-repos/MATLAB/TianfaYao-ADRC


text:
function write_nitf_dessubhdr(obj)
%WRITE_NITF_DESSUBHDR Writes the NITF Data Extension Segment sub-header
%
% It's not expected that the typical user would call this function, rather
% it is called as part of a bigger NITF/SICD writer.
%
% Inputs:
%       obj:       This is the 'hidden' object-specific handle.  It is
%                  similar to 'this' in Java.
%
% References:
%       MIL-STD-2500C, Department of Defense Interface Standard
%                      National Imagery Transmission Format version 2.1
%                      01 May 2006
%
% Written by: Tom Krauss, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

% Always writing version 2.1 NITF

fwriten(obj.FID, 'DE', 2);   % DE
fwriten(obj.FID, 'XML_DATA_CONTENT', 25);  % DESID
fwriten(obj.FID, '01', 2);   % DESVER
obj.write_nitf_security_tags(); % Security tags go here
%fwriten(obj.FID, '0000', 4);   % DESSHL - Complete inclusion of all User-defined Subheader Subfield was not required prior to SICD version 1.0
fwriten(obj.FID, '0773', 4);   % DESSHL
fwriten(obj.FID, '99999', 5); % DESCRC - CRC not computed
fwriten(obj.FID, 'XML', 8); % DESSHFT
if isfield(obj.sicdmeta,'ImageCreation') && isfield(obj.sicdmeta.ImageCreation, 'DateTime')
    fwriten(obj.FID, datestr(obj.sicdmeta.ImageCreation.DateTime,'yyyy-mm-ddTHH:MM:SSZ'), 20); % DESSHDT
else
    fwriten(obj.FID, '', 20); % DESSHDT
end
fwriten(obj.FID, '', 40); % DESSHRP
fwriten(obj.FID, 'SICD Volume 1 Design & Implementation Description Document', 60); % DESSHSI
fwriten(obj.FID, '1.0', 10); % DESSHSV
fwriten(obj.FID, '2011-09-28T00:00:00Z', 20); % DESSHSD
fwriten(obj.FID, 'urn:SICD:1.0.0', 120); % DESSHTN
if isfield(obj.sicdmeta,'GeoData') && isfield(obj.sicdmeta.GeoData,'ImageCorners') && ...
    isfield(obj.sicdmeta.GeoData.ImageCorners,'ICP')
    DESSHLPG = sprintf('%+012.8f%+013.8f%+012.8f%+013.8f%+012.8f%+013.8f%+012.8f%+013.8f',...
        obj.sicdmeta.GeoData.ImageCorners.ICP.FRFC.Lat,...
        obj.sicdmeta.GeoData.ImageCorners.ICP.FRFC.Lon,...
        obj.sicdmeta.GeoData.ImageCorners.ICP.FRLC.Lat,...
        obj.sicdmeta.GeoData.ImageCorners.ICP.FRLC.Lon,...
        obj.sicdmeta.GeoData.ImageCorners.ICP.LRLC.Lat,...
        obj.sicdmeta.GeoData.ImageCorners.ICP.LRLC.Lon,...
        obj.sicdmeta.GeoData.ImageCorners.ICP.LRFC.Lat,...
        obj.sicdmeta.GeoData.ImageCorners.ICP.LRFC.Lon,...
        obj.sicdmeta.GeoData.ImageCorners.ICP.FRFC.Lat,...
        obj.sicdmeta.GeoData.ImageCorners.ICP.FRFC.Lon);
else
    DESSHLPG = '';
end
fwriten(obj.FID, DESSHLPG, 125); % DESSHLPG
fwriten(obj.FID, '', 25); % DESSHLPT
fwriten(obj.FID, '', 20); % DESSHLI
fwriten(obj.FID, '', 120); % DESSHLIN
fwriten(obj.FID, '', 200); % DESSHABS

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
% -*- INTERNAL UNDOCUMENTED FUNCTION -*-
%
% Copyright (C) 2010 Carlo de Falco
% Copyright (C) 2015 Rafael Vazquez
%
% 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 Octave; see the file COPYING.  If not, see
% <http://www.gnu.org/licenses/>.

function [Jinv, geo_inv_der2] = geopdes_inv_der2__ (geo_map_jac, geo_map_der2)

vsize = size (geo_map_jac);
vsize(end+1:4) = 1;

rdim = vsize(1);
ndim = vsize(2);

Jinv = geopdes_inv__ (geo_map_jac);
Jinv = reshape (Jinv, [ndim, rdim, vsize(3:end)]);
Jsize = size (Jinv);
Jsize(end+1:4) = 1;

Jinv1 = reshape (Jinv, [1, Jsize]);
Jinv2 = permute (Jinv1, [1 3 2 4:numel(Jsize)+1]);
Jinv3 = permute (Jinv1, [3 2 1 4:numel(Jsize)+1]);

geo_inv_der2 = zeros ([ndim, rdim, rdim, vsize(3:end)]);

for alp = 1:ndim
  for idim = 1:rdim
    for jdim = 1:rdim
      aux1 = sum (bsxfun (@times, geo_map_der2, Jinv1(:,:,idim,:,:)), 2);
      aux2 = sum (bsxfun (@times, aux1, Jinv2(:,jdim,:,:,:)), 3);
      aux3 = sum (bsxfun (@times, aux2, Jinv3(:,alp,:,:,:)), 1);

      geo_inv_der2(alp,idim,jdim,:,:) = -reshape (aux3, vsize(3:end));
    end
  end
end

end

%!test
%!
%! %% x = u^2 + u + 1
%! %% y = v^2 + v + 1
%!
%! xufun = @(u, v) (2*u + 1);
%! xvfun = @(u, v) (0*u);
%! yufun = @(u, v) (0*u);
%! yvfun = @(u, v) (2*v + 1);
%!
%! xuufun = @(u, v) (2+0*u);
%! xuvfun = @(u, v) (0*u);
%! xvvfun = @(u, v) (0*u);
%!
%! yuufun = @(u, v) (0*u);
%! yuvfun = @(u, v) (0*u);
%! yvvfun = @(u, v) (2+0*v);
%!
%! uxfun = @(u, v) (1./(2*u+1));
%! vyfun = @(u, v) (1./(2*v+1));
%! uxxfun = @(u, v) (-2./(2*u+1).^3);
%! vyyfun = @(u, v) (-2./(2*v+1).^3);
%!
%! u = linspace (0,1, 10);
%! v = linspace (0,1, 10);
%! %u = rand (1, 10);
%! %v = rand (1, 10);
%! geo_map_jac(1,1,:) = xufun(u,v);
%! geo_map_jac(1,2,:) = xvfun(u,v);
%! geo_map_jac(2,1,:) = yufun(u,v);
%! geo_map_jac(2,2,:) = yvfun(u,v);
%!
%! geo_map_der2(1,1,1,:) = xuufun(u,v);
%! geo_map_der2(1,1,2,:) = xuvfun(u,v);
%! geo_map_der2(1,2,1,:) = xuvfun(u,v);
%! geo_map_der2(1,2,2,:) = xvvfun(u,v);
%! geo_map_der2(2,1,1,:) = yuufun(u,v);
%! geo_map_der2(2,1,2,:) = yuvfun(u,v);
%! geo_map_der2(2,2,1,:) = yuvfun(u,v);
%! geo_map_der2(2,2,2,:) = yvvfun(u,v);
%!
%! [Jinv, geo_inv_der2] = geopdes_inv_der2__ (geo_map_jac, geo_map_der2);
%! assert (uxfun (u, v), reshape (Jinv(1,1,:), 1, numel(u)), 1e-14)
%! assert (vyfun (u, v), reshape (Jinv(2,2,:), 1, numel(u)), 1e-14)
%! assert (uxxfun (u, v), reshape (geo_inv_der2(1,1,1,:), 1, numel(u)), 1e-14)
%! assert (vyyfun (u, v), reshape (geo_inv_der2(2,2,2,:), 1, numel(u)), 1e-14)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function write_cell(fid,cell)
% write cell string to text file

for i = 1:length(cell)
  fprintf(fid,'%s\n',cell{i});
end

author: mtex-toolbox

repo: mtex

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


text:
function list = GenesReactions(genelist,rxns)
% Returns genes in model and the associative reactions.
%
% USAGE:
%
%    list = GenesReactions(genelist, rxns)
%
% INPUTS:
%    genelist:    list of genes
%    rxns:        reactions
%
% OUTPUT:
%    list:        genes in model and the associative reactions
%
% .. Author: - Stefan G. Thorleifsson 2011
%
% .. rBioNet is published under GNU GENERAL PUBLIC LICENSE 3.0+
% .. Thorleifsson, S. G., Thiele, I., rBioNet: A COBRA toolbox extension for
% .. reconstructing high-quality biochemical networks, Bioinformatics, Accepted.
% .. rbionet@systemsbiology.is

GPRs = rxns(:,6);
Sg = size(genelist);
list = cell(Sg(1),2);
list(:,1) = genelist;
S = size(rxns);

for i = 1:S(1)
    genes = GPR2Genes({GPRs{i}});
    if ~isempty(genes)
        SG = size(genes);
        for k = 1:SG(2)
            line = find(strcmp(genes{k},genelist));
            if isempty(line)
                msgbox(['Gene ' genes{k} ' is not in genelist....' rxns{i,2}],...
                    'My bad....the programmer.','help');
            else
                if isempty(list{line,2})
                    list{line,2} = rxns{i,2};
                else
                    rxn = list{line,2};
                    list{line,2} = [rxn ', ' rxns{i,2}];
                end
            end
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
close all;
clearvars;
clc;
rng default;
M = 8;
N = 16;
p = randperm(N);
% select some rows randomly
row_pics = sort(p(1:M));
% make sure that the first row is there
row_pics(1) = 1;
col_perm = randperm(N);
%col_perm = 1:N;
Dict = spx.dict.PartialDFT(row_pics, col_perm);
x = ones(N, 1)
y  = Dict.apply(x)
z = Dict.adjoint(y)

x(1) = x(1) + 2
y  = Dict.apply(x)
z = Dict.adjoint(y)

DictMtx = double(Dict);
B = dftmtx(N) / sqrt(N);
row_pics
DictMtx(:, 1:4)
B(:, 1:4)

author: indigits

repo: sparse-plex


text:
% 计算差商
function [p, q] = d_d(x,y)
%
% 参数说明：
% 输入参数： x 为节点，y 为函数值（注意：向量 x 与向量 y 的长度必须一致）
% 输出参数：
%    p 包含所有差商，即差商表中的所有值
%    q 为差商表中对角线的值

m = length(x);
x = x(:);
p = zeros(m, m+1);
p(:,1) = x;
p(:,2) = y(:);
for k = 3 : m+1
    p(k-1:m, k) = diff(p(k-2:m, k-1)) ./ ( x(k-1:m) - x(1:m+2-k) );
end
q = diag(p(1:m,2:m+1));

repo: NumericalAnalysis


text:
function model=combineMetaCycKEGGModels(metacycModel,keggModel)
% combineMetaCycKEGGModels
%	Combine MetaCyc and KEGG draft models into one model structure.
%
%	Input:
%	metacycModel	the reconstructed model from MetaCyc
%	keggModel       the reconstructed model from KEGG
%
%	Output:
%	model           a model structure generated by integrating information
%                   from draft models reconstructed using MetaCyc and KEGG
%                   databases. The 'rxnFrom/metFrom/geneFrom' fields are
%                   included to indicate the source.
%
%	Usage: model=combineMetaCycKEGGModels(metacycModel,keggModel)

%Just return the model
if nargin<2
    disp('Missing input model');
    return;
end

%Add MetaCyc model as template
model=metacycModel;
model.id='COMBINED';
model.name='Combined model from MetaCyc and KEGG draft models';

%Use MetaCyc model as template
model.rxnFrom=cell(numel(model.rxns),1);
model.rxnFrom(:)={'MetaCyc'};
model.metFrom=cell(numel(model.mets),1);
model.metFrom(:)={'MetaCyc'};
model.geneFrom=cell(numel(model.genes),1);
model.geneFrom(:)={'MetaCyc'};
model.grRulesKEGG=cell(numel(model.rxns),1);
model.grRulesKEGG(:)={''};

%Remove the S and rxnGeneMat matrix
if isfield(model,'S')
    model=rmfield(model,'S');
end
if isfield(model,'rxnGeneMat')
    model=rmfield(model,'rxnGeneMat');
end

%Arrange the genes related fields

%Identify the shared genes between KEGG and MetaCyc
sharedGenes=intersect(keggModel.genes,model.genes);
[~, b]=ismember(sharedGenes,model.genes);
model.geneFrom(b)={'Both'};

%Add unique genes from KEGG and update geneFrom
C=setdiff(keggModel.genes,model.genes);
model.genes=[model.genes;C];
geneFrom=cell(numel(C),1);
geneFrom(:)={'KEGG'};
model.geneFrom=[model.geneFrom;geneFrom];

%Prepare for matching grRules
if isfield(keggModel,'grRules')
    keggModel.grRules=strrep(keggModel.grRules,'(','');
    keggModel.grRules=strrep(keggModel.grRules,')','');
end

%Replace rxns in KEGG model with the corresponding ones in MetaCyc by
%updating values rxns field to MetaCyc version whenever possible

%Collect the ones found in MetaCyc as mappedRxns
mappedRxns=[];
rxnsToMove=[];
grRulesToMove={}; %Storing the grRules to be moved to MetaCyc model

%Read in dbLinks, note that linkMetaCycKEGGRxns need to be run in advance
load('metaCycRxns.mat');
num=0;
numToMove=0;
%Loop through the KEGG model and replace the rxn id from KEGG to MetaCyc
%version
for i=1:numel(keggModel.rxns)
    [~, b]=ismember(rxnLinks.kegg,keggModel.rxns{i});
    I=find(b);

    if numel(I)==1
        %Find out the corresponding MetaCyc reactions
        num=num+1;
        %Record the mapped rxn index in KEGG model and remove them later
        mappedRxns=[mappedRxns,i];
        %keggModel.rxns{i}=rxnLinks.metacyc{I(1)};

        [c, d]=ismember(rxnLinks.metacyc{I(1)},model.rxns);
        if c
            model.rxnFrom{d}='Both';
            %Combine the grRule info Check if KEGG grRules equals or a
            %subset of MetaCyc grRules If not, save the KEGG grRules to a
            %new field grRulesKEGG for manual curation
            k=strfind(model.grRules{d},keggModel.grRules{i});
            if isempty(k)
                if isempty(model.grRulesKEGG{d})
                    model.grRulesKEGG{d}=keggModel.grRules{i};
                else
                    model.grRulesKEGG{d}=strcat(model.grRulesKEGG{d},{' or '},keggModel.grRules{i});
                end
            end
        else
            %For the KEGG rxns have MetaCyc version but not in MetaCyc
            %draft model
            [~, y]=ismember(rxnLinks.metacyc{I(1)},metaCycRxns.rxns);

            %Check if this reaction is repitetive and save the grRules
            [Repeat, Index]=ismember(y,rxnsToMove);
            if ~Repeat
                rxnsToMove=[rxnsToMove;y];  %Record rxns to be moved from KEGG to MetaCyc model
                if isempty(keggModel.grRules{i})
                    %Rxns may have empty grRules (e.g. spontaneous) that
                    %are identified and added back here
                    grRulesToMove=[grRulesToMove;{''}];
                else
                    grRulesToMove=[grRulesToMove;keggModel.grRules{i}];
                end
                numToMove=numToMove+1;
            else
                %If this reaction has been recorded, append the grRules if
                %it is different
                k=strfind(grRulesToMove{Index},keggModel.grRules{i});
                if isempty(k)
                    grRulesToMove{Index}=strcat(grRulesToMove{Index},{' or '},keggModel.grRules{i});
                end
            end
        end
    else
        %Here are the case of one KEGG rxn id link to several MetaCyc ids
        %for j=2:numel(I) Ignore this issue for now, and resolve it later,
        %because it never happened This case can be solved by adding above
        %section into a loop below disp(I(j)); end
    end
end
fprintf('NOTE: A total of %d reactions in the KEGG model were mapped to MetaCyc.\n',num);
fprintf('NOTE: %d reactions already in MetaCyc model, %d will be combined.\n',num-numToMove,numToMove);

%Append mapped MetaCyc rxns in KEGG model that are absent from MetaCyc
%model
rxnFrom=cell(numel(rxnsToMove),1);
rxnFrom(:)={'KEGG'};
model.rxnFrom=[model.rxnFrom;rxnFrom];
model.rxns=[model.rxns;metaCycRxns.rxns(rxnsToMove)];

rxnNames=metaCycRxns.rxnNames(rxnsToMove);
%Replace rxnNames with those from metaCycEnzymes
load('metaCycEnzymes.mat');
for i=1:numel(rxnsToMove)
    [a, b]=ismember(metaCycRxns.rxns{rxnsToMove(i)},metaCycEnzymes.rxns);
    if a
        rxnNames{i}=metaCycEnzymes.rxnNames{b};
    end
end
model.rxnNames=[model.rxnNames;rxnNames];

model.eccodes=[model.eccodes;metaCycRxns.eccodes(rxnsToMove)];
model.subSystems=[model.subSystems;metaCycRxns.subSystems(rxnsToMove)];
model.rxnMiriams=[model.rxnMiriams;metaCycRxns.rxnMiriams(rxnsToMove)];
model.rxnReferences=[model.rxnReferences;metaCycRxns.rxnReferences(rxnsToMove)];
model.rev=[model.rev;metaCycRxns.rev(rxnsToMove)];
model.equations=[model.equations;metaCycRxns.equations(rxnsToMove)];
model.lb=[model.lb;metaCycRxns.lb(rxnsToMove)];
model.ub=[model.ub;metaCycRxns.ub(rxnsToMove)];
model.c=[model.c;metaCycRxns.c(rxnsToMove)];
model.grRules=[model.grRules;grRulesToMove];
%model.grRules=[model.grRules;keggModel.grRules(rxnsToMove(:,2))];
%model.grRulesKEGG=[model.grRulesKEGG;keggModel.grRules(rxnsToMove(:,2))];
rxnFields=cell(numel(rxnsToMove),1);
rxnFields(:)={''};
%model.rxnNotes=[model.rxnNotes;rxnFields];
model.grRulesKEGG=[model.grRulesKEGG;rxnFields];

%Remove all mapped reactions from KEGG model
pureKeggModel=removeReactions(keggModel,mappedRxns,true,true);
fprintf('NOTE: A global KEGG model with %d reactions and %d metabolites was obtained.\n',numel(pureKeggModel.rxns),numel(pureKeggModel.mets));

%Replace mets information in KEGG model with the corresponding ones in
%MetaCyc This includes updating all metabolite-related fields, except S
%matrix.
load('metaCycMets.mat');

if ~isfield(pureKeggModel,'metCharges')
    pureKeggModel.metCharges=NaN(numel(pureKeggModel.mets),1);
end

num=0;
for i=1:numel(pureKeggModel.mets)
    [a, b]=ismember(pureKeggModel.mets{i},metaCycMets.keggid);
    if a
        %For now only replace the met ids
        num=num+1;
        pureKeggModel.mets{i}=metaCycMets.mets{b};
        %Record the replaced mets into a vector for later
    end
end
fprintf('NOTE: A total of %d metabolites from the global KEGG model will be re-mapped to MetaCyc.\n',num);

%Generate equations for pureKeggModel
equationStrings=constructEquations(pureKeggModel,'',false,false,false,true);

%Add the pure KEGG reactions to MetaCyc model
rxnFrom=cell(numel(pureKeggModel.rxns),1);
rxnFrom(:)={'KEGG'};
model.rxnFrom=[model.rxnFrom;rxnFrom];
model.rxns=[model.rxns;pureKeggModel.rxns];
model.equations=[model.equations;equationStrings];
model.rxnNames=[model.rxnNames;pureKeggModel.rxnNames];
model.eccodes=[model.eccodes;pureKeggModel.eccodes];
model.subSystems=[model.subSystems;pureKeggModel.subSystems];
model.rxnMiriams=[model.rxnMiriams;pureKeggModel.rxnMiriams];
model.rev=[model.rev;pureKeggModel.rev];
model.lb=[model.lb;pureKeggModel.lb];
model.ub=[model.ub;pureKeggModel.ub];
model.c=[model.c;pureKeggModel.c];
model.grRules=[model.grRules;pureKeggModel.grRules];
rxnFields=cell(numel(pureKeggModel.rxns),1);
rxnFields(:)={''};
model.grRulesKEGG=[model.grRulesKEGG;rxnFields];
model.rxnReferences=[model.rxnReferences;rxnFields];
%model.rxnNotes=[model.rxnNotes;rxnFields];

%It could also be that the metabolite and reaction names are empty for some
%reasons. In that case, use the ID instead
%I=cellfun(@isempty,model.metNames); model.metNames(I)=model.mets(I);
I=cellfun(@isempty,model.rxnNames);
model.rxnNames(I)=model.rxns(I);

%Generate S matrix and mets
[S, mets, ~]=constructS(model.equations);
model.S=S;
model.mets=mets;

%Remap metabolite information
if isfield(model,'metFrom')
    model=rmfield(model,'metFrom');
end
model.metFrom=cell(numel(model.mets),1);
model.metFrom(:)={''};

if isfield(model,'metNames')
    model=rmfield(model,'metNames');
end
model.metNames=cell(numel(model.mets),1);
model.metNames(:)={''};

if isfield(model,'metFormulas')
    model=rmfield(model,'metFormulas');
end
model.metFormulas=cell(numel(model.mets),1);
model.metFormulas(:)={''};

if ~isfield(model,'metCharges')
    model=rmfield(model,'metCharges');
end
model.metCharges=NaN(numel(model.mets),1);

if ~isfield(model,'inchis')
    model=rmfield(model,'inchis');
end
model.inchis=cell(numel(model.mets),1);
model.inchis(:)={''};

if ~isfield(model,'metMiriams')
    model=rmfield(model,'metMiriams');
end
model.metMiriams=cell(numel(model.mets),1);
model.metMiriams(:)={''};

keggMets=getMetsFromKEGG();
for i=1:numel(model.mets)
    [a, b]=ismember(model.mets{i},metaCycMets.mets);
    if a
        %Replace the met related info
        model.mets{i}=metaCycMets.mets{b};
        model.metNames{i}=metaCycMets.metNames{b};
        model.metFormulas{i}=metaCycMets.metFormulas{b};
        model.metMiriams{i}=metaCycMets.metMiriams{b};
        model.inchis{i}=metaCycMets.inchis{b};
        model.metCharges(i,1)=metaCycMets.metCharges(b,1);
        num=num+1;
    else
        [c, d]=ismember(model.mets{i},keggMets.mets);
        if c
            model.mets{i}=keggMets.mets{d};
            model.metNames{i}=keggMets.metNames{d};
            model.metFormulas{i}=keggMets.metFormulas{d};
            model.metMiriams{i}=keggMets.metMiriams{d};
            model.inchis{i}=keggMets.inchis{d};
            num=num+1;
        end
    end

end

%It could also be that the metabolite and reaction names are empty for some
%reasons. In that case, use the ID instead
I=cellfun(@isempty,model.metNames);
model.metNames(I)=model.mets(I);

%Add confidenceScores to model in a rough way
if isfield(model,'rxnConfidenceScores')
    model=rmfield(model,'rxnConfidenceScores');
end
model.rxnConfidenceScores=NaN(numel(model.rxns),1);
model.rxnConfidenceScores(:)=2;

%Put all metabolites in one compartment called 's' (for system). This is
%done just to be more compatible with the rest of the code
model.comps={'s'};
model.compNames={'System'};
model.metComps=ones(numel(model.mets),1);
model.b=zeros(numel(model.mets),1);

%Create the new rxnGene matrix, get back to this later
model.rxnGeneMat=sparse(numel(model.rxns),numel(model.genes));

end

author: SysBioChalmers

repo: RAVEN


text:
function computeAssociations_Textures_Old(pathUnivariate,pathFeatures,cohorts,outcomes,timeToEvent,featType,textType,textName,nPatient,scale_mat,algo_cell,Ng_mat)

% THIS FUNCTION WAS MODIFIED ON SEPTEMBER 29 TO COMPUTE FUSED TEXTURE
% RESULTS. PREVIOUS VERSION (WITHOUT FUSED) IS WITH COMMENTED LINES.

startpath = pwd;
cd(pathUnivariate)
fid = fopen('Texture_SignificanceProportion.txt','w');

nameOutcomes = fieldnames(outcomes.(cohorts{1})); nOutcomes = numel(nameOutcomes); nCohort = numel(cohorts);
nFeatType = numel(featType); nTextType = numel(textType);
nText = 0;
for t = 1:nTextType
    nText = nText + numel(textName{t});
end

% featureNames = {'Best PET texture - GTVp'; 'Best PET texture - GTVtot'; 'Best CT texture - GTVp'; 'Best CT texture - GTVtot'}; nFeature = numel(featureNames);
featureNames = {'Best PET texture - GTVp'; 'Best PET texture - GTVtot'; 'Best CT texture - GTVp'; 'Best CT texture - GTVtot'; 'Best PET/CT texture - GTVp'; 'Best PET/CT texture - GTVtot'}; nFeature = numel(featureNames);
rsMat_Binary = zeros(nFeature,nOutcomes); pMat_Binary = zeros(nFeature,nOutcomes); pCell_Binary = cell(nFeature,nOutcomes); stringCell_Binary = cell(nFeature,nOutcomes); nameTextCell_Binary = cell(nFeature,nOutcomes);
rsMat_Time = zeros(nFeature,nOutcomes); pMat_Time = zeros(nFeature,nOutcomes); pCell_Time = cell(nFeature,nOutcomes); stringCell_Time = cell(nFeature,nOutcomes); nameTextCell_Time = cell(nFeature,nOutcomes);
% scans = {'PT','CT'}; nScans = numel(scans); nParams = [numel(scale_mat)*numel(algo_cell)*numel(Ng_mat),numel(scale_mat)*numel(algo_cell)*numel(Ng_mat)]; paramSizes = {[numel(scale_mat),numel(algo_cell),numel(Ng_mat)],[numel(scale_mat),numel(algo_cell),numel(Ng_mat)]};
scans = {'PT','CT','PTCT'}; nScans = numel(scans); nParams = [numel(scale_mat)*numel(algo_cell)*numel(Ng_mat),numel(scale_mat)*numel(algo_cell)*numel(Ng_mat),numel(CTweight_mat)*numel(scale_mat)*numel(algo_cell)*numel(Ng_mat)]; paramSizes = {[numel(scale_mat),numel(algo_cell),numel(Ng_mat)],[numel(scale_mat),numel(algo_cell),numel(Ng_mat)],[numel(CTweight_mat),numel(scale_mat),numel(algo_cell),numel(Ng_mat)]};
nPat = 0;
for c = 1:nCohort
    nPat = nPat + nPatient.(cohorts{c});
end
i = 0;
for scan = 1:nScans
    nParam = nParams(scan);
    paramSize = paramSizes{scan};
    textParamName = cell(1,nText*nParam); count = 0;
    for t = 1:numel(textType)
        for tStr = 1:numel(textName{t})
            for p = 1:nParam
                count = count + 1;
                if numel(paramSize) == 3
                    [aa,bb,cc] = ind2sub(paramSize,p);
                    textParamName{count} = [textType{t},'/',textName{t}{tStr},' -- ','Scale=',num2str(scale_mat(aa)),',Quant.algo=',algo_cell{bb},',Ng=',num2str(Ng_mat(cc))];
                elseif numel(paramSize) == 4 % Remove this if fused is not used.
                    [aa,bb,cc,dd] = ind2sub(paramSize,p);
                    textParamName{count} = [textType{t},'/',textName{t}{tStr},' -- ','CTweight=',num2str(CTweight_mat(aa),'%.2f'),',Scale=',num2str(scale_mat(bb)),',Quant.algo=',algo_cell{cc},',Ng=',num2str(Ng_mat(dd))];
                end
            end
        end
    end
    for type = 1:nFeatType
        i = i + 1; textMat = zeros(nPat,nText*nParam); countPat = 0;
        outcomeMat = zeros(nPat,nOutcomes);
        timeMat = zeros(nPat,nOutcomes);
        cd(pathFeatures)
        for c = 1:nCohort
            cohort = cohorts{c};
            text = load(['text_',cohort,'_',scans{scan},'_',featType{type}]); text = struct2cell(text); text = text{1};
            patientVect = (countPat+1):(countPat+nPatient.(cohort)); countPat = countPat + nPatient.(cohort);
            count = 0;
            for t = 1:numel(textType)
                for tStr = 1:numel(textName{t})
                    for p = 1:nParam
                        count = count + 1;
                        if numel(paramSize) == 3
                            [aa,bb,cc] = ind2sub(paramSize,p);
                            textMat(patientVect,count) = text{aa,bb,cc}.(textType{t}).(textName{t}{tStr}).Data;
                        elseif numel(paramSize) == 4 % Remove this if fused is not used.
                            [aa,bb,cc,dd] = ind2sub(paramSize,p);
                            textMat(patientVect,count) = text{aa,bb,cc,dd}.(textType{t}).(textName{t}{tStr}).Data;
                        end
                    end
                end
            end
            for j = 1:nOutcomes
                outcomeMat(patientVect,j) = outcomes.(cohort).(nameOutcomes{j});
                timeMat(patientVect,j) = timeToEvent.(cohort).(nameOutcomes{j});
            end
        end
        cd(pathUnivariate)
        for j = 1:nOutcomes
            outcome = outcomeMat(:,j);
            time = timeMat(:,j);
            [rsTemp,pTemp] = corr(textMat,outcome,'type','Spearman','rows','pairwise'); temp = abs(rsTemp);
            [significance] = benjamini_hochberg(pTemp,0.10); proportion = sum(significance)/numel(significance);
            if proportion > 0
                starB = '*';
            else
                starB = '';
            end
            fprintf(fid,['Proportion of significant textures: ',scans{scan},', ',featType{type},', ',nameOutcomes{j},', Binary: ',num2str(proportion),'\n']);
            [~,indMax] = max(temp);
            rsMat_Binary(i,j) = rsTemp(indMax); pMat_Binary(i,j) = pTemp(indMax); nameTextCell_Binary{i,j} = textParamName{indMax};
            pCell_Binary{i,j} = pTemp;
            [rsTemp,pTemp] = corr(textMat,time,'type','Spearman','rows','pairwise'); temp = abs(rsTemp);
            [significance] = benjamini_hochberg(pTemp,0.10); proportion = sum(significance)/numel(significance);
            if proportion > 0
                starT = '*';
            else
                starT = '';
            end
            fprintf(fid,['Proportion of significant textures: ',scans{scan},', ',featType{type},', ',nameOutcomes{j},', Time: ',num2str(proportion),'\n']);
            [~,indMax] = max(temp);
            rsMat_Time(i,j) = rsTemp(indMax); pMat_Time(i,j) = pTemp(indMax); nameTextCell_Time{i,j} = textParamName{indMax};
            pCell_Time{i,j} = pTemp;
            if pMat_Binary(i,j) < 0.01
                stringCell_Binary{i,j} = [starB,'rs = ',num2str(rsMat_Binary(i,j),'%.2f'),', p = ',num2str(pMat_Binary(i,j),'%.2i'),starB];
            else
                stringCell_Binary{i,j} = [starB,'rs = ',num2str(rsMat_Binary(i,j),'%.2f'),', p = ',num2str(pMat_Binary(i,j),'%.2f'),starB];
            end
            if pMat_Time(i,j) < 0.01
                stringCell_Time{i,j} = [starT,'rs = ',num2str(rsMat_Time(i,j),'%.2f'),', p = ',num2str(pMat_Time(i,j),'%.2i'),starT];
            else
                stringCell_Time{i,j} = [starT,'rs = ',num2str(rsMat_Time(i,j),'%.2f'),', p = ',num2str(pMat_Time(i,j),'%.2f'),starT];
            end
        end
    end
end
cd(pathUnivariate), fclose(fid);
Locoregional = stringCell_Binary(:,1); Distant = stringCell_Binary(:,2); Death = stringCell_Binary(:,3);
results.tableCorr = table(Locoregional,Distant,Death,'RowNames',featureNames);
Locoregional = nameTextCell_Binary(:,1); Distant = nameTextCell_Binary(:,2); Death = nameTextCell_Binary(:,3);
results.tableTextName = table(Locoregional,Distant,Death,'RowNames',featureNames);
results.rsMat = rsMat_Binary; results.pMat = pMat_Binary; results.pCell = pCell_Binary;
save('texturesBest_UniV_Binary','results'), clear results
Locoregional = stringCell_Time(:,1); Distant = stringCell_Time(:,2); Death = stringCell_Time(:,3);
results.tableCorr = table(Locoregional,Distant,Death,'RowNames',featureNames);
Locoregional = nameTextCell_Time(:,1); Distant = nameTextCell_Time(:,2); Death = nameTextCell_Time(:,3);
results.tableTextName = table(Locoregional,Distant,Death,'RowNames',featureNames);
results.rsMat = rsMat_Time; results.pMat = pMat_Time; results.pCell = pCell_Time;
save('texturesBest_UniV_Time','results'), clear results

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
addpath('./lib')
run('./SE-Sync/manopt/importmanopt.m')
addpath('./SE-Sync/lib')

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function [S] = co2_sw(x)
% co2_sw.m
% Downloaded from the personal page of Professor Sigurd Skogestad
% http://www.nt.ntnu.no/users/skoge/diplom/prosjekt03/wold/MatLab/
[S] = Helmholtz_anonymous_23829804(x);

function [S] = Helmholtz_anonymous_23829804(x)

[S] = StandardState_anonymous_23817360(x);

function [S] = StandardState_anonymous_23817360(x)

[S]      = EquationOfState_anonymous_23031408(x);
i        = [3];
[S_1]    = MuT_cp_dippr_23809380(x(1));
S.g(1)   = S.g(1) + S_1.dmudT'*x(i);
S.g(i)   = S.g(i) + S_1.mu;
S.H(1,1) = S.H(1,1) + S_1.d2mudTdT'*x(i);
S.H(i,1) = S.H(i,1) + S_1.dmudT;
S.H(1,i) = S.H(i,1)';

function [S] = EquationOfState_anonymous_23031408(x)

[S]   = EquationOfState_anonymous_22908308(x);
[S_1] = ModTVN_ideal_idealgas_23028984(x);
S.g   = S.g + S_1.g;
S.H   = S.H + S_1.H;

function [S] = EquationOfState_anonymous_22908308(x)

[S] = ModTVN_sw_0_5_7_0_23818404(x);

function [S] = ModTVN_sw_0_5_7_0_23818404(x)

R               = 8.314511984;
T_c             = 304.1282;
rho_c           = 10624.90627;
tau             = T_c/x(1);
delta           = x(3)/x(2)/rho_c;
N               = x(3);
a_1             = [0.89875108;-2.1281985;-0.06819032;0.076355306;0.00022053253];
t_1             = [0.25;1.25;1.5;0.25;0.875];
d_1             = [1.0;1.0;1.0;3.0;7.0];
a_2             = [0.41541823;0.71335657;0.00030354234;-0.36643143;-0.0014407781;-0.089166707;-0.023699887];
t_2             = [2.375;2.0;2.125;3.5;6.5;4.75;12.5];
d_2             = [1.0;2.0;5.0;1.0;1.0;4.0;2.0];
p_2             = [1.0;1.0;1.0;2.0;2.0;2.0;3.0];
u_2             = d_2 - p_2.*delta.^p_2;
v_2             = exp(-delta.^p_2);
phir            = a_1'*(delta.^d_1.*tau.^t_1) + a_2'*(v_2.*delta.^d_2.*tau.^t_2);
phi_taur        = a_1'*(delta.^d_1.*t_1.*tau.^(t_1 - 1)) + a_2'*(v_2.*delta.^d_2.*t_2.*tau.^(t_2 - 1));
phi_deltar      = a_1'*(d_1.*delta.^(d_1 - 1).*tau.^t_1) + a_2'*(v_2.*delta.^(d_2 - 1).*u_2.*tau.^t_2);
phi_tautaur     = a_1'*(delta.^d_1.*t_1.*(t_1 - 1).*tau.^(t_1 - 2)) + a_2'*(v_2.*delta.^d_2.*t_2.*(t_2 - 1).*tau.^(t_2 - 2));
phi_deltadeltar = a_1'*(d_1.*(d_1 - 1).*delta.^(d_1 - 2).*tau.^t_1) + a_2'*(v_2.*delta.^(d_2 - 2).*(u_2.*(u_2 - 1) - p_2.^2.*delta.^p_2).*tau.^t_2);
phi_deltataur   = a_1'*(d_1.*delta.^(d_1 - 1).*t_1.*tau.^(t_1 - 1)) + a_2'*(v_2.*delta.^(d_2 - 1).*t_2.*u_2.*tau.^(t_2 - 1));
g_1             = N*R*(phir - tau*phi_taur);
g_2             = -R*T_c*rho_c*delta^2*phi_deltar/tau;
g_i             = R*T_c*(phir + delta*phi_deltar)/tau;
H_11            = N*R*tau^3*phi_tautaur/T_c;
H_21            = -R*delta^2*rho_c*(phi_deltar - tau*phi_deltataur);
H_i1            = R*(phir - tau*phi_taur + delta*(phi_deltar - tau*phi_deltataur));
H_22            = R*T_c*rho_c^2*delta^3*(delta*phi_deltadeltar + 2*phi_deltar)/(tau*N);
H_i2            = -R*T_c*rho_c*delta^2*(delta*phi_deltadeltar + 2*phi_deltar)/(tau*N);
H_ii            = R*T_c*delta*(delta*phi_deltadeltar + 2*phi_deltar)/(tau*N);
S.g             = [g_1;g_2;g_i];
S.H             = [H_11,H_21,H_i1;H_21,H_22,H_i2;H_i1,H_i2,H_ii];

function [S] = ModTVN_ideal_idealgas_23028984(x)

R     = 8.314511984;
pcirc = [101325.0];
xcirc = [1.0];
T     = x(1);
V     = x(2);
i     = [3];
n     = x(i);
NR    = sum(n)*R;
NRT   = NR*T;
e     = [1];
p     = R*log(R*T*(n./pcirc)/V);
g_1   = n'*p;
g_2   = -NRT/V;
g_i   = T*p;
H_11  = NR/T;
H_21  = -NR/V;
H_i1  = p + R*e;
H_22  = NRT/V^2;
H_i2  = -R*(T/V)*e;
H_ii  = R*T*diag(e./n);
S.g   = [g_1;g_2;g_i];
S.H   = [H_11,H_21',H_i1';H_21,H_22,H_i2';H_i1,H_i2,H_ii];

function [S] = MuT_cp_dippr_23809380(T)

t_0        = [298.15];
t_min      = [50.0];
t_max      = [5000.0];
C          = [29.37,34.54,1428.0,26.4,588.0];
c_1        = C(:,1);
c_2        = C(:,2);
c_3        = C(:,3);
c_4        = C(:,4);
c_5        = C(:,5);
c_p        = c_1 + c_2.*c_3.^2./(T^2*sinh(c_3/T).^2) + c_4.*c_5.^2./(T^2*cosh(c_5/T).^2);
h          = c_1*T + c_2.*c_3.*cosh(c_3/T)./sinh(c_3/T) - c_4.*c_5.*sinh(c_5/T)./cosh(c_5/T) - c_1.*t_0 - c_2.*c_3.*cosh(c_3./t_0)./sinh(c_3./t_0) + c_4.*c_5.*sinh(c_5./t_0)./cosh(c_5./t_0);
s          = c_1*log(T) + c_4.*(log(cosh(c_5/T)) - c_5.*sinh(c_5/T)./cosh(c_5/T)/T) + c_2.*(c_3.*cosh(c_3/T)./sinh(c_3/T)/T - log(sinh(c_3/T))) - c_1.*log(t_0) - c_4.*(log(cosh(c_5./t_0)) - c_5.*sinh(c_5./t_0)./cosh(c_5./t_0)./t_0) - c_2.*(c_3.*cosh(c_3./t_0)./sinh(c_3./t_0)./t_0 - log(sinh(c_3./t_0)));
S.mu       = h - T*s;
S.dmudT    = -s;
S.d2mudTdT = -(c_p/T);
i          = [1];
[S_1]      = MuT_hs_h0_23796420(T);
S.mu(i)    = S.mu(i) + S_1.mu;
S.dmudT(i) = S.dmudT(i) + S_1.dmudT;

function [S] = MuT_hs_h0_23796420(T)

hcirc   = [-393510.0];
[S]     = MuT_hs_s0_23790264(T);
S.mu    = S.mu + hcirc;
S.dmudT = S.dmudT + [0];

function [S] = MuT_hs_s0_23790264(T)

scirc   = [213.677];
S.mu    = -(T*scirc);
S.dmudT = -scirc;

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function plotSDCon(sdcone,data)
%PLOTSDCON Plot Semidefinite Constraints on the current figure
%   plotSDCon(sdcone)

%   Copyright (C) 2013 Jonathan Currie (I2C2)

xl = xlim; yl = ylim;
hold on;

%Colour
dkr = [179/255 0.0 70/255];

%Determine number of quad constraints
if(iscell(sdcone))
    no = length(sdcone);
else
    no = 1;
end
%Generate Constraint Surface Points
[x1,x2] = meshgrid(linspace(xl(1),xl(2),data.npts),linspace(yl(1),yl(2),data.npts));

%For each semidefinite constraint, plot
for i = 1:no
    %Get Constraint Variables & Form Constraint Function
    if(iscell(sdcone))
        m = sqrt(size(sdcone{i},1));
        C = reshape(sdcone{i}(:,1),m,m); A0 = reshape(sdcone{i}(:,2),m,m);
        switch(data.ndec)
            case 1
                con = @(x) min(eig(A0*x(1) - C));
            case 2
                A1 = reshape(sdcone{i}(:,3),m,m);
                con = @(x) min(eig(A0*x(1) + A1*x(2) - C));
            case 3
                A1 = reshape(sdcone{i}(:,3),m,m); A2 = reshape(sdcone{i}(:,4),m,m);
                con = @(x) min(eig(A0*x(1) + A1*x(2) + A2*x(3) - C));
            case 4
                A1 = reshape(sdcone{i}(:,3),m,m); A2 = reshape(sdcone{i}(:,4),m,m);
                A3 = reshape(sdcone{i}(:,5),m,m);
                con = @(x) min(eig(A0*x(1) + A1*x(2) + A2*x(3) + A3*x(4) - C));
            case 5
                A1 = reshape(sdcone{i}(:,3),m,m); A2 = reshape(sdcone{i}(:,4),m,m);
                A3 = reshape(sdcone{i}(:,5),m,m); A4 = reshape(sdcone{i}(:,6),m,m);
                con = @(x) min(eig(A0*x(1) + A1*x(2) + A2*x(3) + A3*x(4) + A4*x(5) - C));
            otherwise
                error('Plotting Semidefinite constraints up to 5D is only supported');
        end
    else
        m = sqrt(size(sdcone,1));
        C = reshape(sdcone(:,1),m,m);
        A0 = reshape(sdcone(:,2),m,m);
        switch(data.ndec)
            case 1
                con = @(x) min(eig(A0*x(1) - C));
            case 2
                A1 = reshape(sdcone(:,3),m,m);
                con = @(x) min(eig(A0*x(1) + A1*x(2) - C));
            case 3
                A1 = reshape(sdcone(:,3),m,m);
                A2 = reshape(sdcone(:,4),m,m);
                con = @(x) min(eig(A0*x(1) + A1*x(2) + A2*x(3) - C));
            case 4
                A1 = reshape(sdcone(:,3),m,m);
                A2 = reshape(sdcone(:,4),m,m);
                A3 = reshape(sdcone(:,5),m,m);
                con = @(x) min(eig(A0*x(1) + A1*x(2) + A2*x(3) + A3*x(4) - C));
            case 5
                A1 = reshape(sdcone(:,3),m,m);
                A2 = reshape(sdcone(:,4),m,m);
                A3 = reshape(sdcone(:,5),m,m);
                A4 = reshape(sdcone(:,6),m,m);
                con = @(x) min(eig(A0*x(1) + A1*x(2) + A2*x(3) + A3*x(4) + A4*x(5) - C));
            otherwise
                error('Plotting Semidefinite constraints up to 5D is only supported');
        end
    end
    %Plot Each Quad Con as General Nonlinear Constraint
    plotNLCon(con,[],1e-6,Inf,x1,x2,dkr,data);
end
hold off;

% OLD CODE
% %Plot Semidefinite Constraints (Inefficient.. ideas appreciated!)
% [x1,x2] = meshgrid(linspace(xl(1),xl(2),npts),linspace(yl(1),yl(2),npts));
% nox = size(x1);
% noy = size(x2);
% obj = zeros(nox(1),noy(2));
% if(iscell(sdcone))
%     no = length(sdcone);
% else
%     no = 1;
% end
% for i = 1:no
%     %get vars
%     if(iscell(sdcone))
%         C = sdcone{i}(:,1); A0 = sdcone{i}(:,2); A1 = sdcone{i}(:,3);
%     else
%         m = sqrt(size(sdcone,1));
%         C = reshape(sdcone(:,1),m,m);
%         A0 = reshape(sdcone(:,2),m,m);
%         A1 = reshape(sdcone(:,3),m,m);
%     end
%     % create surface
%     for n = 1:nox(1)
%         for m = 1:noy(2)
%             obj(n,m) = min(eig(A0*x1(n,m) + A1*x2(n,m) - C));
%         end
%     end
%     c = contour(x1,x2,obj,'color',dkr,'levellist',0);
%     %Plot Hatch
%     if(~isempty(c))
%         %See if we have multiple contours (non-convex or sd)
%         len = size(c,2)-1;
%         if(c(2,1) ~= len)
%             %Build contour array
%             cstrt = 2; cend = []; n = 2; ind = 1;
%             while(ind <= len)
%                 ind = ind + c(2,ind) + 1;
%                 cend(n-1) = ind-1; %#ok<AGROW>
%                 cstrt(n) = ind+1;  %#ok<AGROW>
%                 n = n + 1;
%             end
%         else
%             cstrt = 2;
%             cend = len;
%         end
%         %Plot each contour hatch
%         for n = 1:length(cend)
%             %Get contour vectors
%             vecx = diff(c(1,cstrt(n):cend(n)));
%             vecy = diff(c(2,cstrt(n):cend(n)));
%             if(isempty(vecx) || isempty(vecy))
%                 continue;
%             end
%             %Rotate hatch lines based on infeasible region
%             xt = [c(1,cstrt(n))+vecy(1) c(2,cstrt(n))-vecx(1)]'; %check rotated -90
%             fval = all(eig(A0*xt(1) + A1*xt(2) - C) >= 1e-6);
%             if(fval) %rotate 90
%                 hvecx = -vecy;
%                 hvecy = vecx;
%             else %rotate -90
%                 hvecx = vecy;
%                 hvecy = -vecx;
%             end
%             %Normalize
%             av = mean(sqrt(hvecx.^2 + hvecy.^2));
%             dirs = atan2(hvecy,hvecx);
%             hvecx = av*cos(dirs);
%             hvecy = av*sin(dirs);
%             %Shift origin
%             hvecx = c(1,cstrt(n):cend(n)-1) + hvecx;
%             hvecy = c(2,cstrt(n):cend(n)-1) + hvecy;
%             %Plot
%             line([c(1,cstrt(n):cend(n)-1)' hvecx']',[c(2,cstrt(n):cend(n)-1)' hvecy']','Color',dkr)
%         end
%     else
%         optiwarn('opti:plot','Cannot plot semidefinite constraint as contour data is empty!');
%     end
% end
%
% hold off;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [cov_ewma,corr_ewma,vola_ewma] = ewma_covariance(data,lambda)

% calculates the RiskMetrics "Technical Document" (1996) exponentially
% weighted covariance matrix (p.179), correlation and volatilities.
%
% Input:
% data - needs to be in format T x k with T = # observations, k = # assets
% lambda = decay factor

[r,c]         = size(data);
data_mwb      = data-repmat(mean(data,1),r,1);
lambdavec     = lambda.^(0:1:r-1)';
data_tilde    = repmat(sqrt(lambdavec),1,c).*data_mwb;

cov_ewma      = 1/sum(lambdavec)*(data_tilde'*data_tilde);
corr_ewma     = zeros(c);
vola_ewma     = zeros(c,1);

for i = 1:c
    for j = 1:c
        corr_ewma(i,j) = cov_ewma(i,j)/sqrt(cov_ewma(i,i)*cov_ewma(j,j));
    end
    vola_ewma(i) = sqrt(cov_ewma(i,i));
end

author: Sable

repo: mcbench-benchmarks


text:
function a = hlp_wrapresults(f,varargin)
% Wraps all outputs produced by function f for the given arguments into a cell array.
% Results = hlp_wrapresults(Function, Arguments...)
%
% In:
%   Function     : some function handle to execute
%
%   Arguments... : list of arguments to pass to the function
%
% Out:
%   Results : cell array of all function results
%
%
% Notes:
%   It is not (currently) possible to efficiently determine the number of out-args for a
%   varargout function; in this case, at most 10 outputs are supported.
%
% Examples:
%   % wrap both outputs of a particualr sort() call into a cell array
%   results = hlp_wrapresults(@sort,[1 4 2 1 3 6 0])
%
% See also:
%   hlp_getresult
%
%					         	 Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-03-28

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

% find out how many arguments f can maximally produce
len = nargout(f);
% for varargout functions, we assume at most 10
if len < 0
    len = 10; end
% ... then invoke the function
while len >= 1
    try
        % get the appropriate number of results from f()
        [a{1:len}] = f(varargin{:});
        return;
    catch e
        % got an exception, check if it is outarg-related
        if ~any(strcmp(e.identifier,{'MATLAB:TooManyOutputs','MATLAB:maxlhs','MATLAB:unassignedOutputs'}))
            % it isn't: propagate the error
            settings = dbstatus;
            if any(strcmp({settings.cond},'error'))
                % if in dbstop if error mode, we're re-running the line to get the debugger to stop at the right place
                fprintf('hlp_wrapresults: caught error while in "dbstop if error" mode; re-running offending code to get break point...\n');
                fprintf('  error traceback was: %s\n',hlp_handleerror(e));
                [a{1:len}] = f(varargin{:});
            else
                rethrow(e);
            end
        end
        % but if it was, we need to retry with fewer out-arguments
        len = len-1;
    end
end

% len = 0: f produces no outputs
f(varargin{:});
a = {};

repo: REST


text:
function [cfg] = ft_multiplotTFR(cfg, data)

% FT_MULTIPLOTTFR plots the time-frequency representations of power or coherence
% in a topographical layout. The plots of the indivual sensors are arranged
% according to their location specified in the layout.
%
% Use as
%   ft_multiplotTFR(cfg, data)
%
% The data can be a time-frequency representation of power or coherence
% that was computed using the FT_FREQANALYSIS or FT_FREQDESCRIPTIVES
% functions.
%
% The configuration can have the following parameters:
%   cfg.parameter        = field to be represented as color, for example 'powspctrm' or 'cohspctrm' (default depends on data.dimord)
%   cfg.maskparameter    = field in the data to be used for masking of data, can be logical (e.g. significant data points) or numerical (e.g. t-values).
%                        (not possible for mean over multiple channels, or when input contains multiple subjects
%                        or trials)
%   cfg.maskstyle        = style used to masking, 'opacity', 'saturation', or 'outline' (default = 'opacity')
%                        'outline' can only be used with a logical cfg.maskparameter
%                        use 'saturation' or 'outline' when saving to vector-format (like *.eps) to avoid all sorts of image-problems
%   cfg.maskalpha        = alpha value between 0 (transparent) and 1 (opaque) used for masking areas dictated by cfg.maskparameter (default = 1)
%                        (will be ignored in case of numeric cfg.maskparameter or if cfg.maskstyle = 'outline')
%   cfg.masknans         = 'yes' or 'no' (default = 'yes')
%   cfg.xlim             = 'maxmin', 'maxabs', 'zeromax', 'minzero', or [xmin xmax] (default = 'maxmin')
%   cfg.ylim             = 'maxmin', 'maxabs', 'zeromax', 'minzero', or [ymin ymax] (default = 'maxmin')
%   cfg.zlim             = plotting limits for color dimension, 'maxmin', 'maxabs', 'zeromax', 'minzero', or [zmin zmax] (default = 'maxmin')
%   cfg.gradscale        = number, scaling to apply to the MEG gradiometer channels prior to display
%   cfg.magscale         = number, scaling to apply to the MEG magnetometer channels prior to display
%   cfg.channel          = Nx1 cell-array with selection of channels (default = 'all'), see FT_CHANNELSELECTION for details
%   cfg.refchannel       = name of reference channel for visualising connectivity, can be 'gui'
%   cfg.baseline         = 'yes', 'no' or [time1 time2] (default = 'no'), see FT_FREQBASELINE
%   cfg.baselinetype     = 'absolute', 'relative', 'relchange', 'normchange', 'db' or 'zscore' (default = 'absolute')
%   cfg.trials           = 'all' or a selection given as a 1xN vector (default = 'all')
%   cfg.box              = 'yes', 'no', whether to draw a box around each graph (default = 'no', if maskparameter is given default = 'yes')
%   cfg.hotkeys          = enables hotkeys (up/down arrows) for dynamic colorbar adjustment
%   cfg.colorbar         = 'yes', 'no' (default = 'no')
%   cfg.colorbartext     = string indicating the text next to colorbar
%   cfg.colormap         = string, or Nx3 matrix, see FT_COLORMAP
%   cfg.showlabels       = 'yes', 'no' (default = 'no')
%   cfg.showoutline      = 'yes', 'no' (default = 'no')
%   cfg.showscale        = 'yes', 'no' (default = 'yes')
%   cfg.showcomment      = 'yes', 'no' (default = 'yes')
%   cfg.comment          = string of text (default = date + limits)
%                          Add 'comment' to graph (according to COMNT in the layout)
%   cfg.limittext        = add user-defined text instead of cfg.comment, (default = cfg.comment)
%   cfg.fontsize         = font size of comment and labels (if present) (default = 8)
%   cfg.fontweight       = font weight of comment and labels (if present)
%   cfg.interactive      = Interactive plot 'yes' or 'no' (default = 'yes')
%                          In a interactive plot you can select areas and produce a new
%                          interactive plot when a selected area is clicked. Multiple areas
%                          can be selected by holding down the SHIFT key.
%   cfg.figure           = 'yes' or 'no', whether to open a new figure. You can also specify a figure handle from FIGURE, GCF or SUBPLOT. (default = 'yes')
%   cfg.position         = location and size of the figure, specified as [left bottom width height] (default is automatic)
%   cfg.renderer         = string, 'opengl', 'zbuffer', 'painters', see RENDERERINFO (default is automatic, try 'painters' when it crashes)
%   cfg.directionality   = '', 'inflow' or 'outflow' specifies for
%                          connectivity measures whether the inflow into a
%                          node, or the outflow from a node is plotted. The
%                          (default) behavior of this option depends on the dimor
%                          of the input data (see below).
%   cfg.layout           = specify the channel layout for plotting using one of
%                          the supported ways (see below).
%
% The following options for the scaling of the EEG, EOG, ECG, EMG, MEG and NIRS channels
% is optional and can be used to bring the absolute numbers of the different
% channel types in the same range (e.g. fT and uV). The channel types are determined
% from the input data using FT_CHANNELSELECTION.
%   cfg.eegscale         = number, scaling to apply to the EEG channels prior to display
%   cfg.eogscale         = number, scaling to apply to the EOG channels prior to display
%   cfg.ecgscale         = number, scaling to apply to the ECG channels prior to display
%   cfg.emgscale         = number, scaling to apply to the EMG channels prior to display
%   cfg.megscale         = number, scaling to apply to the MEG channels prior to display
%   cfg.gradscale        = number, scaling to apply to the MEG gradiometer channels prior to display (in addition to the cfg.megscale factor)
%   cfg.magscale         = number, scaling to apply to the MEG magnetometer channels prior to display (in addition to the cfg.megscale factor)
%   cfg.nirsscale        = number, scaling to apply to the NIRS channels prior to display
%   cfg.mychanscale      = number, scaling to apply to the channels specified in cfg.mychan
%   cfg.mychan           = Nx1 cell-array with selection of channels
%   cfg.chanscale        = Nx1 vector with scaling factors, one per channel specified in cfg.channel
%
% For the plotting of directional connectivity data the cfg.directionality
% option determines what is plotted. The default value and the supported
% functionality depend on the dimord of the input data. If the input data
% is of dimord 'chan_chan_XXX', the value of directionality determines
% whether, given the reference channel(s), the columns (inflow), or rows
% (outflow) are selected for plotting. In this situation the default is
% 'inflow'. Note that for undirected measures, inflow and outflow should
% give the same output. If the input data is of dimord 'chancmb_XXX', the
% value of directionality determines whether the rows in data.labelcmb are
% selected. With 'inflow' the rows are selected if the refchannel(s) occur in
% the right column, with 'outflow' the rows are selected if the
% refchannel(s) occur in the left column of the labelcmb-field. Default in
% this case is '', which means that all rows are selected in which the
% refchannel(s) occur. This is to robustly support linearly indexed
% undirected connectivity metrics. In the situation where undirected
% connectivity measures are linearly indexed, specifying 'inflow' or
% 'outflow' can result in unexpected behavior.
%
% The layout defines how the channels are arranged and what the size of each
% subplot is. You can specify the layout in a variety of ways:
%  - you can provide a pre-computed layout structure, see FT_PREPARE_LAYOUT
%  - you can give the name of an ASCII layout file with extension *.lay
%  - you can give the name of an electrode file
%  - you can give an electrode definition, i.e. "elec" structure
%  - you can give a gradiometer definition, i.e. "grad" structure
% If you do not specify any of these and the data structure contains an
% electrode or gradiometer structure (common for MEG data, since the header
% of the MEG datafile contains the gradiometer information), that will be
% used for creating a layout. If you want to have more fine-grained control
% over the layout of the subplots, you should create your own layout file.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
% If you specify this option the input data will be read from a *.mat
% file on disk. This mat files should contain only a single variable named 'data',
% corresponding to the input structure. For this particular function, the
% data should be provided as a cell-array.
%
% See also:
%   FT_MULTIPLOTER, FT_SINGLEPLOTER, FT_SINGLEPLOTTFR, FT_TOPOPLOTER, FT_TOPOPLOTTFR,
%   FT_PREPARE_LAYOUT

% Undocumented local options:
% cfg.channel
% cfg.layoutname
% cfg.orient = landscape/portrait

% Copyright (C) 2003-2006, Ole Jensen
% Copyright (C) 2007-2011, Roemer van der Meij & Jan-Mathijs Schoffelen
% Copyright (C) 2012-2017, F.C. Donders Centre
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEVELOPERS NOTE: This code is organized in a similar fashion for multiplot/singleplot/topoplot
% and for ER/TFR and should remain consistent over those 6 functions.
% Section 1: general cfg handling that is independent from the data
% Section 2: data handling, this also includes converting bivariate (chan_chan and chancmb) into univariate data
% Section 3: select the data to be plotted and determine min/max range
% Section 4: do the actual plotting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Section 1: general cfg handling that is independent from the data

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar data
ft_preamble provenance data

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input data is valid for this function
data = ft_checkdata(data, 'datatype', 'freq');

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'channels', 'trial'}); % prevent accidental typos, see issue 1729
cfg = ft_checkconfig(cfg, 'renamed',    {'cohrefchannel', 'refchannel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'matrixside', 'directionality'});
cfg = ft_checkconfig(cfg, 'renamed',    {'zparam', 'parameter'});
cfg = ft_checkconfig(cfg, 'renamedval', {'directionality', 'feedback', 'inflow'});
cfg = ft_checkconfig(cfg, 'renamedval', {'directionality', 'feedforward', 'outflow'});
cfg = ft_checkconfig(cfg, 'renamedval', {'zlim', 'absmax', 'maxabs'});
cfg = ft_checkconfig(cfg, 'unused',     {'cohtargetchannel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'newfigure', 'figure'});

% set the defaults
cfg.parameter      = ft_getopt(cfg, 'parameter',    'powspctrm');
cfg.baseline       = ft_getopt(cfg, 'baseline',     'no');
cfg.baselinetype   = ft_getopt(cfg, 'baselinetype', 'absolute');
cfg.trials         = ft_getopt(cfg, 'trials', 'all', 1);
cfg.xlim           = ft_getopt(cfg, 'xlim',         'maxmin');
cfg.ylim           = ft_getopt(cfg, 'ylim',         'maxmin');
cfg.zlim           = ft_getopt(cfg, 'zlim',         'maxmin');
cfg.colorbar       = ft_getopt(cfg, 'colorbar',     'no');
cfg.colormap       = ft_getopt(cfg, 'colormap',      'default');
cfg.colorbartext   = ft_getopt(cfg, 'colorbartext', '');
cfg.comment        = ft_getopt(cfg, 'comment',       date);
cfg.limittext      = ft_getopt(cfg, 'limittext',    'default');
cfg.showlabels     = ft_getopt(cfg, 'showlabels',   'no');
cfg.showoutline    = ft_getopt(cfg, 'showoutline',  'no');
cfg.showscale      = ft_getopt(cfg, 'showscale',    'yes');
cfg.showcomment    = ft_getopt(cfg, 'showcomment',  'yes');
cfg.channel        = ft_getopt(cfg, 'channel',      'all');
cfg.fontsize       = ft_getopt(cfg, 'fontsize',      8);
cfg.fontweight     = ft_getopt(cfg, 'fontweight');
cfg.interactive    = ft_getopt(cfg, 'interactive',  'yes');
cfg.hotkeys        = ft_getopt(cfg, 'hotkeys',      'yes');
cfg.orient         = ft_getopt(cfg, 'orient',       'landscape');
cfg.maskalpha      = ft_getopt(cfg, 'maskalpha',     1);
cfg.masknans       = ft_getopt(cfg, 'masknans',     'yes');
cfg.maskparameter  = ft_getopt(cfg, 'maskparameter');
cfg.maskstyle      = ft_getopt(cfg, 'maskstyle',    'opacity');
cfg.directionality = ft_getopt(cfg, 'directionality', '');
cfg.figurename     = ft_getopt(cfg, 'figurename');
cfg.commentpos     = ft_getopt(cfg, 'commentpos',   'layout');
cfg.scalepos       = ft_getopt(cfg, 'scalepos',     'layout');
cfg.renderer       = ft_getopt(cfg, 'renderer'); % let MATLAB decide on the default

if ~isfield(cfg, 'box')
  if ~isempty(cfg.maskparameter)
    cfg.box = 'yes';
  else
    cfg.box = 'no';
  end
end

% check if the colormap is in the proper format
if ~isequal(cfg.colormap, 'default')
  if ischar(cfg.colormap)
    cfg.colormap = ft_colormap(cfg.colormap);
  elseif iscell(cfg.colormap)
    cfg.colormap = ft_colormap(cfg.colormap{:});
  elseif isnumeric(cfg.colormap) && size(cfg.colormap,2)~=3
    ft_error('cfg.colormap must be Nx3');
  end
  % the actual colormap will be set below
end

% this is needed for the figure title
if isfield(cfg, 'dataname') && ~isempty(cfg.dataname)
  dataname = cfg.dataname;
elseif isfield(cfg, 'inputfile') && ~isempty(cfg.inputfile)
  dataname = cfg.inputfile;
elseif nargin>1
  dataname = arrayfun(@inputname, 2:nargin, 'UniformOutput', false);
else
  dataname = {};
end

%% Section 2: data handling, this also includes converting bivariate (chan_chan and chancmb) into univariate data

hastime = isfield(data, 'time');
hasfreq = isfield(data, 'freq');

assert((hastime && hasfreq), 'please use ft_multiplotTFR for time-frequency data only');

xparam = ft_getopt(cfg, 'xparam', 'time');
yparam = ft_getopt(cfg, 'yparam', 'freq');

% check whether rpt/subj is present and remove if necessary
dimord = getdimord(data, cfg.parameter);
dimtok = tokenize(dimord, '_');
hasrpt = any(ismember(dimtok, {'rpt' 'subj'}));

if ~hasrpt
  assert(isequal(cfg.trials, 'all') || isequal(cfg.trials, 1), 'incorrect specification of cfg.trials for data without repetitions');
else
  assert(~isempty(cfg.trials), 'empty specification of cfg.trials for data with repetitions');
end

% parse cfg.channel
if isfield(cfg, 'channel') && isfield(data, 'label')
  cfg.channel = ft_channelselection(cfg.channel, data.label);
elseif isfield(cfg, 'channel') && isfield(data, 'labelcmb')
  cfg.channel = ft_channelselection(cfg.channel, unique(data.labelcmb(:)));
end

% Apply baseline correction:
if ~strcmp(cfg.baseline, 'no')
  tmpcfg = removefields(cfg, {'inputfile', 'reproducescript'});
  % keep mask-parameter if it is set
  if ~isempty(cfg.maskparameter)
    tempmask = data.(cfg.maskparameter);
  end
  data = ft_freqbaseline(tmpcfg, data);
  % put mask-parameter back if it is set
  if ~isempty(cfg.maskparameter)
    data.(cfg.maskparameter) = tempmask;
  end
end

% channels SHOULD be selected here, as no interactive action produces a new multiplot
tmpcfg = keepfields(cfg, {'channel', 'trials', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
if hasrpt
  tmpcfg.avgoverrpt = 'yes';
else
  tmpcfg.avgoverrpt = 'no';
end
tmpvar = data;
[data] = ft_selectdata(tmpcfg, data);
% restore the provenance information
[cfg, data] = rollback_provenance(cfg, data);

if isfield(tmpvar, cfg.maskparameter) && ~isfield(data, cfg.maskparameter)
  % the mask parameter is not present after ft_selectdata, because it is
  % not included in all input arguments. Make the same selection and copy
  % it over
  tmpvar = ft_selectdata(tmpcfg, tmpvar);
  data.(cfg.maskparameter) = tmpvar.(cfg.maskparameter);
end

clear tmpvar tmpcfg dimord dimtok hastime hasfreq hasrpt

% ensure that the preproc specific options are located in the cfg.preproc
% substructure, but also ensure that the field 'refchannel' remains at the
% highest level in the structure. This is a little hack by JM because the field
% refchannel can relate to connectivity or to an EEG reference.

if isfield(cfg, 'refchannel'), refchannelincfg = cfg.refchannel; cfg = rmfield(cfg, 'refchannel'); end
cfg = ft_checkconfig(cfg, 'createsubcfg',  {'preproc'});
if exist('refchannelincfg', 'var'), cfg.refchannel  = refchannelincfg; end

if ~isempty(cfg.preproc)
  % preprocess the data, i.e. apply filtering, baselinecorrection, etc.
  fprintf('applying preprocessing options\n');
  if ~isfield(cfg.preproc, 'feedback')
    cfg.preproc.feedback = cfg.interactive;
  end
  data = ft_preprocessing(cfg.preproc, data);
end

% Handle the bivariate case
dimord = getdimord(data, cfg.parameter);
if startsWith(dimord, 'chan_chan_') || startsWith(dimord, 'chancmb_')
  % convert the bivariate data to univariate and call this plotting function with univariate input
  cfg.originalfunction = 'ft_multiplotTFR';
  cfg.trials = 'all'; % trial selection has been taken care off
  bivariate_common(cfg, data);
  return
end

% Apply channel-type specific scaling
fn = fieldnames(cfg);
fn = setdiff(fn, {'skipscale', 'showscale', 'gridscale'}); % these are for the layout and plotting, not for CHANSCALE_COMMON
fn = fn(endsWith(fn, 'scale') | startsWith(fn, 'mychan') | strcmp(fn, 'channel') | strcmp(fn, 'parameter'));
tmpcfg = keepfields(cfg, fn);
if ~isempty(tmpcfg)
  data = chanscale_common(tmpcfg, data);
  % remove the scaling fields from the configuration, to prevent them from being called again in interactive mode
  % but keep the parameter and channel field
  cfg = removefields(cfg, setdiff(fn, {'parameter', 'channel'}));
else
  % do nothing
end

%% Section 3: select the data to be plotted and determine min/max range

% Read or create the layout that will be used for plotting
tmpcfg = keepfields(cfg, {'layout', 'channel', 'rows', 'columns', 'commentpos', 'skipcomnt', 'scalepos', 'skipscale', 'projection', 'viewpoint', 'rotate', 'width', 'height', 'elec', 'grad', 'opto', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
cfg.layout = ft_prepare_layout(tmpcfg, data);

% Take the subselection of channels that is contained in the layout, this is the same in all datasets
[selchan, sellay] = match_str(data.label, cfg.layout.label);

% Get physical min/max range of x, i.e. time
if ~isnumeric(cfg.xlim)
  xmin = nanmin(data.(xparam));
  xmax = nanmax(data.(xparam));
  switch cfg.xlim
    case 'maxmin'
      % keep them as they are
    case 'maxabs'
      xmax = max(abs(xmax), abs(xmin));
      xmin = -xmax;
    case 'zeromax'
      xmin = 0;
    case 'minzero'
      xmax = 0;
    otherwise
      ft_error('invalid specification of cfg.xlim');
  end % switch
else
  xmin = cfg.xlim(1);
  xmax = cfg.xlim(2);
end

% Get the index of the nearest bin, this is the same in all datasets
xminindx = nearest(data.(xparam), xmin);
xmaxindx = nearest(data.(xparam), xmax);
xmin = data.(xparam)(xminindx);
xmax = data.(xparam)(xmaxindx);
selx = xminindx:xmaxindx;
xval = data.(xparam)(selx);

% Get physical min/max range of y, i.e. freq
if ~isnumeric(cfg.ylim)
  ymin = nanmin(data.(yparam));
  ymax = nanmax(data.(yparam));
  switch cfg.ylim
    case 'maxmin'
      % keep them as they are
    case 'maxabs'
      ymax = max(abs(ymax), abs(ymin));
      ymin = -ymax;
    case 'zeromax'
      ymin = 0;
    case 'minzero'
      ymax = 0;
    otherwise
      ft_error('invalid specification of cfg.ylim');
  end % switch
else
  ymin = cfg.ylim(1);
  ymax = cfg.ylim(2);
end

% Get the index of the nearest bin, this is the same in all datasets
yminindx = nearest(data.(yparam), ymin);
ymaxindx = nearest(data.(yparam), ymax);
ymin = data.(yparam)(yminindx);
ymax = data.(yparam)(ymaxindx);
sely = yminindx:ymaxindx;
yval = data.(yparam)(sely);

% test if X and Y are linearly spaced (to within 10^-12): % FROM UIMAGE
dx = min(diff(xval));  % smallest interval for X
dy = min(diff(yval));  % smallest interval for Y
evenx = all(abs(diff(xval)/dx-1)<1e-12);     % true if X is linearly spaced
eveny = all(abs(diff(yval)/dy-1)<1e-12);     % true if Y is linearly spaced

if ~evenx || ~eveny
  ft_warning('(one of the) axis is/are not evenly spaced, but plots are made as if axis are linear')
end

% masking is only possible for evenly spaced axis
if strcmp(cfg.masknans, 'yes') && (~evenx || ~eveny)
  ft_warning('(one of the) axis are not evenly spaced -> nans cannot be masked out -> cfg.masknans is set to ''no'';')
  cfg.masknans = 'no';
end

% the usual data is chan_freq_time, but other dimords should also work
dimtok = tokenize(dimord, '_');
datamatrix = data.(cfg.parameter);
[c, ia, ib] = intersect({'chan', yparam, xparam}, dimtok, 'stable');
datamatrix = permute(datamatrix, ib);
datamatrix = datamatrix(selchan, sely, selx);

if ~isempty(cfg.maskparameter)
  % one value for each channel-freq-time point
  maskmatrix = data.(cfg.maskparameter)(selchan, sely, selx);
  if islogical(maskmatrix) && any(strcmp(cfg.maskstyle, {'saturation', 'opacity'}))
    maskmatrix = double(maskmatrix);
    maskmatrix(~maskmatrix) = cfg.maskalpha;
  elseif isnumeric(maskmatrix)
    if strcmp(cfg.maskstyle, 'outline')
      error('Outline masking with a numeric cfg.maskparameter is not supported. Please use a logical mask instead.')
    end
    if cfg.maskalpha ~= 1
      ft_warning('Using field "%s" for masking, cfg.maskalpha is ignored.', cfg.maskparameter)
    end
    % scale mask between 0 and 1
    minval = min(maskmatrix(:));
    maxval = max(maskmatrix(:));
    maskmatrix = (maskmatrix - minval) / (maxval-minval);
  end
else
  % create an Nx0x0 matrix
  maskmatrix = zeros(length(selchan), 0, 0);
end

chanX      = cfg.layout.pos(sellay, 1);
chanY      = cfg.layout.pos(sellay, 2);
chanWidth  = cfg.layout.width(sellay);
chanHeight = cfg.layout.height(sellay);

%% Section 4: do the actual plotting

% open a new figure, or add it to the existing one
% note that in general adding a TFR to an existing one does not make sense, since they will overlap
open_figure(keepfields(cfg, {'figure', 'position', 'visible', 'renderer', 'figurename', 'title'}));
hold on

% Get physical z-axis range (color axis):
if strcmp(cfg.zlim, 'maxmin')
  zmin = nanmin(datamatrix(:));
  zmax = nanmax(datamatrix(:));
elseif strcmp(cfg.zlim, 'maxabs')
  zmin = -nanmax(abs(datamatrix(:)));
  zmax =  nanmax(abs(datamatrix(:)));
elseif strcmp(cfg.zlim, 'zeromax')
  zmin = 0;
  zmax = nanmax(datamatrix(:));
elseif strcmp(cfg.zlim, 'minzero')
  zmin = nanmin(datamatrix(:));
  zmax = 0;
else
  zmin = cfg.zlim(1);
  zmax = cfg.zlim(2);
end

% Plot channels
for k=1:length(selchan)
  cdata = shiftdim(datamatrix(k, :, :));
  if ~isempty(cfg.maskparameter)
    mdata = shiftdim(maskmatrix(k, :, :));
  end

  % Draw plot (and mask Nan's with maskfield if requested)
  if isequal(cfg.masknans, 'yes') && isempty(cfg.maskparameter)
    nans_mask = ~isnan(cdata);
    mask = double(nans_mask);
    ft_plot_matrix(cdata, 'clim', [zmin zmax], 'tag', 'cip', 'highlightstyle', cfg.maskstyle, 'highlight', mask, 'hpos', chanX(k), 'vpos', chanY(k), 'width', chanWidth(k), 'height', chanHeight(k))
  elseif isequal(cfg.masknans, 'yes') && ~isempty(cfg.maskparameter)
    nans_mask = ~isnan(cdata);
    mask = nans_mask .* mdata;
    mask = double(mask);
    ft_plot_matrix(cdata, 'clim', [zmin zmax], 'tag', 'cip', 'highlightstyle', cfg.maskstyle, 'highlight', mask, 'hpos', chanX(k), 'vpos', chanY(k), 'width', chanWidth(k), 'height', chanHeight(k))
  elseif isequal(cfg.masknans, 'no') && ~isempty(cfg.maskparameter)
    mask = mdata;
    mask = double(mask);
    ft_plot_matrix(cdata, 'clim', [zmin zmax], 'tag', 'cip', 'highlightstyle', cfg.maskstyle, 'highlight', mask, 'hpos', chanX(k), 'vpos', chanY(k), 'width', chanWidth(k), 'height', chanHeight(k))
  else
    ft_plot_matrix(cdata, 'clim', [zmin zmax], 'tag', 'cip', 'hpos', chanX(k), 'vpos', chanY(k), 'width', chanWidth(k), 'height', chanHeight(k))
  end

  % Currently the handle isn't being used below, this is here for possible use in the future
  h = findobj('tag', 'cip');
end % plot channels

% plot the layout, labels and outline
ft_plot_layout(cfg.layout, 'box', cfg.box, 'label', cfg.showlabels, 'outline', cfg.showoutline, 'point', 'no', 'mask', 'no', 'fontsize', cfg.fontsize, 'labelyoffset', 1.4*median(cfg.layout.height/2), 'labelalignh', 'center', 'chanindx', find(~ismember(cfg.layout.label, {'COMNT', 'SCALE'})) );

% apply the colormap
if ~isempty(cfg.colormap)
  set(gcf,  'colormap', cfg.colormap);
end

% Construct comment
if ~strcmp(cfg.comment, 'no')
  if ~strcmp(cfg.limittext, 'default')
    comment = cfg.limittext;
  else
    comment = cfg.comment;
    comment = sprintf('%0s\nxlim=[%.3g %.3g]', comment, xmin, xmax);
    comment = sprintf('%0s\nylim=[%.3g %.3g]', comment, ymin, ymax);
    comment = sprintf('%0s\nzlim=[%.3g %.3g]', comment, zmin, zmax);
  end
end

% Write comment
if strcmp(cfg.comment, 'no')
  comment_handle = [];
elseif strcmp(cfg.commentpos, 'title')
  comment_handle = title(comment, 'FontSize', cfg.fontsize);
elseif ~isempty(strcmp(cfg.layout.label, 'COMNT'))
  x_comment = cfg.layout.pos(strcmp(cfg.layout.label, 'COMNT'), 1);
  y_comment = cfg.layout.pos(strcmp(cfg.layout.label, 'COMNT'), 2);
  % 'HorizontalAlignment', 'left', 'VerticalAlignment', 'bottom',
  comment_handle = ft_plot_text(x_comment, y_comment, comment, 'FontSize', cfg.fontsize, 'FontWeight', cfg.fontweight);
else
  comment_handle = [];
end

% show scale
if istrue(cfg.showscale)
  k = cellstrmatch('SCALE', cfg.layout.label);
  if ~isempty(k)
    % Get average cdata across channels:
    cdata = shiftdim(mean(datamatrix, 1));

    % Draw plot (and mask Nan's with maskfield if requested)
    if isequal(cfg.masknans, 'yes') && isempty(cfg.maskparameter)
      mask = ~isnan(cdata);
      mask = double(mask);
      ft_plot_matrix(cdata, 'clim', [zmin zmax], 'tag', 'cip', 'highlightstyle', cfg.maskstyle, 'highlight', mask, 'hpos', cfg.layout.pos(k, 1), 'vpos', cfg.layout.pos(k, 2), 'width', cfg.layout.width(k), 'height', cfg.layout.height(k))
    elseif isequal(cfg.masknans, 'yes') && ~isempty(cfg.maskparameter)
      mask = ~isnan(cdata);
      mask = mask .* mdata;
      mask = double(mask);
      ft_plot_matrix(cdata, 'clim', [zmin zmax], 'tag', 'cip', 'highlightstyle', cfg.maskstyle, 'highlight', mask, 'hpos', cfg.layout.pos(k, 1), 'vpos', cfg.layout.pos(k, 2), 'width', cfg.layout.width(k), 'height', cfg.layout.height(k))
    elseif isequal(cfg.masknans, 'no') && ~isempty(cfg.maskparameter)
      mask = mdata;
      mask = double(mask);
      ft_plot_matrix(cdata, 'clim', [zmin zmax], 'tag', 'cip', 'highlightstyle', cfg.maskstyle, 'highlight', mask, 'hpos', cfg.layout.pos(k, 1), 'vpos', cfg.layout.pos(k, 2), 'width', cfg.layout.width(k), 'height', cfg.layout.height(k))
    else
      ft_plot_matrix(cdata, 'clim', [zmin zmax], 'tag', 'cip', 'hpos', cfg.layout.pos(k, 1), 'vpos', cfg.layout.pos(k, 2), 'width', cfg.layout.width(k), 'height', cfg.layout.height(k))
    end
  end
end % show scale

% show colorbar
if isfield(cfg, 'colorbar') && (strcmp(cfg.colorbar, 'yes'))
  c = colorbar;
  ylabel(c, cfg.colorbartext);
end

% Set colour axis
caxis([zmin zmax]);
if strcmp('yes', cfg.hotkeys)
  %  Attach data and cfg to figure and attach a key listener to the figure
  set(gcf, 'KeyPressFcn', {@key_sub, zmin, zmax})
end

axis tight
axis off

% Make the axis a little wider when boxes are shown
if strcmp(cfg.box, 'yes')
  abc = axis;
  axis(abc + [-1 +1 -1 +1]*mean(abs(abc))/10)
end

% Set orientation for printing if specified
if ~isempty(cfg.orient)
  orient(gcf, cfg.orient);
end

% set the figure window title
if ~isempty(dataname)
  set(gcf, 'Name', sprintf('%d: %s: %s', double(gcf), mfilename, join_str(', ', dataname)));
else
  set(gcf, 'Name', sprintf('%d: %s', double(gcf), mfilename));
end
set(gcf, 'NumberTitle', 'off');

% Make the figure interactive:
if strcmp(cfg.interactive, 'yes')
  % add the cfg/data/channel information to the figure under identifier linked to this axis
  ident                 = ['axh' num2str(round(sum(clock.*1e6)))]; % unique identifier for this axis
  set(gca,'tag',ident);
  info                  = guidata(gcf);
  info.(ident).x        = cfg.layout.pos(:, 1);
  info.(ident).y        = cfg.layout.pos(:, 2);
  info.(ident).label    = cfg.layout.label;
  info.(ident).dataname = dataname;
  info.(ident).cfg      = cfg;
  info.(ident).data     = data;
  info.(ident).commenth = comment_handle;
  guidata(gcf, info);
  set(gcf, 'WindowButtonUpFcn', {@ft_select_channel, 'multiple', true, 'callback', {@select_singleplotTFR}, 'event', 'WindowButtonUpFcn'});
  set(gcf, 'WindowButtonDownFcn', {@ft_select_channel, 'multiple', true, 'callback', {@select_singleplotTFR}, 'event', 'WindowButtonDownFcn'});
  set(gcf, 'WindowButtonMotionFcn', {@ft_select_channel, 'multiple', true, 'callback', {@select_singleplotTFR}, 'event', 'WindowButtonMotionFcn'});
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous data
ft_postamble provenance
ft_postamble savefig

% add a menu to the figure, but only if the current figure does not have subplots
menu_fieldtrip(gcf, cfg, false);

if ~ft_nargout
  % don't return anything
  clear cfg
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function l = cellstrmatch(str, strlist)
l = [];
for k=1:length(strlist)
  if strcmp(char(str), char(strlist(k)))
    l = [l k];
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION which is called after selecting channels in case of cfg.interactive='yes'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function select_singleplotTFR(label, varargin)
% fetch cfg/data based on axis indentifier given as tag
ident = get(gca,'tag');
info  = guidata(gcf);
cfg   = info.(ident).cfg;
data  = info.(ident).data;
if ~isempty(label)
  cfg = removefields(cfg, 'inputfile');   % the reading has already been done and varargin contains the data
  cfg.baseline = 'no';                    % make sure the next function does not apply a baseline correction again
  cfg.channel = label;
  cfg.dataname = info.(ident).dataname;   % put data name in here, this cannot be resolved by other means
  cfg.trials = 'all';                     % trial selection has already been taken care of
  fprintf('selected cfg.channel = {%s}\n', join_str(', ', cfg.channel));
  % ensure that the new figure appears at the same position
  cfg.figure = 'yes';
  cfg.position = get(gcf, 'Position');
  ft_singleplotTFR(cfg, data);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION which handles hot keys in the current plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function key_sub(handle, eventdata, varargin)
ident       = get(gca, 'tag');
info        = guidata(gcf);

climits = caxis;
incr_c  = abs(climits(2) - climits(1)) /10;

newz = climits;
if length(eventdata.Modifier) == 1 && strcmp(eventdata.Modifier{:}, 'control')
  % TRANSLATE by 10%
  switch eventdata.Key
    case 'pageup'
      newz = [climits(1)+incr_c climits(2)+incr_c];
    case 'pagedown'
      newz = [climits(1)-incr_c climits(2)-incr_c];
  end % switch
else
  % ZOOM by 10%
  switch eventdata.Key
    case 'pageup'
      newz = [climits(1)-incr_c climits(2)+incr_c];
    case 'pagedown'
      newz = [climits(1)+incr_c climits(2)-incr_c];
    case 'm'
      newz = [varargin{1} varargin{2}];
  end % switch
end % if

% update the color axis
caxis(newz);

if ~isempty(ident) && isfield(info.(ident), 'commenth') && ~isempty(info.(ident).commenth)
  commentstr = get(info.(ident).commenth, 'string');
  sel        = contains(commentstr, 'zlim');
  if any(sel)
    commentstr{sel} = sprintf('%0s=[%.3g %.3g]', 'zlim', newz(1), newz(2));
    set(info.(ident).commenth, 'string', commentstr);
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip