text:
function spm_htw_from_fit(varargin)
% When fitting a 1st-level SPM model with multiple basis functions,
% reconstructs fitted response at each voxel beta images and basis functions, and
% estimates height (amplitude), time-to-peak, and response width (duration and half-max)
% For each condition. Saves maps of these across voxels, and contrasts
% across height, time-to-peak, and width maps if contrasts are specified.
% This provides contrast images to take to 2nd-level analyses for group statistics when
% using multiple basis functions at the 1st level.
%
% :Usage:
% ::
%
%     function spm_htw_from_fit(varargin)
%
% :NOTES:
%
% This function loads the SPM.mat file in the current directory and uses
% the basis set specified in the loaded SPM structure.
%
% : Required Inputs:
% 	None - run within an SPM first-level analysis directory
%
% :Optional Inputs:
%
%   **'amplitudes'**:
%       do amplitudes only; no contrasts
%
%   **'contrasts'**:
%       do contrasts only; no amplitudes
%
%   **'noamplitudes'**
%        skip create amping images (combination of betas
%        across basis functions)
%
%   **'nocontrasts'**
%        skip creation of contrast images
%        (running contrast images assumes that amplitude images are already
%        created)
%
%   **'all'**
%        will run both the amplitudes and contrasts sections
%
%          - The second way uses CANlab HTW code to estimate height, time to peak,
%            width, and area under the curve (see Lindquist & Wager 2007 for
%            simulations using a version of this code).
%            It requires SCANlab specific functions, in SCN_Core_Support
%            (unlike the deriv. boost).
%            To turn this OFF, enter 'nohtw' as an optional argument
%
%   **'startend'**
%        followed by starting and ending values in seconds for amplitude
%        estimation window (for HTW estimation only).
%        If you do not enter this, it will show you a plot and ask you to pick
%        these values.
%        If you enter them here as inputs, you can skip the interactive step and
%        loop over subjects.
%
%   **'condition_numbers'**
%        followed by which index numbers in your list should
%        be used to calculate h, t, w from.  You should use this if you
%        are entering regressors of no interest, besides the intercepts.
%
% :Important for Contrasts:
%    disp('Using contrasts entered as F-contrasts. Assuming the first contrast vector in each F-contrast '
%
%    disp('is a contrast of interest across the CANONICAL basis function regressors.')
%
% :Output:
% ::
% Reconstructed amplitude, time-to-peak, and duration (width) images for each event type
% e.g., for Event type (condition) 001:
% htw_amplitude_001.nii
% htw_time_to_peak_001.nii
% htw_width_001.nii
% htw_area_under_curve_001.nii
%
% Contrasts across these, using names stored in SPM.xCon.name:
% e.g.,
% con_htw_ampl_targetvsstandards5.nii
% con_htw_time_targetvsstandards5.nii
% con_htw_widt_targetvsstandards5.nii
% con_htw_area_targetvsstandards5.nii
%
% :Examples:
% ::
%
%    % RUN THIS IN COMMAND WINDOW TO BATCH
%    subj = dir('06*')
%    for i = 1:length(subj), cd(subj(i).name), spm_htw_from_fit, cd('..'); end
%
%    % ANOTHER BATCH EXAMPLE:
%    d = dir('remi*'); d = d(cat(2, d.isdir)); [mydirs{1:length(d)}] = deal(d.name)
%    for i = 1:length(mydirs), cd(mydirs{i}), spm_htw_from_fit('all', 'nodb', 'startend', [4 15]), cd('..'); end
%
%    %An example for an event-related design, specifying condition numbers to get HTW from:
%    spm_htw_from_fit('all', 'contrasts', 'condition_numbers', 1:14, 'startend', [4 10]);
%
%    % CALCULATE CONTRASTS ONLY ON ALREADY-ESTIMATED HTW IMAGES
%    spm_htw_from_fit('contrasts','condition_numbers',1:14);
%
% :References:
% ::
% Lindquist, M. A. & Wager, T. D. (2007). Validity and power in hemodynamic response modeling:
% a comparison study and a new approach. Human Brain Mapping. 8:764-84.
%
% Lindquist, M. A., Meh Loh, J., Atlas, L. Y. & Wager T. D. (2009). Modeling the hemodynamic
% response function in fMRI: efficiency, bias and mis-modeling. Neuroimage. 45:187-98.

% Notes:
% This function is adapted from apply_derivative_boost.m, which was
% originally intended to implement Calhoun derivative boost, but this is
% now deprecated.
% Original documentation:
%        In addition, 'amplitudes' now has two separate parts:
%          - The first uses Vince Calhoun's derivative boost (Calhoun, 2004) to
%             estimate amplitudes.  NOTE: *We have not worked out the scaling yet, so
%             I'm not sure this is working right*
%             To turn this OFF, enter 'nodb' as an optional argument
%

    spmname = fullfile(pwd, 'SPM.mat');     % SETUP INPUTS
    if ~exist(spmname, 'file'), error('You must be in an SPM 1st-level results directory with SPM5 SPM.mat file.'); end

    load(spmname);

    nbf = size(SPM.xBF.bf, 2);

    docontrasts = 1;
    doamps = 1;
    nodb = 1;
    nohtw = 0;
    condition_numbers = [];
    do_downsample = [];          % downsample; default = 1 sec if units are seconds (1/dt)

    for i = 1:length(varargin)
        if ischar(varargin{i})
            switch varargin{i}
                % reserved keywords
                case 'all', docontrasts = 1; doamps = 1;

                case 'amplitudes', doamps = 1; docontrasts = 0;

                case 'contrasts', docontrasts = 1; doamps = 0;

                case 'noamplitudes', doamps = 0;

                case 'nocontrasts', docontrasts = 0;

                case 'db', nodb = 0; % original DB estimation (legacy, deprecated)

                case 'nohtw', nohtw = 1; % skip HTW estimation

                case 'startend', startend = varargin{i + 1};  % starting and ending values in seconds for amplitude estimation window (for HTW estimation only).

                case 'condition_numbers',  condition_numbers = varargin{i + 1};

                case 'nodownsample', do_downsample = 0;

                case 'downsample', do_downsample = varargin{i + 1};

                otherwise, warning(['Unknown input string option:' varargin{i}]);
            end
        end
    end

    if isempty(do_downsample)
        % default downsampling
        do_downsample = round(1 ./ SPM.xBF.dt);
    end

    if ~(docontrasts || doamps)
        disp('Nothing to do! Enter ''contrasts'' ''amplitudes'' or ''all'' as input argument.');
        return
    end

    % ---------------------------------------------
    % FILE NAMES and imgtype
    % ---------------------------------------------

    imgs = dir(sprintf(['beta*img'])); imgs = char(imgs.name);

    if isempty(imgs)
        imgs = dir(sprintf(['beta*nii'])); imgs = char(imgs.name);
        imgtype = '.nii';
    else
        imgtype = '.img';
    end

    if isempty(imgs)
        error('Cannot find beta*img or beta*nii files in current folder')
    end

    n = size(imgs, 1);

    if doamps
        % ---------------------------------------------
        % ---------------------------------------------

        % ESTIMATE AMPLITUDES

        % ---------------------------------------------
        % ---------------------------------------------

        fprintf('Found %3.0f beta images', n); fprintf('\n');

        %load(spmname);
        nsess = length(SPM.Sess);
        fprintf('I think there are %3.0f sessions (runs)', nsess); fprintf('\n');

        if ~isempty(condition_numbers)
            wh_intercept = true(1, size(imgs, 1));  % exclude these
            wh_intercept(condition_numbers) = 0;
            wh_intercept = find(wh_intercept);

            fprintf('\nIncluding only these images: ');
            fprintf('%3.0f ', condition_numbers);
            fprintf('\n');

        else
            wh_intercept = SPM.xX.iB;
            fprintf('\nI think these images are intercepts, and am not using them: ');
            fprintf('%3.0f ', wh_intercept);
            fprintf('\n');
        end

        imgs(wh_intercept, :) = [];

        %not used; only when using image_eval_function
        %mask_img = './mask.img';

        n = size(imgs, 1);

        if n / nbf ~= round(n / nbf), error('Error!  Wrong number of images for the specified number of basis functions.'); end

        if nodb
            % skip Derivative Boost estimation and go straight to HTW
        else
            % DO DB

            switch nbf
                case 2
                    derivative_case = 'timeonly';
                case 3
                    derivative_case = 'timedispersion';

                otherwise
                    warning('Deriv. Boost only works for SPM canonical hrf with time or time + dispersion derivatives.  This SPM.mat doesn''t match those specs.');
            end

            % Not used; only for image eval function
            % boost = @(b) sign(b(1)) .* sqrt(sum(b .^ 2));

            % ---------------------------------------------
            % CALCULATE
            % ---------------------------------------------
            cond_indx = 1;

            for i = 1: nbf : (n - nbf + 1)

                imgs_cond = imgs(i : i+nbf - 1, :);

                disp('Working on :')
                disp(imgs_cond)

                out_name = sprintf(['db_amplitude_%03d',imgtype ], cond_indx);

                % This code uses SCN lab tools to create images
                % ---------------------------------------------
                % %     y = image_eval_function(imgs_cond, boost, 'mask', mask_img, ...
                % %         'outimagelabels', {out_name});
                % ---------------------------------------------

                % This code uses SPM instead
                % ---------------------------------------------
                switch derivative_case
                    case 'timeonly'
                        spm_imcalc(imgs_cond, out_name, 'sign(i1) .* sqrt(i1.^2 + i2.^2)');
                    case 'timedispersion'
                        spm_imcalc(imgs_cond, out_name, 'sign(i1) .* sqrt(i1.^2 + i2.^2 + i3.^2)');
                    otherwise
                        error('Basis set is incompatible with DB estimation!');
                end

                fprintf('Created %s\n', out_name);
                % ---------------------------------------------

                cond_indx = cond_indx + 1;
            end

            % Get and save names
            db_amp_names = [];
            for i = 1:nsess
                sessnames = char(SPM.Sess(i).Fc.name);
                sessnames = [repmat(sprintf('Sess%02d_', i), size(sessnames, 1), 1) sessnames];
                db_amp_names = strvcat(db_amp_names, sessnames);
            end

            save db_amplitude_names db_amp_names
            disp(db_amp_names);
            disp(' ')
            disp('Saved DB amplitude condition names for each image in db_amplitude_names.mat');

            fprintf('\n*-----------------------------*\nApplied DB successfully\n*-----------------------------*\n')

        end

        if nohtw
            % skip this

        else

            % ---------------------------------------------
            % Estimated amplitude from fit: HTW
            % ---------------------------------------------
            % This code uses SCN lab tools to create images

            disp(' ')
            disp('Next: Estimating amplitude, time to peak, width, and area-under-curve images from fitted response using SCN lab code.')
            disp(' ')

            % downsample bf, if requested
            if do_downsample
                mytimeres = SPM.xBF.dt * do_downsample;
                SPM.xBF.bf = downsample(SPM.xBF.bf, do_downsample);

            else
                mytimeres = SPM.xBF.dt;
            end

            if exist('startend', 'var')
                % just check, and use input values
                if length(startend) ~= 2, error('Startend input must have two values, a starting and ending value in seconds for the amp. estimate window'), end

            else
                % Set range in sec
                htw_from_fit(SPM.xBF.bf, ones(size(SPM.xBF.bf, 2), 1), mytimeres, 'plot', 'verbose');

                disp(' ')
                disp('Enter the range in seconds within which to estimate peak amplitude.')
                disp('Example: type [3 12] and press return for a typical event-related setup.');
                disp('More sustained responses, like pain responses, may require a longer window.');
                disp('This estimates the amplitude of the IMPULSE RESPONSE, before convolution with the stimulus function')
                disp('so if you have an epoch design, a typical window of [3 12] sec is still appropriate.');
                disp('Also note: AUC images are calculated as the area under the curve within the window you specify.')
                disp(' ')
                disp('In the future, you can input ''startend'', [3 12] (for example) to skip interactive selection')
                startend = input('Enter your choice in [ ] and press return: ');
            end

            % Test your choice by showing you a plot
            htwfunction = @(b) htw_from_fit(SPM.xBF.bf, b, mytimeres, 'startval', startend(1), 'endval', startend(2), 'plot', 'verbose');
            htwfunction(ones(size(SPM.xBF.bf, 2), 1))
            drawnow

            % Create without plot option for loop through brain.
            htwfunction = @(b) htw_from_fit(SPM.xBF.bf, b, mytimeres, 'startval', startend(1), 'endval', startend(2));

            disp('Check the screen for a plot of your choice of window.')
            disp(' ')

            % ---------------------------------------------
            % CALCULATE
            % ---------------------------------------------
            cond_indx = 1;

            for i = 1: nbf : (n - nbf + 1)

                imgs_cond = imgs(i : i+nbf - 1, :);

                disp('Working on :')
                disp(imgs_cond)

                clear out_name
                out_name{1} = sprintf(['htw_amplitude_%03d',imgtype], cond_indx);
                out_name{2} = sprintf(['htw_time_to_peak_%03d',imgtype], cond_indx);
                out_name{3} = sprintf(['htw_width_%03d',imgtype], cond_indx);
                out_name{4} = sprintf(['htw_area_under_curve_%03d',imgtype], cond_indx);

                % ---------------------------------------------
                [h, t, w, auc] = image_eval_function(imgs_cond, htwfunction, 'mask', fullfile(pwd, ['mask', imgtype]), ...
                    'outimagelabels', out_name);

                h;t;w;auc;  % we need the outputs above to tell it to write 4 images.
                % ---------------------------------------------

                cond_indx = cond_indx + 1;
            end

            % Get and save names
            htw_amp_names = [];
            for i = 1:nsess
                sessnames = char(SPM.Sess(i).Fc.name);
                sessnames = [repmat(sprintf('Sess%02d_', i), size(sessnames, 1), 1) sessnames];
                htw_amp_names = strvcat(htw_amp_names, sessnames);
            end

            if ~exist(fullfile(pwd, 'db_amplitude_names.mat'), 'file')
                save db_amplitude_names htw_amp_names
            else
                save db_amplitude_names -append htw_amp_names
            end
            disp(htw_amp_names);
            disp(' ')
            disp('Saved HTW amplitude condition names for each image in db_amplitude_names.mat');

            fprintf('\n*-----------------------------*\nApplied HTW estimation successfully\n*-----------------------------*\n')

        end

    end  % amplitudes

    % ---------------------------------------------
    % ---------------------------------------------

    % CREATE CONTRAST FOR THIS SUBJECT

    % ---------------------------------------------
    % ---------------------------------------------

    if docontrasts

        % Load contrast vectors

        % -----------------------------------------
        %load(spmname);
        % nsess = length(SPM.Sess);

        % Define which indices to exclude from contrasts
        if ~isempty(condition_numbers)
            nconvals = length(SPM.xCon(1).c(:, 1));  % test contrast
            wh_intercept = true(1, nconvals);  % exclude these
            wh_intercept(condition_numbers) = 0;
            wh_intercept = find(wh_intercept);

            fprintf('\nIncluding only these images: ');
            fprintf('%3.0f ', condition_numbers);
            fprintf('\n');

        else
            wh_intercept = SPM.xX.iB;
            fprintf('\nI think these images are intercepts, and am not using them: ');
            fprintf('%3.0f ', wh_intercept);
            fprintf('\n');
        end

        if ~isfield(SPM, 'xCon')
            error('Enter F-contrasts first, with the first contrast vector in each F-contrast the contrast across the CANONICAL basis function.');
        else
            disp('Using contrasts entered as F-contrasts. Assuming the first contrast vector in each F-contrast ')
            disp('is a contrast of interest across the CANONICAL basis function regressors.')
        end

        wh_F = strmatch('F', char(SPM.xCon.STAT), 'exact');

        wh_T = strmatch('T', char(SPM.xCon.STAT), 'exact');

        wh_F = [wh_F wh_T];

        if isempty(wh_F) && isempty(wh_T), error('No F-contrasts or T-contrasts entered yet.'); end

        % All sets of contrast images
        % ------------------------------------------
        if ~nodb

            ampimgs_name = sprintf(['db_amplitude_*',imgtype]);

            ampimgs = dir(ampimgs_name); ampimgs = char(ampimgs.name);

            if ~isempty(ampimgs)
                contrast_image_names_dbamp = calc_contrasts(ampimgs_name, ampimgs, wh_F, SPM, wh_intercept, nbf, imgtype);
            end

        end

        % HTW amplitude
        % ------------------------------------------
        ampimgs_name = sprintf(['htw_amplitude_*',imgtype]);

        ampimgs = dir(ampimgs_name); ampimgs = char(ampimgs.name);

        if ~isempty(ampimgs)
            contrast_image_names_htwamp = calc_contrasts(ampimgs_name, ampimgs, wh_F, SPM, wh_intercept, nbf, imgtype);
        else
            disp(['Checked for but did not find: ' ampimgs_name]);
        end

        % HTW time
        % ------------------------------------------
        ampimgs_name = sprintf(['htw_time_to_peak_*',imgtype]);

        ampimgs = dir(ampimgs_name); ampimgs = char(ampimgs.name);

        if ~isempty(ampimgs)
            contrast_image_names_htwtime = calc_contrasts(ampimgs_name, ampimgs, wh_F, SPM, wh_intercept, nbf, imgtype);
        else
            disp(['Checked for but did not find: ' ampimgs_name]);
        end

        % HTW width
        % ------------------------------------------
        ampimgs_name = sprintf(['htw_width_*',imgtype]);

        ampimgs = dir(ampimgs_name); ampimgs = char(ampimgs.name);

        if ~isempty(ampimgs)
            contrast_image_names_htwwid = calc_contrasts(ampimgs_name, ampimgs, wh_F, SPM, wh_intercept, nbf, imgtype);
        else
            disp(['Checked for but did not find: ' ampimgs_name]);
        end

        % HTW area
        % ------------------------------------------
        ampimgs_name = sprintf(['htw_area_under_curve_*',imgtype]);

        ampimgs = dir(ampimgs_name); ampimgs = char(ampimgs.name);

        if ~isempty(ampimgs)
            contrast_image_names_htwarea = calc_contrasts(ampimgs_name, ampimgs, wh_F, SPM, wh_intercept, nbf, imgtype);
        else
            disp(['Checked for but did not find: ' ampimgs_name]);
        end

        check_it = whos('contrast_image_names*');

        if isempty(check_it)
            disp('No valid images found to create contrasts on');

        else
            save('db_amplitude_names', '-append', 'contrast_image_names*');
            disp('Saved lists of contrast image names in db_amplitude_names.mat');
        end

    end % end contrasts

    %% INLINE

end  % main function

function contrast_image_names = calc_contrasts(ampimgs_name, ampimgs, wh_F, SPM, wh_intercept, nbf, imgtype)

    n = size(ampimgs, 1);

    fprintf('Found %3.0f images:\n', n);
    disp(ampimgs)

    disp('Reading image data.')
    V = spm_vol(ampimgs);
    vols = spm_read_vols(V);

    % spm_check_registration(ampimgs);
    % colormap jet

    contrast_image_names = [];

    for i = 1:length(wh_F)

        name = SPM.xCon(wh_F(i)).name;
        disp(['Calculating contrast on: ' name])
        original_name = name;

        name = deblank(name);
        wh_bad = (name == ' ' | name == ',' | name == '.' | name == '^' | name == '~' | name == '''' | name == ':' | name == '*' | name == '%' | name == ';' | name == '@' | name == '&');
        name(wh_bad) = [];

        name = ['con_', ampimgs_name(1:8), '_', name, imgtype];

        c = SPM.xCon(wh_F(i)).c(:, 1);
        c(wh_intercept) = [];
        c = c(1 : nbf : end);

        if length(c) ~= n, error('Contrast is wrong length for some reason!  Coding error in this function?  Or wrong number of db_amplitude images.'); end

        fprintf('Contrast values: ')
        fprintf('%01d ', c)
        fprintf('\n')

        % calculate and save

        contrast_image_calc(name, c, vols, V, original_name)

        disp(['Written: ' name]);

        contrast_image_names = strvcat(contrast_image_names, name);

        disp(' ');

    end

    fprintf('\n*-----------------------------*\nContrasts Done successfully!\n*-----------------------------*\n')
    spm_check_registration(contrast_image_names);

end

function contrast_image_calc(Q, myc, vols, V, original_name)

    % FROM:
    % function contrast_image(Q, myc)
    %
    % Tor Wager
    %
    % Creates a contrast image called Q (do not include path)
    % Given a list of img files P (spm format, with path)
    % and a contrast vector myc
    % In the directory of 1st image in P.

    if ~(length(myc) == size(vols,4))
        error('Contrast vector length is not equal to number of image files.')
    end

    myc2 = zeros(size(vols));

    for i = 1:length(myc)
        myc2(:,:,:,i) = myc(i);
    end

    cvol = vols .* myc2;
    cvol = sum(cvol,4);

    % -------------------------
    % write
    % -------------------------
    dd = fileparts(V(1).fname);
    Q = fullfile(dd, Q);
    Vo = V(1);
    Vo.fname = Q;
    Vo.descrip = ['Contrast ' original_name];

    spm_write_vol(Vo,cvol);

end


text:
%% Run this script in the off directory to convert all the .ply files to .mat files with the same name
clear all; close all
addpath(genpath('./../../../Tools/'))
d = dir('./*.ply')
N_models = length(d);

for i=1:N_models
    model_filename = d(i).name; tokens = split(model_filename,'.'); model_name = tokens{1};
    shape = plyread(model_filename);
    VERT = [shape.vertex.x,shape.vertex.y,shape.vertex.z];
    TRIV = cell2mat(shape.face.vertex_indices) + 1; %+1 due to zero based indexing in the .ply file
    n = size(VERT,1);
    m = size(TRIV,1);
    savename = [model_name,'.mat'];
    save(savename,'VERT','TRIV','n','m')
end

author: OshriHalimi


function [l, u] = ComputeDistanceExtremes(X, a, b, M)
% [l, u] = ComputeDistanceExtremes(X, a, b, M)
%
% Computes sample histogram of the distances between rows of X and returns
% the value of these distances at the a^th and b^th percentils.  This
% method is used to determine the upper and lower bounds for
% similarity / dissimilarity constraints.
%
% X: (n x m) data matrix
% a: lower bound percentile between 1 and 100
% b: upper bound percentile between 1 and 100
% M: Mahalanobis matrix to compute distances
%
% Returns l: distance corresponding to a^th percentile
% u: distance corresponding the b^th percentile

if (a < 1 || a > 100),
    error('a must be between 1 and 100')
end
if (b < 1 || b > 100),
    error('b must be between 1 and 100')
end

n = size(X, 1);

num_trials = min(100, n*(n-1)/2);

% we will sample with replacement
dists = zeros(num_trials, 1);
for (i=1:num_trials),
    j1 = ceil(rand(1)*n);
    j2 = ceil(rand(1)*n);
    dists(i) = (X(j1,:) - X(j2,:))*M*(X(j1,:) - X(j2,:))';
end

[f, c] = hist(dists, 100);
l = c(floor(a));
u = c(floor(b));


text:
% VECTDATA - vector data interpolation with optional moving
%              average.
%
% Usage:
%   >> [interparray timesout] = vectdata( array, timesin, 'key', 'val', ... );
%
% Inputs:
%   array      - 1-D or 2-D float array. If 2-D, the second dimension
%                only is interpolated.
%   timesin    - [float vector] time point indices. Same dimension as
%                the interpolated dimension in array.
%
% Optional inputs
%   'timesout' - [float vector] time point indices for interpolating
%                data.
%   'method'   - method for interpolation
%        'linear'    -> Triangle-based linear interpolation (default).
%        'cubic'     -> Triangle-based cubic interpolation.
%        'nearest'   -> Nearest neighbor interpolation.
%        'v4'        -> MATLAB 4 griddata method.
%   'average'  - [real] moving average in the dimension of timesin
%                note that extreme values might be inaccurate (see 'borders').
%                Default none or [].
%   'avgtype'  - ['const'|'gauss'] use a const value when averaging (array of
%                ones) or a gaussian window. Default is 'const'.
%   'border'   - ['on'|'off'] correct border effect when smoothing.
%                default is 'off'.
%
% Outputs:
%   interparray - interpolated array
%   timesout    - output time points
%
% Author: Arnaud Delorme, CNL / Salk Institute, 20 Oct 2002
%
% See also: GRIDDATA

% Copyright (C) 2002 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [interparray, timesout] = vectdata( array, timevect, varargin );

if nargin < 3
    help vectdata;
    return;
end

g = finputcheck( varargin, { 'timesout'   'real'  []           [];
                             'average'    'real'  []           [];
                             'gauss'      'real'  []           [];
                             'border'     'string' { 'on','off' } 'off';
                             'avgtype'    'string' { 'const','gauss' } 'const';
                             'method'     'string' { 'linear','cubic','nearest','v4' } 'linear'});
if ischar(g), error(g); end

if size(array,2) == 1
    array = transpose(array);
end

if ~isempty(g.average)
    timediff = timevect(2:end) -timevect(1:end-1);
    if any( (timediff - mean(timediff)) > 1e-8 ) % not uniform values
        fprintf('Data has to be interpolated uniformly for moving average\n');
        minspace = mean(timediff);
        newtimevect = linspace(timevect(1), timevect(end), ceil((timevect(end)-timevect(1))/minspace));
        array = interpolate( array, timevect, newtimevect, g.method);
        timevect = newtimevect;
    end
    oldavg = g.average;
    g.average = round(g.average/(timevect(2)-timevect(1)));
    if oldavg ~= g.average
        fprintf('Moving average updated from %3.2f to %3.2f (=%d points)\n', ...
                oldavg, g.average*(timevect(2)-timevect(1)), g.average);
    end
    if strcmpi(g.border, 'on')
        if strcmpi(g.avgtype, 'const')
            array = convolve(array, ones(1, g.average));
        else
            convolution = gauss2d(1,g.average,1,round(0.15*g.average));
            array = convolve(array, convolution);
        end
    else
        if strcmpi(g.avgtype, 'const')
            array = conv2(array, ones(1, g.average)/g.average, 'same');
        else
            convolution = gauss2d(1,g.average,1,round(0.15*g.average));
            array = conv2(array, convolution/sum(convolution), 'same');
        end
    end
end

interparray = interpolate( array, timevect, g.timesout, g.method);
timesout = g.timesout;

% interpolation function
% ----------------------
function [interparray] = interpolate( array, timesin, timesout, method);
    interparray = zeros(size(array,1), length(timesout));
    for index = 1:size(array,1)
        tmpa = [array(index,:) ; array(index,:)];

        [Xi,Yi,Zi] = griddata(timesin, [1 2]', tmpa, timesout, [1 2]', method); % Interpolate data
        interparray(index,:) = Zi(1,:);
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function blas0_test01 ( )

%******************************************************************************/
%
%% BLAS0_TEST01 tests R4_ABS.
%
%  Licensing:
%
%   This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 June 2008
%
%  Author:
%
%    John Burkardt
%
  r4_hi = 5.0;
  r4_lo = -3.0;
  test_num = 10;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BLAS0_TEST01\n' );
  fprintf ( 1, '  R4_ABS returns the absolute value of an R4.\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num
    [ r4, seed ] = r4_uniform_ab ( r4_lo, r4_hi, seed );
    r4_absolute = r4_abs ( r4 );
    fprintf ( 1, '  %10.6f  %10.6f\n', r4, r4_absolute );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [dLdp,iCpY,L] = mci_approach_deriv (P,M,U,Y)
% Gradient of log-likelihood for approach model
% FORMAT [dLdp,iCpY,L] = mci_approach_deriv (P,M,U,Y)
%
% dLdp      gradient of log joint
% iCpY      curvature (Fisher Information)
% L         log joint
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: mci_approach_deriv.m 6548 2015-09-11 12:39:47Z will $

G = mci_approach_gen (P,M,U);
if isstruct(Y)
    e = Y.y-G;
else
    e = Y-G;
end

V=exp(P(1));
tau=exp(P(2));
t=U.X;

y=-60+V*(1-exp(-t/tau));

dydp = [V*(1-exp(-t/tau)),-V*exp(-t/tau).*(t/tau)];
dLdp = dydp'*M.iCe*e;
iCpY = dydp'*M.iCe*dydp;

if nargout > 2
    L = mci_approach_like (P,M,U,Y);
end

author: spm


text:
%   Matlab script to input initial parameters for genas and genascum
%
%                                                R. Zuniga, 4/94

report_this_filefun(mfilename('fullpath'));

minmg = minma;
maxmg = maxma;
magstep = 0.5
%
if isempty(newcat), newcat = a; end   % verify whether to start with
% original catalogue
% make the interface for input
%
figure_w_normalized_uicontrolunits(mess)
clf
set(gcf,'Name','Genas Control Panel');
%  set(gcf,'visible','off');
set(gca,'visible','off');
set(gcf,'pos',[ 0.02  0.6 0.35 0.35])
%par1 = ((teb-t0b)*365)/100;
%
freq_field2=uicontrol('Style','edit',...
    'Position',[.80 .60 .15 .10],...
    'Units','normalized','String',num2str(minmg),...
    'Callback','minmg=str2double(get(freq_field2,''String'')); set(freq_field2,''String'',num2str(minmg));');

freq_field3=uicontrol('Style','edit',...
    'Position',[.80 .45 .15 .10],...
    'Units','normalized','String',num2str(maxmg),...
    'Callback','maxmg=str2double(get(freq_field3,''String'')); set(freq_field3,''String'',num2str(maxmg));');

freq_field=uicontrol('Style','edit',...
    'Position',[.80 .30 .15 .10],...
    'Units','normalized','String',num2str(magstep),...
    'Callback','magstep=str2double(get(freq_field,''String'')); set(freq_field,''String'',num2str(magstep));');

uicontrol('Style','Pushbutton',...
    'Position',[.80 .05 .15 .15 ],...
    'Units','normalized','Callback','web([''file:'' hodi ''/zmapwww/chap8.htm#997197''])','String','Info');

close_button=uicontrol('Style','Pushbutton',...
    'Position',[.60 .05 .15 .15 ],...
    'Units','normalized','Callback','welcome','String','Cancel');

genas_button1=uicontrol('Style','Pushbutton',...
    'Position',[.20 .05 .15 .15 ],...
    'Units','normalized',...
    'Callback','welcome,gogenas',...
    'String','Genas');

% genas_button2=uicontrol('Style','Pushbutton',...
%     'Position',[.40 .05 .15 .15 ],...
%     'Units','normalized',...
%     'Callback','welcome,inmakegr',...
%     'String','GenGrid');

txt5 = text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0. 0.65 0 ],...
    'Rotation',0 ,...
    'FontSize',ZmapGlobal.Data.fontsz.m ,...
    'Fontweight','bold',...
    'String','Minimum magnitude:');

txt4 = text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0. 0.48 0 ],...
    'Rotation',0 ,...
    'Fontweight','bold',...
    'FontSize',ZmapGlobal.Data.fontsz.m ,...
    'String','Maximum magnitude:');

txt1 = text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0. 0.31 0 ],...
    'Rotation',0 ,...
    'Fontweight','bold',...
    'FontSize',ZmapGlobal.Data.fontsz.m ,...
    'String','Magnitude increment (step):');

set(gcf,'visible','on');

author: CelsoReyes


text:
function [ar,xi,kappa,ehat] = lpcana(x,M)
%  lpcana --> Linear prediction analysis.
%
%    [ar,xi,kappa,ehat] = lpcana(x,M)
%
%    The function performs autocorrelation based LP analysis on the
%    signal vector x using the Levinson-Durbin recursion. Thus, the
%    function finds the coefficients, ar=[1 -a(1) ... -a(M)], of an
%    M'th order forward linear predictor
%
%      xhat(n) = a(1)*x(n-1) + a(2)*x(n-2) + ... + a(M)*x(n-M)
%
%    such that the sum of the squares of the prediction errors
%
%      ehat(n) = x(n) - xhat(n)
%
%    is minimized. The reflection coefficients are returned in the
%    vector kappa, and the prediction error energies for the 0'th to
%    the M'th order solution are returned in the vector xi. Finally,
%    the residual signal, ehat, is obtained by applying the inverse
%    filter A(z) to the signal frame.

% Short-term autocorrelation.
[rx,eta] = xcorr(x,M,'biased');

% LP analysis based on Levinson-Durbin recursion.
[a,xi,kappa] = durbin(rx(M+1:2*M+1),M);
ar = [1; -a];

% Prediction error signal obtained by inverse filtering.
ehat = filter(ar,1,x);

author: Sable

repo: mcbench-benchmarks


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Q = INVERSEKINEMATIC_IRB1200(robot, T)
%   Solves the inverse kinematic problem for the ABB IRB 1200 robot
%   where:
%   robot stores the robot parameters.
%   T is an homogeneous transform that specifies the position/orientation
%   of the end effector.
%
%   A call to Q=INVERSEKINEMATIC_IRB1200 returns 8 possible solutions, thus,
%   Q is a 6x8 matrix where each column stores 6 feasible joint values.
%
%
%   Example code:
%
%   abb=load_robot('abb', 'IRB1200');
%   q = [0 0 0 0 0 0];
%   T = directkinematic(abb, q);
%   %Call the inversekinematic for this robot
%   qinv = inversekinematic(abb, T);
%   check that all of them are feasible solutions!
%   and every Ti equals T
%   for i=1:8,
%        Ti = directkinematic(abb, qinv(:,i))
%   end
%	See also DIRECTKINEMATIC.
%
%   Author: Arturo Gil Aparicio
%           Universitas Miguel Hernandez, SPAIN.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function q = inversekinematic_irb1200(robot, T)

%initialize q,
%eight possible solutions are generally feasible
q=zeros(6,8);

% %Evaluate the parameters
% theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
L6=d(6);

%T= [ nx ox ax Px;
%     ny oy ay Py;
%     nz oz az Pz];
Px=T(1,4);
Py=T(2,4);
Pz=T(3,4);

%Compute the position of the wrist, being W the Z component of the end effector's system
W = T(1:3,3);

% Pm: wrist position
Pm = [Px Py Pz]' - L6*W;

%first joint, two possible solutions admited:
% if q(1) is a solution, then q(1) + pi is also a solution
q1=atan2(Pm(2), Pm(1));

%solve for q2
q2_1=solve_for_theta2(robot, [q1 0 0 0 0 0 0], Pm);
%the other possible solution is q1 + pi
q2_2=solve_for_theta2(robot, [q1+pi 0 0 0 0 0 0], Pm);

%solve for q3
q3_1=solve_for_theta3(robot, [q1 0 0 0 0 0 0], Pm);
%solver for q3 for both cases
q3_2=solve_for_theta3(robot, [q1+pi 0 0 0 0 0 0], Pm);

%the next matrix doubles each column. For each two columns, two different
%configurations for theta4, theta5 and theta6 will be computed. These
%configurations are generally referred as wrist up and wrist down solution
q = [q1         q1         q1        q1       q1+pi   q1+pi   q1+pi   q1+pi;
     q2_1(1)    q2_1(1)    q2_1(2)   q2_1(2)  q2_2(1) q2_2(1) q2_2(2) q2_2(2);
     q3_1(1)    q3_1(1)    q3_1(2)   q3_1(2)  q3_2(1) q3_2(1) q3_2(2) q3_2(2);
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0];

%leave only the real part of the solutions
q=real(q);

%Note that in this robot, the joint q3 has a non-simmetrical range. In this
%case, the joint ranges from 60� to -219�, thus, the typical normalizing
%step is avoided in this angle (the next line is commented). When solving
%for the orientation, the solutions are normalized to the [-pi, pi] range
%only for the theta4, theta5 and theta6 joints.

%normalize q to [-pi, pi]
q(1,:) = normalize(q(1,:));
q(2,:) = normalize(q(2,:));

% solve for the last three joints
% for any of the possible combinations (theta1, theta2, theta3)
for i=1:2:size(q,2),
    qtemp = solve_spherical_wrist2(robot, q(:,i), T, 1,'algebraic'); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i)=qtemp;

    qtemp = solve_spherical_wrist2(robot, q(:,i), T, -1, 'algebraic'); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i+1)=qtemp;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for second joint theta2, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q2 = solve_for_theta2(robot, q, Pm)

%Evaluate the parameters
%theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
%alpha = eval(robot.DH.alpha);

%See geometry
L2=a(2);
L3=d(4);
A2=a(3);

%See geometry of the robot
%compute L4
L4 = sqrt(A2^2 + L3^2);

%The inverse kinematic problem can be solved as in the IRB 140 (for example)

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

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

beta = atan2(-p1(2), p1(1));
gamma = real(acos((L2^2+r^2-L4^2)/(2*r*L2)));

%return two possible solutions
%elbow up and elbow down
%the order here is important and is coordinated with the function
%solve_for_theta3
q2(1) = pi/2 - beta - gamma; %elbow up
q2(2) = pi/2 - beta + gamma; %elbow down

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for third joint theta3, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q3 = solve_for_theta3(robot, q, Pm)

%Evaluate the parameters
%theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
%alpha = eval(robot.DH.alpha);

%See geometry
L2=a(2);
L3=d(4);
A2=a(3);

%See geometry of the robot
%compute L4
L4 = sqrt(A2^2 + L3^2);

%the angle phi is fixed
phi=acos((A2^2+L4^2-L3^2)/(2*A2*L4));

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

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

eta = real(acos((L2^2 + L4^2 - r^2)/(2*L2*L4)));

%return two possible solutions
%elbow up and elbow down solutions
%the order here is important
q3(1) = pi - phi- eta;
q3(2) = pi - phi + eta;

repo: ARTE


text:
function nlogL = DN_nlogL_norm(y)
% DN_nlogL_norm     Negative log likelihood of data coming from a Gaussian distribution.
%
% Fits a Gaussian distribution to the data using the normfit function in
% Matlab's Statistics Toolbox and returns the negative log likelihood of the
% data coming from a Gaussian distribution using the normlike function.
%
%---INPUT:
% y, a vector of data

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function 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/>.
% ------------------------------------------------------------------------------

[muhat, sigmahat] = normfit(y);
nlogL = normlike([muhat, sigmahat],y)/length(y);

% ** Somehow this just scales with length, regardless of the distribution of y?

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [x,y]=limits(a)
% LIMITS returns min & max values of matrix; else scalar value.
%
%   [lo,hi]=LIMITS(a) returns LOw and HIgh values respectively.
%
%   lim=LIMITS(a) returns 1x2 result, where lim = [lo hi] values

% Copyright 2003-2010 The MathWorks, Inc.

if nargin~=1 | nargout>2 %bogus syntax
  error('usage: [lo,hi]=limits(a)')
end

siz=size(a);

if prod(siz)==1 %scalar
  result=a;                         % value
else %matrix
  result=[min(a(:)) max(a(:))];     % limits
end

if nargout==1 %composite result
  x=result;                         % 1x2 vector
elseif nargout==2 %separate results
  x=result(1);                      % two scalars
  y=result(2);
else %no result
  ans=result                        % display answer
end


text:
function ins=ins_time_updata(ins)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INS time update
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2020-2025, by Kai Chen, All rights reserved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ins.Phi=update_trans_mat(ins);

G=zeros(15,15);
G(1:3,1:3)=-ins.Cnb;
G(4:6,4:6)=ins.Cnb;
G(10:12,10:12)=eye(3);
G(13:15,13:15)=eye(3);

Q0=G*ins.Q*G';
P0=ins.P+0.5*Q0;
ins.P=ins.Phi*P0*ins.Phi'+0.5*Q0;

% ins.P=ins.Phi*ins.P*ins.Phi'+ins.Q;

return

repo: GINav


text:
function test_bug1956

% WALLTIME 02:00:00
% MEM 3gb
% DEPENDENCY ft_volumesegment ft_prepare_sourcemodel volumesmooth

mri = ft_read_mri(dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/Subject01.mri'));
mri.coordsys = 'ctf'; % this can also be determined with ft_determine_coordsys

cfg = [];
tpm = ft_volumesegment(cfg,mri);  % gray, white, csf tissue prob. map

% there are some spm functions which change their own input

cfg = [];
cfg.mri = tpm;
grid01 = ft_prepare_sourcemodel(cfg);

cfg = [];
cfg.mri = tpm;
grid02 = ft_prepare_sourcemodel(cfg);

cfg = [];
cfg.mri = tpm;
grid03 = ft_prepare_sourcemodel(cfg);

grid01 = rmfield(grid01, 'cfg');
grid02 = rmfield(grid02, 'cfg');
grid03 = rmfield(grid03, 'cfg');

% make sure that they are the same, and not super-smoohted
assert(isequal(grid01,grid02))
assert(isequal(grid01,grid03))

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


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

%*****************************************************************************80
%
%% SCHUR_BLOCK returns the SCHUR_BLOCK matrix.
%
%  Formula:
%
%    if ( i == j )
%      a(i,j) = x( (i+1)/2 )
%    else ( mod ( i, 2 ) == 1 & j == i + 1 )
%      a(i,j) = y( (i+1)/2 )
%    else ( mod ( i, 2 ) == 0 & j == i - 1 )
%      a(i,j) = -y( (i+1)/2 )
%    else
%      a(i,j) = 0.0
%
%  Example:
%
%    N = 5, X = ( 1, 2, 3 ), Y = ( 4, 5 )
%
%    1  4    0    0    0
%   -4  1    0    0    0
%    0  0    2    5    0
%    0  0   -5    2    0
%    0  0    0    0    3
%
%  Properties:
%
%    A is generally not symmetric: A' /= A.
%
%    A is block diagonal, with the blocks being 2 by 2 or 1 by 1 in size.
%
%    A is in real Schur form.
%
%    The eigenvalues of A are X(I) +/- sqrt ( - 1 ) * Y(I)
%
%    A is tridiagonal.
%
%    A is banded, with bandwidth 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Francoise Chatelin,
%    Section 4.2.7,
%    Eigenvalues of Matrices,
%    John Wiley, 1993.
%
%    Francoise Chatelin, Valerie Fraysse,
%    Qualitative computing: Elements of a theory for finite precision
%    computation, Lecture notes,
%    CERFACS, Toulouse, France and THOMSON-CSF, Orsay, France, June 1993.
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real X( (N+1)/2 ), specifies the diagonal elements
%    of A.
%
%    Input, real Y( N/2 ), specifies the off-diagonal elements
%    of the Schur blocks.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    for j = 1 : n

      if ( i == j )
        a(i,j) = x( floor ( (i+1)/2 ) );
      elseif ( mod ( i, 2 ) == 1 & j == i + 1 )
        a(i,j) = y( floor ( (i+1)/2 ) );
      elseif ( mod ( i, 2 ) == 0 & j == i - 1 )
        a(i,j) = - y( floor ( (i+1)/2 ) );
      else
        a(i,j) = 0.0;
      end

    end
  end

  return
end


text:
function y = nnz(x)
% Octave treats ranges specially to save memory: https://www.gnu.org/software/octave/doc/v4.2.0/Ranges.html
% When applying `nnz` to a range an error occurs. It seems to be solved by linearizing the input (that is,
% applying index `(:)`). That probably turns the range into an actual column vector
y = builtin('nnz', x(:));
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


text:
function [r, s] = randInt(rangeNums, dimensions)
    start = rangeNums(1);
    last = rangeNums(2);
    r = start + (last - start) * rand(dimensions);
    r = round(r);
    s = sum(sum(r));
end

author: anishLearnsToCode

repo: introduction-to-programming-with-matlab


text:
function [ n_data, x, fx ] = sqrt_values ( n_data )

%*****************************************************************************80
%
%% SQRT_VALUES returns some values of the square root function.
%
%  Discussion:
%
%    SQRT(X) = positive real number Y such that Y * Y = X.
%
%    In Mathematica, the function can be evaluated by:
%
%      Sqrt[x]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output real FX, the value of the function.
%
  n_max = 14;

  fx_vec = [ ...
     0.0000000000000000E+00, ...
     0.9000000040950000E-04, ...
     0.3000000000000000E+00, ...
     0.3162277660168379E+00, ...
     0.6324555320336759E+00, ...
     0.1000000000000000E+01, ...
     0.1414213562373095E+01, ...
     0.1732050807568877E+01, ...
     0.1772453850905516E+01, ...
     0.4358898943540674E+01, ...
     0.5385164807134504E+01, ...
     0.8426149773176359E+01, ...
     0.9848857801796105E+01, ...
     0.1111111106055556E+05 ];

  x_vec = [ ...
     0.0000000000000000E+00, ...
     0.8100000073710001E-08, ...
     0.9000000000000000E-01, ...
     0.1000000000000000E+00, ...
     0.4000000000000000E+00, ...
     0.1000000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.3000000000000000E+01, ...
     0.3141592653589793E+01, ...
     0.1900000000000000E+02, ...
     0.2900000000000000E+02, ...
     0.7100000000000000E+02, ...
     0.9700000000000000E+02, ...
     0.1234567890000000E+09 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function [obX,obY,newlp,choice] = mrSelInplane(numofanats,obXM,obYM,lp,prevChoice)
%NAME:   [obX,obY,newlp,choice] = mrSelInplane(numofanats,obXM,obYM,lp,prevChoice)
%AUTHOR:  Poirson
%DATE:	  08.04.96
%PURPOSE: One of a set of routines that allows user
%         to set and select a set of oblique planes in saggital slice.
%         Returns the choice of inplane.
%	  The routines are mrTransInplanes.m, mrRotInplanes.m,
%	  mrClipInplanes.m, mrSelInplane.m, mrSetupInplanes.m,
%         mrSpreadInplanes.m
%HISTORY: Started with mrGetOblPlane from G. Boynton 4/6/96
%NOTES:   Modified 2/25/97 by SPG to return the inplane number so
%         that can use inplane angle as a degree of freedom in the rotation
%         control window

global sagwin
figure(sagwin)

newlp = lp;

% One more pair of points for the perpendicular line
nPtPairs = numofanats + 1;

% See if the user has already been working on some inPlanes
if size(obXM,1) == 0
	disp('You must first create a candidate set of Inplanes');
	disp('Choose Set_Up_Inplanes first');
	return
end

%prevChoice = -99;

xlim=get(gca,'XLim');
ylim=get(gca,'YLim');
xt = diff(xlim)*(-1.0) * 0.8;
yt=ylim(1)+diff(xlim)*[0.05,0.125,0.20,0.275];
msg(1) = text(xt,yt(1),'Select Inplane:');
msg(2) = text(xt,yt(2),'Left = Select Line');
msg(3) = text(xt,yt(3),'Right= Show Selection');

button = 0;
while(button~=3)
	[tempx,tempy,button]=mrGinput(1,'arrow');
	if (button~=3)
		%find the closest line and give it a different color
		for i=1:numofanats
			% p(1) = slope, p(2) = offset
			p = polyfit(obXM(i,:),obYM(i,:),1);
			d(i) = abs((p(1)*tempx)+(-1*tempy)+p(2))/((p(1)^2 + 1)^0.5);
		end
		%find the minimum distance
		choice = find(d == min(d));
		% remove previous choice
		if prevChoice > 0
			% delete(lp(prevChoice))
			newlp(prevChoice)=line(obXM(prevChoice,:),obYM(prevChoice,:),'Color',[0 0 1]);
		end
		% re-draw choice in red
		% delete(lp(choice));
		newlp(choice)=line(obXM(choice,:),obYM(choice,:),'Color',[1 0 0]);
	end
end

% Return this choice
obX = obXM(choice,:);
obY = obYM(choice,:);

for i=1:3
	delete(msg(i));
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function quadrule_test06 ( )

%*****************************************************************************80
%
%% QUADRULE_TEST06 tests CHEBYSHEV2_COMPUTE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2009
%
%  Author:
%
%    John Burkardt
%
  order_max = 4;

  nfunc = func_set ( 'COUNT', 'DUMMY' );

  a = - 1.0;
  b =   1.0;
  nsub = 1;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUADRULE_TEST06\n' );
  fprintf ( 1, '  CHEBYSHEV2_COMPUTE sets up Gauss-Chebyshev type 2 quadrature;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Integration interval is [%f, %f]\n', a, b );
  fprintf ( 1, '  Number of subintervals is %d\n', nsub );
  fprintf ( 1, '  Quadrature order will vary.\n' );
  fprintf ( 1, '  Integrand will vary.\n' );
  fprintf ( 1, '  The weight function is sqrt ( 1 - X**2 )\n' );
  fprintf ( 1, '\n' );

  for ilo = 1 : 5 : nfunc

    ihi = min ( ilo + 4, nfunc );

    fprintf ( 1, '\n' );
    fprintf ( 1, '    ' );
    for i = ilo : ihi
      fprintf ( '%14s', fname(i) );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '\n' );

    for order = 1 : order_max

      fprintf ( 1, '  %2d', order );

      for i = ilo : ihi

        func_set ( 'SET', i );

        [ xtab, weight ] = chebyshev2_compute ( order );

        result = 0.0;
        for j = 1 : order
          result = result + weight(j) * func ( xtab(j) );
        end

        fprintf ( 1, '  %12f', result );

      end

      fprintf ( 1, '\n' );

    end

  end

  return
end


text:
function [ t_start, y_start ] = p11_start ( neqn )

%*****************************************************************************80
%
%% P11_START returns the starting point for problem p11.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_START, Y_START(NEQN,1), the initial data.
%
  y_start = zeros ( neqn, 1 );

  t_start = 0.0;
  y_start(1) = 1.0;
  y_start(2:10) = 0.0;

  return
end


text:
function [x, infos] = online_mu_nmf(V, rank, in_options)
% Online non-negative matrix factorization (Online-NMF) algorithm.
%
% Inputs:
%       matrix      V
%       rank        rank
%       in_options  options
% Output:
%       w           solution of w
%       infos       information
%
% References:
%       S. S. Bucak, B. Gunsel,
%       "Incremental Subspace Learning via Non-negative Matrix Factorization,"
%       Pattern Recognition, 2009.
%
%
% This file is part of NMFLibrary.
%
% Created by H.Kasai and H.Sakai on Feb. 12, 2017
%
% Change log:
%
%       Oct. 27, 2017 (Hiroyuki Kasai): Fixed algorithm.
%
%       May. 20, 2019 (Hiroyuki Kasai): Added initialization module.
%
%       Jul. 12, 2022 (Hiroyuki Kasai): Modified code structures.
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % initialize factors
    init_options = options;
    [init_factors, ~] = generate_init_factors(V, rank, init_options);
    Wt = init_factors.W;
    H = init_factors.H;
    R = init_factors.R;

    % initialize
    method_name = 'Online-MU';
    epoch = 0;
    grad_calc_count = 0;

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    % store initial info
    clear infos;
    [infos, f_val, optgap] = store_nmf_info(V, Wt, H, R, options, [], epoch, grad_calc_count, 0);

    if options.verbose > 1
        fprintf('%s: Epoch = 0000, cost = %.16e, optgap = %.4e\n', method_name, f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main outer loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag] = check_stop_condition(epoch, infos, options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        % Reset sufficient statistic
        At = zeros(m, rank);
        Bt = zeros(rank, rank);

        % main inner loop
        for t = 1 : options.batch_size : n - 1

            % Retrieve vt and ht
            vt = V(:, t:t+options.batch_size-1);
            ht = H(:, t:t+options.batch_size-1);

            % uddate ht
            ht = ht .* (Wt.' * vt) ./ (Wt.' * (Wt * ht));
            ht = ht + (ht<eps) .* eps;

            % update sufficient statistics
            At = At + vt * ht';
            Bt = Bt + ht * ht';

            % update W
            Wt = Wt .* At ./ (Wt * Bt);
            Wt = Wt + (Wt<eps) .* eps;

            % store new h
            H(:,t:t+options.batch_size-1) = ht;

            grad_calc_count = grad_calc_count + m * options.batch_size;
        end

        % measure elapsed time
        elapsed_time = toc(start_time);

        % update epoch
        epoch = epoch + 1;

        % store info
        infos = store_nmf_info(V, Wt, H, R, options, infos, epoch, grad_calc_count, elapsed_time);

        % display info
        display_info(method_name, epoch, infos, options);

    end

    x.W = Wt;
    x.H = H;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
%Chromosome Segregation
%
% @wholeCellModelID Process_ChromosomeSegregation
% @name             Chromosome Segregation
% @description
%   Biology
%   ======================
%   Chromosome segregation in Mycoplasma genitalium has not been well described.
%   Chromosome segregation in M. genitalium is believed to result from two
%   factors
%   - entropically favorable segregation throughout DNA replication
%   - four chromosome segregation proteins which act to decatenate the
%     chromosomes following the completion of replication
%   Little is known about the times at which these two factors contribute to
%   chromosome segregation.
%
%   Bacterial chromosome segregation is believed to be entropically favorable
%   and occur during replication (Bloom et al. 2010, Jun et al. 2006). That is,
%   its believed to be entropically unfavorable for two chromosomes to be
%   located near each other, and so as the replication fork moves, the
%   replicated DNA molecules migrate away from each other toward the poles. We
%   assume that by the end of replication the chromosomes have already migrated
%   toward the poles.
%
%   Four proteins are believed to be involved in chromosome decatenation in M.
%   genitalium:
%   - MG_470_MONOMER  nucleotide binding domain protein CobQ/CobB/MinD/ParA
%   - MG_221_OCTAMER  MraZ
%   - MG_387_MONOMER  GTP binding protein Era
%   - MG_384_MONOMER  GTPase Obg
%   The specific function, kinetics, and metabolic costs of these proteins are
%   not known. Furthemore, because M. genitalium contains a reduced complement
%   of segregation proteins, it is difficult to infer their functions from
%   studies of other bacterial species.
%
%   Knowledge Base
%   ======================
%   The knowledge base contains the value of the gtpCost parameter.
%
%   Representation
%   ======================
%   The substrates and enzymes properties represent the counts of available
%   metabolites and segregation enzymes. gtpCost represents the energetic
%   (GTP) cost of chromosome segregation. The segregation status of the
%   chromosomes is represented as a boolean properties (segregated) of the
%   Chromosome class.
%
%   Initialization
%   ======================
%   Chromosome initializes to a state with 1 chromosome (which obviously
%   hasn't yet segregated from the second chromosome which will be later
%   produced).
%
%   Simulation
%   ======================
%   Because little is known about the molecular biology of chromosome
%   segregation, we have chosen to implement this process as a simple boolean
%   rule: chromosome segregation can occur if:
%   - the chromosome is replicated
%   - the chromosome is properly supercoiled
%   - there is at least one free molecule of each segregation proteins, and
%   - there is at least gtpCost free GTP molecules
%
%   Note, although the Glass et al. gene essentiality study suggests that the
%   cobQ/cobB/minD/parA gene is non-essential, we model this gene as essential
%   because we don't know its specific function and how the other segregation
%   proteins compensate in its absence.
%
%   References
%   ======================
%   1. Bloom, K., Joglekar, A. (2010). Towards building a chromosome segregation
%      machine. Nature 463: 446-456.
%   2. Jun, S., Mulder, B. (2006). Entropy-driven special organization of highly
%      confined polymers: Lessons for the bacterial chromosome. PNAS 103:
%      12388-12393.
%   3. Glass, J.I., Assad-Garcia, N., Alperovich, N., Yooseph, S., Lewis, M.R.,
%      Maruf, M., Hutchison III, C.A., Smith, H.O., Venter, J.C. (2006).
%      Essential genes of a minimal bacterium. PNAS 103: 425-430.
%
% Author: Jayodita Sanghvi, jayodita@stanford.edu
% Author: Jonathan Karr, jkarr@stanford.edu
% Author: Jared Jacobs, jmjacobs@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 7/7/2010
classdef  ChromosomeSegregation < edu.stanford.covert.cell.sim.Process & edu.stanford.covert.cell.sim.ChromosomeProcessAspect
    %property annotations
    properties (Constant)
        optionNames__              = {}; %names of option properties
        fixedConstantNames__       = {   %names of fixed constant properties
            'gtpCost'};
        fittedConstantNames__      = {}; %names of fitted constant properties
        localStateNames__          = {}; %names of simulation state properties redundant with timecourses in this or other processes or the simulation
    end

    %IDs, names, and local indices
    properties
        stimuliWholeCellModelIDs   = {}; %whole cell model IDs of stimuli

        substrateWholeCellModelIDs = {   %whole cell model IDs of substrates
            'GTP';'GDP';'H';'H2O';'PI'};
        substrateIndexs_gtp       = 1; %index of ATP within substrateWholeCellModelIDs
        substrateIndexs_gdp       = 2; %index of ADP within substrateWholeCellModelIDs
        substrateIndexs_hydrogen  = 3; %index of hydrogen within substrateWholeCellModelIDs
        substrateIndexs_water     = 4; %index of water within substrateWholeCellModelIDs
        substrateIndexs_phosphate = 5; %index of inorganic phosphate within substrateWholeCellModelIDs

        enzymeWholeCellModelIDs = {      %whole cell model IDs of enzymes
            'MG_470_MONOMER'       %CobQ/CobB/MinD/ParA nucleotide binding domain
            'MG_221_OCTAMER'       %mraZ protein
            'MG_387_MONOMER'       %GTP-binding protein Era
            'MG_384_MONOMER'       %GTPase1 Obg
            'MG_203_204_TETRAMER'  %DNA topoisomerase IV
            };
        enzymeIndexs_cobQ   = 1;
        enzymeIndexs_mraZ   = 2;
        enzymeIndexs_obg    = 3;
        enzymeIndexs_era    = 4;
        enzymeIndexs_topoIV = 5;
    end

    %fixed biological constants
    properties
        gtpCost %number of GTP required for chromosome segregation (1); true value unknown
    end

    %constructor
    methods
        function this = ChromosomeSegregation(wholeCellModelID, name)
            this = this@edu.stanford.covert.cell.sim.Process(wholeCellModelID, name);
        end
    end

    %communication between process/simulation
    methods
        %set references to state objects
        function storeObjectReferences(this, simulation)
            this.storeObjectReferences@edu.stanford.covert.cell.sim.Process(simulation);
            this.storeObjectReferences@edu.stanford.covert.cell.sim.ChromosomeProcessAspect(simulation);
        end

        %initialize constants
        function initializeConstants(this, knowledgeBase, simulation, varargin)
            this.initializeConstants@edu.stanford.covert.cell.sim.Process(...
                knowledgeBase, simulation, varargin{:});
            this.initializeConstants@edu.stanford.covert.cell.sim.ChromosomeProcessAspect(...
                knowledgeBase, simulation, varargin{:});
        end
    end

    %model
    methods
        %Calculate
        %- contribution to FBA objective
        %- minimum expression consistent with cell cycle length
        function [bmProd, byProd, minEnzExp, maxEnzExp] = calcResourceRequirements_LifeCycle(this, ~, ~)
            %% initialize
            bmProd = zeros(size(this.substrateWholeCellModelIDs));
            byProd = zeros(size(this.substrateWholeCellModelIDs));
            minEnzExp = zeros(size(this.enzymeWholeCellModelIDs));
            maxEnzExp = Inf(size(this.enzymeWholeCellModelIDs));

            %% substrate and byproducts: GTP required to decatenate chromosomes, and its hydrolysis (GTP + H2O ==> GDP + PI + H)
            bmProd(this.substrateIndexs_gtp)       = this.gtpCost;
            bmProd(this.substrateIndexs_water)     = this.gtpCost;
            byProd(this.substrateIndexs_gdp)       = this.gtpCost;
            byProd(this.substrateIndexs_phosphate) = this.gtpCost;
            byProd(this.substrateIndexs_hydrogen)  = this.gtpCost;

            %% Segregation requires at least 1 copy of every enzyme
            minEnzExp(:) = 2;
        end

        %initialization: chromosome state initializes chromosome to state with 1
        %chromosome (which trivially is not segregated)
        function initializeState(~)
        end

        %resource requirements
        function result = calcResourceRequirements_Current(this)
            c = this.chromosome;
            result = zeros(size(this.substrates));
            if ...
                    ~c.segregated && ...
                    collapse(c.polymerizedRegions) == c.nCompartments * c.sequenceLen && ...
                    all(this.enzymes)

                result(this.substrateIndexs_gtp) = this.gtpCost;
                result(this.substrateIndexs_water) = result(this.substrateIndexs_gtp);
            end
        end

        %simulation
        function evolveState(this)
            c = this.chromosome;

            if ...
                    ~c.segregated && ...
                    collapse(c.polymerizedRegions) == c.nCompartments * c.sequenceLen && ...
                    collapse(c.supercoiled) == c.nCompartments && ...
                    all(this.enzymes) && ...
                    this.substrates(this.substrateIndexs_gtp) >= this.gtpCost && ...
                    this.substrates(this.substrateIndexs_water) >= this.gtpCost

                c.segregated = true;

                this.substrates(this.substrateIndexs_gtp)       = this.substrates(this.substrateIndexs_gtp)       - this.gtpCost;
                this.substrates(this.substrateIndexs_water)     = this.substrates(this.substrateIndexs_water)     - this.gtpCost;
                this.substrates(this.substrateIndexs_gdp)       = this.substrates(this.substrateIndexs_gdp)       + this.gtpCost;
                this.substrates(this.substrateIndexs_phosphate) = this.substrates(this.substrateIndexs_phosphate) + this.gtpCost;
                this.substrates(this.substrateIndexs_hydrogen)  = this.substrates(this.substrateIndexs_hydrogen)  + this.gtpCost;
            end
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
classdef metricTree < handle
%%METRICTREE A metric tree class for performing searches in a radius
%            around a given point. If a C++ class interface for the metric
%            tree has been compiled, then its functions will be called in
%            place of the Matlab routines, since the C++ implementation
%            can be significantly faster.
%
%The metric tree data structure is implemented as described in [1] and [2].
%
%This implementation builds the tree from a batch of data all at once. As
%long as there are not numerous points equidistance from a given point, the
%tree will be balanced. The Euclidean distance is used as the metric,
%though the code could be changed to replace it with almost any distance
%metric that satisfies the triangle inequality.
%
%Note that unlike certain metric tree implementations, not all of the nodes
%are held in the leaves.
%
%Modification of the CPPData member of this class or of the error checking
%code in the members can potentially lead to Matlab crashing as CPPData is
%a pointer to data in the C++ implementation.
%
%REFERENCES:
%[1] J. K. Uhlmann. (1991, Nov.) Implementing metric trees to satisfy
%    general proximity/similarity queries. [Online]. Available:
%    http://people.cs.missouri.edu/ uhlmannj/ImplementGH.pdf
%[2] J. K. Uhlmann, "Satisfying general proximity/similarity queries with
%    metric trees," Information Processing Letters, vol. 40, no. 4, pp.
%    175-179, 25 Nov. 1991.
%
%December 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    properties
        DATAIDX%The index of the data at a node.
        innerChild
        outerChild
        innerRadii
        outerRadii
        data%A matrix of the data points.

        CPPData%Only used if an interface to a C++ implementation exists.
    end

    methods
        function newTree=metricTree(k,N)
        %%METRICTREE Construct a new metric tree with space to hold a given
        %            number of nodes.
        %
        %INPUTS: k The dimensionality of the nodes that will be placed in
        %          the tree.
        %        N The number of nodes that the tree will hold.
        %
        %OUTPUTS: newTree A new metricTree instance with the proper amount
        %                 of space.
        %
        %Once a tree has been allocated, it can be initialized using the
        %buildTreeFromBatch method. The size of the data batch given with
        %that method should match the size of the tree allocated here.
            if(exist('metricTreeCPPInt','file'))
                newTree.CPPData=metricTreeCPPInt('metricTreeCPP',k,N);
            else
                newTree.DATAIDX=zeros(N,1);
                newTree.innerChild=zeros(N,1);
                newTree.outerChild=zeros(N,1);
                newTree.innerRadii=zeros(N,1);
                newTree.outerRadii=zeros(N,1);
                newTree.data=[];
            end
        end

        function buildTreeFromBatch(theTree,dataBatch)
        %BUILDTREEFROMBATCH Build a metric tree from a batch of data.
        %
        %INPUTS: theTree The implicitly passed metricTree object.
        %      dataBatch A kXN array of points that are to be stored in
        %                the metric tree. k is the dimensionality of the
        %                points and N is the number of points.
        %
        %The tree is constructed roughly as described in [1].
        %
        %REFERENCES:
        %[1] J. K. Uhlmann. (1991, Nov.) Implementing metric trees to
        %    satisfy general proximity/similarity queries. [Online].
        %    Available: http://people.cs.missouri.edu/ uhlmannj/ImplementGH.pdf

            N=size(dataBatch,2);
            if(exist('metricTreeCPPInt','file'))
                N=metricTreeCPPInt('getN',theTree.CPPData);
                k=metricTreeCPPInt('getk',theTree.CPPData);

                if(N~=size(dataBatch,2)||k~=size(dataBatch,1))
                   error('The data batch given does not match the preallocated size');
                end

                metricTreeCPPInt('buildTreeFromBatch',theTree.CPPData,dataBatch);
            else
                theTree.data=dataBatch;

                %Preallocate an adjacency matrix that will be used in the
                %recursion.
                adjMat=zeros(N,1);
                idx=(1:N)';
                theTree.treeGrow(adjMat,dataBatch,idx,1);
            end
        end

        function [retSet,distSet]=searchRadius(theTree,point,radius)
        %%SEARCHRADIUS Return all of the points in the metric tree that are
        %              located within a given search radius of within a
        %              given radius about a given point.
        %
        %INPUTS: theTree The implicitly passed metricTree object.
        %          point A kXm matrix of k-dimensional points about which
        %                range searches will be performed.
        %         radius An mX1 vector of the search radii about each
        %                point.
        %
        %OUTPUTS: retSet An instance of the ClusterSet class containing all
        %                of the items within the given radius for each
        %                point. If no points in the tree fall into a given
        %                region, then the corresponding clusterSize element
        %                in the ClusterSet will be zero.
        %        distSet An instance of the clusterSet class containing all
        %                of the distances associated with the found items
        %                in retSet

            m1=size(point,2);
            m2=size(radius,1);

            if(m1~=m2)
                error('Inconsistent numbers of points and radii given.');
            end

            if(exist('metricTreeCPPInt','file'))
                [retSet,distSet]=metricTreeCPPInt('searchRadius',theTree.CPPData,point,radius);
                %Convert indices to Matlab indices
                retSet.clusterEls=retSet.clusterEls+1;
            else
                %Create the new ClusterSet instances with m1 empty clusters.
                retSet=ClusterSet([],zeros(m1,1),zeros(m1,1));
                distSet=ClusterSet([],zeros(m1,1),zeros(m1,1));

                %Add the clusters
                for curPoint=1:m1
                    [idxRange,distList]=theTree.searchRadRecur(point(:,curPoint),radius(curPoint),1);

                    retSet.clusterSizes(curPoint)=size(idxRange,1);
                    retSet.clusterEls=[retSet.clusterEls;idxRange];
                    distSet.clusterSizes(curPoint)=size(idxRange,1);
                    distSet.clusterEls=[distSet.clusterEls;distList];
                    if(curPoint<m1)
                        retSet.offsetArray(curPoint+1)=retSet.offsetArray(curPoint)+retSet.clusterSizes(curPoint);
                        distSet.offsetArray(curPoint+1)=distSet.offsetArray(curPoint)+distSet.clusterSizes(curPoint);
                    end
                end
            end
        end

        function display(theTree)
        %%DISPLAY Display information about the tree, including whether the
        %         C++ implementation (wrapped by a Matlab class) or the
        %         Matlab-only implementation is being used.

            if(exist('metricTreeCPPInt','file'))
                display('A metric tree as implemented via a mex wrapper around a C++ class.')
            else
                display('A metric tree as implemented as a Matlab class.')
            end
            theTree.disp();
        end

        function copyDataFromCPP(theTree)
    %%COPYDATAFROMCPP Copy the data that makes up the structure of the tree
    %                 out of C++ into the variables in the Matlab kdTree
    %                 class. This can be useful when debugging the C++
    %                 implementation. This function just raises an error
    %                 when executed when using the Matlab implementation.

            if(exist('metricTreeCPPInt','file'))
                [theTree.DATAIDX,theTree.innerChild,theTree.outerChild,theTree.innerRadii,theTree.outerRadii,theTree.data]=metricTreeCPPInt('getAllData',theTree.CPPData);
            else
                error('This metric tree class instance does not stem from a C++ class');
            end
        end

        function delete(theTree)
        %%DELETE The destructor method. This method is used when the metric
        %        tree is implemented as a C++ class. This method prevents a
        %        memory leak.

            if(exist('metricTreeCPPInt','file'))
                metricTreeCPPInt('~metricTreeCPP',theTree.CPPData);
            end
        end
    end

    methods(Access=private)
        function [idxRange,distList]=searchRadRecur(theTree,point,r2,curNode)
        %%SEARCHRADRECUR A recursion function for performing a search of a
        %                particular radius about a point.

            distCur=dist(point,theTree.data(:,theTree.DATAIDX(curNode)));
            if(distCur<=r2)
                idxRange=theTree.DATAIDX(curNode);
                distList=distCur;
            else
                idxRange=[];
                distList=[];
            end

            if(distCur+r2>=theTree.outerRadii(curNode)&&theTree.outerChild(curNode)~=-1)
                [idxRest,distSqRest]=theTree.searchRadRecur(point,r2,theTree.outerChild(curNode));
                idxRange=[idxRange;idxRest];
                distList=[distList;distSqRest];
            end

            if(distCur-r2<=theTree.innerRadii(curNode)&&theTree.innerChild(curNode)~=-1)
                [idxRest,distSqRest]=theTree.searchRadRecur(point,r2,theTree.innerChild(curNode));
                idxRange=[idxRange;idxRest];
                distList=[distList;distSqRest];
            end
        end

        function nextFreeNode=treeGrow(theTree,adjMat,dataBatch,idx,curNode)
        %%TREEGROW A recursion function for creating a new metric tree.

            NSubTree=size(dataBatch,2);

            %The first index in idx is added to the tree.
            theTree.DATAIDX(curNode)=idx(1);

            %If this is a leaf node.
            if(NSubTree==1)
                theTree.innerChild(curNode)=-1;
                theTree.outerChild(curNode)=-1;
                theTree.innerRadii(curNode)=-Inf;
                theTree.outerRadii(curNode)=Inf;
                nextFreeNode=curNode+1;
                return
            end

            %Fill the NCurX1 subvector of adjMat with all pairwise
            %distances from the current point to the other points in idx.
            NSubTree=NSubTree-1;
            for curPoint=1:NSubTree
                adjMat(curPoint)=dist(dataBatch(:,1),dataBatch(:,curPoint+1));
            end

            %Next, sort the points
            [adjMat(1:NSubTree), sortIdx]=sort(adjMat(1:NSubTree));

            idx=idx(2:end);
            idx=idx(sortIdx);

            dataBatch=dataBatch(:,2:end);
            dataBatch=dataBatch(:,sortIdx);

            %Find the first occurance of the median distance.
            midIdx=findFirstMax(adjMat(1:ceil((NSubTree+1)/2)));

            %The inner and outer radii of how things are split in terms of
            %distances from the given point.
            theTree.outerRadii(curNode)=adjMat(midIdx);
            if(midIdx>1)
                theTree.innerRadii(curNode)=adjMat(midIdx-1);
            else
                theTree.innerRadii(curNode)=-Inf;
            end

            %Continue the recursion.
            theTree.outerChild(curNode)=curNode+1;
            nextFreeNode=theTree.treeGrow(adjMat,dataBatch(:,midIdx:end),idx(midIdx:end),curNode+1);

            %If a full partitioning of the nodes took place
            if(midIdx>1)
                theTree.innerChild(curNode)=nextFreeNode;
                nextFreeNode=theTree.treeGrow(adjMat,dataBatch(:,1:(midIdx-1)),idx(1:(midIdx-1)),nextFreeNode);
            else
                theTree.innerChild(curNode)=-1;
            end
        end
    end
end

function val= dist(a,b)
%%DIST Evaluate the Euclidean distance between two vectors.

    %The square root must be used to obey the triangle inequality.
    diff=a-b;
    val=sqrt(diff'*diff);
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [imdbs, roidbs] = subsample_images(imdbs, roidbs, max_num_neg_images, seed)

if ~exist('seed', 'var')
  seed = 6;
end

% class_num = cellfun(@(x) length(x.class_ids), imdbs, 'UniformOutput', true);
% assert(length(unique(class_num)) == 1);
% class_num = unique(class_num);

rois = cellfun(@(x) x.rois(:), roidbs, 'UniformOutput', false);
rois_combine = cell2mat(rois(:));

% fix the random seed for repeatability
prev_rng = seed_rand(seed);
inds = randperm(length(rois_combine), max_num_neg_images);
inds = sort(inds);

img_idx_start = 1;
for i = 1:length(imdbs)
    imdb_img_num = length(imdbs{i}.image_ids);
    img_idx_end = img_idx_start + imdb_img_num - 1;
    inds_start = find(inds >= img_idx_start, 1, 'first');
    inds_end = find(inds <= img_idx_end, 1, 'last');

    inds_sub = inds(inds_start:inds_end);
    inds_sub = inds_sub - img_idx_start + 1;

    imdbs{i}.image_ids = imdbs{i}.image_ids(inds_sub);
    imdbs{i}.sizes = imdbs{i}.sizes(inds_sub, :);
    if isfield(imdbs{i}, 'image_dir')
        imdbs{i}.image_at = @(x) ...
            sprintf('%s/%s.%s', imdbs{i}.image_dir, imdbs{i}.image_ids{x}, imdbs{i}.extension);
    else
       imdbs{i}.image_at = @(x) ...
            sprintf('%s/%s.%s', imdbs{i}.imagedir, imdbs{i}.image_ids{x}, imdbs{i}.extension);
    end
    roidbs{i}.rois = roidbs{i}.rois(inds_sub);

    img_idx_start = img_idx_start + imdb_img_num;
end

% restore previous rng
rng(prev_rng);

end

author: ShaoqingRen


text:
function test_errors
%TEST_ERRORS tests error handling for the factorize object methods
%
% Example
%   test_errors
%
% See also test_all, factorize1, factorize.

% Copyright 2009, Timothy A. Davis, University of Florida

% the matrix A must be square for factorize1
fprintf ('\nTesting error-handling, error messages expected:\n') ;
ok = true ;
A = rand (3,2) ;
try
    F = factorize1 (A) ;                                                %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

% the matrix A must be 2D for factorize
A = ones (2,2,2) ;
try
    F = factorize (A) ;                                                 %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

% the matrix A must be full-rank
A = ones (3) ;
try
    F = factorize (A) ;                                                 %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

try
    F = factorize1 (A) ;                                                %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

A = ones (3,2) ;
try
    F = factorize (A) ;                                                 %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

A = ones (2,3) ;
try
    F = factorize (A) ;                                                 %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

% cannot use cell indexing
A = rand (3,2) ;
F = factorize (A) ;
try
    C = F {1} ;                                                         %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

% cannot use linear indexing of the inverse
S = inverse (A) ;
try
    C = S (1) ;                                                         %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

% invalid indexing
try
    C = F (1,1).L ;                                                     %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

% invalid indexing
try
    C = F.L (1,1).stuff ;                                               %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

% non-existent field
try
    C = F.junk ;                                                        %#ok
    ok = false ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

% can only update/downdate a dense Cholesky factorization
A = rand (2) ;
F = factorize (A) ;
w = rand (2,1) ;
try
    F = F + w ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end
try
    F = F - w ;
catch me
    fprintf ('Expected error: [%s]\n', me.message) ;
end

if (~ok)
    error ('error-handling failed') ;
end
fprintf ('\nAll error-handing tests passed\n') ;


text:
function [ p2, dp2, p1 ] = gegenbauer_recur ( x, n, alpha, c )

%*****************************************************************************80
%
%% GEGENBAUER_RECUR finds the value and derivative of a Gegenbauer polynomial.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2008
%
%  Author:
%
%    Original FORTRAN77 version by Arthur Stroud, Don Secrest,
%    MATLAB version by John Burkardt,
%
%  Reference:
%
%    Arthur Stroud, Don Secrest,
%    Gaussian Quadrature Formulas,
%    Prentice Hall, 1966,
%    LC: QA299.4G3S7.
%
%  Parameters:
%
%    Input, real X, the point at which polynomials are evaluated.
%
%    Input, integer N, the order of the polynomial to be computed.
%
%    Input, real ALPHA, the exponent of (1-X^2) in the quadrature rule.
%
%    Input, real C(N), the recursion coefficients.
%
%    Output, real P2, the value of J(N)(X).
%
%    Output, real DP2, the value of J'(N)(X).
%
%    Output, real P1, the value of J(N-1)(X).
%
  p1 = 1.0;
  dp1 = 0.0;

  p2 = x;
  dp2 = 1.0;

  for i = 2 : n

    p0 = p1;
    dp0 = dp1;

    p1 = p2;
    dp1 = dp2;

    p2 = x * p1 - c(i) * p0;
    dp2 = x * dp1 + p1 - c(i) * dp0;

  end

  return
end


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

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

t = 10;                                      % maturity
f0 = 0.03;                                  % spot value

figname = 'Risk Neutral Density - CEV';

x = 0:0.001:2;                            % range

sigma_base = 0.2;                           % volatility base scenario
beta_base = 0.5;                            % CEV exponent base scenario

y_base = pcev(t,x,f0,sigma_base, beta_base);
legendname_base = 'Base';
%% Changing sigma_base
sigma_low = 0.15;
sigma_high = 0.25;

y_low = pcev(t,x,f0,sigma_low, beta_base);
y_high = pcev(t,x,f0,sigma_high, beta_base);

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

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

%% Changing beta_base
beta_low = 0.3;
beta_high = 0.7;

y_low = pcev(t,x,f0,sigma_base, beta_low);
y_high = pcev(t,x,f0,sigma_base, beta_high);

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

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


text:
function test_pull360

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY

mri				= ft_read_mri('single_subj_T1_1mm.nii', 'datatype', 'nifti');
mri.transform	= [-1,0,0,95.7;0,1,0,-122.5;0,0,1,-132.5;0,0,0,1]; % this will require intervention by getinside in ft_prepare_headmodels

cfg                         = [];
cfg.mri                     = mri;
cfg.sourcemodel.resolution         = 10;       % in mm
cfg.sourcemodel.unit               = 'mm';
ft_prepare_sourcemodel(cfg);


text:
function [varargout] = likGauss(hyp, y, mu, s2, inf, i)

% likGauss - Gaussian likelihood function for regression. The expression for the
% likelihood is
%   likGauss(t) = exp(-(t-y)^2/2*sn^2) / sqrt(2*pi*sn^2),
% where y is the mean and sn is the standard deviation.
%
% The hyperparameters are:
%
% hyp = [  log(sn)  ]
%
% Several modes are provided, for computing likelihoods, derivatives and moments
% respectively, see likFunctions.m for the details. In general, care is taken
% to avoid numerical issues when the arguments are extreme.
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch, 2014-03-04.
%                                      File automatically generated using noweb.
%
% See also LIKFUNCTIONS.M.

if nargin<3, varargout = {'1'}; return; end   % report number of hyperparameters

sn2 = exp(2*hyp);

if nargin<5                              % prediction mode if inf is not present
  if numel(y)==0,  y = zeros(size(mu)); end
  s2zero = 1; if nargin>3, if norm(s2)>0, s2zero = 0; end, end         % s2==0 ?
  if s2zero                                                    % log probability
    lp = -(y-mu).^2./sn2/2-log(2*pi*sn2)/2; s2 = 0;
  else
    lp = likGauss(hyp, y, mu, s2, 'infEP');                         % prediction
  end
  ymu = {}; ys2 = {};
  if nargout>1
    ymu = mu;                                                   % first y moment
    if nargout>2
      ys2 = s2 + sn2;                                          % second y moment
    end
  end
  varargout = {lp,ymu,ys2};
else
  switch inf
  case 'infLaplace'
    if nargin<6                                             % no derivative mode
      if numel(y)==0, y=0; end
      ymmu = y-mu; dlp = {}; d2lp = {}; d3lp = {};
      lp = -ymmu.^2/(2*sn2) - log(2*pi*sn2)/2;
      if nargout>1
        dlp = ymmu/sn2;                      % dlp, derivative of log likelihood
        if nargout>2                    % d2lp, 2nd derivative of log likelihood
          d2lp = -ones(size(ymmu))/sn2;
          if nargout>3                  % d3lp, 3rd derivative of log likelihood
            d3lp = zeros(size(ymmu));
          end
        end
      end
      varargout = {lp,dlp,d2lp,d3lp};
    else                                                       % derivative mode
      lp_dhyp = (y-mu).^2/sn2 - 1;  % derivative of log likelihood w.r.t. hypers
      dlp_dhyp = 2*(mu-y)/sn2;                               % first derivative,
      d2lp_dhyp = 2*ones(size(mu))/sn2;   % and also of the second mu derivative
      varargout = {lp_dhyp,dlp_dhyp,d2lp_dhyp};
    end

  case 'infEP'
    if nargin<6                                             % no derivative mode
      lZ = -(y-mu).^2./(sn2+s2)/2 - log(2*pi*(sn2+s2))/2;    % log part function
      dlZ = {}; d2lZ = {};
      if nargout>1
        dlZ  = (y-mu)./(sn2+s2);                    % 1st derivative w.r.t. mean
        if nargout>2
          d2lZ = -1./(sn2+s2);                      % 2nd derivative w.r.t. mean
        end
      end
      varargout = {lZ,dlZ,d2lZ};
    else                                                       % derivative mode
      dlZhyp = ((y-mu).^2./(sn2+s2)-1) ./ (1+s2./sn2);   % deriv. w.r.t. hyp.lik
      varargout = {dlZhyp};
    end

  case 'infVB'
    % variational lower site bound
    % t(s) = exp(-(y-s)^2/2sn2)/sqrt(2*pi*sn2)
    % the bound has the form: (b+z/ga)*f - f.^2/(2*ga) - h(ga)/2
    n = numel(s2); b = zeros(n,1); y = y.*ones(n,1); z = y;
    varargout = {b,z};
  end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO