text:
%ProcessTestCase
% Base class for whole cell process test classes. This purposes of this
% class are to provide functions for
% - loading the process test fixture prior to the execution of each test.
%   This is automatically called by the setUp function.
% - testing for run-time errors in these process methods:
%   - allocateMemoryForState
%   - computeBiomassProductionByproducts
%   - computeMinimumEnzymeExpression
%   - resourceRequirement
%   - initializeState
%   - evolveState
% - testing that the process obeys mass balance
% - testing gene essentiality with regard to the process's primary function
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 9/10/2010
classdef ProcessTestCase < TestCase
    properties
        process               %process object
        processClass          %full class name of process
        processClassName      %short class name of process
        processBaseLocation   %base directory of process
    end

    %constructor
    methods
        function this = ProcessTestCase(name)
            %parent class constructor
            this = this@TestCase(name);

            %get class name of process
            this.processClass     = edu.stanford.covert.util.substr(class(this), 1, -5);
            this.processClassName = this.processClass(find(this.processClass == '.', 1, 'last') + 1:end);

            %find parent folder of process -- used to store test fixture and
            %expected outputs
            locationParts = strfind(this.Location, filesep);
            if this.Location(locationParts(end - 1) + 1) == '@'
                this.processBaseLocation = this.Location(1:locationParts(end-1));
            else
                this.processBaseLocation = this.Location(1:locationParts(end));
            end
        end
    end

    %setUp and tearDown
    methods
        %Called prior to the execution of each test case.
        %Performs three functions:
        % - instantiate process
        % - initialize process using test fixture
        % - seed random number generator
        function setUp(this)
            constructor = str2func(this.processClass);
            this.process = constructor([], this.processClassName);
            this.loadTestFixture();
            this.process.verbosity = 0;
            this.process.seed = 1;
            this.process.seedRandStream();
        end

        %Called following the execution of each test case.
        function tearDown(this)
            this.process = []; %free memory
        end

        %Loads knowledge base test fixture into process. Called by setUp prior
        %to the execution of each test. You may wish to override this method to
        %set additional properties.
        function loadTestFixture(this)
            this.process = edu.stanford.covert.cell.sim.ProcessFixture.load(this.process);
        end
    end

    %Suggested tests for subclasses to implement:
    % - Demonstrate main purpose of process using minimal, hard-coded fake data.
    % - Withhold one or all stimuli.
    % - Withhold one or all substrates.
    % - Withhold one or all enzymes.
    % - Predict which genes are essential to the process (testGeneEssentiality).

    %tests for absence of run-time errors
    methods (Sealed)
        function testAllocateMemoryForState(this)
            this.process.allocateMemoryForState(1);  %no run-time errors
        end

        function testCalcResourceRequirements_LifeCycle(this)
            %references
            sim = edu.stanford.covert.cell.sim.SimulationFixture.load([], true);
            m = sim.process(this.process.wholeCellModelID(9:end));
            g = sim.gene;
            t = sim.state('Time');
            r = sim.state('Rna');
            pm = sim.state('ProteinMonomer');
            pc = sim.state('ProteinComplex');

            %process local state
            m.stimuli = this.process.stimuli;
            m.substrates = this.process.substrates;
            m.enzymes = this.process.enzymes;
            m.boundEnzymes = this.process.boundEnzymes;
            this.process = m;

            %constants
            constants = edu.stanford.covert.util.StructUtil.catstruct(...
                sim.getFixedConstants(), sim.getFittedConstants());
            constants.processes.DNADamage.substrateGlobalIndexs = sim.process('DNADamage').substrateGlobalIndexs;
            constants.processes.DNADamage.substrateStimulusGlobalIndexs = sim.process('DNADamage').substrateStimulusGlobalIndexs;
            constants.processes.DNADamage.substrateStimulusCompartmentIndexs = sim.process('DNADamage').substrateStimulusCompartmentIndexs;
            constants.processes.DNADamage.substrateStimulusGlobalCompartmentIndexs = sim.process('DNADamage').substrateStimulusGlobalCompartmentIndexs;
            constants.processes.DNADamage.substrateStimulusLocalIndexs = sim.process('DNADamage').substrateStimulusLocalIndexs;

            %RNA and protein
            rnas = sum( r.counts( r.matureIndexs, :), 2) + ...
                sum(r.counts(r.boundIndexs, :), 2) + ...
                sum(r.counts(r.aminoacylatedIndexs, :), 2) + ...
                sum(r.counts(r.damagedIndexs, :), 2);
            mons = ...
                sum(pm.counts(pm.matureIndexs, :), 2) + ...
                sum(pm.counts(pm.boundIndexs, :), 2) + ...
                sum(pm.counts(pm.damagedIndexs, :), 2) + ...
                sum(pm.counts(pm.inactivatedIndexs, :), 2);
            cpxs = sum(pc.counts(pc.matureIndexs, :), 2) + ...
                sum(pc.counts(pc.boundIndexs, :), 2) + ...
                sum(pc.counts(pc.damagedIndexs, :), 2) + ...
                sum(pc.counts(pc.inactivatedIndexs, :), 2);

            rnaInCpxs = zeros(size(rnas));
            rnaInCpxs(setdiff(1:end, r.matureMRNAIndexs)) = ...
                sum(pc.proteinComplexComposition(setdiff(1:end, g.mRNAIndexs), :, :), 3) * cpxs;
            monInCpxs = sum(pc.proteinComplexComposition(g.mRNAIndexs, :, :), 3) * cpxs;

            cpxDecays = cpxs .* pc.decayRates(pc.matureIndexs) * t.cellCycleLength / log(2);
            rnaDecays = rnas .*  r.decayRates( r.matureIndexs) * t.cellCycleLength / log(2);
            rnaDecays(setdiff(1:end, r.matureMRNAIndexs)) = ...
                rnaDecays(setdiff(1:end, r.matureMRNAIndexs)) + ...
                sum(pc.proteinComplexComposition(setdiff(1:end, g.mRNAIndexs), :, :), 3) * cpxDecays;
            monDecays = mons .* pm.decayRates(pm.matureIndexs) * t.cellCycleLength / log(2) + ...
                sum(pc.proteinComplexComposition(g.mRNAIndexs, :, :), 3) * cpxDecays;

            states = struct;

            states.rnas0               = rnas;
            states.rnas                = rnas * t.cellCycleLength / log(2);
            states.rnaDecays           = rnaDecays;
            states.rnaDecays0          = rnaDecays * log(2) / t.cellCycleLength;
            states.rnaProductions      = rnas + rnaInCpxs + rnaDecays;
            states.rnaProductions0     = (rnas + rnaInCpxs) * t.cellCycleLength / log(2) + (rnas .* r.decayRates(r.matureIndexs));
            states.rnaProductions0(setdiff(1:end, r.matureMRNAIndexs)) = ...
                states.rnaProductions0(setdiff(1:end, r.matureMRNAIndexs)) + ...
                sum(pc.proteinComplexComposition(setdiff(1:end, g.mRNAIndexs), :, :), 3) * (cpxs .* pc.decayRates(pc.matureIndexs));

            states.monomers0           = mons;
            states.monomers            = mons * t.cellCycleLength / log(2);
            states.monomerDecays       = monDecays;
            states.monomerDecays0      = monDecays * log(2) / t.cellCycleLength;
            states.monomerProductions  = mons + monInCpxs + monDecays;
            states.monomerProductions0 = (mons + monInCpxs) * t.cellCycleLength / log(2) + (mons .* pm.decayRates(pm.matureIndexs));
            states.monomerProductions0 = ...
                states.monomerProductions0 + ...
                sum(pc.proteinComplexComposition(g.mRNAIndexs, :, :), 3) * (cpxs .* pc.decayRates(pc.matureIndexs));

            states.complexs0           = cpxs;
            states.complexs            = cpxs * t.cellCycleLength / log(2);
            states.complexDecays       = cpxDecays;
            states.complexDecays0      = cpxDecays * log(2) / t.cellCycleLength;
            states.complexProductions  = cpxs + cpxDecays;
            states.complexProductions0 = cpxs * t.cellCycleLength / log(2) + (cpxs .* pc.decayRates(pc.matureIndexs));

            %calculate life cycle resource requirements
            [bmProd, byProd, minEnzExp, maxEnzExp] = m.calcResourceRequirements_LifeCycle(constants, states);

            %assertions
            assertEqual(size(m.substrateWholeCellModelIDs), size(bmProd));
            assertEqual(size(m.substrateWholeCellModelIDs), size(byProd));
            assertEqual(size(m.enzymeWholeCellModelIDs), size(minEnzExp));
            assertEqual(size(m.enzymeWholeCellModelIDs), size(maxEnzExp));
        end

        function testCalcResourceRequirements_Current(this)
            this.process.calcResourceRequirements_Current();  %no run-time errors
        end

        function testInitializeState(this)
            this.process.initializeState();  %no run-time errors
        end

        function testEvolveState(this)
            m = this.process;
            m.evolveState();  %no run-time errors
            validateattributes(m.substrates(setdiff(1:end, m.substrateStimulusLocalIndexs), :), {'numeric'}, {'real', 'nonnegative', 'integer'});
            validateattributes(m.enzymes, {'numeric'}, {'real', 'nonnegative', 'integer'});
            validateattributes(m.boundEnzymes, {'numeric'}, {'real', 'nonnegative', 'integer'});
        end
    end

    %other tests provided by this class
    methods (Sealed)
        %Test that process obeys mass balance
        function testMassBalance(this)
            import edu.stanford.covert.util.ConstantUtil;
            m = this.process;
            initialWeight = sum(m.dryWeight);

            m.evolveState();
            assertElementsAlmostEqual(...
                initialWeight * ConstantUtil.nAvogadro, sum(m.dryWeight) * ConstantUtil.nAvogadro, ...
                'relative', 1e-8);
        end

        %Test that process uses rand streams correctly
        %- only changes state of process (and possibly state) rand
        %  stream(s)
        %- doesn't change state of global rand stream
        function testRandStreamUsage(this)
            %% hold onto state of default stream
            defaultStream = RandStream.getDefaultStream;
            defaultStreamType = defaultStream.Type;
            defaultStreamSeed = defaultStream.Seed;
            defaultStreamNumStreams = defaultStream.NumStreams;
            defaultStreamStreamIndex = defaultStream.StreamIndex;
            defaultStreamState = defaultStream.State;
            defaultStreamSubstram = defaultStream.Substream;
            defaultStreamRandnAlg = defaultStream.RandnAlg;
            defaultStreamAntithetic = defaultStream.Antithetic;
            defaultStreamFullPrecision = defaultStream.FullPrecision;

            %% simulate
            %run-1
            this.setUp();
            m = this.process;
            m.evolveState();
            streams1 = cell(1+numel(m.states), 1);
            streams1{1} = this.process.randStream;
            for i = 1:numel(m.states)
                streams1{i+1} = m.states{i}.randStream;
            end
            this.tearDown();

            %run-2
            this.setUp();
            m = this.process;
            m.evolveState();
            streams2 = cell(1+numel(m.states), 1);
            streams2{1} = this.process.randStream;
            for i = 1:numel(m.states)
                streams2{i+1} = m.states{i}.randStream;
            end
            this.tearDown();

            %% assert state of process and state random streams updated the same in both runs
            for i = 1:numel(streams1)
                assertTrue(streams1{i} == streams2{i});
            end

            %% assert state of default stream unchanged
            defaultStream = RandStream.getDefaultStream;
            assertEqual(defaultStreamType, defaultStream.Type);
            assertEqual(defaultStreamSeed, defaultStream.Seed);
            assertEqual(defaultStreamNumStreams, defaultStream.NumStreams);
            assertEqual(defaultStreamStreamIndex, defaultStream.StreamIndex);
            assertEqual(defaultStreamState, defaultStream.State);
            assertEqual(defaultStreamSubstram, defaultStream.Substream);
            assertEqual(defaultStreamRandnAlg, defaultStream.RandnAlg);
            assertEqual(defaultStreamAntithetic, defaultStream.Antithetic);
            assertEqual(defaultStreamFullPrecision, defaultStream.FullPrecision);
        end
    end

    %test helpers
    methods (Sealed)
        function helpTestGeneEssentiality(...
                this, expectedEssentialGenes, isFunctioningOkay, options)
            %process options
            if ~exist('options','var')
                options = struct;
            end
            lengthSec = this.optionVal(options, 'lengthSec', 1);
            knockoutStimuli = this.optionVal(options, 'knockoutStimuli', 1);
            knockoutSubstrates = this.optionVal(options, 'knockoutSubstrates', 1);
            knockoutEnzymes = this.optionVal(options, 'knockoutEnzymes', 1);

            m = this.process;

            base_globalStates = struct;
            for i = 1:numel(m.states)
                base_globalStates.(m.states{i}.wholeCellModelID) = struct;
                for j = 1:numel(m.states{i}.stateNames)
                    base_globalStates.(m.states{i}.wholeCellModelID).(m.states{i}.stateNames{j}) = ...
                        m.states{i}.(m.states{i}.stateNames{j});
                end
            end

            base_localStates = struct;
            for i = 1:numel(m.localStateNames)
                base_localStates.(m.localStateNames{i}) = m.(m.localStateNames{i});
            end

            baseStates = edu.stanford.covert.util.StructUtil.catstruct(base_globalStates, base_localStates);

            %compute gene composition of stimuli, substrates, enzymes
            geneWholeCellModelIDs    = m.gene.wholeCellModelIDs;
            stimuliGeneComposition   = m.stimuliGeneComposition();
            substrateGeneComposition = m.substrateGeneComposition();
            enzymeGeneComposition    = m.enzymeGeneComposition();

            includedGenes = false(size(geneWholeCellModelIDs));
            if knockoutStimuli
                includedGenes = includedGenes | any(stimuliGeneComposition,2);
            end
            if knockoutSubstrates
                includedGenes = includedGenes | any(substrateGeneComposition,2);
            end
            if knockoutEnzymes
                includedGenes = includedGenes | any(enzymeGeneComposition,2);
            end
            includedGenes = find(includedGenes);

            geneWholeCellModelIDs    = geneWholeCellModelIDs(includedGenes);
            stimuliGeneComposition   = stimuliGeneComposition(includedGenes, :);
            substrateGeneComposition = substrateGeneComposition(includedGenes, :);
            enzymeGeneComposition    = enzymeGeneComposition(includedGenes, :);

            %verify all expected essential genes included in test
            assertTrue(...
                all(ismember(expectedEssentialGenes, geneWholeCellModelIDs)), ...
                'not all expected essential genes are included in this test');

            %verify first that the process functions properly with all genes intact
            for j = 1:lengthSec
                m.evolveState();
                if isFunctioningOkay(m, baseStates)
                    break;
                end
            end
            assertTrue(...
                isFunctioningOkay(m, baseStates),...
                'process is not functioning properly with all genes intact');

            mispredictedAsEssential = false(size(geneWholeCellModelIDs));
            mispredictedAsNotEssential = false(size(geneWholeCellModelIDs));

            %for each gene
            for i = 1:length(geneWholeCellModelIDs)
                %reset process state
                m.seedRandStream();

                for j = 1:numel(m.states)
                    for k = 1:numel(m.states{j}.stateNames)
                        m.states{j}.(m.states{j}.stateNames{k}) = ...
                            base_globalStates.(m.states{j}.wholeCellModelID).(m.states{j}.stateNames{k});
                    end
                end

                for j = 1:numel(m.localStateNames)
                    m.(m.localStateNames{j}) = base_localStates.(m.localStateNames{j});
                end

                %knock out gene
                m.stimuli(     stimuliGeneComposition(i,:) >0, :) = 0;
                m.substrates(substrateGeneComposition(i,:) >0, :) = 0;
                m.enzymes(      enzymeGeneComposition(i,:) >0, :) = 0;
                m.boundEnzymes( enzymeGeneComposition(i,:) >0, :) = 0;

                this.helpKnockoutGene(geneWholeCellModelIDs{i});

                initial_localStates = base_localStates;
                initial_localStates.stimuli      = m.stimuli;
                initial_localStates.substrates   = m.substrates;
                initial_localStates.enzymes      = m.enzymes;
                initial_localStates.boundEnzymes = m.boundEnzymes;

                initialStates = edu.stanford.covert.util.StructUtil.catstruct(base_globalStates, initial_localStates);

                expectedEssential = ...
                    any(strcmp(geneWholeCellModelIDs{i}, expectedEssentialGenes));

                %execute process
                for j = 1:lengthSec
                    m.evolveState();
                    if isFunctioningOkay(m, initialStates)
                        if expectedEssential
                            mispredictedAsEssential(i) = true;
                        end
                        break;
                    end
                end

                %determine whether process correctly predicted gene essentiality
                if isFunctioningOkay(m, initialStates)
                    if expectedEssential
                        mispredictedAsEssential(i) = true;
                    end
                else
                    if ~expectedEssential
                        mispredictedAsNotEssential(i) = true;
                    end
                end
            end

            assertTrue(all(~mispredictedAsEssential), ...
                ['genes are not essential: '...
                strjoin(', ', geneWholeCellModelIDs{mispredictedAsEssential})]);
            assertTrue(all(~mispredictedAsNotEssential), ...
                ['genes are essential: '...
                strjoin(', ', geneWholeCellModelIDs{mispredictedAsNotEssential})]);
        end

        function result = optionVal(~, options, name, defaultValue)
            if isfield(options, name)
                result = options.(name);
            else
                result = defaultValue;
            end
        end
    end

    %helper methods optionally overridden by subclasses
    methods
        function helpKnockoutGene(this, geneWholeCellModelID) %#ok<MANU,INUSD>
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function results = run_DeepSTRCF(seq, res_path, bSaveImage)

% Feature specific parameters
hog_params.cell_size = 4;
hog_params.compressed_dim = 10;
hog_params.nDim = 31;

grayscale_params.colorspace='gray';
grayscale_params.cell_size = 4;

cn_params.tablename = 'CNnorm';
cn_params.useForGray = false;
cn_params.cell_size = 4;
cn_params.nDim = 10;

cnn_params.nn_name = 'imagenet-vgg-m-2048.mat'; % Name of the network
cnn_params.output_layer = [10];           % Which layers to use
cnn_params.downsample_factor = [1];           % How much to downsample each output layer
cnn_params.input_size_mode = 'adaptive';        % How to choose the sample size
cnn_params.input_size_scale = 1;                % Extra scale factor of the input samples to the network (1 is no scaling)

% Which features to include
params.t_features = {
    struct('getFeature',@get_cnn_layers, 'fparams',cnn_params),...
    struct('getFeature',@get_colorspace, 'fparams',grayscale_params),...
    struct('getFeature',@get_fhog,'fparams',hog_params),...
    struct('getFeature',@get_table_feature, 'fparams',cn_params),...
};

% Global feature parameters
params.t_global.cell_size = 4;                  % Feature cell size

% Image sample parameters
params.search_area_shape = 'square';    % The shape of the samples
params.search_area_scale = 4.5;         % The scaling of the target size to get the search area
params.min_image_sample_size = 200^2;   % Minimum area of image samples
params.max_image_sample_size = 250^2;   % Maximum area of image samples

% Spatial regularization window_parameters
params.feature_downsample_ratio = [4, 14]; %  Feature downsample ratio (We found that decreasing the downsampling ratios of CNN layer may benefit the performance)
params.reg_window_max = 1e5;           % The maximum value of the regularization window
params.reg_window_min = 1e-3;           % The minimum value of the regularization window

% Detection parameters
params.refinement_iterations = 1;       % Number of iterations used to refine the resulting position in a frame
params.newton_iterations = 5;           % The number of Newton iterations used for optimizing the detection score
params.clamp_position = false;          % Clamp the target position to be inside the image

% Learning parameters
params.output_sigma_factor = 1/16;		% Label function sigma
params.temporal_regularization_factor = [15 15]; % The temporal regularization parameters

% ADMM parameters
params.max_iterations = [2 2];
params.init_penalty_factor = [1 1];
params.max_penalty_factor = [0.1, 0.1];
params.penalty_scale_step = [10, 10];

% Scale parameters for the translation model
params.number_of_scales = 5;            % Number of scales to run the detector
params.scale_step = 1.01;               % The scale factor

% Visualization
params.visualization = 1;               % Visualiza tracking and detection scores

% GPU
params.use_gpu = true;                 % Enable GPU or not
params.gpu_id = [];                     % Set the GPU id, or leave empty to use default

% Initialize
params.seq = seq;

% Run tracker
results = tracker(params);

repo: STRCF


text:
function is=lpcrf2is(rf)
%LPCRF2IS Convert reflection coefficients to inverse sines IS=(RF)

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: lpcrf2is.m,v 1.4 2007/05/04 07:01:39 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

is=asin(rf)*2/pi;

author: decouples


text:
function [ a, seed ] = wathen_ge ( nx, ny, n, seed )

%*****************************************************************************80
%
%% WATHEN_GE returns the Wathen matrix, using general (GE) storage.
%
%  Discussion:
%
%    The Wathen matrix is a finite element matrix which is sparse.
%
%    The entries of the matrix depend in part on a physical quantity
%    related to density.  That density is here assigned random values between
%    0 and 100.
%
%    The matrix order N is determined by the input quantities NX and NY,
%    which would usually be the number of elements in the X and Y directions.
%    The value of N is
%
%      N = 3*NX*NY + 2*NX + 2*NY + 1,
%
%    The matrix is the consistent mass matrix for a regular NX by NY grid
%    of 8 node serendipity elements.
%
%    Here is an illustration for NX = 3, NY = 2:
%
%     23-24-25-26-27-28-29
%      |     |     |     |
%     19    20    21    22
%      |     |     |     |
%     12-13-14-15-16-17-18
%      |     |     |     |
%      8     9    10    11
%      |     |     |     |
%      1--2--3--4--5--6--7
%
%    For this example, the total number of nodes is, as expected,
%
%      N = 3 * 3 * 2 + 2 * 2 + 2 * 3 + 1 = 29
%
%    The matrix is symmetric positive definite for any positive values of the
%    density RHO(X,Y).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Nicholas Higham,
%    Algorithm 694: A Collection of Test Matrices in MATLAB,
%    ACM Transactions on Mathematical Software,
%    Volume 17, Number 3, September 1991, pages 289-305.
%
%    Andrew Wathen,
%    Realistic eigenvalue bounds for the Galerkin mass matrix,
%    IMA Journal of Numerical Analysis,
%    Volume 7, Number 4, October 1987, pages 449-457.
%
%  Parameters:
%
%    Input, integer NX, NY, values which determine the size of the matrix.
%
%    Input, integer N, the number of variables.
%
%    Input/output, integer SEED, the random number seed.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  em = [
     6.0, -6.0,  2.0, -8.0,  3.0, -8.0,  2.0, -6.0;
    -6.0, 32.0, -6.0, 20.0, -8.0, 16.0, -8.0, 20.0;
     2.0, -6.0,  6.0, -6.0,  2.0, -8.0,  3.0, -8.0;
    -8.0, 20.0, -6.0, 32.0, -6.0, 20.0, -8.0, 16.0;
     3.0, -8.0,  2.0, -6.0,  6.0, -6.0,  2.0, -8.0;
    -8.0, 16.0, -8.0, 20.0, -6.0, 32.0, -6.0, 20.0;
     2.0, -8.0,  3.0, -8.0,  2.0, -6.0,  6.0, -6.0;
    -6.0, 20.0, -8.0, 16.0, -8.0, 20.0, -6.0, 32.0 ]';

  [ rho, seed ] = r8mat_uniform_01 ( nx, ny, seed );
  rho = 100.0 * rho;

  node = zeros(8,1);

  for j = 1 : ny

    for i = 1 : nx
%
%  For the element (I,J), determine the indices of the 8 nodes.
%
      node(1) = 3 * j * nx + 2 * j + 2 * i + 1;
      node(2) = node(1) - 1;
      node(3) = node(1) - 2;

      node(4) = ( 3 * j - 1 ) * nx + 2 * j + i - 1;
      node(8) = node(4) + 1;

      node(5) = ( 3 * j - 3 ) * nx + 2 * j + 2 * i - 3;
      node(6) = node(5) + 1;
      node(7) = node(5) + 2;

      for krow = 1 : 8
        for kcol = 1 : 8
          a(node(krow),node(kcol)) = a(node(krow),node(kcol)) ...
            + rho(i,j) * em(krow,kcol);
        end
      end

    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function proposals= calcRigorForIm( input, rigorconfig )
%RIGOR Summary of this function goes here
%   Detailed explanation goes here

%RIGOR doesnt suppoert .mat as an input yet
if ismac
    error('rigor not supported on macOS');
end
im = input;
	[masks]=rigor_obj_segments(im,'force_recompute',true);

	if(isfield((rigorconfig.opts),'numProposals'))
		numProposals=rigorconfig.opts.numProposals;
	        if(size(masks,3)>=numProposals)
        	        masks=masks(:,:,1:numProposals);
        	else
                	fprintf('Only %d proposals were generated for image: %s\n',size(masks,3),imname);
        	end
	end
	boxes=zeros(size(masks,3),4);
	for j=1:size(masks,3)
		boxes(j,:)=mask2box(masks(:,:,j));
	end
	proposals.boxes= boxes;
	proposals.regions=masks;
end

author: Cloud-CV

repo: object-proposals


text:
function plot_tree_structure2(P,Bal,Segs,SChi,fig,ms,BO,segind)

% Plots the branch-segmented tree point cloud so that each branching order
% has its own color Blue = trunk, green = 1st-order branches,
% red = 2nd-order branches, etc.
%
% Inputs
% P         Point cloud
% Bal       Cover sets, Bal = cover.bal
% Segs      Segments, Segs = segment.segments
% SChi      Child segments, SChi = segment.ChildSegment
% fig       Figure number
% ms        Marker size
% BO        How many branching orders are plotted. 0 = all orders
% segind    Index of the segment where the plotting of tree structure
%                   starts. If segnum = 1 and BO = 0, then plots the whole
%                   tree. If segnum = 1 and B0 = 2, then plots the stem and
%                   the 1st-order branches. If segnum = 2 and BO = 0, then
%                   plots the branch whose index is 2 and all its sub-branches.

col = [
	0.00  0.00  1.00
	0.00  0.50  0.00
	1.00  0.00  0.00
	0.00  0.75  0.75
	0.75  0.00  0.75
	0.75  0.75  0.00
	0.25  0.25  0.25
	0.75  0.25  0.25
	0.95  0.95  0.00
	0.25  0.25  0.75
	0.75  0.75  0.75
	0.00  1.00  0.00
	0.76  0.57  0.17
	0.54  0.63  0.22
	0.34  0.57  0.92
	1.00  0.10  0.60
	0.88  0.75  0.73
	0.10  0.49  0.47
	0.66  0.34  0.65
	0.99  0.41  0.23];
col = repmat(col,[1000,1]);

if iscell(Segs{1})
    n = max(size(Segs));
    Seg = cell(n,1);
    for i = 1:n
        m = size(Segs{i},1);
        S = zeros(0);
        for j = 1:m
            s = Segs{i}(j);
            s = s{:};
            S = [S; s];
        end
        Seg{i} = S;
    end
else
    Seg = Segs;
end

if BO == 0
    BO = 1000;
end

S = vertcat(Bal{Seg{segind}});
figure(fig)
plot3(P(S,1),P(S,2),P(S,3),'.','Color',col(1,:),'Markersize',ms)
axis equal
forb = S;
if BO > 1
    %pause
    hold on
    c = SChi{segind};
    i = 2;
    while (i <= BO) && (~isempty(c))
        C = vertcat(Bal{unique(vertcat(Seg{c}))});
        C = setdiff(C,forb);
        figure(fig)
        plot3(P(C,1),P(C,2),P(C,3),'.','Color',col(i,:),'Markersize',ms)
        axis equal
        c = unique(vertcat(SChi{c}));
        i = i+1;
        forb = union(forb,C);
        if i <= BO
            %pause
        end
    end
    hold off
end

author: InverseTampere

repo: TreeQSM


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function MU_funcContourf(Temp,Event,handles)
handles = guidata(handles.MU_matrix_display);

figure;
contourf(handles.BMatrix);
colormap(handles.V.Color_map);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function yout = naninsert(nanvec, y)
% Insert NaNs back into data series from which they were removed
%
% :Usage:
% ::
%
%     yout = naninsert(nanvec, y)
%
% nanvec is logical indicator for NaNs, of size size(yout).  y is data.
%
% If y is a matrix, length(nanvec) should equal number of rows in y.
% NaNs are inserted for true values in nanvec in all columns of y.
%
% ..
%    Tor Wager, July 2007; Bug fix, March 2008;
%    Updated functionality to include all columns in matrix, April 2015
%  ..
%
% :See also:
% nanremove.m

    wh = find(nanvec);

    if isempty(wh), yout = y; return, end

    yout = NaN * ones(length(nanvec), size(y, 2));

    yout(~nanvec, :) = y;

%     yout(1:wh(1) - 1) = y(1:wh(1) - 1);
%     for i = 2:length(wh)
%         ystart = wh(i-1) + 2 - i; % NaN index value - num previous removed; wh(i-1) + 1 - i + 1;
%         yend = wh(i) - i; % wh(i) - 1 - i + 1
%
%         yout(wh(i-1) + 1 : wh(i) - 1, :) = y(ystart : yend, :);
%     end
%     % last segment
%         i = i+1;
%         if isempty(i), i = 2; end
%         ystart = wh(end) + 2 - i;
%         yout(wh(end)+1:end, :) = y(ystart:size(y, 1), :);
%
%     yout(wh, :) = NaN;

end  % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function value = r8_gamma_01_pdf ( alpha, rval )

%*****************************************************************************80
%
%% R8_GAMMA_01_PDF evaluates the PDF of a standard gamma distribution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2013
%
%  Author:
%
%    John Burkardt.
%
%  Parameters:
%
%    Input, real ALPHA, the shape parameter.
%    0.0 < ALPHA.
%
%    Input, real RVAL, the point where the PDF is evaluated.
%
%    Output, real VALUE, the value of the PDF at RVAL.
%
  if ( alpha <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_GAMMA_01_PDF - Fatal error!\n' );
    fprintf ( 1, '  Parameter ALPHA is not positive.\n' );
    error ( 'R8_GAMMA_01_PDF - Fatal error!' );
  end if

  if ( rval <= 0.0 )

    value = 0.0;

  else

    temp = ( alpha - 1.0 ) * log ( rval ) - rval - r8_gamma_log ( alpha );

    value = exp ( temp );

  end

  return
end


text:
function varargout = size(f, varargin)
%SIZE   Size of a DELTAFUN.
%   [S1, S2] = SIZE(F) returns the size of the funPart of F.
%
%   S = SIZE(F) returns the same as above in a 1x2 vector, S = [S1, S2].
%
% See also LENGTH.

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

% The size of a DELTAFUN is the size of its funPart.
varargout{1:nargout} = size(f.funPart, varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef AMICO_WATERFREE

properties
    id, name                % id and name of the model
    dPar                    % parallel diffusivity of the tensors [units of mm^2/s]
    dPer                    % perpendicular diffusivities of the tensors [units of mm^2/s]
    dIso                    % isotropic diffusivities [units of mm^2/s]
    OUTPUT_names            % suffix of the output maps
    OUTPUT_descriptions     % description of the output maps
end

methods

    % =================================
    % Setup the parameters of the model
    % =================================
	function obj = AMICO_WATERFREE()
        global CONFIG

        % set the parameters of the model
        obj.id        = 'WATERFREE';
        obj.name      = 'Water free';
        obj.dPar      = 1.7 * 1E-3;
        obj.dIso      = [2.0 3.0] * 1E-3;
        obj.dPer      = linspace(0.1,1.0,10) * 1E-3;

        obj.OUTPUT_names        = { 'ICVF', 'ISOVF' };
        obj.OUTPUT_descriptions = {'Intra-cellular volume fraction', 'Isotropic volume fraction'};

        % set the parameters to fit it
        CONFIG.OPTIMIZATION.SPAMS_param.mode    = 2;
        CONFIG.OPTIMIZATION.SPAMS_param.pos     = true;
        CONFIG.OPTIMIZATION.SPAMS_param.lambda  = 0;    % l1 regularization
        CONFIG.OPTIMIZATION.SPAMS_param.lambda2 = 1e-3; % l2 regularization
    end

    % ==================================================================
    % Generate high-resolution kernels and rotate them in harmonic space
    % ==================================================================
    function GenerateKernels( obj, ATOMS_path, schemeHR, AUX, idx_IN, idx_OUT )
        global CONFIG AMICO_data_path

        % Tensor compartments
        % ===================
        idx = 1;
        for i = 1:numel(obj.dPer)
            TIME = tic();
            fprintf( '\t\t- A_%03d... ', idx );

            % generate
            D = diag( [obj.dPer(i) obj.dPer(i) obj.dPar] );
            signal = obj.TensorSignal( D, schemeHR.camino );

            % rotate and save
            lm = AMICO_RotateKernel( signal, AUX, idx_IN, idx_OUT, false );
            save( fullfile( ATOMS_path, sprintf('A_%03d.mat',idx) ), '-v6', 'lm' )
            idx = idx + 1;

            fprintf( '[%.1f seconds]\n', toc(TIME) );
        end

        % Isotropic compartments
        % ======================
        for i = 1:numel(obj.dIso)
            TIME = tic();
            fprintf( '\t\t- A_%03d... ', idx );

            % generate
            D = diag( [obj.dIso(i) obj.dIso(i) obj.dIso(i)] );
            signal = obj.TensorSignal( D, schemeHR.camino );

            % resample and save
            lm = AMICO_RotateKernel( signal, AUX, idx_IN, idx_OUT, true );
            save( fullfile( ATOMS_path, sprintf('A_%03d.mat',idx) ), '-v6', 'lm' )
            idx = idx + 1;

            fprintf( '[%.1f seconds]\n', toc(TIME) );
        end

    end

    % ==============================================
    % Project kernels from harmonic to subject space
    % ==============================================
    function ResampleKernels( obj, ATOMS_path, idx_OUT, Ylm_OUT )
        global CONFIG AMICO_data_path KERNELS

        % Setup the KERNELS structure
        % ===========================
        n1 = numel(obj.dPer);
        n2 = numel(obj.dIso);
        KERNELS = {};
        KERNELS.nS       = CONFIG.scheme.nS;
        KERNELS.nA       = n1 + n2; % number of atoms
        KERNELS.A1       = zeros( [KERNELS.nS n1 181 181], 'single' );
        KERNELS.A2       = zeros( [KERNELS.nS n2], 'single' );

        % Tensors
        % =======
        idx = 1;
        for i = 1:n1
            TIME = tic();
            fprintf( '\t- A_%03d...  ', idx );

            load( fullfile( ATOMS_path, sprintf('A_%03d.mat',idx) ), 'lm' );
            KERNELS.A1(:,i,:,:) = AMICO_ResampleKernel( lm, idx_OUT, Ylm_OUT, false );
            idx = idx + 1;

            fprintf( '[%.1f seconds]\n', toc(TIME) );
        end

        % Isotropic
        % =========
        for i = 1:n2
            TIME = tic();
            fprintf( '\t- A_%03d...  ', idx );

            load( fullfile( ATOMS_path, sprintf('A_%03d.mat',idx) ), 'lm' );
            KERNELS.A2(:,i,:,:) = AMICO_ResampleKernel( lm, idx_OUT, Ylm_OUT, true );
            idx = idx + 1;

            fprintf( '[%.1f seconds]\n', toc(TIME) );
        end

    end

    % ===========================
    % Fit the model to each voxel
    % ===========================
    function [DIRs, MAPs] = Fit( obj )
        global CONFIG
        global niiSIGNAL niiMASK
        global KERNELS bMATRIX

        % setup the output files
        MAPs         = zeros( [CONFIG.dim(1:3) numel(obj.OUTPUT_names)], 'single' );
        DIRs         = zeros( [CONFIG.dim(1:3) 3], 'single' );
        niiY = niiSIGNAL;

        n1 = numel(obj.dPer);
        n2 = numel(obj.dIso);

        fprintf( '\n-> Fitting "%s" model to data:\n', obj.name );
        TIME = tic();
        for iz = 1:niiSIGNAL.hdr.dime.dim(4)
        for iy = 1:niiSIGNAL.hdr.dime.dim(3)
        for ix = 1:niiSIGNAL.hdr.dime.dim(2)
            if niiMASK.img(ix,iy,iz)==0, continue, end

            % read the signal
            b0 = mean( squeeze( niiSIGNAL.img(ix,iy,iz,CONFIG.scheme.b0_idx) ) );
            if ( b0 < 1e-3 ), continue, end
            y = double( squeeze( niiSIGNAL.img(ix,iy,iz,:) ) ./ ( b0 + eps ) );
            y( y < 0 ) = 0; % [NOTE] this should not happen!

            % build the DICTIONARY
            [ i1, i2 ] = AMICO_Dir2idx( Vt );
            A = double( [ KERNELS.A1(CONFIG.scheme.dwi_idx,:,i1,i2) KERNELS.A2(CONFIG.scheme.dwi_idx,:) ] );

            % fit AMICO
            y = y(CONFIG.scheme.dwi_idx);
            yy = [ 1 ; y ];
            AA = [ ones(1,size(A,2)) ; A ];

            % estimate CSF partial volume and remove it
            x = full( mexLasso( yy, AA, CONFIG.OPTIMIZATION.SPAMS_param ) );

            % STORE results
            DIRs(ix,iy,iz,:) = Vt; % fiber direction

            MAPs(ix,iy,iz,1) = sum( x(1:n1) ) / ( sum(x) + eps ); % intracellular volume fraction

            MAPs(ix,iy,iz,2) = 1 - MAPs(ix,iy,iz,1); % isotropic volume fraction

            x(n1+1:end) = 0;
            niiY.img(ix,iy,iz,CONFIG.scheme.dwi_idx) = AA(2:end,:)*x*b0;

        end
        end
        end
        save_untouch_nii(niiY, 'dwi_fw_corrected.nii');
        TIME = toc(TIME);
        fprintf( '   [ %.0fh %.0fm %.0fs ]\n', floor(TIME/3600), floor(mod(TIME/60,60)), mod(TIME,60) )

        % compute MAPS
        n1 = numel(obj.dPer);
        MAPs    = zeros( [1 numel(obj.OUTPUT_names)], 'single' );
        MAPs(1) = sum( x(1:n1) ) / ( sum(x) + eps ); % intracellular volume fraction
        MAPs(2) = 1 - MAPs(1);                       % isotropic volume fraction
    end

    % ================================================================
    % Simulate signal according to tensor model (1 fiber along z-axis)
    % ================================================================
    function [ signal ] = TensorSignal( obj, D, XYZB )
        nDIR   = size( XYZB, 1 );
        signal = zeros( nDIR, 1 );
        for d = 1:nDIR
            signal(d) = exp(-XYZB(d,4) * XYZB(d,1:3) * D * XYZB(d,1:3)');
        end
    end

end

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function give_Me_Jelly()

L = 3;           % Length of Computational Domain
ds = L/(2*600);  % Makes Lagrangian Spaced based on 1024x1024 finest grid

frac = 0.5;

% Makes Geometry for EQUILIBRIUM State %
%
a1 = 0.63;  %Semi-major axis length for position ONE (Expanded State)
b1 = 0.45;  %Semi-minor axis length for position ONE (Expanded State)
a1 = frac*a1; b1 = frac*b1;
a1
b1
[X1 Y1 Ninfo arcL] = make_Geometry_Position_ONE(ds,a1,b1);

% Makes Geometry for OVER-CONTRACTED State for RIGHT side %
%
a2 = 0.325; %Semi-major axis length for position TWO (Contracted State)
b2 = 0.6;   %GUESS: Semi-minor axis length for position TWO (Contracted State)
a2 = frac*a2; b2 = frac*b2;
[X2 Y2 b2] = make_Geometry_Position_TWO(a2,b2,b1,Ninfo,arcL);

% Makes Geometry for LESS-CONTRACTED State for LEFT side %
%
a3 = 0.4;   %Semi-major axis length for position TWO (Contracted State)
b3 = 0.7;   %GUESS: Semi-minor axis length for position TWO (Contracted State)
a3 = frac*a3; b3 = frac*b3;
a3
b3
[X3 Y3 b3] = make_Geometry_Position_TWO(a3,b3,b1,Ninfo,arcL);

% Makes Geometry for OVER-EXPANDED State for LEFT SIDE %
%
a4 = 0.75;  %Semi-major axis length for position TWO (Contracted State)
b4 = 0.3;   %GUESS: Semi-minor axis length for position TWO (Contracted State)
a4 = frac*a4; b4 = frac*b4;
[X4 Y4 b4] = make_Geometry_Position_TWO(a4,b4,b1,Ninfo,arcL);

% Translate Points from [-L/2,L/2]x[-L/2,L/2] -> [0,L]x[0,L]
X1 = X1 + L/2; Y1 = Y1 + L/2;
X2 = X2 + L/2; Y2 = Y2 + L/2;
X3 = X3 + L/2; Y3 = Y3 + L/2;
X4 = X4 + L/2; Y4 = Y4 + L/2;

nArm = (length(X1)-1)/2;

% Plots Geometry for Both Phases %
val = frac*0.75;
plot(X1,Y1,'*',X2(1:nArm),Y2(1:nArm),'r*',X3(nArm+2:end),Y3(nArm+2:end),'g*',X4(nArm+2:end),Y4(nArm+2:end),'m*'); hold on;
axis([0 L 0 L]);
legend('1','2','3','4');

figure(2)
plot(X2,Y2,'r*'); hold on;
plot(X1,Y1,'b*'); hold on;

% Prints Useful Info about Geometry %
fprintf('\n\n');
fprintf('Number of Total Geometry Pts: %d\n',Ninfo(2));
fprintf('Number of Geometry Pts in EACH Arm: %d\n',Ninfo(1));
fprintf('Arc length of bell (one side): %d\n',arcL);
fprintf('\n');
fprintf('(X1,Y1): Equilibrium State\n');
fprintf('(X2,Y2): OVER-Contracted State for RIGHT Side\n');
fprintf('(X3,Y3): Contracted State for LEFT Side (slightly less than right side)\n');
fprintf('(X4,Y4): OVER-EXPANDED State for LEFT Side\n');
fprintf('\n');
fprintf('Idea: \n');
fprintf('RHS contracts further than LHS, rests & returns to equil. state after LHS begins expanding\n');
fprintf('LHS contracts less than RHS, then expands further than equil. state and then returns to equil. state\n');
fprintf('\n');
fprintf('PUT FOLLOWING IN UPDATE_TARGET_POINTS.C:\n');
fprintf('numPts (# of Total Geometry Pts): %d\n',Ninfo(2));
fprintf('NL_Start (FIRST Point in LEFT Arm): %d\n\n\n',Ninfo(1)+1);

% Prints INPUT files %
print_Vertex_Pts(X1,Y1,Ninfo);
target_force = 5e6;
print_Target_Pts(target_force,Ninfo);

% Prints all vectors to .TXT files %
%print_Text_Files(Ninfo,X1,Y1,X2,Y2,X3,Y3,X4,Y4)

% Prints 2-Position Right/Left Side vectors to .TXT files %
print_2_Position_Text_File(Ninfo,X1,Y1,X2,Y2,X3,Y3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Function: Prints 2 Position Phases:
%           XY(:,1) = original position x-Values
%           XY(:,3) = original position y-Values
%           XY(2:N+1,2) = 2nd position x-Values for RIGHT bell
%           XY(N+2:end,2) = 2nd position x-Values for LEFT bell
%           XY(2:N+1,4) = 2nd position y-Values for RIGHT bell
%           XY(N+2:end,4) = 2nd position y-Values for LEFT bell
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_2_Position_Text_File(Ninfo,X1,Y1,X2,Y2,X3,Y3)

N = Ninfo(1);    % Number of pts. in ONE arm
NTot = Ninfo(2); % Total # of Lag Pts. on Jelly Bell

% Initialization
XY = zeros(NTot,4);

% Construct Matrices For Storing Points
XY(:,1) = X1(1:NTot); XY(:,2) = X1(1:NTot);
XY(:,3) = Y1(1:NTot); XY(:,4) = Y1(1:NTot);

% 2nd Phase Right Arm
XY(2:N+1,2) = X2(2:N+1);
XY(2:N+1,4) = Y2(2:N+1);

% 2nd Phase Left Arm
XY(N+2:end,2) = X2(N+2:end);
XY(N+2:end,4) = Y2(N+2:end);

% PLOT TEST! %
%figure(3)
%plot(XY(:,1),XY(:,3),'k*'); hold on;
%plot(XY(:,2),XY(:,4),'g*'); hold on;

% PRINT x/y-Information %

fileID_XY = fopen('XY_2Pos.txt','w');

for i=1:length( XY(:,1) )
    fprintf(fileID_XY,'%1.16e %1.16e %1.16e %1.16e\n',XY(i,1),XY(i,2),XY(i,3),XY(i,4));
end

fclose(fileID_XY);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Function to PRINT all TEXT files for both PHASES
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Text_Files(Ninfo,X1,Y1,X2,Y2,X3,Y3,X4,Y4)

N = Ninfo(1); %Number of pts. in ONE arm

% x-Information %

fileID_x1 = fopen('x1.txt','w');
fprintf(fileID_x1,'%1.16e\n',X1);
fclose(fileID_x1);

fileID_x2 = fopen('x2.txt','w');
fprintf(fileID_x2,'%1.16e\n',X2);
fclose(fileID_x2);

fileID_x3 = fopen('x3.txt','w');
fprintf(fileID_x3,'%1.16e\n',X3);
fclose(fileID_x3);

fileID_x4 = fopen('x4.txt','w');
fprintf(fileID_x4,'%1.16e\n',X4);
fclose(fileID_x4);

% y-Information %

fileID_y1 = fopen('y1.txt','w');
fprintf(fileID_y1,'%1.16e\n',Y1);
fclose(fileID_y1);

fileID_y2 = fopen('y2.txt','w');
fprintf(fileID_y2,'%1.16e\n',Y2);
fclose(fileID_y2);

fileID_y3 = fopen('y3.txt','w');
fprintf(fileID_y3,'%1.16e\n',Y3);
fclose(fileID_y3);

fileID_y4 = fopen('y4.txt','w');
fprintf(fileID_y4,'%1.16e\n',Y4);
fclose(fileID_y4);

% xR1 = X1(2:2+N-1); % x-Pts. on RHS of center point for Position ONE
% xL1 = X1(2+N:end); % x-Pts. on LHS of center point for Position ONE
%
% xR2 = X2(2:2+N-1); % x-Pts. on RHS of center point for Position TWO
% xL2 = X2(2+N:end); % x-Pts. on LHS of center point for Position TWO
%
% yR1 = Y1(2:2+N-1); % y-Pts. on RHS of center point for Position ONE
% yL1 = Y1(2+N:end); % y-Pts. on RHS of center point for Position ONE
%
% yR2 = Y2(2:2+N-1); % y-Pts. on LHS of center point for Position TWO
% yL2 = Y2(2+N:end); % y-Pts. on LHS of center point for Position TWO
%
%
% fileID_xR1 = fopen('xR_1.txt','w');
% fprintf(fileID_xR1,'%1.16e\n',xR1);
% fclose(fileID_xR1);
%
% fileID_xR2 = fopen('xR_2.txt','w');
% fprintf(fileID_xR2,'%1.16e\n',xR2);
% fclose(fileID_xR2);
%
% fileID_xL1 = fopen('xL_1.txt','w');
% fprintf(fileID_xL1,'%1.16e\n',xL1);
% fclose(fileID_xL1);
%
% fileID_xL2 = fopen('xL_2.txt','w');
% fprintf(fileID_xL2,'%1.16e\n',xL2);
% fclose(fileID_xL2);
%
%
% fileID_yR1 = fopen('yR_1.txt','w');
% fprintf(fileID_yR1,'%1.16e\n',yR1);
% fclose(fileID_yR1);
%
% fileID_yR2 = fopen('yR_2.txt','w');
% fprintf(fileID_yR2,'%1.16e\n',yR2);
% fclose(fileID_yR2);
%
% fileID_yL1 = fopen('yL_1.txt','w');
% fprintf(fileID_yL1,'%1.16e\n',yL1);
% fclose(fileID_yL1);
%
% fileID_yL2 = fopen('yL_2.txt','w');
% fprintf(fileID_yL2,'%1.16e\n',yL2);
% fclose(fileID_yL2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Makes the Oblate Geometry (via ellipses)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [X1 Y1 Ninfo arcL] = make_Geometry_Position_ONE(ds,a,b)

%a: Semi-major axis (along horizontal)
%b: Semi-minor axis (along vertical)

%Finding elliptical arc length
h = ( (a-b)/(a+b) )^2;
arcL = (1/4)* pi*(a+b)*(1 + 3*h/ (10 + sqrt(4-3*h) ) ); % (1/4)th of the entire ELLIPTICAL ARC

n = 1;        %counting for array
X(n) = 0;     %first x-value at (0,-b)
Y(n) = -b;    %first y-value at (0,-b)
tprev = -pi/2;  %initial angle

tol = ds / 10; %error tolerance for root-finding algorithm

maxy = max(a,b);

    while tprev < 0

        n = n+1;
        tnext = tprev + ds/maxy; %Guess for next angle value

        if tprev < -pi/4
            tfar =  0.0;%abs(2*tprev);      %Far guess
        else
            tfar = pi/4;
        end

        %initiating guess for bisection-algorithm
        xn = a*cos(tnext);
        yn = b*sin(tnext);
        errSign = ( ds - sqrt( (xn-X(n-1))^2 + (yn-Y(n-1))^2 ) );
        err = abs(errSign);

        %Bisection algorithm to make points equally spaced
        while ( err > tol )

            if errSign < 0
                tfar = tnext;
                tnext = (tnext+tprev)/2;
            elseif errSign > 0
                tprev = tnext;
                tnext = (tnext+tfar)/2;
            end

            if tnext<tprev
                fprintf('NOT CONVERGING AT ANGLE %d\n',tprev)
                break;
            end

            xn = a*cos(tnext);
            yn = b*sin(tnext);
            errSign = ( ds - sqrt( (xn-X(n-1))^2 + (yn-Y(n-1))^2 ) );
            err = abs(errSign);
        end

        X(n) = xn;   %Store X-value
        Y(n) = yn;   %Store Y-value
        tprev = tnext; %Update previous angle

        %plot(xn,yn,'*'); hold on;

    end

 N = length(X);
 X(N+1) = a; %Add one more pt.
 Y(N+1) = 0; %Add one more pt.

 X1 = [X -X(2:end)]; %Reflects geometry
 Y1 = -[Y  Y(2:end)]; %Y-values for reflected X-value points

 Ninfo(1) = N;          % # of pts. along one arm (NOT counting center pt.)
 Ninfo(2) = length(X1); % Total # of pts. in geometry

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Makes the Oblate Geometry (via ellipses)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [X2 Y2 b] = make_Geometry_Position_TWO(a,bguess,b1,Ninfo,L)

%Find correct b to give same length arm
b = Newtons_Give_Me_Semi_Minor(a,bguess,L);

%Finding elliptical arc length
h = ( (a-b)/(a+b) )^2;
arcL = (1/4)* pi*(a+b)*(1 + 3*h/ (10 + sqrt(4-3*h) ) ); % (1/4)th of the entire ELLIPTICAL ARC

ds2 = arcL / Ninfo(1); %Ensures same number of pts on arm for position ONE and TWO

[X Y Ninfo] = make_Elliptic_Geometry_Position_TWO(ds2,a,b,Ninfo);

dV = b1-b;   %Makes sure first pt is the same as in position ONE

Y2 = -(Y-dV);%Minus sign to put jelly into right orientation
X2 = X;

%Ninfo
%plot(X2,Y2,'r*'); hold on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Makes the Oblate Geometry (via ellipses)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function b = Newtons_Give_Me_Semi_Minor(a,bn,L)

tol = 1e-6;

bprev = 0.2*bn; %"Previous Guess" (left side of guess)
bnext = bn;     % First "Guess" at b
bfar =  1.5*bn; %"Far Guess" (right side of guess)

%Finding elliptical arc length
h = ( (a-bn)/(a+bn) )^2;
arcL = (1/4)* pi*(a+bn)*(1 + 3*h/ (10 + sqrt(4-3*h) ) ); % (1/4)th of the entire ELLIPTICAL ARC

errSign = L - arcL;
err = abs( errSign );

%Bisection algorithm to make points equally spaced
while ( err > tol )

    if errSign < 0
        bfar = bnext;
        bnext = (bnext+bprev)/2;
    elseif errSign > 0
        bprev = bnext;
        bnext = (bnext+bfar)/2;
    end

    %if bnext<bprev
    %    fprintf('NOT CONVERGING AT ANGLE %d\n',tprev)
    %    break;
    %end

    %Finding elliptical arc length
    h = ( (a-bnext)/(a+bnext) )^2;
    arcL = (1/4)* pi*(a+bnext)*(1 + 3*h/ (10 + sqrt(4-3*h) ) ); % (1/4)th of the entire ELLIPTICAL ARC

    errSign = ( L - arcL );
    err = abs(errSign);
end

b = bnext; %Assign value

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Makes the Oblate Geometry (via ellipses)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [X1 Y1 Ninfo] = make_Elliptic_Geometry_Position_TWO(ds,a,b,Ninfo)

%a: Semi-major axis (along horizontal)
%b: Semi-minor axis (along vertical)

n = 1;        %counting for array
X(n) = 0;     %first x-value at (0,-b)
Y(n) = -b;    %first y-value at (0,-b)
tprev = -pi/2;  %initial angle

tol = ds / 10; %error tolerance for root-finding algorithm

maxy = max(a,b);

    while n < Ninfo(1)+1

        n = n+1;
        tnext = tprev + ds/maxy; %Guess for next angle value

        if tprev < -pi/4
            tfar =  0.0;%abs(2*tprev);      %Far guess
        else
            tfar = pi/4;
        end

        %initiating guess for bisection-algorithm
        xn = a*cos(tnext);
        yn = b*sin(tnext);
        errSign = ( ds - sqrt( (xn-X(n-1))^2 + (yn-Y(n-1))^2 ) );
        err = abs(errSign);

        %Bisection algorithm to make points equally spaced
        while ( err > tol )

            if errSign < 0
                tfar = tnext;
                tnext = (tnext+tprev)/2;
            elseif errSign > 0
                tprev = tnext;
                tnext = (tnext+tfar)/2;
            end

            if tnext<tprev
                fprintf('NOT CONVERGING AT ANGLE %d\n',tprev)
                break;
            end

            xn = a*cos(tnext);
            yn = b*sin(tnext);
            errSign = ( ds - sqrt( (xn-X(n-1))^2 + (yn-Y(n-1))^2 ) );
            err = abs(errSign);
        end

        X(n) = xn;   %Store X-value
        Y(n) = yn;   %Store Y-value
        tprev = tnext; %Update previous angle

        %plot(xn,yn,'*'); hold on;

    end

 N = length(X)-1;
 %X(N+1) = a; %Add one more pt.
 %Y(N+1) = 0; %Add one more pt.

 %plot(X,Y,'*'); hold on;

 X1 = [X -X(2:end)]; %Reflects geometry
 Y1 = [Y  Y(2:end)]; %Y-values for reflected X-value points

 Ninfo(1) = N;          % # of pts. along one arm (NOT counting center pt.)
 Ninfo(2) = length(X1); % Total # of pts. in geometry

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to print VERTEX points
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Vertex_Pts(X,Y,Ninfo)

%Ninfo(1): # of pts. on one arm on Jellyfish
%Ninfo(2): # of total pts. in Jellyfish Geometry

Ntot = Ninfo(2);

vertex_fid = fopen('jellyfish.vertex', 'w');

% first line is the number of vertices in the file
fprintf(vertex_fid, '%d\n', Ntot);

for i=1:Ntot
    fprintf(vertex_fid, '%1.16e %1.16e\n', X(i), Y(i) );
end

fclose(vertex_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to print TARGET pts.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Target_Pts(target_force,Ninfo)

%Ninfo(1): # of pts. on one arm on Jellyfish
%Ninfo(2): # of total pts. in Jellyfish Geometry

Ntot = Ninfo(2);

% Write out the target point information
target_fid = fopen('jellyfish.target', 'w');

% First line is the number of target pts. in the file
fprintf(target_fid, '%d\n', Ntot);

for s=1:Ntot
        fprintf(target_fid, '%d %1.16e\n', s , target_force);
end

fclose(target_fid);

author: nickabattista


text:
clear;
clc;

addpath('../../utils/');

lambda1 = 1;
lambda2 = 1;

opts.tol = 10^-4;
opts.max_iter = 1000;
opts.tFlag = 1;
opts.init = 0;

sample_size_arr = [500 500 500 500 500];
task_num = length(sample_size_arr);
dimension = 200;

X = cell(task_num, 1);
Y = cell(task_num, 1);
W = rand(dimension, task_num);
for t = 1 :task_num
    % sampling
    X{t} = rand(sample_size_arr(t), dimension);
    %Y{t} = X{t} * W(:, t) + rand(sample_size_arr(t), 1) * 0.01;
    Y{t} = X{t} * W(:, t);
end

tic;
[W_1inf,fun_val,iter_1inf] = Least_Dirty(X,Y,lambda1,lambda2, opts);
toc;

plot(fun_val);

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function taille = imgsize(adresse)

M=imread(adresse);
taille = size(M);

if length(taille) ==2
    taille(3) = 1;
end

author: Sable

repo: mcbench-benchmarks


text:
function [mu, con, B] = clvm_Mstep(w, Y, YY, YTY, X, XX, XY, varargin)
% MSTEP_CLG Compute ML/MAP estimates for a conditional linear Von Mises
% [mu, Con, B] = clvm_Mstep(w, Y, YY, YTY, X, XX, XY, varargin)

% This currently only accounts for conditions where there are no cts parents.
% Dsc-->Cts is the structure that this can handle.

% INPUTS:
% w(i) = sum_t w(i,t) = responsibilities for each mixture component
%  If there is only one mixture component (i.e., Q does not exist),
%  then w(i) = N = nsamples,  and
%  all references to i can be replaced by 1.
% Y(:,i) = sum_t w(i,t) cos(y(:,t)) = weighted observations
% YY(:,:,i) = sum_t w(i,t) sin(y(:,t)) = weighted outer product

[cov_type, tied_cov, ...
 clamped_cov, clamped_mean, clamped_weights,  cov_prior, ...
 xs, ys, post] = ...
    process_options(varargin, ...
		    'cov_type', 'full', 'tied_cov', 0,  'clamped_cov', [], 'clamped_mean', [], ...
		    'clamped_weights', [], 'cov_prior', [], ...
		    'xs', [], 'ys', [], 'post', []);

[Ysz Q] = size(Y);

if isempty(X) % no regression
    B2 = zeros(Ysz, 1, Q);
  for i=1:Q
    B(:,:,i) = B2(:,1:0,i); % make an empty array of size Ysz x 0 x Q
  end
  [mu, con] = mixvonMises_Mstep(w, Y, YY, YTY, varargin{:});
  return;
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [xd yd] = mapping_homography(xs, ys, mode, Theta)

% [xd yd] = mapping_homography(xs, ys, mode, Theta)
%
% DESC:
% Apply the homography mapping to the source points (xs, ys)
%
% AUTHOR
% Marco Zuliani - marco.zuliani@gmail.com
%
% VERSION:
% 1.0.0
%
% INPUT:
% xs, ys            = source points
% mode              = true for forward mapping, false for backward mapping
% Theta             = parameters of the mapping function
%
% OUTPUT:
% xd, yd            = destination points
%
% SEE ALSO: imcanvas

% HISTORY
% 1.0.0             - 11/28/08 - Initial version

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% check for the input parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% preliminary operations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H = reshape(Theta, 3, 3);

xs = transpose(xs(:));
ys = transpose(ys(:));
Xs = [xs; ys];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% mapping
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if mode
    % forward mapping
    Xd = homo2cart(H*cart2homo(Xs));
else
    % backward mapping
    Xd = homo2cart(H\cart2homo(Xs));
end;

% make sure the points are returned as row vectors
xd = Xd(1,:);
yd = Xd(2,:);

return;

author: RANSAC

repo: RANSAC-Toolbox


text:
function [ n_data, n, c ] = omega_values ( n_data )

%*****************************************************************************80
%
%% OMEGA_VALUES returns some values of the OMEGA function.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by
%
%      Length [ FactorInteger [ n ] ]
%
%  First values:
%
%     N   OMEGA(N)
%
%     1    0
%     2    1
%     3    1
%     4    1
%     5    1
%     6    2
%     7    1
%     8    1
%     9    1
%    10    2
%    11    1
%    12    2
%    13    1
%    14    2
%    15    2
%    16    1
%    17    1
%    18    2
%    19    1
%    20    2
%
%  Formula:
%
%    If N = 1, then
%
%      OMEGA(N) = 0
%
%    else if the prime factorization of N is
%
%      N = P1^E1 * P2^E2 * ... * PM^EM,
%
%    then
%
%      OMEGA(N) = M
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 April 2013
%
%  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, integer N, the argument of the OMEGA function.
%
%    Output, integer C, the value of the OMEGA function.
%
  n_max = 23;

  c_vec = [ ...
      0,   1,   1,   1,   1, ...
      2,   1,   1,   1,   2, ...
      3,   1,   4,   4,   3, ...
      1,   5,   2,   2,   1, ...
      6,   7,   8 ];

  n_vec = [ ...
           1, ...
           2, ...
           3, ...
           4, ...
           5, ...
           6, ...
           7, ...
           8, ...
           9, ...
          10, ...
          30, ...
         101, ...
         210, ...
        1320, ...
        1764, ...
        2003, ...
        2310, ...
        2827, ...
        8717, ...
       12553, ...
       30030, ...
      510510, ...
     9699690 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    c = 0;
  else
    n = n_vec(n_data);
    c = c_vec(n_data);
  end

  return
end


text:
function spp_demo()
clear mex;

spp_model_file = '.\data\spp_model\VOC2007\spp_model.mat';
if ~exist(spp_model_file, 'file')
  error('%s not exist ! \n', spp_model_file);
end
try
    load(spp_model_file);
catch err
    fprintf('load spp_model_file : %s\n', err.message);
end
caffe_net_file     = fullfile(pwd, 'data\cnn_model\Zeiler_conv5\Zeiler_conv5');
caffe_net_def_file = fullfile(pwd, 'data\cnn_model\Zeiler_conv5\Zeiler_spm_scale224_test_conv5.prototxt');

use_gpu = true;
gpu_id = 1;
if use_gpu
    gpuDevice(gpu_id);
end

caffe('init', caffe_net_def_file, caffe_net_file);
caffe('set_phase_test');
if use_gpu
    spp_model.cnn.layers = spp_layers_in_gpu(spp_model.cnn.layers);
    caffe('set_mode_gpu');
else
    caffe('set_mode_cpu');
end

spm_im_size = [480 576 688 874 1200];
% spm_im_size = [ 688 ];

im = imread('.\datasets\VOCdevkit2007\VOC2007\JPEGImages\000015.jpg');

dets = spp_detect(im, spp_model, spm_im_size, use_gpu);

classes = spp_model.classes;
boxes = cell(length(classes), 1);
thres = -0.5;
for i = 1:length(boxes)
    I = dets{i}(:, 5) >= thres;
    boxes{i} = dets{i}(I, :);
end
showboxes_new(im, boxes, classes);

caffe('release');

if use_gpu
    gpuDevice([]);
end

author: ShaoqingRen


text:
function [cube, metadata] = matRad_readCube(filename)
% matRad Cube read wrapper
% determines the extension and assigns the appropriate reader to it
%
% call
%   matRad_readCube(filename)
%
% input
%   filename:   full path of the file
%
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist(filename,'file')
    error(['File ' filename ' does not exist!']);
end

[pathstr,name,ext] = fileparts(filename);

switch ext
    case {'.nrrd','.NRRD'}
        disp(['Reading NRRD: ' filename '...']);
        [cube, metadata] = matRad_readNRRD(filename);
        disp('Done!');
    otherwise
        error(['Extension ' ext ' not (yet) supported!']);
end
metadata.name = name;
metadata.path = pathstr;

end

repo: matRad


text:
function index_test02 ( )

%*****************************************************************************80
%
%% INDEX_TEST02 tests INDEX01, INDEX10, INDEX12 and INDEX21.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 November 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'INDEX_TEST02\n' );
  fprintf ( 1, '  For a 2D array,\n' );
  fprintf ( 1, '  INDEX01 column indexes with zero base,\n' );
  fprintf ( 1, '  INDEX10 row indexes with zero base,\n' );
  fprintf ( 1, '  INDEX12 column indexes with unit base,\n' );
  fprintf ( 1, '  INDEX21 row indexes with unit base.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '                Min   Index     Max\n' );
  fprintf ( 1, '\n' );

  i_min = 1;
  i = 3;
  i_max = 5;
  j_min = 1;
  j = 2;
  j_max = 4;
  fprintf ( 1, '  2D Index:  %3d%3d  %3d%3d  %3d%3d\n', ...
    i_min, j_min, i, j, i_max, j_max );

  value = index01 ( i_min, i, i_max, j_min, j, j_max );
  index_min = 0;
  index_max = index_min + ( i_max - i_min + 1 ) * ( j_max - j_min + 1 ) - 1;
  fprintf ( 1, '  INDEX01:   %6d  %6d  %6d\n', ...
    index_min, value, index_max );

  value = index10 ( i_min, i, i_max, j_min, j, j_max );
  index_min = 0;
  index_max = index_min + ( i_max - i_min + 1 ) * ( j_max - j_min + 1 ) - 1;
  fprintf ( 1, '  INDEX10:   %6d  %6d  %6d\n', ...
    index_min, value, index_max );

  value = index12 ( i_min, i, i_max, j_min, j, j_max );
  index_min = 1;
  index_max = index_min + ( i_max - i_min + 1 ) * ( j_max - j_min + 1 ) - 1;
  fprintf ( 1, '  INDEX12:   %6d  %6d  %6d\n', ...
    index_min, value, index_max );

  value = index21 ( i_min, i, i_max, j_min, j, j_max );
  index_min = 1;
  index_max = index_min + ( i_max - i_min + 1 ) * ( j_max - j_min + 1 ) - 1;
  fprintf ( 1, '  INDEX21:   %6d  %6d  %6d\n', ...
    index_min, value, index_max );

  return
end


text:
function dd_process_scp(dbname,scpfile,varargin)

%DD_PROCESS_SCP Process correlations by station and channel
%
% DD_PROCESS_SCP(DBNAME,SCP_FILE) cross correlates arrivals for the same
% phase and station. SCP_FILE follows the format written by the program
% DD_MAKE_SCP. For each line in an SCP_FILE, DD_PROCESS_SCP reads all
% available waveforms, performs cross correlations and stores the results
% is a temporary directory. If no directory name is specified the files are
% written into TMP_MATDD_CORR.
%
% DD_PROCESS_SCP(DBNAME,SCP_FILE,DIRNAME) same as above, except
% that the correlation files are written into a directory named DIRNAME.

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

% CHECK FOR VALID DATABASE
try
    db = dbopen(dbname,'r');
    dbclose(db);
catch
    error(['Could not open database: ' dbname]);
end

% CHECK INPUTS
if nargin==3
    directoryname = varargin{1};
elseif nargin==2
    directoryname = 'TMP_MATDD_CORR';
else
    error('Incorrect number of inputs');
end

% READ AND PARSE SCP FILE
try
    fid = fopen(scpfile);
catch
   error(['Unable to open scp file: ' scpfile]);
end
% check for header file
C = textscan(fid,'%s %s %s %s %s %s %s %s %s %s',1);
test1 = C{4};
test2 = C{8};
test1 = str2num(test1{1});
test2 = str2num(test2{1});
if ~isempty(test1) || ~isempty(test2)
    error('SCP file is formatted incorrectly or missing a header line');
end
C = textscan(fid,'%s %s %s %n %n %n %n %n %n %n');
fclose(fid);
sta = C{1};
chan = C{2};
iphase = C{3};
npick = C{4};
pretrig = C{5};
posttrig = C{6};
minCorr = C{7};
maxLag = C{8};
hpf = C{9};
lpf = C{10};
clear C test1 test2

% CREATE OUTPUT DIRECTORY IF NEEDED
tmp = exist(directoryname,'dir');
if (tmp~=7)
    disp(['Creating directory: ' directoryname ' ...']);
    mkdir(directoryname);
end

% CYCLE THROUGH ALL SCP LINES
for n =1:length(sta)
    if npick(n)>3
        [W,tracesFound] = load_by_orid(dbname,sta{n},chan{n},iphase{n},pretrig(n),posttrig(n));
        if tracesFound && ~isempty(W)
            [linenum,linestr,C,minimumCorrelation,maximumLag] = correlate(W,pretrig(n),posttrig(n),hpf(n),lpf(n),sta{n},iphase{n},minCorr(n),maxLag(n));
            save([directoryname '/corr_' sta{n} '_' chan{n} '_' iphase{n}],'C','linenum','linestr','minimumCorrelation','maximumLag');
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LOAD DATA FOR AN INDIVIDUAL STATION, CHANNEL, IPHASE TRIPLET

function [W,tracesFound] = load_by_orid(dbname,STA,CHAN,IPHASE,pretrig,posttrig)

W = [];

% LOAD VALUES FROM DATABASE
try
    db = dbopen(dbname,'r');
    db = dblookup_table(db,'origin');
    db1 = dblookup_table(db,'event');
    db = dbjoin(db,db1);
    db = dbsubset(db,'orid==prefor');
    db1 = dblookup_table(db,'assoc');
    db = dbjoin(db,db1);
    db1 = dblookup_table(db,'arrival');
    db = dbjoin(db,db1);
catch
    error(['Could not open all required database tables in: ' dbname]);
end

% subset required rows
db = dbsubset(db,['sta=="' STA '"']);
db = dbsubset(db,['iphase=="' IPHASE '"']);
tracesFound= dbquery(db,'dbRECORD_COUNT')
if tracesFound == 0
    warning(['No arrivals found. Skipping ' STA '_' CHAN ' ' IPHASE']);
    return
end

% nrecords = dbquery(db,'dbRECORD_COUNT')
[or_time,orid,ar_time] = dbgetv(db,'origin.time','orid','arrival.time');
ar_time = epoch2datenum(ar_time);
or_time = epoch2datenum(or_time);
dbclose(db)

% LOAD WAVEFORMS
disp(['Loading ' STA '_' CHAN ' ' IPHASE ' phase (' num2str(numel(ar_time)) ' arrivals) ...']);
ds = datasource('antelope',dbname);
for n = 1:numel(ar_time)
    scnl = scnlobject(STA,CHAN);
    startTime = ar_time(n)+(pretrig-3)/86400;
    endTime = ar_time(n)+(posttrig+3)/86400;
   try
      w = waveform(ds,scnl,startTime,endTime);
      w = addfield(w,'ORID',orid(n));
      w = addfield(w,'ORIGIN_TIME',or_time(n));
      w = addfield(w,'ARRIVAL_TIME',ar_time(n));
      W = cat(1,W,w);
   catch
      disp(['    Not able to load arrival at ' datestr(ar_time(n),31) '   (no. ' num2str(n) ')' ]);
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CORRELATE WAVEFORMS

function [linenum,linestr,C,minimumCorrelation,maximumLag] = correlate(W,pretrig,posttrig,hpf,lpf,sta,iphase,minCorr,maxLag)

orid = get(W,'ORID');
ar_time = get(W,'ARRIVAL_TIME');
C = correlation(W,ar_time);
C = butter(C,[hpf lpf]);
C = xcorr(C,[pretrig posttrig],'interp');
corr = get(C,'CORR');
lag = get(C,'LAG');

% PRE-ALLOCATE ARRAYS
numberOfRows = sum(1:numel(orid)-1);
linestr = repmat({'placeHolder'},[numberOfRows 1]);
linenum = zeros(numberOfRows,1);
minimumCorrelation = -999 * ones(numberOfRows,1);
maximumLag = -999 * ones(numberOfRows,1);

% WRITE STRINGS OF CORRELATION VALUES
rowNumber = 0;
fprintf('   writing cross correlation pairs:       ');
for n = 1:length(orid)
    for m = n+1:length(orid)
        rowNumber = rowNumber+1;
        params = sprintf('%-8s %6.3f %6.3f %-3s' , sta , lag(n,m) , (corr(n,m))^2 , iphase );
        linenum(rowNumber,1:2) = [orid(n) orid(m)];
        linestr(rowNumber) = {params};
        minimumCorrelation(rowNumber) = minCorr;
        maximumLag(rowNumber) = maxLag;
    end
    fprintf('\b\b\b\b\b\b%5.0f%%',rowNumber/numberOfRows*100);
end
fprintf('\n');

author: geoscience-community-codes

repo: GISMO


text:
function Y = sladdvec(X, v, d)
%SLADDVEC adds a vector to columns or rows of a matrix
%
% $ Syntax $
%   - Y = sladdvec(X, v, d)
%   - Y = sladdvec(X, v)
%
% $ Arguments $
%   - X:        The original matrix
%   - v:        The addend vector
%   - d:        The dimension along which the vector is to add
%   - Y:        The resultant matrix
%
% $ Description $
%   - Y = sladdvec(X, v, d) selects the most efficienct way to add a
%     vector v to every column/row of X. If d == 1, then v should be
%     a column vector, and is added to each column of X, if d == 2,
%     then v should be a row vector, and is added to each row of X.
%
%   - Y = sladdvec(X, v) will automatically determine d according to
%     the shape of v.
%
% $ Remarks $
%   - The implementation simply wraps the mex function vecop_core.
%
% $ History $
%   - Created by Dahua Lin, on Sep 10, 2006
%

if nargin < 3
    if size(v, 2) == 1
        d = 1;
    else
        d = 2;
    end
end

Y = vecop_core(X, v, d, 1);  % 1 is the opcode of addition in vecop_core

author: lmthang


text:
% RICOH MEG Reader toolbox for MATLAB
% Version 1.0.3 13-Sep-2018
%
% Files
%   getRData          - Get measured data
%   getRHdrSystem     - Get header of the system information
%   getRHdrChannel    - Get header of the system information
%   getRHdrAcqCond    - Get header of the system information
%   getRHdrEvent      - Get header of the trigger event
%   getRHdrCoregist   - Get header of the system information
%   getRHdrDigitize   - Get header of the digitization information
%   getRHdrSubject    - Get header of the subject information
%   getRHdrAnnotation - Get header of annotation
%   getRHdrSource     - Get header of the source information
%   getRMriHdr        - Get header of mri
%   getRVersion       - Get toolbox version
%
%------------------------------------------------------------
%                          History
%------------------------------------------------------------
% R1.0.3 : 2018.09.13 - support Axial2ndOrderGradiometer
% R1.0.2 : 2018.06.04 - modify comment(fear of misinterpretation)
% R1.0.1 : 2018.03.14 - modify conflict 3rd party file
% R1.0.0 : 2017.09.27 - First release
%
%-----------------------------------------------------------------------
% Copyright (C) 2017-2018     Ricoh Company Ltd. All Rights Reserved.
% All rights reserved.

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [label, Theta, w, llh] = mixGaussGb(X, opt)
% Collapsed Gibbs sampling for Dirichlet process (infinite) Gaussian mixture model (a.k.a. DPGM).
% This is a wrapper function which calls underlying Dirichlet process mixture model.
% Input:
%   X: d x n data matrix
%   opt(optional): prior parameters
% Output:
%   label: 1 x n cluster label
%   Theta: 1 x k structure of trained Gaussian components
%   w: 1 x k component weight vector
%   llh: loglikelihood
% Written by Mo Chen (sth4nth@gmail.com).
[d,n] = size(X);
mu = mean(X,2);
Xo = bsxfun(@minus,X,mu);
s = sum(Xo(:).^2)/(d*n);
if nargin == 1
    kappa0 = 1;
    m0 = mean(X,2);
    nu0 = d;
    S0 = s*eye(d);
    alpha0 = 1;
else
    kappa0 = opt.kappa;
    m0 = opt.m;
    nu0 = opt.nu;
    S0 = opt.S;
    alpha0 = opt.alpha;
end
prior = GaussWishart(kappa0,m0,nu0,S0);
[label, Theta, w, llh] = mixDpGb(X,alpha0,prior);

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function g = p40_g ( n, x )

%*****************************************************************************80
%
%% P40_G evaluates the gradient for problem 40.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real G(N), the gradient of the objective function.
%
  g = zeros ( n, 1 );

  g(1) = 2.0 * x(1) + 0.9 * pi * sin ( 3.0 * pi * x(1) );
  g(2) = 4.0 * x(2) - 4.0 * pi * sin ( 4.0 * pi * x(2) );

  return
end

author: johannesgerer


text:
function x = modhs52(p)
  n=4;

  x(1)=4.0*p(1)-p(2);
  x(2)=p(2)+p(3)-2.0;
  x(3)=p(4)-1.0;
  x(4)=p(5)-1.0;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
classdef matRad_EUD < DoseObjectives.matRad_DoseObjective
% matRad_EUD Implements a penalized equivalent uniform dose objective
%   See matRad_DoseObjective for interface description
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2020 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    properties (Constant)
        name = 'EUD';
        parameterNames = {'EUD^{ref}', 'k'};
        parameterTypes = {'dose','numeric'};
    end

    properties
        parameters = {0, 3.5};
        penalty = 1;
    end

    methods
        function obj = matRad_EUD(penalty,eudRef, eudExponent)
            %If we have a struct in first argument
            if nargin == 1 && isstruct(penalty)
                inputStruct = penalty;
                initFromStruct = true;
            else
                initFromStruct = false;
                inputStruct = [];
            end

            %Call Superclass Constructor (for struct initialization)
            obj@DoseObjectives.matRad_DoseObjective(inputStruct);

            %now handle initialization from other parameters
            if ~initFromStruct
                if nargin >= 3 && isscalar(eudExponent)
                    obj.parameters{2} = eudExponent;
                end

                if nargin >= 2 && isscalar(eudRef)
                    obj.parameters{1} = eudRef;
                end

                if nargin >= 1 && isscalar(penalty)
                    obj.penalty = penalty;
                end
            end
        end

        %% Calculates the Objective Function value
        function fDose = computeDoseObjectiveFunction(obj,dose)
            % get exponent for EUD
            k = obj.parameters{2};

            % calculate power sum
            powersum = sum(dose.^k);

            %Calculate objective

            %This check is not needed since dose is always positive
            %if powersum > 0
            fDose = obj.penalty * (nthroot(powersum/numel(dose),k) - obj.parameters{1})^2;
            %end
        end

        %% Calculates the Objective Function gradient
        function fDoseGrad  = computeDoseObjectiveGradient(obj,dose)
            % get exponent for EUD
            k = obj.parameters{2};

            %numerical stability
            dose(dose == 0) = 0.001;

            % calculate power sum
            powersum = sum(dose.^k);

            %This check is not needed since dose is always positive
            %if powersum > 0

            %derivatives = nthroot(1/numel(dose),k) * powersum^((1-k)/k) * (dose.^(k-1));
            fDoseGrad = 2 * obj.penalty * nthroot(1/numel(dose),k) * powersum^((1-k)/k) * (dose.^(k-1)) .* (nthroot(powersum/numel(dose),k) - obj.parameters{1});
            %end
            if any(~isfinite(fDoseGrad)) % check for inf and nan for numerical stability
                error(['EUD computation failed. Reduce exponent to resolve numerical problems.']);
            end
        end
    end

end


text:
classdef PathVertexesToBoundaryComputer < handle

    properties (Access = public)

    end

    properties (Access = private)
       mesh
       singularityCoord
    end

    properties (Access = private)
       boundaryPointCoord
       pathVector
       closestVertex
       boundaryNodes
       connectedVertex
       pathVertexes
       pathCells
    end

    methods (Access = public)

        function obj = PathVertexesToBoundaryComputer(cParams)
            obj.init(cParams)
        end

        function v = compute(obj)
            obj.computeBoundaryPointCoord();
            obj.createStraightPathVector();
            obj.computeClosestVertex();
            obj.computeBoundaryNodes();
            obj.computeMeshEdges();
            obj.computePath();
            v = obj.pathVertexes;
        end

        function plot(obj)
            figure()
            hold on
            obj.plotMesh();
            obj.plotStraightPath();
            obj.plotClosestVertex();
            obj.plotVerticesPath();
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.mesh               = cParams.mesh;
            obj.singularityCoord   = cParams.singularityCoord;
        end

        function computeBoundaryPointCoord(obj)
        %    obj.computeToyBoundaryPoint();
            obj.computeBenchmarkBoundaryPoint();
        end

        function computeToyBoundaryPoint(obj)
            obj.boundaryPointCoord = [0 150];
        end

        function computeBenchmarkBoundaryPoint(obj)
           obj.boundaryPointCoord(:,2) = obj.singularityCoord(:,2);
           obj.boundaryPointCoord(:,1) = max(obj.mesh.coord(:,1));
        end

        function createStraightPathVector(obj)
            OA = obj.singularityCoord;
            OB = obj.boundaryPointCoord;
            AB = OB - OA;
            obj.pathVector = AB ;
        end

        function computeClosestVertex(obj)
            xy  = obj.mesh.coord;
            xyS = obj.singularityCoord;
            xD  = xy(:,1) - xyS(:,1);
            yD  = xy(:,2) - xyS(:,2);
            dist = sqrt(xD.^2 + yD.^2);
            [~,iD] = min(dist);
            obj.closestVertex = iD;
        end

        function computeBoundaryNodes(obj)
            coord  = obj.mesh.coord;
            nodesB = boundary(coord);
            obj.boundaryNodes = nodesB;
        end

        function computeMeshEdges(obj)
            obj.mesh.computeEdges();
        end

        function itIs = isInBoundary(obj,node)
            nodesB = obj.boundaryNodes;
            itIs = any(node == nodesB);
        end

        function computePath(obj)
            i = 1;
            vertex       = obj.closestVertex;
            pVertexes    = vertex;
            while ~obj.isInBoundary(vertex)
                otherVertex = obj.mesh.computeConnectedVertex(vertex);
                iD          = obj.computeOptimalVertex(vertex,otherVertex);
                newVertex   = otherVertex(iD);
                i = i + 1;
                vertex       = newVertex;
                pVertexes(i) = vertex;
            end
            obj.pathVertexes = pVertexes;
        end

        function iD = computeOptimalVertex(obj,currentVertex,trialVertexes)
            s.currentVertexCoord  = obj.mesh.coord(currentVertex,:);
            s.trialVertexesCoord  = obj.mesh.coord(trialVertexes,:);
            s.boundaryVertexCoord = obj.boundaryPointCoord;
            s.lineVector          = obj.pathVector;
            n = NextVertexToBoundaryComputer(s);
            iD = n.compute();
        end

        function plotMesh(obj)
            obj.mesh.plot();
        end

        function plotStraightPath(obj)
            xP = [obj.singularityCoord(:,1),obj.boundaryPointCoord(:,1)];
            yP = [obj.singularityCoord(:,2),obj.boundaryPointCoord(:,2)];
            plot(xP,yP,'+-')
        end

        function plotClosestVertex(obj)
            cV = obj.closestVertex;
            x = obj.mesh.coord(cV,1);
            y = obj.mesh.coord(cV,2);
            plot(x,y,'r+')
        end

        function plotVerticesPath(obj)
            cV = obj.pathVertexes;
            x = obj.mesh.coord(cV,1);
            y = obj.mesh.coord(cV,2);
            plot(x,y,'g-','LineWidth',5)
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [g1, g2] = simwhiteXrbfinfwhiteKernGradient(simKern, rbfKern, t1, varargin)

% SIMWHITEXRBFINFWHITEKERNGRADIENT Compute a cross gradient between a
% SIM-WHITE kernel and an RBF-WHITE kernel (with integration limits between
% minus infinity and infinity).
% FORMAT
% DESC computes cross gradient of parameters of a cross kernel between a
% SIM-WHITE and an RBF-WHITE kernels for the multiple output kernel.
% ARG simKern : the kernel structure associated with the SIM-WHITE kernel.
% ARG rbfKern : the kernel structure associated with the RBF-WHITE kernel.
% ARG t1 : inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of the parameters of the SIM-WHITE kernel, for
% ordering see simwhiteKernExtractParam.
% RETURN g2 : gradient of the parameters of the RBF-WHITE kernel, for
% ordering see rbfwhiteKernExtractParam.
%
% FORMAT
% DESC computes cross kernel terms between a SIM-WHITE kernel and an
% RBF-WHITE kernel (with integration limits between minus infinity and
% infinity) for the multiple output kernel.
% ARG simKern : the kernel structure associated with the SIM-WHITE kernel.
% ARG rbfKern : the kernel structure associated with the RBF-WHITE kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of the parameters of the SIM-WHITE kernel, for
% ordering see simwhiteKernExtractParam.
% RETURN g2 : gradient of the parameters of the RBF-WHITE kernel, for
% ordering see rbfwhiteKernExtractParam.
%
% SEEALSO : multiKernParamInit, multiKernCompute, simwhiteKernParamInit,
% rbfinfwhiteKernParamInit, simwhiteKernExtractParam, rbfinfwhiteKernExtractParam
%
% COPYRIGHT : David Luengo, 2009

% KERN

if nargin < 5
    t2 = t1;
else
    t2 = varargin{1};
end
covGrad = varargin{end};

if size(t1, 2) > 1 | size(t2, 2) > 1
  error('Input can only have one column');
end
if simKern.variance ~= rbfKern.variance
  error('Kernels cannot be cross combined if they have different variances.')
end

g1 = zeros(1, simKern.nParams);
g2 = zeros(1, rbfKern.nParams);

T1 = repmat(t1, 1, size(t2, 1));
T2 = repmat(t2.', size(t1, 1), 1);
deltaT = T1 - T2;

% Parameters required for further computations
isStationary = simKern.isStationary;
variance = simKern.variance;
decay = simKern.decay;
sensitivity = simKern.sensitivity;
invWidth = rbfKern.inverseWidth;

% Setting normalised parameters for computation of K
simKern.variance = 1;
rbfKern.variance = 1;
simKern.sensitivity = 1;
K = simwhiteXrbfinfwhiteKernCompute(simKern, rbfKern, t1, t2);

% Gradient w.r.t. the decay (simKern)
g1(1) = variance * sensitivity * sum(sum( ...
    ((decay/invWidth-deltaT) .* K  + 1/sqrt(2*pi*invWidth) ...
        * (exp(-(decay*T1+0.5*invWidth*(T2.^2))) * (~isStationary) ...
            - exp(-0.5*invWidth*(deltaT.^2)))) ...
    .* covGrad));
%       * exp(0.5*(decay^2)/invWidth-decay*deltaT) ...
%         .* (exp(-0.5*invWidth*((T2+decay/invWidth).^2)) * (~isStationary) ...
%           - exp(-0.5*invWidth*((deltaT-decay/invWidth).^2)))) ...

% Gradient w.r.t. the inverse width (rbfKern)
g2(1) = 0.5 * variance * sensitivity * sum(sum( ...
    ((-(decay/invWidth)^2) * K  + 1/sqrt(2*pi*invWidth) ...
        * ((T2-decay/invWidth) .* exp(-(decay*T1+0.5*invWidth*(T2.^2))) * (~isStationary) ...
            + (deltaT+decay/invWidth) .* exp(-0.5*invWidth*(deltaT.^2)))) ...
    .* covGrad));
%       * exp(0.5*(decay^2)/invWidth - decay*deltaT) ...
%         .* ((T2-decay/invWidth) .* exp(-0.5*invWidth*((T2+decay/invWidth).^2)) * (~isStationary) ...
%           + (deltaT+decay/invWidth) .* exp(-0.5*invWidth*((deltaT-decay/invWidth).^2)))) ...

% Gradient w.r.t. sigma_r^2
g1(2) = sensitivity * sum(sum(K .* covGrad));
g2(2) = 0; % Otherwise it is counted twice

% Gradient w.r.t. sensitivity (only simKern)
g1(3) = variance * sum(sum(K .* covGrad));

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function str = parSimp(str)
%PARSIMP   Remove unnecessary parentheses from string inputs.
%   STROUT = PARSIMP(STRIN) returns the string STROUT, obtained by doing some
%   basic simplifications of the string STRIN, attempting to remove unnecessary
%   parenthesis, zeros, and consecutive +/- pairs.

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

% We can't really expect to simplify a one character string...
if ( length(str) < 2 )
    return
end

% Simplify until we can't simplify more!
oldlength = 0;
while ( length(str) ~= oldlength )
    oldlength = length(str);
    str = parSimpMain(str);
end

end

function str = parSimpMain(str)
%PARSIMMAIN   Do the actual simplification

% Find all locations of ( and ):
leftParLoc = strfind(str, '(');
rightParLoc = strfind(str, ')');

% Find weighted vector of operators:
opVec = (str == '-') + (str == '+') + 2*((str == '*') + (str == '/')) + ...
    3*(str == '^');

% Operator locations
opLoc = find(opVec);

% Location of characters and numbers
charLoc = regexp(str, '[A-Za-z0-9@]');
ltgtLoc = regexp(str, '[\<\>]');

% Error if the number of ( and ) are not equal
if ( length(leftParLoc) ~= length(rightParLoc) )
    error('CHEBFUN:STRINGPARSER:parSimp:parenthSimplify', ...
        'Incorrect number of parenthesis.');
end

% Store number of parenthesis
numOfPars = length(leftParLoc);
% Pair them together
pairsLoc = zeros(numOfPars, 2);
pairsVec = zeros(1, length(str));
parCounter = 1;
leftStack = [];
for strIndex = 1:length(str)
    if ( str(strIndex) == '(' ) % Push to the stack
        leftStack = [leftStack, strIndex]; %#ok<AGROW>
    elseif ( str(strIndex) == ')' ) % Pop from the stack
        % Write information to pairs
        pairsLoc(parCounter, 1) = leftStack(end);
        pairsLoc(parCounter, 2) = strIndex;
        pairsVec([leftStack(end) strIndex]) = parCounter;
        % Update stack and counter
        leftStack(end) = [];
        parCounter = parCounter + 1;
    end
end

for k = 1:numel(ltgtLoc)
    idx = find((pairsLoc(:, 1) < ltgtLoc(k)) & (pairsLoc(:, 2) > ltgtLoc(k)));
    [ignored, idx2] = max(pairsLoc(idx, 1));
    charLoc = sort([charLoc pairsLoc(idx(idx2), :)]);
    pairsLoc(idx(idx2),:) = [];
    numOfPars = numOfPars - 1;
end

leftParsLoc = pairsLoc(:, 1);
rightParsLoc = pairsLoc(:, 2);

    function removeExtras(varargin)
        % Remove things other than parentheses, like consecutive +/- pairs
        % and hanging zeros. If there's an input argument, then we don't
        % attempt to remove zeros.

        % Remove leading + signs
        if ( strcmp(str(1), '+') )
            str(1) = [];
            opVec(1) = [];
            leftParsLoc = leftParsLoc - 1;
            rightParsLoc = rightParsLoc - 1;
            opLoc = opLoc - 1;
            charLoc = charLoc - 1;
        end

        % Remove consecutive +/- pairs
        opVec1 = ( opVec == 1 );
        dOV = ~[abs(diff(opVec1)) 1] & opVec1; % Find adjacent opVec == 1's
        pm = find(dOV, 1); % Get the first
        if ( ~isempty(pm) ) % If +/- pair exists, remove the first occurence
            if ( str(pm) == str(pm+1) )
                str(pm) = '+';
            else
                str(pm) = '-';
            end
            str(pm + 1) = [];
            % Need to update indices and operators
            opVec(pm+1) = [];
            leftParsLoc(leftParsLoc > pm) = leftParsLoc(leftParsLoc > pm) - 1;
            rightParsLoc(rightParsLoc > pm) = ...
                rightParsLoc(rightParsLoc > pm) - 1;
            opLoc(opLoc > pm) = opLoc(opLoc > pm) - 1;
            charLoc(charLoc > pm) = charLoc(charLoc > pm) - 1;
            removeExtras(1) % recurse
        end

        % +/- pairs recurse locally
        if ( nargin > 0 )
            return
        end

    end

removeExtras()

% Loop over the pairs
for pIndex = 1:numOfPars

    removeExtras()

    pLeft = leftParsLoc(pIndex);
    pRight = rightParsLoc(pIndex);
    nextOpLeft = max(opLoc(opLoc < pLeft));
    nextOpRight = min(opLoc(opLoc > pRight));

    if ( nextOpLeft == 0 )
        nextOpLeft = [];
    end

    if ( nextOpRight == 0 )
        nextOpRight = [];
    end

    if ( any(charLoc == pLeft - 1) )
        % If the character in str next to the left of the left parenthesis is a
        % letter or a number (e.g. from sin or log2), we have a function from
        % which we cannot remove the () pair.
        continue
    end

    % Just mask the interior of these parentheses
    mask = zeros(1, length(str));
    mask(pLeft+1:pRight-1) = true;
    idx = find((leftParsLoc > pLeft) & (rightParsLoc < pRight));
    for k = 1:numel(idx)
        mask(leftParsLoc(idx(k)):rightParsLoc(idx(k))) = false;
    end
    interiorOpVec = opVec.*mask;
    interiorOpVec(interiorOpVec==0) = inf;
    minOpInside = min(interiorOpVec(pLeft:pRight));

    % If not, we perform a check to see whether we can remove them.
    if ( (isempty(nextOpLeft) || (minOpInside >= opVec(nextOpLeft))) && ...
            (isempty(nextOpRight) || minOpInside >= opVec(nextOpRight)) )

        if ( (minOpInside == 1) && ~isempty(nextOpLeft) && ...
                strcmp(str(pLeft - 1), '-') )
            % Deal with the case of minus (switch interior pluses and minuses
%            continue % uncomment if we don't want to simplify -(a+b) = -a-b
            interiorPM = interiorOpVec == 1;
            interiorPM(pLeft+1) = 0;
            m = (str == '-') & interiorPM;
            p = (str == '+') & interiorPM;
            % But we have to be careful that we don't change 1e-2 to 1e+2!
            mPos = find(m);
            pPos = find(p);
            mThrow = str(mPos-1) == 'e';
            pThrow = str(pPos-1) == 'e';
            m(mPos(mThrow)) = 0;
            p(pPos(pThrow)) = 0;
            str(m) = '+';
            str(p) = '-';
        end

        if ( (minOpInside == 2) && ~isempty(nextOpLeft) && ...
                strcmp(str(pLeft - 1), '/') )
            % Deal with division.  Unlike subtraction, we don't intend to
            % generally change a/(b/c) into a*c/b, so just move on.
            % (See GitHub issue #2357.)
            continue;
        end

        % Remove parenthesis pairs
        str(pLeft) = [];
        str(pRight-1) = [];

        % Update indices and operators
        opVec(pLeft) = [];
        opVec(pRight-1) = [];
        leftParsLoc(pIndex) = NaN;
        rightParsLoc(pIndex) = NaN;
        leftParsLoc = updateLoc(leftParsLoc, pLeft, pRight);
        rightParsLoc = updateLoc(rightParsLoc, pLeft, pRight);
        opLoc = updateLoc(opLoc, pLeft, pRight);
        charLoc = updateLoc(charLoc, pLeft, pRight);
    end

end

% Remove leading + signs
if ( strcmp(str(1), '+') )
    str(1) = [];
end

end

function out = updateLoc(in, pLeft, pRight)
out = in - (in > pLeft) - (in > pRight);
end


text:
function [xn, xmean, xstd] = normdata(x,xmean,xstd)
%NORMDATA  Normalize input to zero mean and unit variance
%
%  Description
%    [XN, XMEAN, XSTD] = NORMDATA(X) normalizes X to zero mean and
%    unit variance. Returns normalized XN, mean of X XMEAN and
%    standard deviance of X XSTD.
%
%    XN = NORMDATA(X,XMEAN,XSTD) normalizes X using precomputed
%    XMEAN and XSTD.
%
%  See also DENORMDATA
%

% Copyright (c) 2010 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  if nargin<=1
    xmean=nanmean(x);
    xstd=nanstd(x);
  end
  xn=bsxfun(@rdivide,bsxfun(@minus,x,xmean),xstd);


text:
% EEG_FINDBOUNDARIES - return indices of boundary events
%
% Usage:
%   >> boundaryIndices = eeg_findboundaries(INEEG);
%   >> boundaryIndices = eeg_findboundaries(INEEG.event);
%
% Inputs:
%   INEEG         - input EEG dataset structure. Used to check if event
%                   types are string or numerical.
%
% Outputs:
%   boundaryIndices - indices of boundary events
%
% Author: Arnaud Delorme, 2022
%
% see also: EEGLAB

% Copyright (C) 2022 Arnaud Delorme
%
% 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 boundaries = eeg_findboundaries(EEG)

    if nargin < 1
        help eeg_findboundaries;
        return
    end

    boundaries = [];
    if isempty(EEG)
        return;
    end

    if isfield(EEG, 'event') && isfield(EEG, 'setname')
        tmpevent = EEG.event;
    else
        tmpevent = EEG;
    end

    if ~isfield(tmpevent, 'type')
        return;
    end

    % type of boundary event
    eeglab_options;
    if ischar(tmpevent(1).type)
        boundaries = strmatch('boundary', { tmpevent.type });
    elseif option_boundary99
        boundaries = find([ tmpevent.type ] == -99);
    else
        boundaries = [];
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab