text:
function out = regress(dat, varargin)
% Multiple regression with an fmri_data object (dat), predicting brain data with a design matrix stored in dat.X (or vice versa)
%
% Regress dat.X on dat.dat at each voxel, and return voxel-wise statistic
% images. Each column of dat.X is a predictor in a multiple regression,
% and the intercept is the last column. Intercept will automatically be
% added if not detected unless 'nointercept' is specified.
%
% Key pointers:
% - Output structure (regression_results_ols) contains beta and t images in statistic_image objects
%
% - If contrasts are entered, output structure contains contrast images and contrast t images in statistic_image objects
%
% - T images are thresholded, but beta/contrast images are not. These can be re-thresholded with the threshold( ) method
%
% - statistic_image objects can be visuallized with any CANlab display methods, e.g., surface( ), montage( ), and orthviews( )
%
% - Use the get_wh_image( ) method to select specific beta, t, contrast, contrast t images from statistic_image objects
%
% - Enter condition names and contrast names as cell vectors of strings Labels are saved in the .image_labels property in output objects
%
% - Some diagnostics, including VIFs and leverages, are automatically checked and returned
%
% - There is an option for robust regression, specified using the 'robust' flag
%
% - For first-level models (only), the grandmeanscale option is recommended to increase homogeneity in scale across participants
%
% - The regress( ) method does not use covariates field of fmri_data(). You must include covariates
%    manually in dat.X.
%
% - The regress( ) method can also create a map of brain regions that predict the dat.Y vector using the 'brainony' option.
%   This is essentially a univariate version of the 'predict' command.  Warning: this is very slow as it loops
%   through all voxels.
%
% Creates thresholded plot by default
%
% :Usage:
% ::
%
%    out = regress(dat, varargin)
%
% :Inputs:
%  **dat:**
%        should be an fmri_data object with X field defined.
%        dat.X can be a design_matrix() object.
%
% :Optional Inputs:
%  **[threshold, 'unc']:**
%        p-value threshold string indicating threshold type
%        (see help statistic_image.threshold for options)
%
%  **robust:**
%        Run a robust regression (default is OLS).  Robust is considerably
%        slower than OLS
%
%  **grandmeanscale:**
%        Scale overall grand mean to a value of 100. Intended to reduce inter-subject
%        variability in 1st-level analysis (single-subject) when doing
%        multi-subject group analyses on resulting contrasts/beta images.
%        Assumes mask and overall brain size are consistent across replicates (e.g., participants, in a first-level analysis)
%        Do not use for 2nd-level or standard single-level analyses.
%
%  **C**
%       Followed by a contrast matrix, each column is a contrast across conditions/events
%       [k x c] matrix, where k = size(X, 2) and c is number of contrasts
%
%       Note: k must be the number including the intercept. If your input
%       model X does not include an intercept, contrast matrix must have
%       one more element than input X has rows, to account for the added
%       intercept.
%
%  **nointercept:**
%        Do not add intercept to model
%
%  **display, display_results:**
%        Show thresholded results usin orthviews (if < 10 regressors)
%
%  **nodisplay:**
%        Do not plot thresholded results using orthviews
%
%  **brainony:**
%        univariate approach to predict obj.Y from brain data
%
%  **residual:**
%        Output residual as fmri_data() object
%
%  **noverbose:**
%        Suppress verbose outputs
%
%  **variable_names:** (or 'names')
%       Followed by a cell array of variable/regressor names, for
%       non-intercept regressors
%
%  **contrast_names:**
%       Followed by a cell array of contrast names
%
%  **analysis_name:**
%       Followed by a string with a name/description for this analysis.
%
%
% :Outputs:
%
%  **out:**
%        A structure containing stats_img and fmri_data objects.
%        In addition to the main outputs below, the out structure also has
%        fields for input_parameters, the design matrix (X), variable
%        names, and warnings.
%
%  **out.b:**
%        stats_img object of beta values estimated from regression
%
%  **out.t:**
%        stats_img object of t-values with input threshold
%
%  **out.df:**
%        fmri_data object of degrees of freedom
%
%  **out.sigma:**
%        fmri_data object of variance of residual
%
%  **out.residual:**
%        fmri_data object of residual data after model has been regressed out (optional).
%
%  **out.diagnostics:***
%        A structure containing VIFs and leverage values for the design matrix
%        out.diagnostics.Variance_inflation_factors = VIFs
%        out.diagnostics.Leverages = leverage values
%
% :Examples:
% ::
%
%    % Run regression with liberal threshold
%    out = regress(dat, .05, 'unc');
%
%    % Run regression with conservative threshold and save residual
%    out = regress(dat, .001, 'unc', 'residual);
%
%    % Run robust regression with fdr threshold
%    out = regress(dat, .05, 'fdr','robust');
%
%    % Run a regression predicting behavior from brain at liberal threshold
%    out  = regress(data_comb, .05, 'unc', 'brainony')
%
%    % Re-threshold at different values
%    out.t = threshold(out.t, .05, 'fdr');
%    out.t = threshold(out.t, .001, 'unc');
%
%    % Re-display results of thresholding
%    orthviews(out.t);
%
%    % Write out beta image to current directory
%    out.b.fullpath = fullfile(pwd,'beta.nii');
%    write(out)
%
%    % Plot diagnostics
%    figure; subplot(1,2,1); title('VIFs')
%    plot(regression_results.diagnostics.Variance_inflation_factors);
%    subplot(1,2,2); title('Leverage of each observation')
%    plot(regression_results.diagnostics.Leverages);
%
%   % Run with options:
%   regression_results = regress(dat, 'variable_names', names, 'analysis_name', 'Pinel localizer 1st-level GLM', 'noverbose');
%
%   regression_results_ols = regress(dat, 'variable_names', names, ...
%     'analysis_name', 'Pinel localizer 1st-level GLM', ...
%     'C', C.weights, 'contrast_names', C.names);
%
% Show some results for selected individual conditions
% % Thresholded t map for the first condition
% t = get_wh_image(regression_results_ols.t, 1);
%
% create_figure('montage'); axis off
% display_obj = montage(t);
%
% figure; surface(t);
%
% % Montages for thresholded t images, condition [2 4 6 8]
%
% t = get_wh_image(regression_results_ols.t, [2 4 6 8]);
%
% create_figure('montage2'); axis off
% display_obj = montage(t);
% Show  results for contrasts
% t = get_wh_image(regression_results_ols.con_t, 1:3);
%
% create_figure('montage3'); axis off
% display_obj = montage(t);
%
% t = get_wh_image(regression_results_ols.con_t, 4:6);
%
% create_figure('montage4'); axis off
% display_obj = montage(t);

% ..
%    Copyright (c) 2015 Tor Wager & Luke Chang
%
%    Permission is hereby granted, free of charge, to any person obtaining a
%    copy of this software and associated documentation files (the "Software"),
%    to deal in the Software without restriction, including without limitation
%    the rights to use, copy, modify, merge, publish, distribute, sublicense,
%    and/or sell copies of the Software, and to permit persons to whom the
%    Software is furnished to do so, subject to the following conditions:
%
%    The above copyright notice and this permission notice shall be included
%    in all copies or substantial portions of the Software.
%
%    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
%    OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
%    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
%    THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
%    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
%    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
%    DEALINGS IN THE SOFTWARE.
% ..

% ..
%    Programmers' Notes:
%    c Tor Wager, Dec 2010
%    Edited by Luke Chang, 9/27/2012 to add optional input to reverse X & Y (i.e., create a map of voxels that predict the behavioral variable)
%    Edited by Luke Chang, 9/28/2012 to add optional input to run robust regression for brainony
%    Edited by Luke Chang, 10/24/2012 to save residuals (i.e., out.r), which is helpful for denoising an image
%    Edited by Luke Chang, 3/26/2013 to add optional input to not add an intercept - allows for more flexible modeling options
%    Code completely refactored by Luke Chang 2/24/25
%    Verbose option updated by Tor, 7/2015
%    Help, some outputs, robust % update: added by Tor, 5/2021
% ..

% ..
%    ---------------------------------------------------------------------
%    Defaults
%    ---------------------------------------------------------------------
% ..
inputargs = {.001, 'uncorrected'}; % default options for thresholding
do_display = false;
brain_is_outcome = 1; %else do brain on Y
do_robust = 0;
do_intercept = 1;
do_resid = 0;
doverbose = true;
variable_names = {};
contrast_names = {};
analysis_name = '';
C = [];
grandmeanscale = false;

% ---------------------------------------------------------------------
% Parse Inputs
% ---------------------------------------------------------------------
for varg = 1:length(varargin)

    if ischar(varargin{varg})
        % reserved keywords
        if strcmpi('unc',varargin{varg})
            inputargs = {varargin{varg-1}, 'uncorrected'};
            varargin{varg} = {}; varargin{varg - 1} = {};
        end
        if strcmpi('fdr',varargin{varg})
            inputargs = {varargin{varg-1}, 'fdr'};
            varargin{varg} = {}; varargin{varg - 1} = {};
        end

        if strcmpi('nodisplay',varargin{varg})
            do_display = 0;
            varargin{varg} = {};
        end

        if strcmpi('display',varargin{varg}) | strcmpi('display_results',varargin{varg})
            do_display = 0;
            varargin{varg} = {};
        end

        if strcmpi('robust',varargin{varg})
            do_robust = 1;
            varargin{varg} = {};
        end

        if strcmpi('brainony',varargin{varg}) | strcmpi('brain_is_predictor',varargin{varg}) %#ok<*OR2>
            brain_is_outcome = 0;
            varargin{varg} = {};
        end
        if strcmpi('nointercept',varargin{varg})
            do_intercept = 0;
            varargin{varg} = {};
        end
        if strcmpi('residual',varargin{varg})
            do_resid = 1;
            varargin{varg} = {};
        end
        if strcmpi('noverbose',varargin{varg})
            doverbose = false;
            varargin{varg} = {};
        end

        if strcmpi('names',varargin{varg}) | strcmpi('variable_names',varargin{varg})
            variable_names = varargin{varg + 1};
            varargin{varg} = {}; varargin{varg + 1} = {};
        end

        if strcmpi('C',varargin{varg}) | strcmpi('contrasts',varargin{varg})
            C = varargin{varg + 1};
            varargin{varg} = {}; varargin{varg + 1} = {};
        end

        if strcmpi('contrast_names',varargin{varg})
            contrast_names = varargin{varg + 1};
            varargin{varg} = {}; varargin{varg + 1} = {};
        end

        if strcmpi('analysis_name',varargin{varg})
            analysis_name = varargin{varg + 1};
            varargin{varg} = {}; varargin{varg + 1} = {};
        end

        if strcmpi('grandmeanscale',varargin{varg})
            grandmeanscale = true;
            varargin{varg} = {};
        end

    end % if ischar
end

if ~doverbose % add to pass into threshold( )
    inputargs{end+1} = 'noverbose';
end

% ---------------------------------------------------------------------
% Check Data and Diagnostics
% ---------------------------------------------------------------------

mywarnings = {};

% Check if fmri_data or image_vector
% ---------------------------------------------------------------------

if ~isa(dat,'fmri_data')
    error('dat input must be fmri_data object')
end

% Check data bit rate
% ---------------------------------------------------------------------

nuniquevals = length(unique(dat.dat(:)));

if nuniquevals < 2^10
    mywarnings{end+1} = sprintf('Number of unique values in dataset is low (%d, %3.2f bits), indicating possible restriction of bit rate. For comparison, Int16 has 65,536 unique values', nuniquevals, log2(nuniquevals));
end

% Check of dat.X exists and is correct format
% ---------------------------------------------------------------------

if brain_is_outcome
    if isempty(dat.X)
        error('Make sure you include a design matrix in dat.X')
    end
    if size(dat.dat, 2) ~= size(dat.X, 1)
        error('dat.dat must have same number of columns as dat.X has rows.')
    end
    if isa(dat.X,'design_matrix')
        dat.X = dat.X.dat;
    end
else % Check if dat.Y exists and is in correct format if running brainony
    if isempty(dat.Y)
        error('Make sure you include a vector in dat.Y.')
    end
    if size(dat.dat, 2) ~= size(dat.Y, 1)
        error('dat.dat must have same number of columns as dat.Y has rows.')
    end
end

% Check if Rank Deficient
% ---------------------------------------------------------------------

if rank(dat.X) < size(dat.X,2)
    mywarnings{end+1} = 'Warning:  dat.X is rank deficient.';
end

intercept_string = 'intercept is last';

% Check if Intercept is in model or specified for x_on_brain default
% ---------------------------------------------------------------------

if do_intercept && brain_is_outcome
    wh_int = intercept(dat.X, 'which');

    if isempty(wh_int)
        % add intercept and update wh_int (used later)
        if doverbose, mywarnings{end+1} = 'No intercept detected, adding intercept to last column of design matrix'; end
        X = intercept(dat.X, 'add');
        variable_names{end + 1} = 'Intercept';
%         wh_int = intercept(X, 'which');

    else
        intercept_string = sprintf('Intercept detected in column %1.0f of dat.X', wh_int);

        if doverbose, mywarnings{end+1} = intercept_string; end
        X = dat.X;

    end

else
    % No intercept or exogenous variables are outcome

    X = dat.X;
end

% Predictor centering
% ---------------------------------------------------------------------

m = mean(X);
wh_int = intercept(X, 'which');
m(wh_int) = [];
non_centered = abs(m) > 100 * eps;

% For effects-coded values [-1 1], ok, we want the intercept to reflect stats at average of 2 groups
iseffectcode = all(abs(X) == 1); % for each column
iseffectcode(wh_int) = [];

if any(non_centered & iseffectcode)

    mywarnings{end+1} = 'Warning:  Group sizes are unequal for effects-coded [1 -1] variable.';

end

if any(non_centered & ~iseffectcode)
    mywarnings{end+1} = 'Warning:  Predictors are not centered -- intercept is not interpretable as stats for average subject';
end

% Variance inflation
% ---------------------------------------------------------------------

vifs = getvif(X);

if any(vifs > 4)

    mywarnings{end+1} = 'Warning!!!  Design multicolinearity. Some regressors have variance inflation factors > 4. Check out.diagnostics';

end

% Leverages
% ---------------------------------------------------------------------

H = X*pinv(X);
%H = X*inv(X'*X)*X'  will be identical if not rank deficient
leverages = diag(H);

if any(abs(zscore(leverages)) >= 3)
     mywarnings{end+1} = 'Warning!!!  Some observations have high leverage values relative to others, regression may be unstable. abs(z(leverage)) > 3';
end

% Names
% ---------------------------------------------------------------------

k = size(X, 2);
if length(variable_names) < k
   if ~isempty(variable_names), mywarnings{end+1} = 'Warning!!!  Too few variable names entered, less than size(X, 2). Names may be inaccurate.'; end % suppress warning if NO names entered

   for i = length(variable_names)+1:k %#ok<*FXUP>
       variable_names{i} = sprintf('R%d', i); %#ok<*AGROW>
   end
end

if length(variable_names) > k
   mywarnings{end+1} = 'Warning!!!  Too many variable names entered, more than size(X, 2). Names may be inaccurate.';

   variable_names = variable_names(1:k);
end

% Check contrasts
% ---------------------------------------------------------------------

if ~isempty(C) && ~(size(C, 1) == size(X, 2))
    % Do this *after* adding intercept to X if needed

    disp('Contrasts entered, but size(C, 1) does not equal size(X, 2).');
    disp('Must have a contrast entry for each column of X (including the intercept)');
    error('Quitting.')

end

% Contrast names
% ---------------------------------------------------------------------

if ~isempty(C)

    kc = size(C, 2);

    if length(contrast_names) < kc
        if ~isempty(contrast_names), mywarnings{end+1} = 'Warning!!!  Too few contrast names entered, less than size(C, 2). Names may be inaccurate.'; end % suppress warning if NO names entered

        for i = length(contrast_names)+1:kc
            contrast_names{i} = sprintf('Con%d', i); %#ok<*AGROW>
        end
    end

    if length(contrast_names) > kc
        mywarnings{end+1} = 'Warning!!!  Too many contrast names entered, more than size(C, 2). Names may be inaccurate.';

        contrast_names = contrast_names(1:kc);
    end

end

% Enforce proper shape
if ~iscolumn(variable_names), variable_names = variable_names'; end
if ~isrow(contrast_names), contrast_names = contrast_names'; end

% Enforce valid names: Eliminate special characters and leading numbers
[variable_names, namewarnings] = format_text_letters_only(variable_names, 'numbers', 'cleanup', 'squeeze', 'underscore_ok');
mywarnings = [mywarnings namewarnings];
if ~isempty(namewarnings), mywarnings{end+1} = 'Enter valid variable_names'; end

if ~isempty(C)

    [contrast_names, namewarnings] = format_text_letters_only(contrast_names, 'numbers', 'cleanup', 'squeeze', 'underscore_ok');
    mywarnings = [mywarnings namewarnings];
    if ~isempty(namewarnings), mywarnings{end+1} = 'Enter valid contrast_names'; end

end

% Data scaling and format
% ---------------------------------------------------------------------

% Enforce double-format (just in case)
dat.dat = double(dat.dat);

if grandmeanscale

    % Scale grand mean to 100; approximates what SPM and other packages do. Assumes mask and overall brain size are consistent across replicates (e.g., participants, in a first-level analysis
    dat.dat = dat.dat .* 100 / nanmean(dat.dat(:));

end

if doverbose

    fprintf('Analysis: %s\n', analysis_name);
    disp('----------------------------------');

    nowarnings = all(cellfun(@isempty, mywarnings));

    disp('Design matrix warnings:');
    disp('----------------------------------');
    if nowarnings
        disp('None')

    else
        disp(char(mywarnings{:}))

    end

    disp(' ');
end

% ---------------------------------------------------------------------
% Run Regression
% ---------------------------------------------------------------------

tic
warning off

if brain_is_outcome
    % default is to regress dat.X on dat.dat (x on brain)
    % ---------------------------------------------------------------------

    % display info about regression
    if doverbose
        linestr = '______________________________________________________';
        disp(linestr);

        fprintf('Running regression: %3.0f voxels. Design: %3.0f obs, %3.0f regressors, %s\n', size(dat.dat, 1), size(X, 1), size(X, 2), intercept_string);
        if brain_is_outcome
            fprintf('\nPredicting exogenous variable(s) in dat.X using brain data as predictors, mass univariate');

        else % default
            fprintf('\nPredicting Brain Activity from dat.X, mass univariate');
        end
    end

    if do_robust
        % Robust - Regress stim/behavior(X) on brain (Y)
        %need to loop through voxels - Slow!

        if doverbose
            fprintf('\nRunning in Robust Mode ___%%');
        end

        v = size(dat.dat, 1);
        [n, k] = size(X);

        % Initialize outputs
        [b, t, stderr] = deal(zeros(k, v));
        p = ones(k, v);
        [dfe, sigma] = deal(zeros(1, v));

        for i = 1:v
            % For each voxel

            [bb,stats] = robustfit(X, dat.dat(i,:)', 'bisquare', [], 'off');

            b(:,i)=bb; %Betas
            t(:,i)=stats.t; %t-values
            p(:,i)=stats.p; %p-values
            dfe(:,i)=stats.dfe; %degrees of freedom
            stderr(:,i)=stats.se; %Standard Error
            sigma(:,i)=stats.robust_s; %robust estimate of sigma. LC not sure this is the best choice can switch to 'OLS_s','MAD_s', or 's'

            if doverbose && mod(i, 100) == 0
                fprintf('\b\b\b\b%03d%%', 100 * round(i/v))
            end
        end
        r = dat.dat' - X*b; %residual

    else
        % OLS - X predicting brain
        % - vectorized - Fast!

        if doverbose, fprintf('\nRunning in OLS Mode'); end

        % Estimate Betas in vector

        b = pinv(X) * dat.dat';

        % Error
        r = dat.dat' - X*b;

        % Residual variance
        [stderr, sigma] = get_std_errors(r, X);

        % Inference
        [t,dfe,p,sigma] = param_t_test(X,b,stderr,sigma);

    end

else
    % Regress brain (X) on stim/behavior (Y) - loops through voxels, slow!
    % ---------------------------------------------------------------------

    if doverbose
        % display info about regression
        fprintf('regression > X: %3.0f voxels. Y: %3.0f obs, %3.0f regressors, %s\n', size(dat.dat, 1), size(dat.Y, 2), intercept_string);
        fprintf('\nPredicting dat.Y from Brain Activity');
    end

    if do_robust %need to loop through voxels - Slow!
        if doverbose
            fprintf('\nRunning in Robust Mode ___%%');
        end

        v = size(dat.dat, 1);
        n = size(dat.dat, 2);
        k = 2; % one predictor (brain) and an intercept

        % Initialize outputs
        [b, t, stderr] = deal(zeros(k, v));
        p = ones(k, v);
        [dfe, sigma] = deal(zeros(1, v));

        for i = 1:v
            % Create X from brain Data
            if do_intercept
                X = intercept(dat.dat(i,:)','add');
            else
                X = dat.dat(i,:)';
            end

            [bb,stats] = robustfit(X, dat.Y, 'bisquare', [], 'off');

            b(:,i)=bb; %Betas
            t(:,i)=stats.t; %t-values
            p(:,i)=stats.p; %p-values
            dfe(:,i)=stats.dfe; %degrees of freedom
            stderr(:,i)=stats.se; %Standard Error
            sigma(:,i)=stats.robust_s; %robust estimate of sigma. LC not sure this is the best choice can switch to 'OLS_s','MAD_s', or 's'
            r(:,i) = dat.Y - X * b(:,i); %residual

            if doverbose && mod(i, 100) == 0
                fprintf('\b\b\b\b%03d%%', 100 * round(i/v))
            end

        end % voxel

    else
        % ---------------------------------------------------------------------
        %OLS -- Regress brain on Y

        if doverbose, fprintf('\nRunning in OLS Mode'); end

        for i = 1:size(dat.dat,1)

            % Create X from brain Data
            if do_intercept
                X = intercept(dat.dat(i,:)','add');
            else
                X = dat.dat(i,:)';
            end

            % Estimate Betas in vector
            b(:,i) = pinv(X) * dat.Y;

            % Error
            r(:,i) = dat.Y - X * b(:,i);

            %             sigma(i) = std(r(:,i)); % wrong
            %             stderr(:,i) = ( diag(inv(X' * X)) .^ .5 ) * sigma(i);  % params x voxels matrix of std. errors
        end

        % Residual variance - can do this at end because X is always the same size
        [stderr, sigma] = get_std_errors(r, X);

        % Inference
        [t,dfe,p,sigma] = param_t_test(X,b,stderr,sigma);

    end % dorobust

end % brain_is_outcome or not

stop = toc;
if doverbose, fprintf('\nModel run in %d minutes and %.2f seconds\n',floor(stop/60),rem(stop,60)); end

% ---------------------------------------------------------------------
% Contrasts
% ---------------------------------------------------------------------

if ~isempty(C)

    con_vals = C' * b;

    % Contrast STE
    con_ste = diag(C' * inv(X' * X) * C) .^ .5 * sigma;

    [con_t, ~, con_p] = param_t_test(X, con_vals, con_ste, sigma);

end

% ---------------------------------------------------------------------
% Create Output
% ---------------------------------------------------------------------

out = struct;

out.analysis_name = analysis_name;

out.input_parameters = struct( ...
    'brain_is_predictor', brain_is_outcome, 'do_robust', do_robust, 'grandmeanscale', grandmeanscale, ...
    'do_intercept', do_intercept, ...
    'do_resid', do_resid, 'doverbose', doverbose, 'do_display', do_display);

out.input_parameters.initial_statistical_threshold = inputargs;

out.input_image_metadata.source_notes = dat.source_notes;
out.input_image_metadata.history = dat.history;
out.input_image_metadata.image_names = dat.image_names;
out.input_image_metadata.fullpath = dat.fullpath;

% design, contrasts, and diagnostics

out.X = X;
out.variable_names = variable_names;
out.C = C;
out.contrast_names = contrast_names;

out.contrast_summary_table = table();

if ~isempty(C)
    % Contrast summary table

    for i = 1:size(C, 2)
        out.contrast_summary_table(:, i) = table(C(:, i));
    end

    out.contrast_summary_table.Properties.RowNames = variable_names;
    out.contrast_summary_table.Properties.VariableNames = contrast_names;

    if doverbose
        fprintf('\nSummary of conditions and contrasts\n%s\n', linestr);
        disp(out.contrast_summary_table);
    end

end

out.diagnostics = struct('Variance_inflation_factors', vifs, 'Leverages', leverages);
out.warnings = mywarnings;

% Create objects
if doverbose
    fprintf('\nCreating beta and t images, thresholding t images\n%s\n', linestr);
end

% Betas
out.b = statistic_image;
out.b.type = 'Beta';
out.b.p = p';
out.b.ste = stderr';
out.b.N = n;
out.b.dat = b';
out.b.dat_descrip = sprintf('Beta Values from regression, intercept is column %d', wh_int);
out.b.volInfo = dat.volInfo;
out.b.removed_voxels = dat.removed_voxels;
out.b.removed_images = false;  % this image does not have the same dims as the original dataset
out.b.image_labels = variable_names;

out.b = enforce_variable_types(out.b);

% - beta and contrast images are unthresholded, t images from both are thresholded
% if doverbose, fprintf('Thresholding b images at %3.6f %s\n', inputargs{1}, inputargs{2}); end
% out.b = threshold(out.b, inputargs{:}, 'noverbose'); % Threshold image

% T stats
out.t = statistic_image;
out.t.type = 'T';
out.t.p = p';
out.t.ste = stderr';
out.t.N = n;
out.t.dat = t';
out.t.dat_descrip = sprintf('t-values from regression, intercept is column %d', wh_int);
out.t.volInfo = dat.volInfo;
out.t.removed_voxels = dat.removed_voxels;
out.t.removed_images = false;  % this image does not have the same dims as the original dataset
out.t.image_labels = variable_names;

out.t = enforce_variable_types(out.t);

if doverbose
    fprintf('Thresholding t images at %3.6f %s\n', inputargs{1}, inputargs{2});
end
out.t = threshold(out.t, inputargs{:}); %Threshold image

% DF as fmri_data
out.df = dat;
out.df.dat = dfe';
out.df.dat_descrip = sprintf('Degrees of Freedom');

% Sigma as fmri_data
out.sigma = dat;
out.sigma.dat = sigma';
out.sigma.dat_descrip = sprintf('Sigma from Regression');

% Residual as fmri_data
if do_resid
    out.resid = dat;
    out.resid.dat = r';
    out.resid.dat_descrip = sprintf('Residual from Regression');
end

if ~isempty(C)

    if doverbose, fprintf('\nCreating contrast and t images and thresholding t images\n%s\n', linestr); end

    % Contrast values
    out.contrast_images = statistic_image;
    out.contrast_images.type = 'Contrast';
    out.contrast_images.p = con_p';
    out.contrast_images.ste = stderr';
    out.contrast_images.N = n;
    out.contrast_images.dat = con_vals';
    out.contrast_images.dat_descrip = 'Contrast Values from regression';
    out.contrast_images.volInfo = dat.volInfo;
    out.contrast_images.removed_voxels = dat.removed_voxels;
    out.contrast_images.removed_images = false;  % this image does not have the same dims as the original dataset
    out.contrast_images.image_labels = contrast_names;

    out.contrast_images = enforce_variable_types(out.contrast_images);

%     if doverbose, fprintf('Thresholding contrast images at %3.6f %s\n', inputargs{1}, inputargs{2}); end
%     out.contrast_images = threshold(out.contrast_images, inputargs{:}, 'noverbose'); % Threshold image

    % T stats
    out.con_t = statistic_image;
    out.con_t.type = 'T';
    out.con_t.p = con_p';
    out.con_t.ste = con_ste';
    out.con_t.N = n;
    out.con_t.dat = con_t';
    out.con_t.dat_descrip = sprintf('t-values from regression, intercept is column %d', wh_int);
    out.con_t.volInfo = dat.volInfo;
    out.con_t.removed_voxels = dat.removed_voxels;
    out.con_t.removed_images = false;  % this image does not have the same dims as the original dataset
    out.con_t.image_labels = contrast_names;

    out.con_t = enforce_variable_types(out.con_t);

    if doverbose
        fprintf('Thresholding t images at %3.6f %s\n', inputargs{1}, inputargs{2});

    end
    out.con_t = threshold(out.con_t, inputargs{:}); %Threshold image

end

% ---------------------------------------------------------------------
% Plot Results
% ---------------------------------------------------------------------
if k < 10 && do_display

    orthviews(out.t);

elseif do_display && doverbose
    disp('Warning: No display because >= 10 images.');

end

warning on

% ---------------------------------------------------------------------
% Subfunctions
% ---------------------------------------------------------------------

    function [stderr, sigma] = get_std_errors(r, X)

        [n, k] = size(X);
        v = size(dat.dat, 1);

        % We want diag(r' * r), but matrix size can be large
        % std(r) does not account for k parameters used, so is incorrect
        % residual std. not sqrt(var(resid)) -- we must account for k params used
        for i = 1:v
            sigma(1, i) = (r(:, i)' * r(:, i) ./ (n - k)) .^ .5;
        end

        stderr = ( diag(inv(X' * X)) .^ .5 ) * sigma;  % params x voxels matrix of std. errors

    end

    function [t,dfe,p,sigma] = param_t_test(X,b,stderr,sigma)
        % test whether parameter is significantly different from zero
        %
        % Inputs:
        % X:        design matrix
        % b:        beta values
        % stderr:   standard error of beta estimate
        % sigma:    standard deviation of residual
        %
        % Returns:
        % t:        t-value
        % dfe:      degrees of freedom
        % p:        p-value
        % sigma:    standard deviation of residual

        [n, k] = size(X);
        t = b ./ stderr;
        dfe = n - k;
        p = 2 * (1 - tcdf(abs(t), dfe));

        sigma(sigma == 0) = Inf;
        t(isnan(t)) = 0;
        p(isnan(p)) = 0;
        dfe = repmat(dfe,1,size(t,2));
    end

end % Main Function


text:
%  anova2rm_cell() - compute F-values in cell array using repeated measure
%                    ANOVA.
%
% Usage:
%    >> [FC FR FI dfc dfr dfi] = anova2rm_cell( data );
%
% Inputs:
%   data       = data consisting of PAIRED arrays to be compared. The last
%                dimension of the data array is used to compute ANOVA.
% Outputs:
%   FC   - F-value for columns.
%   FR   - F-value for rows.
%   FI   - F-value for interaction.
%   dfc  - degree of freedom for columns.
%   dfr  - degree of freedom for rows.
%   dfi  - degree of freedom for interaction.
%
% Note: this function is inspired from rm_anova available at
%       http://www.mathworks.se/matlabcentral/fileexchange/6874-two-way-rep
%       eated-measures-anova
%       It allows for fast computation of about 20 thousands ANOVA per
%       second. It is different from anova2_cell which mimics the ANOVA
%       fonction from the Matlab statistical toolbox. This function
%       computes true repeated measure ANOVA.
%
% Example:
%   a = { rand(1,10) rand(1,10) rand(1,10); rand(1,10) rand(1,10) rand(1,10) }
%   [FC FR FI dfc dfr dfi] = anova2rm_cell(a)
%   signifC = 1-fcdf(FC, dfc(1), dfc(2))
%   signifR = 1-fcdf(FR, dfr(1), dfr(2))
%   signifI = 1-fcdf(FI, dfi(1), dfi(2))
%
%   % for comparison
%   z = zeros(10,1); o = ones(10,1); t = ones(10,1)*2;
%   rm_anova2(  [ a{1,1}';a{1,2}';a{1,3}';a{2,1}';a{2,2}';a{2,3}' ], ...
%               repmat([1:10]', [6 1]), [o;o;o;z;z;z], [z;o;t;z;o;t], {'a','b'})
%
%   c = { rand(200,400,10) rand(200,400,10); ...
%         rand(200,400,10) rand(200,400,10)};
%   [FC FR FI dfc dfr dfi] = anova2rm_cell(c) % computes 200x400 ANOVAs
%
% Author: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 2010

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

function [fA fB fAB dfApair dfBpair dfABpair] = anova2rm_cell(data)

% compute all means and all std
% -----------------------------
a = size(data,1);
b = size(data,2);
nd = myndims( data{1} );
n  = size( data{1} ,nd);

% only for paired stats
% ---------------------
if nd == 1
    AB = zeros(a,b,'single');
    AS = zeros(a,n,'single');
    BS = zeros(b,n,'single');
    sq = single(0);
    for ind1 = 1:a
        for ind2 = 1:b
            AB(ind1,ind2) = sum(data{ind1,ind2});
            AS(ind1,:)    = AS(ind1,:) + data{ind1,ind2}';
            BS(ind2,:)    = BS(ind2,:) + data{ind1,ind2}';
            sq            = sq + sum(data{ind1,ind2}.^2);
        end;
    end;
    dimA = 2;
    dimB = 1;
elseif nd == 2
    AB = zeros(size(data{1},1),a,b,'single');
    AS = zeros(size(data{1},1),a,n,'single');
    BS = zeros(size(data{1},1),b,n,'single');
    sq = zeros(size(data{1},1),1,'single');
    for ind1 = 1:a
        for ind2 = 1:b
            AB(:,ind1,ind2) = sum(data{ind1,ind2},nd);
            AS(:,ind1,:)    = AS(:,ind1,:) + reshape(data{ind1,ind2},size(data{1},1),1,n);
            BS(:,ind2,:)    = BS(:,ind2,:) + reshape(data{ind1,ind2},size(data{1},1),1,n);
            sq              = sq + sum(data{ind1,ind2}.^2,nd);
        end;
    end;
    dimA = 3;
    dimB = 2;
elseif nd == 3
    AB = zeros(size(data{1},1),size(data{1},2),a,b,'single');
    AS = zeros(size(data{1},1),size(data{1},2),a,n,'single');
    BS = zeros(size(data{1},1),size(data{1},2),b,n,'single');
    sq = zeros(size(data{1},1),size(data{1},2),'single');
    for ind1 = 1:a
        for ind2 = 1:b
            AB(:,:,ind1,ind2) = sum(data{ind1,ind2},nd);
            AS(:,:,ind1,:)    = AS(:,:,ind1,:) + reshape(data{ind1,ind2},size(data{1},1),size(data{1},2),1,n);
            BS(:,:,ind2,:)    = BS(:,:,ind2,:) + reshape(data{ind1,ind2},size(data{1},1),size(data{1},2),1,n);
            sq                = sq + sum(data{ind1,ind2}.^2,nd);
        end;
    end;
    dimA = 4;
    dimB = 3;
elseif nd == 4
    AB = zeros(size(data{1},1),size(data{1},2),size(data{1},3),a,b,'single');
    AS = zeros(size(data{1},1),size(data{1},2),size(data{1},3),a,n,'single');
    BS = zeros(size(data{1},1),size(data{1},2),size(data{1},3),b,n,'single');
    sq = zeros(size(data{1},1),size(data{1},2),size(data{1},3),'single');
    for ind1 = 1:a
        for ind2 = 1:b
            AB(:,:,:,ind1,ind2) = sum(data{ind1,ind2},nd);
            AS(:,:,:,ind1,:)    = AS(:,:,:,ind1,:) + reshape(data{ind1,ind2},size(data{1},1),size(data{1},2),size(data{1},3),1,n);
            BS(:,:,:,ind2,:)    = BS(:,:,:,ind2,:) + reshape(data{ind1,ind2},size(data{1},1),size(data{1},2),size(data{1},3),1,n);
            sq                = sq + sum(data{ind1,ind2}.^2,nd);
        end;
    end;
    dimA = 5;
    dimB = 4;
end;

A = sum(AB,dimA); % sum across columns, so result is ax1 column vector
B = sum(AB,dimB); % sum across rows, so result is 1xb row vector
S = sum(AS,dimB); % sum across columns, so result is 1xs row vector
T = sum(sum(A,dimB),dimA); % could sum either A or B or S, choice is arbitrary

% degrees of freedom
dfA = a-1;
dfB = b-1;
dfAB = (a-1)*(b-1);
dfS = n-1;
dfAS = (a-1)*(n-1);
dfBS = (b-1)*(n-1);
dfABS = (a-1)*(b-1)*(n-1);

% bracket terms (expected value)
expA  = sum(A.^2,dimB)./(b*n);
expB  = sum(B.^2,dimA)./(a*n);
expAB = sum(sum(AB.^2,dimA),dimB)./n;
expS  = sum(S.^2,dimA)./(a*b);
expAS = sum(sum(AS.^2,dimB),dimA)./b;
expBS = sum(sum(BS.^2,dimB),dimA)./a;
expY  = sq; %sum(Y.^2);
expT  = T.^2 / (a*b*n);

% sums of squares
ssA   = expA - expT;
ssB   = expB - expT;
ssAB  = expAB - expA - expB + expT;
ssS   = expS - expT;
ssAS  = expAS - expA - expS + expT;
ssBS  = expBS - expB - expS + expT;
ssABS = expY - expAB - expAS - expBS + expA + expB + expS - expT;
ssTot = expY - expT;

% mean squares
msA   = ssA / dfA;
msB   = ssB / dfB;
msAB  = ssAB / dfAB;
msS   = ssS / dfS;
msAS  = ssAS / dfAS;
msBS  = ssBS / dfBS;
msABS = ssABS / dfABS;

% f statistic
fA = msA ./ msAS;
fB = msB ./ msBS;
fAB = msAB ./ msABS;
dfApair  = [dfA dfAS];
dfBpair  = [dfB dfBS];
dfABpair = [dfAB dfABS];

function val = myndims(a)
    if ndims(a) > 2
        val = ndims(a);
    else
        if size(a,1) == 1,
            val = 2;
        elseif size(a,2) == 1,
            val = 1;
        else
            val = 2;
        end;
    end;

author: PatternRecognition

repo: OpenBMI


text:
function meshnd_quality (do_metis)
%MESHND_QUALITY test the ordering quality computed by meshnd.
% The fill-in and flop count for sparse Cholesky factorization using the meshnd
% nested dissection ordering is computed with AMD.  If SuiteSparse is installed
% with METIS, and if requested, then the metis nested dissection ordering is
% also compared.
%
% Example:
%   meshnd_quality          % compare MESHND and AMD
%   meshnd_quality (1)      % also compare with METIS
%
% See also meshnd, meshsparse, nested, amd, metis.

% Copyright 2007, Timothy A. Davis, Univ. of Florida

stencils = [5 9 7 27] ;

if (nargin < 1)
    do_metis = 0 ;
end
if (do_metis)
    if (exist ('metis') ~= 3)                                               %#ok
        % METIS not installed
        do_metis = 0 ;
    end
end

figure (1)
clf

for sk = 1:4

    stencil = stencils (sk) ;

    is3D = (stencil == 7 | stencil == 27) ;     %#ok
    if (is3D)
	s = 2.^(3:.1:7) ;		% mesh size up to 127-by-127-by-127
    else
	s = 2.^(3:.1:10) - 1 ;		% mesh size up to 1023-by-1023
    end
    t = length (s) ;
    lnz = nan * zeros (3,t) ;
    fl  = nan * zeros (3,t) ;

    try

        for t = 1:length (s)

            n = floor (s (t)) ;

            % create the mesh and the matrix, and get nested dissection ordering
            if (is3D)
                fprintf ('3D mesh: %d-by-%d-by-%d, %d-point stencil\n', ...
                    n, n, n, stencil) ;
                [G p] = meshnd (n, n, n) ;
            else
                fprintf ('2D mesh: %d-by-%d, %d-point stencil\n', n, n,stencil);
                [G p] = meshnd (n, n) ;
            end
            A = meshsparse (G, stencil) ;

            % ND results
            c = symbfact (A (p,p)) ;
            lnz (1,t) = sum (c) ;
            fl  (1,t) = sum (c.^2) ;
            fprintf ('    MESHND:            nnz(L) %8.3e  flops %8.3e\n', ...
                lnz (1,t), fl (1,t)) ;
            clear G

            % AMD results
            try
                p = amd (A) ;
            catch
                % assume SuiteSparse is installed
                p = amd2 (A) ;
            end
            c = symbfact (A (p,p)) ;
            lnz (2,t) = sum (c) ;
            fl  (2,t) = sum (c.^2) ;
            fprintf ('    AMD:               nnz(L) %8.3e  flops %8.3e\n', ...
                lnz (2,t), fl (2,t)) ;

            % METIS results (requires SuiteSparse and METIS)
            if (do_metis)
                p = metis (A) ;
                c = symbfact (A (p,p)) ;
                lnz (3,t) = sum (c) ;
                fl  (3,t) = sum (c.^2) ;
                fprintf (...
                    '    METIS:             nnz(L) %8.3e  flops %8.3e\n', ...
                    lnz (3,t), fl (3,t)) ;
            end

            % plot the relative nnz(L) results
            subplot (2, 4, 2*sk - 1) ;
            loglog (s (1:t), lnz (2,1:t) ./ lnz (1,1:t), 'b-') ;
            hold on
            if (do_metis)
                loglog (s (1:t), lnz (3,1:t) ./ lnz (1,1:t), 'r-') ;
            end
            loglog (s (1:t), ones (1,t), 'k-') ;
            if (do_metis)
                ylabel ('nnz(L) for AMD or METIS / nnz(L) for meshnd') ;
                legend ('AMD', 'METIS') ;
            else
                ylabel ('nnz(L) for AMD / nnz(L) for meshnd') ;
            end
            xlabel ('mesh size') ;
            axis ([min(s) max(s) .1 10]) ;
            set (gca, 'YTick', [.1 .25 .5 .8 1 1.25 2 4 10]) ;
            if (is3D)
                set (gca, 'XTick', [1 10 100]) ;
                title (sprintf ('3D mesh, %d-point stencil', stencil)) ;
            else
                set (gca, 'XTick', [1 10 100 1000]) ;
                title (sprintf ('2D mesh, %d-point stencil', stencil)) ;
            end

            % plot the relative flop results
            subplot (2, 4, 2*sk) ;
            loglog (s (1:t), fl (2,1:t) ./ fl (1,1:t), 'b-') ;
            hold on
            if (do_metis)
                loglog (s (1:t), fl (3,1:t) ./ fl (1,1:t), 'r-') ;
            end
            loglog (s (1:t), ones (1,t), 'k-') ;
            ylabel ('flops for AMD or METIS / flops for meshnd') ;
            if (do_metis)
                ylabel ('flops for AMD or METIS / flops for meshnd') ;
                legend ('AMD', 'METIS') ;
            else
                ylabel ('nnz(L) for AMD / nnz(L) for meshnd') ;
            end
            xlabel ('mesh size') ;
            axis ([min(s) max(s) .1 10]) ;
            set (gca, 'YTick', [.1 .25 .5 .8 1 1.25 2 4 10]) ;
            if (is3D)
                set (gca, 'XTick', [1 10 100]) ;
                title (sprintf ('3D mesh, %d-point stencil', stencil)) ;
            else
                set (gca, 'XTick', [1 10 100 1000]) ;
                title (sprintf ('2D mesh, %d-point stencil', stencil)) ;
            end

            drawnow

        end

    catch
        % out-of-memory is OK, other errors are not
        disp (lasterr) ;
        if (isempty (strfind (lasterr, 'Out of memory')))
            error (lasterr) ;                                               %#ok
        else
            fprintf ('test terminated early, but otherwise OK\n') ;
        end
    end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function MeshStat = fem_meshstats(FemFile)
% FEM_MESHSTATS: Computes and display FEM mesh volume (tetrahedral only).
%
% INPUTS:
%    - FemFile : Relative file path to a Braistorm tetrahedral FEM mesh
% OUTPUTS:
%    - MeshStat : Mtlab structure that contains all FEM mesh stqtistics (edge length, elem volum and mesh quality)
%
% DEPENDENCIES:
%    This function require the iso2mesh toolbox

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Takfarinas Medani, Francois Tadel, 2023

% Install/load iso2mesh plugin
isInteractive = 1;
[isInstalled, errInstall] = bst_plugin('Install', 'iso2mesh', isInteractive);
if ~isInstalled
    error('Plugin "iso2mesh" not available.');
end

% Get data in database
bst_progress('start', 'Mesh volume', 'Loading file...');
FemFullFile = file_fullpath(FemFile);
femmat = load(FemFullFile);
% Check type of mesh: accept only tetrahedral
if (size(femmat.Elements,2) ~= 4)
    error('This menu is available for tetrahedral meshes only.');
end
% Display results in figures if no variable in output
isDisplay = (nargout == 0);
hFig = [];

% Convert to millimeter for convenience
femmat.Vertices = 1000 .* femmat.Vertices;

% Loop over the tissues
TissueID = unique(femmat.Tissue);
for iTissue = 1:length(TissueID)
    iTissueID = find(femmat.Tissue == TissueID(iTissue));

    % 1. Edges length
    bst_progress('text', sprintf('Computing edges length...  [%d/%d]', iTissue, length(TissueID)));
    Edges = meshedge(femmat.Elements(iTissueID,:));
    n1 = femmat.Vertices(Edges(:,1),:);
    n2 = femmat.Vertices(Edges(:,2),:);
    EdgeLength = sqrt((n1(:,1)- n2(:,1)).^2 + (n1(:,2)- n2(:,2)).^2 + (n1(:,3)- n2(:,3)).^2);

    tstat.EdgeLengthMax = max(EdgeLength);
    tstat.EdgeLengthMin = min(EdgeLength);
    tstat.EdgeLengthStd = std(EdgeLength);
    tstat.EdgeLengthMean = mean(EdgeLength);
    tstat.EdgeLengthRMS = rms(EdgeLength);

    % 2. Mesh quality: Joe-Liu mesh quality metric (0-1)
    %  quality: a vector of the same length as size(elem,1), with
    %            each element being the Joe-Liu mesh quality metric (0-1) of
    %            the corresponding element. A value close to 1 represents
    %            higher mesh quality (1 means equilateral tetrahedron);
    %            a value close to 0 means nearly degenerated element.
    bst_progress('text', sprintf('Computing mesh quality...  [%d/%d]', iTissue, length(TissueID)));
    quality = 100 .*meshquality(femmat.Vertices, femmat.Elements(iTissueID,:));
    tstat.MeshQualityMax = max(quality);
    tstat.MeshQualityMin = min(quality);
    tstat.MeshQualityStd = std(quality);
    tstat.MeshQualityMean = mean(quality);

    % 3. Volume of elem
    bst_progress('text', sprintf('Computing volume of elements...  [%d/%d]', iTissue, length(TissueID)));
    voli = elemvolume(femmat.Vertices, femmat.Elements(iTissueID,:));
    tstat.MeshVolumeMax = max(voli);
    tstat.MeshVolumeMin = min(voli);
    tstat.MeshVolumeStd = std(voli);
    tstat.MeshVolumeMean = mean(voli);
    tstat.MeshVolumeSum = sum(voli);

    MeshStat.(femmat.TissueLabels{iTissue}) = tstat;

    % Visualization
    if isDisplay
        bst_progress('text', sprintf('Visualisation... [%d/%d]', iTissue, length(TissueID)));
        hFig(end+1) = figure('Name', ['Mesh stat: ' femmat.TissueLabels{iTissue}], 'NumberTitle', 'off');

        nbins = 30;
        subplot(3,1,1)
        histogram(EdgeLength,nbins);
        xlabel(sprintf('Edge length (mm):   mean=%1.2f | std=%1.2f | min=%1.2f | max=%1.2f', tstat.EdgeLengthMean, tstat.EdgeLengthStd, tstat.EdgeLengthMin, tstat.EdgeLengthMax))
        drawnow

        subplot(3,1,2)
        histogram(quality,nbins);
        xlabel(sprintf('Mesh quality (%%):   mean=%1.2f | std=%1.2f | min=%1.2f | max=%1.2f', tstat.MeshQualityMean, tstat.MeshQualityStd, tstat.MeshQualityMin, tstat.MeshQualityMax))
        drawnow

        subplot(3,1,3)
        histogram(voli,nbins);
        xlabel(sprintf('Element volume (mm3):   mean=%1.2f | std=%1.2f | min=%1.2f | max=%1.2f | sum=%1.2f', tstat.MeshVolumeMean, tstat.MeshVolumeStd, tstat.MeshVolumeMin, tstat.MeshVolumeMax, tstat.MeshVolumeSum))
        drawnow
    end
end

% For all the full Model
% 1. Edges length
bst_progress('text', 'Computing edges length...');
Edges = meshedge(femmat.Elements);
n1 = femmat.Vertices(Edges(:,1),:);
n2 = femmat.Vertices(Edges(:,2),:);
EdgeLength = sqrt((n1(:,1)- n2(:,1)).^2 + (n1(:,2)- n2(:,2)).^2 + (n1(:,3)- n2(:,3)).^2);

MeshStat.FullModel.EdgeLengthMax = max(EdgeLength);
MeshStat.FullModel.EdgeLengthMin = min(EdgeLength);
MeshStat.FullModel.EdgeLengthStd = std(EdgeLength);
MeshStat.FullModel.EdgeLengthMean = mean(EdgeLength);
MeshStat.FullModel.EdgeLengthRMS = rms(EdgeLength);

% 2. Mesh quality: Joe-Liu mesh quality metric (0-100)
bst_progress('text', 'Computing mesh quality...');
quality = 100 .* meshquality(femmat.Vertices, femmat.Elements);
MeshStat.FullModel.MeshQualityMax = max(quality);
MeshStat.FullModel.MeshQualityMin = min(quality);
MeshStat.FullModel.MeshQualityStd = std(quality);
MeshStat.FullModel.MeshQualityMean = mean(quality);

% 3. Volume of elem
bst_progress('text', 'Computing volume of elements...');
voli = elemvolume(femmat.Vertices, femmat.Elements);
MeshStat.FullModel.MeshVolumeMax = max(voli);
MeshStat.FullModel.MeshVolumeMin = min(voli);
MeshStat.FullModel.MeshVolumeStd = std(voli);
MeshStat.FullModel.MeshVolumeMean = mean(voli);
MeshStat.FullModel.MeshVolumeSum = sum(voli);

% Visualization
if isDisplay
    bst_progress('text', 'Visualisation...');
    hFig(end+1) = figure('Name', 'Mesh stat: all tissues combined', 'NumberTitle', 'off');

    nbins = 30;
    subplot(3,1,1)
    histogram(EdgeLength,nbins);
    xlabel(sprintf('Edge length (mm):   mean=%1.2f | std=%1.2f | min=%1.2f | max=%1.2f', MeshStat.FullModel.EdgeLengthMean, MeshStat.FullModel.EdgeLengthStd, MeshStat.FullModel.EdgeLengthMin, MeshStat.FullModel.EdgeLengthMax))
    drawnow

    subplot(3,1,2)
    histogram(quality,nbins);
    xlabel(sprintf('Mesh quality (%%):   mean=%1.2f | std=%1.2f | min=%1.2f | max=%1.2f', MeshStat.FullModel.MeshQualityMean, MeshStat.FullModel.MeshQualityStd, MeshStat.FullModel.MeshQualityMin, MeshStat.FullModel.MeshQualityMax))
    drawnow

    subplot(3,1,3)
    histogram(voli,nbins);
    xlabel(sprintf('Element volume (mm3):   mean=%1.2f | std=%1.2f | min=%1.2f | max=%1.2f | sum=%1.2f', MeshStat.FullModel.MeshVolumeMean, MeshStat.FullModel.MeshVolumeStd, MeshStat.FullModel.MeshVolumeMin, MeshStat.FullModel.MeshVolumeMax, MeshStat.FullModel.MeshVolumeSum))

    % Close all the figures at once
    set(hFig, 'DeleteFcn', @(h,ev)delete(setdiff(hFig,h)));
end

bst_progress('stop');

author: brainstorm-tools


text:
function createNonemptySolution_FuncGen(OCP)
    dim = OCP.dim;
    N = OCP.N;
    fileID = fopen('./funcgen/createNonemptySolution.m','w');

    fprintf(fileID, 'function solution = createNonemptySolution()\n');
    fprintf(fileID, '   solution.lambda = zeros(%d,%d);\n',dim.lambda,N);
    if dim.mu ~= 0
        fprintf(fileID, '   solution.mu = zeros(%d,%d);\n',dim.mu,N);
    end
    fprintf(fileID, '   solution.u = zeros(%d,%d);\n',dim.u,N);
    fprintf(fileID, '   solution.x = zeros(%d,%d);\n',dim.x,N);
    if dim.z ~= 0
        fprintf(fileID, '   solution.z = zeros(%d,%d);\n',dim.z,N);
    end
    fprintf(fileID, '   solution.LAMBDA = zeros(%d,%d,%d);\n',dim.x,dim.x,N);
    fprintf(fileID, 'end');

    fclose(fileID);

end

author: deng-haoyang

repo: ParNMPC


text:
function sin_test ( )

%*****************************************************************************80
%
%% SIN_TEST tests R4_SIN and R8_SIN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../test_values' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SIN_TEST:\n' );
  fprintf ( 1, '  Test SIN_VALUES, R4_SIN, R8_SIN.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '             X         SIN(X)\n' );
  fprintf ( 1, '                    R4_SIN(X)         Diff\n' );
  fprintf ( 1, '                    R8_SIN(X)         Diff\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx1 ] = sin_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = r4_sin ( single ( x ) );
    fx3 = r8_sin ( x );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %14.4f  %14.6g\n', x, fx1 );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx2, abs ( fx1 - fx2 ) );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx3, abs ( fx1 - fx3 ) );

  end

  rmpath ( '../test_values' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [strikb,dipb,rakeb,dipdib,ierr] = focal_pl2pl(strika,dipa,rakea)

    % compute strike, dip and rake of a nodal plane from strike, dip and rake of the other one
    %
    %
    %     usage:
    %     call pl2pl(strika,dipa,rakea,strikb,dipb,rakeb,dipdib,ierr)
    %
    %     arguments:
    %     strika         strike angle in degrees of the first nodal plane (INPUT)
    %     dipa           dip angle in degrees of the first nodal plane (INPUT)
    %     rakea          rake angle in degrees of the first nodal plane (INPUT)
    %     strikb         strike angle in degrees of the second nodal plane (OUTPUT)
    %     dipb           dip angle in degrees of the second nodal plane (OUTPUT)
    %     rakeb          rake angle in degrees of the second nodal plane (OUTPUT)
    %     dipdib         dip direction in degrees of the second nodal plane (OUTPUT)
    %     ierr           error indicator (OUTPUT)
    %
    %     errors:
    %     1              input STRIKE angle out of range
    %     2              input DIP angle out of range
    %     4              input RAKE angle out of range
    %     3              1+2
    %     5              1+4
    %     7              1+2+4
    %     8              internal error
    %
    %      implicit none
    %-------------------------------------------------------------------------------
    %%       integer io
    %%       real amistr,amastr,amidip,amadip,amirak,amarak,amitre,amatre
    %%      1,amiplu,amaplu,orttol,ovrtol,tentol,dtor,c360,c90,c0,c1,c2,c3
    %%       common /fpscom/amistr,amastr,amidip,amadip,amirak,amarak,amitre
    %%      1,amatre,amiplu,amaplu,orttol,ovrtol,tentol,dtor,c360,c90,c0,c1,c2
    %%      2,c3,io
    %-------------------------------------------------------------------------------
    %%       real strika,dipa,rakea,anx,any,anz,dx,dy,dz,strikb,dipb,rakeb,
    %%      1dipdib
    %%       integer ierr
    %
    %%      call fpsset
    amistr=-360.;
    amastr=360.;
    amidip=0.;
    amadip=90.;
    amirak=-360.;
    amarak=360.;
    amitre=-360.;
    amatre=360.;
    amiplu=0.;
    amaplu=90.;
    orttol=2.;
    ovrtol=0.001;
    tentol=0.0001;
    dtor=0.017453292519943296;
    c360=360.;
    c90=90.;
    c0=0.;
    c1=1.;
    c2=2.;
    c3=3.;
    io=6;
    ifl=1;

    [anx,any,anz,dx,dy,dz,ierr] = focal_pl2nd(strika,dipa,rakea);
    if (ierr ~= 0)
        disp(['PL2PL: ierr = ' num2str(ierr)]);
        return;
    end
    [strikb,dipb,rakeb,dipdib,ierr] = focal_nd2pl(dx,dy,dz,anx,any,anz);
    if (ierr ~= 0)
        ierr = 8;
        disp(['PL2PL: ierr = ' num2str(ierr)]);
    end
end

author: CelsoReyes


text: %SIGN Datafile overload

author: marianux

repo: ecg-kit

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


text:
function [execError] = getExecutionErrors()
%getExecutionErrors Summary of this function goes here
%   Detailed explanation goes here
    global script_execution_errors;
    execError = script_execution_errors;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [xMin,fMin,exitCode]=quasiNetwonBFGS(f,x0,D0,epsilon,deltaTestDist,delta,lineSearchParams,scaleD,maxIter)
%%QUASINEWTONBFGS Perform unconstrained nonlinear optimization using the
%                 Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm. The
%                 algorithm performs unconstrained minimization of a
%                 nonlinear function without one having to provide a
%                 Hessian matrix. Note that if a minimization over a
%                 least-squares problem is desired, then the
%                 Levenberg-Marquardt algorithm in LSEstLMarquardt is
%                 often preferable. The limited memory version of thie
%                 algorithm, quasiNewtonLBFGS, is more appropriate for use
%                 with very large matrices. The L-BFGS algorithm also
%                 supports an additional L1 norm term in the objective
%                 function. If one wishes to zero a vector (not a scalar),
%                 then NewtonsMethod is more appropriate as this function
%                 assumes the Hessian matrix is symmetric.
%
%INPUTS: f A handle to the function (and its gradient) over which the
%          minimization is to be performed. The function [fVal,gVal]=f(x)
%          takes the NX1 x vector and returns the real scalar function
%          value fVal and gradient gVal at the point x.
%       x0 The NX1-dimensional point from which the minimization starts.
%       D0 An estimate of the inverse Hessian matrix at x0. If omitted or
%          an empty matrix is passed, then the identity matrix is used.
%  epsilon The parameter determining the accuracy of the desired solution
%          in terms of the gradient. The function terminates when
%          norm(g) < epsilon*max([1, norm(x)])
%          where g is the gradient. The default if omitted or an empty
%          matrix is passed is 1e-6.
% deltaTestDist The number of iterations back to use to compute the
%          decrease of the objective function if a delta-based convergence
%          test is performed. If zero, then no delta-based convergence
%          testing is done. The default if omitted or an empty matrix is
%          passed is zero.
%    delta The delta for the delta convergence test. This determines the
%          minimum rate of decrease of the objective function. Convergence
%          is determined if (f'-f)<=delta*f, where f' is the value of the
%          objective function f deltaTestDist iterations ago,and f is the
%          current objective function value. The default if this parameter
%          is omitted or an empty matrix is passed is 0.
% lineSearchParams An optional structure whose members specify tolerances
%          for the line search. The parameters are described as in the
%          lineSearch function.
%   scaleD A boolean parameter indicating whether the inverse Hessian
%          estimate should be scaled as in Equation 1.201 of Section 1.7 of
%          [1]. The default if omitted or an empty matrix is passed is
%          false.
%  maxIter The maximum number of iterations to use for the overall BFGS
%          algorithm. The default if this parameter is omitted or an empty
%          matrix is passed is 1000.
%
%OUTPUTS: xMin The value of x at the minimum point found. If exitCode is
%              negative, then this might be an empty matrix.
%         fMin The cost function value at the minimum point found. If
%              exitCode is negative, then this might be an empty
%              matrix.
%     exitCode A value indicating the termination condition of the
%              algorithm. Nonnegative values indicate success; negative
%              values indicate some type of failure. Possible values are:
%                  0 The algorithm termiated successfully based on the
%                    gradient criterion.
%                  1 The algorithm terminated successfully based on the
%                    accuracy criterion.
%                 -1 The maximum number of overall iterations was reached.
%              Other negative values correspond to a failure in lineSearch
%              and correspond to the exitCode returned by the lineSearch
%              function.
%
%The algorithm is implemented based on the description in Chapter 1.7 of
%[1] with the function lineSearch used to perform the line search.
%
%EXAMPLE 1:
%The first example is that used in the lineSearch file.
% f=@(x)deal((x(1)+x(2)-3)*(x(1)+x(2))^3*(x(1)+x(2)-6)^4,... %The function
%            [(-6+x(1)+x(2))^3*(x(1)+x(2))^2*(54+8*x(1)^2+x(2)*(-45+8*x(2))+x(1)*(-45+16*x(2)));
%            (-6+x(1)+x(2))^3*(x(1)+x(2))^2*(54+8*x(1)^2+x(2)*(-45+8*x(2))+x(1)*(-45+16*x(2)))]);%And the gradient as the second return.
% %Note that the deal function is used to make an anonymous function have
% %two outputs.
% x0=[0.5;0.25];
% [xMin,fMin,exitCode]=quasiNetwonBFGS(f,x0)
%The optimum point found is such that sum(xMin) is approximately
%1.73539450 with a minimum function value of approximately -2.1860756.
%
%EXAMPLE 2:
%The second example requires that the cost function be placed in a separate
%file. The cost function is
% function [fx,g]=objFun(x)
%    n=length(x);
%    fx=0;
%
%    g=zeros(n,1);
%    for i=0:(n-1)
%        if(mod(i,2)==1)
%            continue;
%        end
%        t1=1-x(i+1);
%        t2=10*(x(i+1+1)-x(i+1)^2);
%        g(i+1+1)=20*t2;
%        g(i+1)=-2*(x(i+1)*g(i+1+1)+t1);
%        fx =fx+t1^2+t2^2;
%    end
% end
%It is used as
% n=100;
% x0=zeros(n,1);
% for i=0:(n-1)
%    if(mod(i,2)==1)
%        continue;
%    end
%
%    x0(i+1)=-1.2;
%    x0(i+1+1)=1;
% end
% f=@(x)objFun(x);
% [xMin,fMin,exitCode]=quasiNetwonBFGS(f,x0)
%whereby the optimal solution is all ones with a minimum function value of
%zero. This second example is the same as that provided with the L-BFGS
%library in C.
%
%REFERENCES:
%[1] D. P. Bertsekas, Nonlinear Programming, 2nd ed. Belmont, MA: Athena
%    Science, 1999.
%
%July 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

xDim=size(x0,1);
if(nargin<3||isempty(D0))
    D=eye(xDim,xDim);
else
    D=D0;
end

if(nargin<4||isempty(epsilon))
    epsilon=1e-6;
end

if(nargin<5||isempty(deltaTestDist))
    deltaTestDist=0;
end

if(nargin<6||isempty(delta))
    delta=0;
end

if(nargin<7)
    lineSearchParams=[];
end

if(nargin<8||isempty(scaleD))
   scaleD=false;
end

if(nargin<9||isempty(maxIter))
   maxIter=1000;
end

xPrev=x0;
[fValPrev,gradFPrev]=f(x0);

if(deltaTestDist>0)
    pastFVals=zeros(deltaTestDist,1);
    pastFVals(1)=fValPrev;
end
for curIter=1:maxIter
    %Equation 1.181 for the descent direction.
    d=-D*gradFPrev;

    %Perform a line search in the given descent direction.
    [xCur,fValCur,gradFCur,~,~,exitCode]=lineSearch(f,xPrev,d,[],[],lineSearchParams);

    if(isempty(xCur))
        xMin=[];
        fMin=[];
        return;
    end

    %Check for convergence based on the gradient.
    if(norm(gradFCur)<epsilon*max([1, norm(xCur)]))
        xMin=xCur;
        fMin=fValCur;
        exitCode=0;
        return;
    end

    %Check for convergence based on the actual function value.
    if(deltaTestDist~=0)
        if(pastFVals(end)-fValCur<=delta*fValCur)
            xMin=xCur;
            fMin=fValCur;
            exitCode=1;
            return;
        end
        pastFVals=circshift(pastFVals,[1,0]);
        pastFVals(1)=fValCur;
    end

    %After taking the step, update the inverse Hessian approximation.
    D=updateHessianApprox(D,xCur,xPrev,gradFCur,gradFPrev,scaleD);

    xPrev=xCur;
    gradFPrev=gradFCur;
end

%The maximum number of iterations elapsed without convergence
xMin=xCur;
fMin=fValCur;
exitCode=-1;
end

function D=updateHessianApprox(D,xCur,xPrev,gradFCur,gradFPrev,scaleD)
%This is the update for the inverse Hessian estimate in the BFGS algorithm.
%All equation numbers refer to Section 1.7 of [1].

q=gradFCur-gradFPrev;%Equation 1.184
p=xCur-xPrev;%Equation 1.183

%A test is added so that D does not change if p'*q is very close to zero.
%Otherwise, D would just be full of NaN and Inf values. The extra realmin
%test avoids denormalized numbers.
if(p'*q>max(realmin,max(eps(p),eps(q))))
    if(scaleD==true)
    %If D should be scaled as in Equation 1.201. As noted, this can improve
    %the condition number of D and is often recommended after the first
    %iteration.
        D=(p'*q/(q'*D*q))*D;
    end

    %The BFGS update from Problem 1.7.2 with the correction from the errata
    %for the denominator of the p*p' term.
    D=D+(1+q'*D*q/(p'*q))*(p*p')/(p'*q)-(D*q*p'+p*q'*D)/(p'*q);
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% Assume a target positioned at x = 1, travelling with speed v = 0.1
state1 = [1000;0.1;1000;0];    % Assume state with four dimensions [x_pos, x_vel, y_pos, y_vel]
state2 = [1500;0.1;1500;0];

% Create an instance of a 1D Constant Velocity model
obs = RangeBearing2CartesianX('NumStateDims',4,'MeasurementErrVariance',...
                              [(pi/45)^2,10^2],'Mapping',[1 3]);

% View the transition matrix and process covariance matrices
Q = obs.covar();

% Predict the target's position and velocity after the interval has passed
measurement1  = obs.feval(state1);
measurement2 = obs.feval(state2);

% Generate 50 random noise samples from the dynamic model
noise1 = obs.random(500000);
noise2 = obs.random(500000);

% Add noise to the measurements
Yk1 = noise1 + measurement1;
Xk1 = obs.finv(Yk1);

Yk2 = noise2 + measurement2;
Xk2 = obs.finv(Yk2);

Xk = [Xk1,Xk2];
Yk = [Yk1,Yk2];

figure;

% Left plot
subplot(1,2,1);
[bandwidth,density,X,Y]=kde2d(Yk([1,2],:)');
surf(X,Y,density);
hold on;
shading interp
colormap(jet(3000))
title('Radar measurement noise (Polar)')
xlabel('Bearing (rad)');
ylabel('Range (m)');

% Right plot
subplot(1,2,2);
[bandwidth,density,X,Y]=kde2d(Xk([1,3],:)');
%contour3(X,Y,density,50);
surf(X,Y,density);
hold on;
shading interp
colormap(jet(3000))
title('Radar measurement noise (Cartesian)')
xlabel('X (m)');
ylabel('Y (m)');

% model = PositionalObsModelX(config);%PositionalObsModelX(config);%Polar2CartGaussianModelX(config);
%
% xkm1 = [0.5; 0.3; sqrt(2); sqrt(2)];
% pkm1 = [0 0.5; 0 0.3; sqrt(2) sqrt(2); sqrt(2) sqrt(2)];
% Pkm1 = [0 1 0 1; 0 1 0 1; 0 1 0 1; 1 1 0 1];
% yk = model.obs(1,xkm1);
% Yk = model.sample(1, yk, 500000);
% pk = model.obs_cov();
% Pk = model.eval(1,Yk,xkm1);
% [bandwidth,density,X,Y]=kde2d(Yk');
% figure;
% %contour3(X,Y,density,50);
% surf(X,Y,density);
% hold on;

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function [ y, m, d, f ] = jed_to_ymdf_alexandrian ( jed )

%*****************************************************************************80
%
%% JED_TO_YMDF_ALEXANDRIAN converts a JED to an Alexandrian YMDF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real JED, the Julian Ephemeris Date.
%
%    Output, integer Y, M, D, real F,
%    the YMDF date.
%

%
%  Determine the computational date (Y'/M'/D').
%
  j = floor ( jed + 0.5 );
  f = ( jed + 0.5 ) - j;

  j_prime = j + 124;

  y_prime = floor ( ( 4 * j_prime + 3 ) / 1461 );
  t_prime = floor ( mod ( 4 * j_prime + 3, 1461 ) / 4 );
  m_prime = floor ( t_prime / 30 );
  d_prime = mod ( t_prime, 30 );
%
%  Convert the computational date to a calendar date.
%
  d = d_prime + 1;
  m = mod ( m_prime, 13 ) + 1;
  y = y_prime - 4690 + floor ( ( 13 - m ) / 13 );

  return
end


text:
function obj = sym(obj)
% the symmetric part of the tensors

for i = 1:length(obj.u)
  if  obj.rank(i) > 1
    obj.u{i} = sym(obj.u{i});
  end
end

end

author: mtex-toolbox

repo: mtex

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


text:
function [grRules_new,genes,rxnGeneMat] = replaceGrRules(grRules,idMapping)
%replaceGrRules  Replace grRules with another set of gene IDs.
%
% NOTE: This function is adapted from "translateGrRules" and specifically
%       designed for converting grRules from current gene ID to another
%       gene ID type
%
% NOTE: The function also "cleans" the grRules, meaning duplicate genes,
%       extra parentheses, etc. will be removed from each grRule.
%
% Usage:
%
%   [grRules_new,genes,rxnGeneMat] = replaceGrRules(grRules,idMapping)
%
%
% Inputs:
%
%   grRules      A cell array of gene-reaction rules (e.g., model.grRules).
%
%
%   idMapping    must be an Nx2 cell array, where the first column contains
%                gene IDs corresponding to the current grRules, and the
%                second column contains new gene IDs to which the rules
%                will be translated.
%
% Outputs:
%
%   grRules_new  grRules with the gene IDs converted to the new IDs.
%
%   genes        A list of all genes found in the converted grRules.
%
%   rxnGeneMat   A matrix associating reactions to genes, built from the
%                converted grRules and genes list.
%

% handle input arguments
if nargin < 2
    error('Missing input.');
end

% check input data

% get original list of genes from the grRules
genes_orig = getGenesFromGrRules(grRules);
if ismember('and',genes_orig) || ismember('or',genes_orig)
    error('Problem reading grRules. Verify that all "and" and "or" elements are lowercase and surrounded by spaces.');
end

% remove duplicate genes
genes_orig  = unique(genes_orig,'stable');

% idMapping should be a NX2 cell array
if ~(size(idMapping,2) == 2)
    error('The idMapping data structure must be a NX2 cell array.');
% input genes should be consistent with those retrieved from grRules
elseif all(~ismember(idMapping(:,1), genes_orig))
    error('The input genes are Not consistent with those retrieved from grRules.');
end

% conduct replacement

% begin by "cleaning" the original grRules
% this is not necessary but can speed up the translation if the original
% grRules are not in a simplified format.
rules_orig = cleanGrRules(grRules);

% determine logical operator type and change to "&/|" for easy manipulation
textBooleanType = false;
if any(contains(grRules,{' and ',' or '}))
    textBooleanType = true;
    grRules_new = rules_orig;
    grRules_new = regexprep(grRules_new, ' or ', ' | ');
    grRules_new = regexprep(grRules_new, ' and ', ' & ');
end

% define function convertGeneId to convert gene IDs
% use strjoin to combine ids with 'or' operator if multiple hits are found
convertGeneId = @(g) strjoin(idMapping(ismember(idMapping(:,1),g),2), ' | ');

% replace grRules with another ID type
% The next line identifies gene IDs as collections of characters that
% are not spaces, parentheses, or the symbols '&' or '|'. It then replaces
% those gene IDs with the new gene ID type, which calls the convertGeneId
% inline function to retrieve.
grRules_new = regexprep(grRules_new, '[^&|\(\) ]+', '(${convertGeneId($0)})');

% prepare output

% clean up rules (removes extra parentheses, repeated genes, etc.)
grRules_new = cleanGrRules(grRules_new);

% restore "&" as "and" and "|" as "or"
if textBooleanType
    grRules_new = regexprep(grRules_new, ' \| ', ' or ');
    grRules_new = regexprep(grRules_new, ' & ', ' and ');
end

% generate new rxnGeneMat and gene list based on converted grRules
[genes,rxnGeneMat] = getGenesFromGrRules(grRules_new);

author: SysBioChalmers

repo: Human-GEM


text:
function W = threshold_proportional(W, p)
%THRESHOLD_PROPORTIONAL     Proportional thresholding
%
%   W_thr = threshold_proportional(W, p);
%
%   This function "thresholds" the connectivity matrix by preserving a
%   proportion p (0<p<1) of the strongest weights. All other weights, and
%   all weights on the main diagonal (self-self connections) are set to 0.
%
%   Inputs: W,      weighted or binary connectivity matrix
%           p,      proportion of weights to preserve
%                       range:  p=1 (all weights preserved) to
%                               p=0 (no weights preserved)
%
%   Output: W_thr,  thresholded connectivity matrix
%
%
%   Mika Rubinov, U Cambridge,
%   Roan LaPlante, Martinos Center, MGH
%   Zitong Zhang, Penn Engineering

%   Modification history:
%   2010: Original (MR)
%   2012: Bug fix for symmetric matrices (RLP)
%   2015: Improved symmetricity test (ZZ)

n=size(W,1);                                %number of nodes
W(1:n+1:end)=0;                             %clear diagonal

if max(max(abs(W-W.'))) < 1e-10             %if symmetric matrix
    W=triu(W);                              %ensure symmetry is preserved
    ud=2;                                   %halve number of removed links
else
    ud=1;
end

ind=find(W);                                %find all links
E=sortrows([ind W(ind)], -2);               %sort by magnitude
en=round((n^2-n)*p/ud);                     %number of links to be preserved

W(E(en+1:end,1))=0;                         %apply threshold

if ud==2                                    %if symmetric matrix
    W=W+W.';                                %reconstruct symmetry
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%% logicRemoveInterior
% Below is a demonstration of the features of the |logicRemoveInterior| function

%%
clear; close all; clc;

%% Syntax
% |[Lb]=logicRemoveInterior(L);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function h=plotgeometry(geo,angle)
%PLOTGEOMETRY(GEO,ANGLE) plots a simplified version of the CBCT geometry with the
% given geomerty GEO and angle ANGLE. If angle is not given, 0 will be chosen.
%
% h=PLOTGEOMETRY(...) will return the figure handle
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Ander Biguri
%--------------------------------------------------------------------------
if nargin<2
    angle=0;
end
angle=angle*180/pi;
%% Figure stuff
h=figure('Name','Cone Beam Compute Tomography geometry');
hold on
title('Current CBCT geometry, in scale')
xlabel('X');
ylabel('Y');
zlabel('Z');
set(gcf, 'color', [1 1 1])

%% CUBE/Image

drawCube(geo.offOrigin,geo.sVoxel,'k',0.05,0);

%% Detector
drawCube([-geo.DSD+geo.DSO; geo.offDetector],[1; geo.sDetector],'r',1,angle)
plotCircle3D([0 0 0],-geo.DSD+geo.DSO);
%% source
p=plot3(geo.DSO,0,0,'.','MarkerSize',30);
rotate(p,[0  0  1],angle,[0 0 0]);
plotCircle3D([0 0 0],geo.DSO);

%% Arrows.
arrow=geo.sVoxel;
%XYZ arrows
try
arrow3d([0 arrow(1)],[0 0],[0 0],.90,5,15,'r');
arrow3d([0 0],[0 arrow(2)],[0 0],.90,5,15,'b');
arrow3d([0 0],[0 0],[0 arrow(3)],.90,5,15,'g');
catch e
    error('CBCT:plotgeometry:Arrow','arrow3D not found, make sure its added to the path');
end
%UV arrows

%%
axis equal;
view(128,26)
% grid on
% axis off;

end
function drawCube( origin, size,color,alpha,a)
% From
% http://www.mathworks.com/matlabcentral/newsreader/view_thread/235581

x=([0 1 1 0 0 0;1 1 0 0 1 1;1 1 0 0 1 1;0 1 1 0 0 0]-0.5)*size(1)+origin(1);
y=([0 0 1 1 0 0;0 1 1 0 0 0;0 1 1 0 1 1;0 0 1 1 1 1]-0.5)*size(2)+origin(2);
z=([0 0 0 0 0 1;0 0 0 0 0 1;1 1 1 1 0 1;1 1 1 1 0 1]-0.5)*size(3)+origin(3);
for i=1:6
    h=patch(x(:,i),y(:,i),z(:,i),color);
    set(h,'facealpha',alpha)
    rotate(h,[0 0 1],a,[0 0 0]);
end

end

function plotCircle3D(center,radius)

theta=0:0.1:2*pi;
v=null([0 0 1]);
points=repmat(center',1,size(theta,2))+radius*(v(:,1)*cos(theta)+v(:,2)*sin(theta));
plot3(points(1,:),points(2,:),points(3,:),'k-.');

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
% compute the posterior and likelihood of data on a GMM
function [likelihood, posterior, avgLL] = compLikelihoodGMM(data, prior, mu, invCov, diagCov, useGPU)
nGaussian = length(prior);
[dim,nFr] = size(data);

if useGPU
    likelihood = gpuArray.zeros(nGaussian, nFr);
else
    likelihood = zeros(nGaussian, nFr);
end

for j=1:nGaussian
    if diagCov
        likelihood(j,:) = prior(j) * my_mvnpdf2(data, mu(:,j)', invCov(:,j)');
    else
        likelihood(j,:) = prior(j) * my_mvnpdf2(data, mu(:,j)', invCov(:,:,j));
    end
end

if nargout>=2
    evidence = sum(likelihood,1);
    posterior = bsxfun(@times, likelihood, 1./evidence);
end
if nargout==3
    avgLL = mean(log(evidence));
end

end

author: singaxiong

repo: SignalGraph


text:
function [cmap,cax] = mcColorMap(values,cmap)
%MCCOLORMAP returns the color mappin for values given cmap.
%   MCCOLORMAP uses the color map (CMAP) to define the color mapping for
%   the range in VALUES.
%

% find range of data in values and normalize to [0 1]
maxVal = max(values);
minVal = min(values);
normVal = (values - minVal)./ (maxVal - minVal);

% map to color
noColors = length(cmap);
indx = interp1(0:1/(noColors-1):1,1:noColors,normVal,'nearest');
cmap = cmap(indx,:);
cax = [minVal maxVal];

author: Sable

repo: mcbench-benchmarks


text:
function balanceStructure=getElementalBalance(model,rxns,printUnbalanced,printUnparsable)
% getElementalBalance
%   Checks a model to see if the reactions are elementally balanced
%
%   model             a model structure
%   rxns              either a cell array of reaction IDs, a logical vector
%                     with the same number of elements as reactions in the model,
%                     of a vector of indexes. Only these reactions will be
%                     checked (opt, default model.rxns)
%   printUnbalanced   print warnings about the reactions that were
%                     unbalanced (opt, default false)
%   printUnparsable   print warnings about the reactions that cannot be
%                     parsed (opt, default false)
%
%   balanceStructure
%       balanceStatus   1 if the reaction is balanced, 0 if it's unbalanced,
%                      -1 if it couldn't be balanced due to missing information,
%                      -2 if it couldn't be balanced due to an error
%       elements
%           abbrevs     cell array with abbreviations for all used elements
%           names       cell array with the names for all used elements
%       leftComp        MxN matrix with the sum of coefficients for each of
%                       the elements (N) for the left side of the
%                       reactions (M)
%       rightComp       the corresponding matrix for the right side
%
%   Usage: balanceStructure=getElementalBalance(model,rxns,printUnbalanced,printUnparsable)

if nargin<2
    rxns=[];
elseif ~islogical(rxns) && ~isnumeric(rxns)
    rxns=convertCharArray(rxns);
end

if nargin<3
    printUnbalanced=false;
end

if nargin<4
    printUnparsable=false;
end

if ~isempty(rxns)
    indexes=~getIndexes(model,rxns,'rxns',true);
    model=removeReactions(model,indexes,true);
end

balanceStructure.balanceStatus=nan(numel(model.rxns),1);

%Get the formulas
if isfield(model,'metFormulas')
    [balanceStructure.elements, useMat, exitFlag]=parseFormulas(model.metFormulas, true);
else
    if isfield(model,'inchis')
        [balanceStructure.elements, useMat, exitFlag]=parseFormulas(model.inchis, true,true);
    else
        EM='The model must contain either the "metFormulas" or the "inchis" field in order to test for elemental balancing';
        dispEM(EM);
    end
end

balanceStructure.leftComp=zeros(numel(model.rxns),numel(balanceStructure.elements.names));
balanceStructure.rightComp=balanceStructure.leftComp;

%Look at the right and left side of the reactions separately
S{1}=model.S;
S{2}=model.S;
S{1}(S{1}>0)=0; %Left side
S{2}(S{2}<0)=0; %Right side
S{1}=abs(S{1}); %Both should have positive values

%First do it for left side and then for right
for i=1:2
    for j=1:numel(model.rxns)
        %Get the balance status of the involved mets
        I=exitFlag(S{i}(:,j)~=0);
        if any(I==-1)
            balanceStructure.balanceStatus(j)=-2;
        end
        if any(I==0)
            %Don't change a more serious error to a less serious one
            balanceStructure.balanceStatus(j)=min(-1,balanceStructure.balanceStatus(j));
        end
        %Loop through each element
        for k=1:numel(balanceStructure.elements.names)
            if i==1
                balanceStructure.leftComp(j,k)=sum(S{i}(:,j).*useMat(:,k));
            else
                balanceStructure.rightComp(j,k)=sum(S{i}(:,j).*useMat(:,k));
            end
        end
    end
end

%Now compare the left and right sides to find which are unbalanced. This is
%done even if the reaction as a whole couldn't be balanced
total=abs(balanceStructure.rightComp-balanceStructure.leftComp)>10^-8; %To deal with roundoff error

%Get which reactions are unbalanced. Don't change an error to just
%unbalanced
balanceStructure.balanceStatus(any(total,2))=min(balanceStructure.balanceStatus(any(total,2)),0);

%The remaining ones are all balanced
balanceStructure.balanceStatus(isnan(balanceStructure.balanceStatus))=1;

%Print warnings
toPrint=[];
if printUnbalanced==true
    toPrint=[toPrint;find(balanceStructure.balanceStatus==0)];
end
if printUnparsable==true
    toPrint=[toPrint;find(balanceStructure.balanceStatus<0)];
end

toPrint=sort(toPrint);
for i=1:numel(toPrint)
    if balanceStructure.balanceStatus(toPrint(i))<0
        if balanceStructure.balanceStatus(toPrint(i))==-1
            EM=['The reaction ' model.rxns{toPrint(i)} ' could not be balanced due to missing information'];
            dispEM(EM,false);
        else
            EM=['The reaction ' model.rxns{toPrint(i)} ' could not be balanced due to a parsing error'];
            dispEM(EM,false);
        end
    else
        %Find the compounds that it's not balanced for
        notBalanced=find(total(toPrint(i),:));
        for j=1:numel(notBalanced)
            EM=['The reaction ' model.rxns{toPrint(i)} ' is not balanced with respect to ' balanceStructure.elements.names{notBalanced(j)}];
            dispEM(EM,false);
        end
    end
end

% Re-order the structure entries so they're consistent with the ordering of
% the input reaction indexes
if ~isempty(rxns)
    rxns = getIndexes(model,rxns,'rxns');
    [~,i] = sort(rxns);
    balanceStructure.balanceStatus(i) = balanceStructure.balanceStatus;
    balanceStructure.leftComp(i,:) = balanceStructure.leftComp;
    balanceStructure.rightComp(i,:) = balanceStructure.rightComp;
end

repo: RAVEN


text:
function colorSegmentation ( )

%% COLORSEGMENTATION segments the "fabric.jpg" image by color.
%
%    COLORSEGMENTATION displays a figure for each of the four segmented
%    colors: red, yellow, green, and magenta.
%
%  This function requires the Image Processing Toolbox.
%

%
%  Create a cell array of color names.
%
  colors = { 'red', 'yellow', 'green', 'magenta' };
%
%  Read the image.
%
  image = imread ( 'fabric.png' );

  imagesc ( image )
  title ( 'Original image to be color-filtered' )

  for index = 1: 4
%
% Each iteration will filter a different color
%
    filteredImage = colorFilter ( image, colors{index} );
%
%  Display the extracted segment.
%
    figure;
    imagesc ( filteredImage );
    title_string = sprintf ( 'Image filtered by %s', colors{index} );
    title ( title_string )

  end

  return
end


text:
%PREX_SOFT Simple example of handling soft labels in PRTools
%
% Soft labels are implemented next to the 'crisp' and 'targets' labels.
% Like 'targets' labels they are stored in the target field of a dataset.
% Their values should be between 0 and 1. For every class a soft label
% values should be given. The density based classifiers can handle soft
% labels, interpreting them as class weights for every objects in the
% density estimation.
%
% The posterior probabilities found by classifying objects can be
% interpreted as soft labels. They, however, sum to one (over the classes),
% while this is not necessary for training and test objects.
%
% Note that the routine CLASSSIZES returns the sum of the soft labels over
% the dataset for every class separately. In contrast to crisp labels the
% sum over the classes of the output of CLASSSIZES is not necessarily
% equal to number of objects in the dataset.
%
% The routine SELDATA(A,N) returns the entire dataset in case of a soft
% labeled dataset A for every value of N and not just class N, as all
% objects may participate in all classes.

help prex_soft;
echo on

% Generate artificial soft labeled dataset using posteriors as soft labels
a = gendath([100 100]);

% retrieve a dataset with posteriors to be used for soft labels
labels = a*qdc(a)*classc;

% create a new dataset with soft labels
s = prdataset(+a);

% we just need the values of 'labels'
s = setlabtype(s,'soft',+labels);

% give the classes a name (optional, just to show how this is done)
s = setlablist(s,{'A','B'});

% experiment: % generate train set and test set
[train_s,test_s] = gendat(s,0.5);

% compute classifier that outputs posteriors
w_s = parzenc(train_s)*classc;

% apply classifier on testdata
d_s = test_s*w_s;

% result, by default for soft labeled data
% the 'soft' test type is used in testc
testc(d_s)

% compare with crisp labeling, convert train and test set to crisp labels
train_c = setlabtype(train_s,'crisp');
test_c  = setlabtype(test_s,'crisp');

 % compute classifier
w_c = parzenc(train_c)*classc;

% apply classifier on testdata
d_c = test_c*w_c;

% result, by default for crisp labeled data
% the 'crisp' test type is used in testc
testc(d_c)

echo off


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ============================================
% align_cps.m
%
% Goal: Align P to X using ICP
%
% Li Shen
% 07/29/2003 - customized using Hany Farid's icp code.
% 01/08/2004 - customized again for corresponding point set registration.

function [P,M] = align_cps(P,X)

disp('running cps_align');

M = eye(4);

%%%%%%%%%%%%%%%%%%%%
% Added by Li Shen
%%%%%%%%%%%%%%%%%%%%
% center P and X first
[P,cmP] = center(P);
[X,cmX] = center(X);

Xc = X;

muP   = mean(P);
muXc  = mean(Xc);
N     = length(P);
sig   = zeros(3);
for k = 1 : N
  sig = sig + (P(k,:)' * Xc(k,:));
end
sig   = sig/N;
sig   = sig - (muP' * muXc); % cross-covariance matrix
A     = sig - sig';
delta = [A(2,3) ; A(3,1) ; A(1,2)];
Q     = [trace(sig)  delta' ; delta  sig+sig'-trace(sig)*eye(3)];
[V,D] = eig(Q);
[val,ind] = max(diag(D));
qR    = V(:,ind); % optimal rotation
q0    = qR(1);
q1    = qR(2);
q2    = qR(3);
q3    = qR(4);
R     = [q0^2+q1^2-q2^2-q3^2  2*(q1*q2 - q0*q3)  2*(q1*q3 + q0*q2) ; ...
    2*(q1*q2 + q0*q3)  q0^2+q2^2-q1^2-q3^2  2*(q2*q3 - q0*q1) ; ...
    2*(q1*q3 - q0*q2) 2*(q2*q3 + q0*q1)  q0^2+q3^2-q1^2-q2^2];
qT    = muXc' - R*muP'; % optimal translation

Pnew  = (R*P')';
diff  = mean( sum(((Pnew - P).^2)') );

P = Pnew;

P(:,1) = P(:,1) + qT(1);
P(:,2) = P(:,2) + qT(2);
P(:,3) = P(:,3) + qT(3);

%%%%%%%%%%%%%%%%%%%%
% Added by Li Shen
%%%%%%%%%%%%%%%%%%%%
% recover P and X
P(:,1) = P(:,1) + cmX(1);
P(:,2) = P(:,2) + cmX(2);
P(:,3) = P(:,3) + cmX(3);
X(:,1) = X(:,1) + cmX(1);
X(:,2) = X(:,2) + cmX(2);
X(:,3) = X(:,3) + cmX(3);

Mr = eye(4); Mr(1:3,1:3) = R;
Mt = eye(4); Mt(1:3,4) = qT(1:3)+cmX'-cmP';
M = Mt*Mr*M;

return;

%
% center point set
%

function [P,muP] = center(P)

n = size(P,1);
muP = mean(P);
P = P - muP(ones(1,n),:);

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function rgb = dna_to_rgb_average ( dna )

%*****************************************************************************80
%
%% DNA_TO_RGB_AVERAGE creates an RGB image from DNA.
%
%  Discussion:
%
%    The DNA is subdivided into 32 separate "chromosomes", each of which
%    represents on rectangle.
%
%    For each rectangle, we have 7 "genes", which represent the (x,y) of
%    corner 1, (x,y) of corner 2, and the R, G, and B values for the rectangle.
%
%    Each gene, in turn is represented by 8 bits, because each gene is
%    an integer between 0 and 255.
%
%    This function computes the RGB image by averaging the R, G, and B
%    values from each block.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 January 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Nick Berry,
%    A "Practical" Use for Genetic Programming,
%    http://www.datagenetics.com/blog.html
%
%  Parameters:
%
%    Input, integer DNA(56,32), the 8 bits of the 7 "genes" of the
%    32 "chromosomes" of the DNA.
%
%    Output, uint8 RGB(256,256,3), the RGB information.
%
  c = zeros ( 256, 256 );
  rgb = zeros ( 256, 256, 3 );

  for j = 1 : 32

    x1 = b8_to_i ( dna(1:8,j) );
    y1 = b8_to_i ( dna(9:16,j) );
    x2 = b8_to_i ( dna(17:24,j) );
    y2 = b8_to_i ( dna(25:32,j) );
    r = b8_to_i ( dna(33:40,j) );
    g = b8_to_i ( dna(41:48,j) );
    b = b8_to_i ( dna(49:56,j) );

    xlo = min ( x1, x2 ) + 1;
    xhi = max ( x1, x2 ) + 1;
    ylo = min ( y1, y2 ) + 1;
    yhi = max ( y1, y2 ) + 1;

    c(xlo:xhi,ylo:yhi) = c(xlo:xhi,ylo:yhi) + 1;

    rgb(xlo:xhi,ylo:yhi,1) = rgb(xlo:xhi,ylo:yhi,1) + r;
    rgb(xlo:xhi,ylo:yhi,2) = rgb(xlo:xhi,ylo:yhi,2) + g;
    rgb(xlo:xhi,ylo:yhi,3) = rgb(xlo:xhi,ylo:yhi,3) + b;

  end

  rgb(1:256,1:256,1) = rgb(1:256,1:256,1) ./ c(1:256,1:256);
  rgb(1:256,1:256,2) = rgb(1:256,1:256,2) ./ c(1:256,1:256);
  rgb(1:256,1:256,3) = rgb(1:256,1:256,3) ./ c(1:256,1:256);
  rgb = uint8 ( rgb );

  return
end


text:
function [y, slope] = groundHeight(x)
%
% This function returns the height (y) and slope (Dy) of the ground as a
% function of the horizontal position. An interesting ground profile is
% generated by summing three sine waves. This also makes the derivative
% (slope) easy to compute.
%

%Parameters for each of the three sine ewaves:
    Amp1 = 0.3;
    Freq1 = 1;
    Phase1 = 1.58;

    Amp2 = 0.1;
    Freq2 = 6;
    Phase2 = 1.9;

    Amp3 = 0.2;
    Freq3 = 3.7;
    Phase3 = 1.58;

%Compute the ground height
y = Amp1*sin(Freq1*x+Phase1)+...
    Amp2*sin(Freq2*x+Phase2)+...
    Amp3*sin(Freq3*x+Phase3);

%Compute the ground slope
slope = Amp1*Freq1*cos(Freq1*x+Phase1)+...
        Amp2*Freq2*cos(Freq2*x+Phase2)+...
        Amp3*Freq3*cos(Freq3*x+Phase3);

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [PC,Vmea]=prange(obs,nav,opt,azel,iter) %#ok

global glc
lam=nav.lam(obs.sat,:);
PC=0;Vmea=0;
i=1; j=2;

[sys,~]=satsys(obs.sat);
if sys==0,return;end

% if glc.NFREQ>=3&&(sys==glc.SYS_GAL),j=3;end

if glc.NFREQ<2||lam(i)==0||lam(j)==0,return;end

% test snr mask
if iter>1
end

% pseudorange with code bias correction
[cbias,use_dcb_flag]=getdcb(nav,obs,opt);
if use_dcb_flag==0&&sys~=glc.SYS_GLO
    cbias=gettgd(nav,obs,opt);
end
P1=obs.P(i)-cbias(1);
P2=obs.P(j)-cbias(2);

gamma=lam(j)^2/lam(i)^2;
if opt.ionoopt==glc.IONOOPT_IFLC
    if P1==0||P2==0,return;end
    PC=(gamma*P1-P2)/(gamma-1);
else
    if P1==0,return;end
    PC=P1;
end

Vmea=0.3^2;

return

repo: GINav


text:
%  Scaled least squares solver for A*x = b
%
%   Input:  A -- matrix
%           v -- right-side vector
%           w -- (optional) weight vector
%
%  Output:  (return) -- the least squares solution
%
%  syntax  >> x = ScaledLeastSquares(A,b)
%          >> x = ScaledLeastSquares(A,b,w)

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [logp, yhat, res] = tapas_rs_surprise(r, infStates, ptrans)
% Calculates the log-probability of response speed y (in units of ms^-1) according to the surprise
% model introduced in:
%
% Vossel, S.*, Mathys, C.*, Daunizeau, J., Bauer, M., Driver, J., Friston, K. J., and Stephan, K. E.
% (2013). Spatial Attention, Precision, and Bayesian Inference: A Study of Saccadic Response Speed.
% Cerebral Cortex.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Transform zetas to their native space
ze1v = exp(ptrans(1));
ze1i = exp(ptrans(2));
ze2  = exp(ptrans(3));
ze3  = exp(ptrans(4));

% Initialize returned log-probabilities, predictions,
% and residuals as NaNs so that NaN is returned for all
% irregualar trials
n = size(infStates,1);
logp = NaN(n,1);
yhat = NaN(n,1);
res  = NaN(n,1);

% Weed irregular trials out from inferred states, responses, and inputs
mu1hat = infStates(:,1,1);
mu1hat(r.irr) = [];

y = r.y(:,1);
y(r.irr) = [];

u = r.u(:,1);
u(r.irr) = [];

% Calculate alpha (i.e., attention)
alpha = 1./(1-log2(mu1hat));

% Calculate predicted response speed
rs = u.*(ze1v + ze2*alpha) + (1-u).*(ze1i + ze2*(1-alpha));

% Calculate log-probabilities for non-irregular trials
% Note: 8*atan(1) == 2*pi (this is used to guard against
% errors resulting from having used pi as a variable).
reg = ~ismember(1:n,r.irr)
logp(reg) = -1/2.*log(8*atan(1).*ze3) -(y-rs).^2./(2.*ze3);
yhat(reg) = rs;
res(reg) = y-rs;

return;

author: translationalneuromodeling

repo: tapas