text:
function h=huniform(p,varargin)

%   Copyright (C) 2004-2006 Per-Olof Persson. See COPYRIGHT.TXT for details.

h=ones(size(p,1),1);

repo: ifem


text:
% HOMOGRAPHY2D - computes 2D homography
%
% Usage:   H = homography2d(x1, x2)
%          H = homography2d(x)
%
% Arguments:
%          x1  - 3xN set of homogeneous points
%          x2  - 3xN set of homogeneous points such that x1<->x2
%
%           x  - If a single argument is supplied it is assumed that it
%                is in the form x = [x1; x2]
% Returns:
%          H - the 3x3 homography such that x2 = H*x1
%
% This code follows the normalised direct linear transformation
% algorithm given by Hartley and Zisserman "Multiple View Geometry in
% Computer Vision" p92.

% Copyright (c) 2003-2005 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% pk at csse uwa edu au
% http://www.csse.uwa.edu.au/~pk
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.

% May 2003  - Original version.
% Feb 2004  - Single argument allowed for to enable use with RANSAC.
% Feb 2005  - SVD changed to 'Economy' decomposition (thanks to Paul O'Leary)

function H = homography2d(varargin)

    [x1, x2] = checkargs(varargin(:));

    % Attempt to normalise each set of points so that the origin
    % is at centroid and mean distance from origin is sqrt(2).
    [x1, T1] = normalise2dpts(x1);
    [x2, T2] = normalise2dpts(x2);

    % Note that it may have not been possible to normalise
    % the points if one was at infinity so the following does not
    % assume that scale parameter w = 1.

    Npts = length(x1);
    A = zeros(3*Npts,9);

    O = [0 0 0];
    for n = 1:Npts
	X = x1(:,n)';
	x = x2(1,n); y = x2(2,n); w = x2(3,n);
	A(3*n-2,:) = [  O  -w*X  y*X];
	A(3*n-1,:) = [ w*X   O  -x*X];
	A(3*n  ,:) = [-y*X  x*X   O ];
    end

    [U,D,V] = svd(A,0); % 'Economy' decomposition for speed

    % Extract homography
    H = reshape(V(:,9),3,3)';

    % Denormalise
    H = T2\H*T1;

%--------------------------------------------------------------------------
% Function to check argument values and set defaults

function [x1, x2] = checkargs(arg);

    if length(arg) == 2
	x1 = arg{1};
	x2 = arg{2};
	if ~all(size(x1)==size(x2))
	    error('x1 and x2 must have the same size');
	elseif size(x1,1) ~= 3
	    error('x1 and x2 must be 3xN');
	end

    elseif length(arg) == 1
	if size(arg{1},1) ~= 6
	    error('Single argument x must be 6xN');
	else
	    x1 = arg{1}(1:3,:);
	    x2 = arg{1}(4:6,:);
	end
    else
	error('Wrong number of arguments supplied');
    end

author: SuTanTank

repo: VideoStitchingViaShakinessRemoving


text:
classdef MimCombinedImageDatabaseController < CoreBaseClass
    % MimCombinedImageDatabaseController. Part of the internal framework of the TD MIM Toolkit.
    %
    %     You should not use this class within your own code. It is intended to
    %     be used internally within the framework of the TD MIM Toolkit.
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    events
        ProjectChangedEvent
        PatientChangedEvent
    end

    properties (Access = private)
        GuiCallback
        CurrentProject = MimImageDatabase.LocalDatabaseId;
        MatNatDatabase
    end

    methods
        function obj = MimCombinedImageDatabaseController(gui_callback, matnat_database)
            obj.MatNatDatabase = matnat_database;
            obj.GuiCallback = gui_callback;
        end

        function ProjectClicked(obj, project_id)
            obj.CurrentProject = project_id;
            notify(obj, 'ProjectChangedEvent', CoreEventData(project_id));
        end

        function PatientClicked(obj, patient_id)
            if strcmp(obj.CurrentProject, MimImageDatabase.LocalDatabaseId)
                obj.GuiCallback.LoadPatient(patient_id);
            else
                notify(obj, 'PatientChangedEvent', CoreEventData(patient_id));
            end
        end

        function SeriesClicked(obj, patient_id, series_uid)
            if strcmp(obj.CurrentProject, MimImageDatabase.LocalDatabaseId)
                obj.GuiCallback.LoadFromPatientBrowser(series_uid, patient_id);
            else
                if ~isempty(obj.MatNatDatabase)
                    import_path = obj.MatNatDatabase.downloadScan(obj.CurrentProject, patient_id, series_uid);
                    obj.GuiCallback.ImportFromPath(import_path);
                end
            end
        end

        function AddPatient(obj)
            obj.GuiCallback.ImportMultipleFiles;
        end

        function AddSeries(obj)
            obj.GuiCallback.ImportMultipleFiles;
        end

        function DeletePatient(obj, patient_id)
            obj.GuiCallback.DeletePatient(patient_id);
        end

        function DeleteSeries(obj, series_uid)
            obj.GuiCallback.DeleteDataset(series_uid);
        end

        function RefreshProjects(obj)
        end

        function BringToFront(obj)
            obj.GuiCallback.BringToFront;
        end

        function UnlinkDataset(obj, series_uid)
            obj.GuiCallback.UnlinkDataset(series_uid);
        end

        function project_id = GetCurrentProject(obj)
            project_id = obj.CurrentProject;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function varargout = ihashfind(varargin)
% VL_IHASHFIND
%  SEL = VL_IHASHFIND(ID, NEXT, K, X)
[varargout{1:nargout}] = vl_ihashfind(varargin{:});

author: jianxiongxiao

repo: ProfXkit


text:
function s = structrmfield(s,flds)
% s = structrmfield(s,flds)
% rmfield with removal of nested fields allowed
%
% Elements of flds which are not present in s are silently ignored.
%
% Example:
% flds = {'fld1','fld2.subfld'};
% s = structrmfield(s,flds);

assert(isstruct(s) && isscalar(s));
if ischar(flds)
  flds = cellstr(flds);
end
assert(iscellstr(flds));

for f = flds(:)',f=f{1}; %#ok<FXSET>
  fsplit = regexp(f,'\.','split');
  if isscalar(fsplit)
    fldRm = fsplit{1};
    if isfield(s,fldRm)
      s = rmfield(s,fldRm);
    else
      %warning('structrmfield:noField','No field ''%s''.',fldRm);
    end
  else
    fldTop = fsplit{1};
    fldRemainder = f(numel(fldTop)+2:end);
    valTop = s.(fldTop);
    if ~isstruct(valTop)
      %warning('structrmfield:badField','Expected field ''%s'' to be a struct.',fldTop);
    else
      s.(fldTop) = structrmfield(valTop,fldRemainder);
    end
  end
end

author: kristinbranson

repo: JAABA


text:
function [gal] = km32gal(km3)
% Convert volume from cubic kilometers to US liquid gallons.
% Chad Greene 2012
gal = km3*264172052360;

author: Sable

repo: mcbench-benchmarks


text:
%DNA medium test
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 4/19/2011
classdef DNADamageRepair_Test < TestCase
    properties
        simulation
    end

    methods
        function this = DNADamageRepair_Test(name)
            this = this@TestCase(name);
        end

        function setUp(this)
            sim = edu.stanford.covert.cell.sim.SimulationFixture.load([], {
                'DNADamage'
                'DNARepair'
                });
            sim.applyOptions('verbosity', 0, 'seed', 1);

            this.simulation = sim;
        end

        function testExpectedDamageRates_sampleAccessibleSites(this)
            %references
            sim = this.simulation;
            c = sim.state('Chromosome');

            %test sampleAccessibleSites
            iterMax = 1000;
            cnts = zeros(iterMax, 1);
            prob = 1e-6;
            for i = 1:iterMax
                cnts(i) = size(c.sampleAccessibleSites(prob, Inf, 'T'), 1);
            end
            assertEqual([0; 1], unique(cnts));
            assertElementsAlmostEqual(...
                sum(cnts), ...
                iterMax * prob * 2 * c.sequenceLen * (1-c.sequenceGCContent)/2, ...
                'relative', 0.10);
        end

        function testExpectedDamageRates_setSiteDamaged(this)
            %references
            sim = this.simulation;
            c = sim.state('Chromosome');

            %test setSiteDamaged
            iterMax = 1000;
            cnts = zeros(iterMax, 1);
            prob = 1e-6;
            for i = 1:iterMax
                c.initialize();
                cnts(i) = size(c.setSiteDamaged('damagedBases', 1, prob, Inf, 'T', 'sequence'), 1);
            end
            assertEqual([0; 1], unique(cnts));
            assertElementsAlmostEqual(...
                sum(cnts), ...
                iterMax * prob * 2 * c.sequenceLen * (1-c.sequenceGCContent)/2, ...
                'relative', 0.10);
        end

        function testExpectedDamageRates_evolveState_simpleReactionSet(this)
            %references
            sim = this.simulation;
            c = sim.state('Chromosome');
            d = sim.process('DNADamage');

            %test evolve state with simplified reaction set
            iterMax = 1000;

            d.reactionWholeCellModelIDs = d.reactionWholeCellModelIDs(1:4);
            d.reactionStoichiometryMatrix = d.reactionStoichiometryMatrix(:, 1:4);
            d.reactionDamageTypes = {'damagedBases'; 'damagedBases'; 'intrastrandCrossLinks'; 'intrastrandCrossLinks'};
            d.reactionBounds = [zeros(4, 1) repmat(1e-6, 4, 1)];
            d.reactionRadiation = zeros(4, 1);
            d.reactionDNAProduct = (1:4)';
            d.reactionVulnerableMotifs = {'A'; 'C'; 'G'; 'T'};
            d.reactionVulnerableMotifTypes = repmat({'sequence'}, 4, 1);

            c.initialize();
            expRates = d.calcExpectedReactionRates() * iterMax;

            d.substrates(any(d.reactionSmallMoleculeStoichiometryMatrix,2)) = 1e6;
            for i = 1:iterMax
                d.evolveState();
            end
            [~, vals] = find(c.damagedBases);
            assertElementsAlmostEqual(histc(vals, d.reactionDNAProduct), expRates .* [1; 1; 0; 0], ...
                'relative', 0.10);
            [~, vals] = find(c.intrastrandCrossLinks);
            assertElementsAlmostEqual(histc(vals, d.reactionDNAProduct), expRates .* [0; 0; 1; 1], ...
                'relative', 0.10);
        end

        function testExpectedDamageRates_evolveState_lessSimpleReactionSet(this)
            %references
            sim = this.simulation;
            c = sim.state('Chromosome');
            m = sim.state('Metabolite');
            d = sim.process('DNADamage');

            %simulate damage
            iterMax = 1000;

            d.reactionBounds(d.reactionBounds(:, 2) < 1000, 2) = 1e-6;
            d.substrates(d.reactionRadiation(d.reactionRadiation ~= 0)) = 1;
            d.reactionBounds(~cellfun(@ischar, d.reactionVulnerableMotifs), 2) = 0;

            expRates = d.calcExpectedReactionRates() * iterMax;

            assertFalse(any(d.reactionDNAProduct == m.m6ADIndexs & strcmp(d.reactionDamageTypes, 'damagedBases')));

            d.substrates(any(d.reactionSmallMoleculeStoichiometryMatrix,2)) = 1e6;
            for i = 1:iterMax
                d.evolveState();
            end

            %assert damage rates ~= expected
            this.assertDamageRateIsExpected(c, d, expRates, iterMax, 0.4, 0.5);
        end

        function testExpectedDamageRates_evolveState_fullReactionSet(this)
            %references
            sim = this.simulation;
            c = sim.state('Chromosome');
            m = sim.state('Metabolite');
            s = sim.state('Stimulus');
            d = sim.process('DNADamage');
            r = sim.process('DNARepair');

            %simulate damage
            iterMax = 1000;

            d.reactionBounds(d.reactionBounds(:, 2) < 1000, 2) = 1e-6;
            d.substrates(:) = 1e12;
            d.substrates(d.reactionRadiation(d.reactionRadiation ~= 0)) = 1;
            s.setValues(s.getIndexs(d.substrateWholeCellModelIDs(d.reactionRadiation(d.reactionRadiation ~= 0))), 3) = 1;

            expRates = d.calcExpectedReactionRates() * iterMax;

            assertFalse(any(d.reactionDNAProduct == m.m6ADIndexs & ...
                strcmp(d.reactionDamageTypes, 'damagedBases')));

            d.reactionSmallMoleculeStoichiometryMatrix(:) = 0;

            r.reactionBounds(:, 2) = 0;
            r.enzymeBounds(:, 2) = 0;

            sim.applyOptions('lengthSec', iterMax);
            sim.state('MetabolicReaction').initialGrowthFilterWidth = Inf;
            sim.run();

            %assert damage rates ~= expected
            this.assertDamageRateIsExpected(c, d, expRates, iterMax, 0.1, 10);
        end

        function testAllDamageCanBeRepaired(this)
            sim = this.simulation;
            c = sim.state('Chromosome');
            d = sim.process('DNADamage');
            r = sim.process('DNARepair');

            repairRates = 10000 * r.enzymeBounds;

            for i = 1:numel(d.reactionWholeCellModelIDs)
                c.initialize();
                d.initializeState();
                r.initializeState();

                sim.state('Metabolite').counts(:, sim.compartment.cytosolIndexs) = 1e6;
                sim.state('Stimulus').values(:) = 0;
                sim.state('Stimulus').setValues(:, 3) = 0;

                d.reactionBounds(:, 2) = 0;

                switch d.reactionVulnerableMotifTypes{i}
                    case 'sequence'
                        d.reactionBounds(i, 2) = 3e-5 * 4 ^ (numel(d.reactionVulnerableMotifs{i})-1);
                    case {'damagedSugarPhosphate', 'damagedBases', 'intrastrandCrossLinks'}
                        c.(d.reactionVulnerableMotifTypes{i})((1:5)*10000, 1) = ...
                            d.reactionVulnerableMotifs{i};
                        c.(d.reactionVulnerableMotifTypes{i})((6:10)*10000, 2) = ...
                            d.reactionVulnerableMotifs{i};
                        d.reactionBounds(i, 2) = 1;
                    case {'gapSites', 'abasicSites', 'strandBreaks', 'hollidayJunctions'}
                        c.(d.reactionVulnerableMotifTypes{i})((1:5) * 10000, 1) = 1;
                        c.(d.reactionVulnerableMotifTypes{i})((6:10) * 10000, 2) = 1;
                        d.reactionBounds(i, 2) = 1;
                    otherwise, throw(MException('DNADamageRepair_Test:error', ...
                            'Unimplemented vulnerable motif %s', d.reactionVulnerableMotifTypes{i}));
                end

                if d.reactionRadiation(i)
                    sim.state('Stimulus').setValues(sim.state('Stimulus').setValues(:,1) == ...
                        d.substrateGlobalIndexs(d.reactionRadiation(i)), 3) = 1;
                end

                r.enzymeBounds(:) = 0;

                %stochastically damage
                for j = 1:2
                    sim.evolveState();
                end

                switch d.reactionDamageTypes{i}
                    case {'damagedSugarPhosphate', 'damagedBases', 'intrastrandCrossLinks'}
                        reactionDNAProduct = d.reactionDNAProduct(i);
                    case {'gapSites', 'abasicSites', 'strandBreaks', 'hollidayJunctions'}
                        reactionDNAProduct = 1;
                    otherwise, throw(MException('DNADamageRepair_Test:error', ...
                            'Unimplemented damage type %s', d.reactionDamageTypes{i}));
                end
                assertIn(nnz(c.(d.reactionDamageTypes{i}) == reactionDNAProduct), [2 50], ...
                    sprintf('%d %s %s\n', i, d.reactionWholeCellModelIDs{i}, d.reactionDamageTypes{i}));

                %repair
                d.reactionBounds(:, 2) = 0;
                r.enzymeBounds = repairRates;

                posStrnds = find(any(c.(d.reactionDamageTypes{i}) == reactionDNAProduct, 2));
                nDamagesMax = sum(diff([posStrnds(:, 1); posStrnds(1,1) + c.sequenceLen]) <= max(r.enzymeDNAFootprints));

                initDamagedSites = nnz(c.damagedSites);

                for j = 1:350
                    sim.evolveState();

                    damagedSites = c.damagedSites;
                    damagedSites(find(c.damagedBases)) = 0; %#ok<FNDSB>
                    damagedSites(c.shiftPositionsStrandsBase3Prime(find(c.damagedBases))) = 0; %#ok<FNDSB>
                    damagedSites(c.shiftPositionsStrandsBase5Prime(find(c.damagedBases))) = 0; %#ok<FNDSB>
                    if nnz(damagedSites) <= nDamagesMax
                        break;
                    end
                end

                assertIn(nnz(damagedSites), [0 nDamagesMax], ...
                    sprintf('%d %s %d\n', i, d.reactionWholeCellModelIDs{i}, initDamagedSites));
            end
        end

        function testDNARepairOutweighsDamage(this)
            %% references
            sim = this.simulation;
            c = sim.state('Chromosome');
            d = sim.process('DNADamage');
            r = sim.process('DNARepair');

            c.initialize();
            r.initializeState();

            %% indices of proteins which should never bind
            monGblIdxs = [];
            cpxGblIdxs = [];
            for i = 1:numel(sim.processes)
                p = sim.processes{i};
                monGblIdxs = [monGblIdxs; ...
                    p.enzymeGlobalIndexs(p.enzymeMonomerLocalIndexs)]; %#ok<*AGROW>
                cpxGblIdxs = [cpxGblIdxs; ...
                    p.enzymeGlobalIndexs(p.enzymeComplexLocalIndexs)];
            end
            monGblIdxs = sort(monGblIdxs);
            cpxGblIdxs = sort(cpxGblIdxs);

            monGblIdxs = setdiff(monGblIdxs, r.enzymeGlobalIndexs(r.enzymeIndexs_recombinationStrandExchange));
            cpxGblIdxs = setdiff(cpxGblIdxs, r.enzymeGlobalIndexs(r.enzymeIndexs_DisA));

            %% indices of allowed damage
            dmgSPIdxs = [];
            dmgBaseIdxs = sort([
                d.reactionDNAProduct(ismember(d.reactionDamageTypes, {'damagedBases'}))
                r.substrateGlobalIndexs(r.substrateIndexs_m6AD)
                ]);
            dmgCLIdxs = sort(...
                d.reactionDNAProduct(ismember(d.reactionDamageTypes, {'intrastrandCrossLinks'})));

            %% simulate
            iterMax = 250;

            d.reactionBounds = 1e4 * d.reactionBounds;
            r.enzymeBounds = 1e4 * r.enzymeBounds;

            sim.state('Metabolite').counts(:, sim.compartment.cytosolIndexs) = 1e6;

            nDamages = zeros(iterMax, 1);
            for i = 1:iterMax
                this.assertSecondChromosomeIgnored(c);
                this.assert1ChromsomePolymerized(c);
                this.assertCorrectProteinsBound(c, monGblIdxs, cpxGblIdxs);
                this.assertCorrectDamages(c, dmgSPIdxs, dmgBaseIdxs, dmgCLIdxs);

                sim.evolveState();

                nDamages(i) = nnz(c.damagedSites);
            end

            this.assertDamageDoesntAccumulate(c, r);

            d.reactionBounds(:, 2) = 0;

            for i = 1:20
                sim.evolveState();
            end

            assertEqual(0, nnz(c.damagedSites));
        end
    end

    %helper methods
    methods
        %check that second chromosome is never modified
        function assertSecondChromosomeIgnored(~, c)
            posStrnds = find(c.polymerizedRegions);     assertFalse(any(posStrnds(:, 2) > 2));
            posStrnds = find(c.linkingNumbers);         assertFalse(any(posStrnds(:, 2) > 2));

            posStrnds = find(c.gapSites);               assertFalse(any(posStrnds(:, 2) > 2));
            posStrnds = find(c.abasicSites);            assertFalse(any(posStrnds(:, 2) > 2));
            posStrnds = find(c.damagedSugarPhosphates); assertFalse(any(posStrnds(:, 2) > 2));
            posStrnds = find(c.damagedBases);           assertFalse(any(posStrnds(:, 2) > 2));
            posStrnds = find(c.strandBreaks);           assertFalse(any(posStrnds(:, 2) > 2));
            posStrnds = find(c.intrastrandCrossLinks);  assertFalse(any(posStrnds(:, 2) > 2));
            posStrnds = find(c.hollidayJunctions);      assertFalse(any(posStrnds(:, 2) > 2));

            posStrnds = find(c.monomerBoundSites);      assertFalse(any(posStrnds(:, 2) > 2));
            posStrnds = find(c.complexBoundSites);      assertFalse(any(posStrnds(:, 2) > 2));

            assertFalse(c.segregated);
        end

        %check that only first chromosome is polymerized
        function assert1ChromsomePolymerized(~, c)
            import edu.stanford.covert.util.CircularSparseMat;

            lk0 = c.sequenceLen / c.relaxedBasesPerTurn * (1 + c.equilibriumSuperhelicalDensity);
            assertEqual(CircularSparseMat([1 1; 1 2], [c.sequenceLen; c.sequenceLen], [c.sequenceLen c.nCompartments], 1), c.polymerizedRegions);
            assertEqual(CircularSparseMat([1 1; 1 2], [lk0; lk0], [c.sequenceLen c.nCompartments], 1), c.linkingNumbers);
        end

        %check that only allowed proteins ever bind chromosome
        function assertCorrectProteinsBound(~, c, monGblIdxs, cpxGblIdxs)
            [~, monVals] = find(c.monomerBoundSites);
            [~, cpxVals] = find(c.complexBoundSites);

            assertFalse(any(ismembc(monVals, monGblIdxs)));
            assertFalse(any(ismembc(cpxVals, cpxGblIdxs)));
        end

        %check damage is of allowed types
        function assertCorrectDamages(~, c, dmgSPIdxs, dmgBaseIdxs, dmgCLIdxs)
            [~, vals] = find(c.gapSites);               assertTrue(all(ismembc(vals, [false; true])));
            [~, vals] = find(c.abasicSites);            assertTrue(all(ismembc(vals, [false; true])));
            [~, vals] = find(c.damagedSugarPhosphates); assertTrue(all(ismembc(vals, dmgSPIdxs)));
            [~, vals] = find(c.damagedBases);           assertTrue(all(ismembc(vals, dmgBaseIdxs)));
            [~, vals] = find(c.strandBreaks);           assertTrue(all(ismembc(vals, [false; true])));
            [~, vals] = find(c.intrastrandCrossLinks);  assertTrue(all(ismembc(vals, dmgCLIdxs)));
            [~, vals] = find(c.hollidayJunctions);      assertTrue(all(ismembc(vals, [false; true])));
        end

        %check damage doesn't accumulate -- damage that occurs is repaired quickly
        function assertDamageDoesntAccumulate(~, c, r)
            assertIn(nnz(c.gapSites), [0 10 + r.NER_UvrABC_IncisionMargin3 + r.NER_UvrABC_IncisionMargin5]);
            assertIn(nnz(c.abasicSites), [0 10 + r.NER_UvrABC_IncisionMargin3 + r.NER_UvrABC_IncisionMargin5]);
            assertIn(nnz(c.damagedSugarPhosphates), [0 10]);
            assertIn(nnz(c.damagedBases), [0 10 + 2 * size(r.RM_MunI_RecognitionSites, 1)]);
            assertIn(nnz(c.strandBreaks), [0 10]);
            assertIn(nnz(c.intrastrandCrossLinks), [0 20]);
            assertIn(nnz(c.hollidayJunctions), [0 20]);
        end

        %check that damage rate follow expectations
        function assertDamageRateIsExpected(~, c, d, expRates, iterMax, relTol, absTol)
            [types, ~, idxs] = unique(d.reactionDamageTypes);
            reactionDNAProduct = d.reactionDNAProduct;
            reactionDNAProduct(ismember(d.reactionDamageTypes, ...
                {'strandBreaks','abasicSites','gapSites','hollidayJunctions'})) = 1;
            dmgTypesProds = unique([idxs reactionDNAProduct], 'rows');
            calcRates = zeros(size(dmgTypesProds, 1), 1);
            obsRates = zeros(size(dmgTypesProds, 1), 1);
            for i = 1:size(dmgTypesProds, 1)
                rxnTfs = ...
                    ismember(d.reactionDamageTypes, types(dmgTypesProds(i,1))) & ...
                    reactionDNAProduct == dmgTypesProds(i, 2);
                calcRates(i) = sum(expRates(rxnTfs));
                obsRates(i) = nnz(c.(types{dmgTypesProds(i, 1)}) == dmgTypesProds(i, 2));
            end

            reactionVulnerableMotifs = d.reactionVulnerableMotifs;
            reactionVulnerableMotifs(cellfun(@ischar, reactionVulnerableMotifs)) = {0};
            reactionVulnerableMotifs(cellfun(@islogical, reactionVulnerableMotifs)) = {1};
            reactionVulnerableMotifs = cell2mat(reactionVulnerableMotifs);
            rates = d.reactionBounds(:, 2);
            rates(d.reactionRadiation~=0) = ...
                rates(d.reactionRadiation~=0) .* ...
                d.substrates(d.reactionRadiation(d.reactionRadiation~=0));
            for i = 1:size(dmgTypesProds, 1)
                j = find(...
                    rates > 0 & ...
                    ismember(reactionVulnerableMotifs, dmgTypesProds(i, 2)) & ...
                    ismember(d.reactionVulnerableMotifTypes, types(dmgTypesProds(i, 1))));

                [~, k] = ismember([idxs(j) reactionDNAProduct(j)], dmgTypesProds, 'rows');

                if isempty(j)
                    continue;
                end

                calcRates(k) = calcRates(k) + min(1, iterMax * rates(j)) * calcRates(i);
                if all(dmgTypesProds(k,1) == dmgTypesProds(i,1))
                    calcRates(i) = max(0, 1 - iterMax * sum(rates(j))) * calcRates(i);
                end
            end

            %metIds = d.metabolite.wholeCellModelIDs(dmgTypesProds(:,2));
            %tfs =ismember(types(dmgTypesProds(:,1)), {'strandBreaks','abasicSites','gapSites','hollidayJunctions'});
            %metIds(tfs) = cell(sum(tfs), 1);
            %[types(dmgTypesProds(:,1)) metIds num2cell(calcRates) num2cell(obsRates)]

            assertElementsAlmostEqual(calcRates', obsRates', 'relative', relTol, absTol);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function MatingPool = MatingSelection(PopObj,Range,N)
% The mating selection of RSEA

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

	%% Calculate the convergence of each solution
    PopObj = (PopObj-repmat(Range(1,:),size(PopObj,1),1))./repmat(Range(2,:)-Range(1,:),size(PopObj,1),1);
    Con = sum(PopObj.^2,2).^0.5;

    %% Calculate the radar grid of each solution
    [Site,~] = RadarGrid(PopObj,ceil(sqrt(size(PopObj,1))));
    temp     = tabulate(Site);
    CrowdG   = temp(:,2);

    %% Binary tournament selection
    MatingPool = zeros(1,ceil(N/2)*2);
    grids      = TournamentSelection(2,length(MatingPool),CrowdG);
    for i = 1 : length(MatingPool)
        current = find(Site==grids(i));
        if isempty(current)
             MatingPool(i) = randi(size(PopObj,1),1);
        else
            parents       = current(randi(length(current),1,2));
            [~,best]      = min(Con(parents));
            MatingPool(i) = parents(best);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function DrawTextureHead(vertex, tri, color, keyPoints)
n = size(vertex, 2);

% draw
trisurf(tri', vertex(1, :), vertex(2, :), vertex(3, :), 1 : n, 'edgecolor', 'none');

if nargin == 4
    hold on
    plot3(keyPoints(1, :), keyPoints(2, :), keyPoints(3, :), 'g.', 'MarkerSize', 20)
end

% color
colormap(color');

shading interp
axis equal
%axis([1 250 1 250])
axis vis3d
title('fitted head');
xlabel('x');
ylabel('y');
zlabel('z');
view([0,90]);

author: XgTu


text:
% GSPBOX - Embeddings
%
%  Availlable embeddings
%    gsp_lle                       -  Local Linear Embedding
%    gsp_laplacian_eigenmaps       -  Laplacian Eigenmaps
%    gsp_isomap                    -  Isomap
%    gsp_eigenspace_estimation     -  Fast Eigenspace Approximation using Random Signals (FEARS)
%
%  Utils
%    gsp_weight2distance           -  Distance matrix from weight matrix
%    gsp_compute_coordinates       -  Compute new coordinates for a graph
%
%  For help, bug reports, suggestions etc. please send email to
%  gspbox 'dash' support 'at' groupes 'dot' epfl 'dot' ch
%

repo: gspbox


text:
function se3mat = VecTose3(V)
% *** CHAPTER 3: RIGID-BODY MOTIONS ***
% Takes a 6-vector (representing a spatial velocity).
% Returns the corresponding 4x4 se(3) matrix.
% Example Input:
%
% clear; clc;
% V = [1; 2; 3; 4; 5; 6];
% se3mat = VecTose3(V)
%
% Output:
% se3mat =
%     0    -3     2     4
%     3     0    -1     5
%    -2     1     0     6
%     0     0     0     0

se3mat = [VecToso3(V(1: 3)), V(4: 6); 0, 0, 0, 0];
end

author: ShuoYangRobotics

repo: QuadrupedSim


text:
function [] = visPatchLearning(RBM,iE,jB);
%-----------------------------------------
%  [] = visPatchLearning(RBM,jB,iE,rotFlag);
%-----------------------------------------
% DES

if RBM.useGPU
	RBM = gpuGather(RBM);
end

if notDefined('iE'),iE = numel(RBM.e); end

if isfield(RBM.auxVars,'invXForm');
	invXForm = RBM.auxVars.invXForm;
else
	invXForm = eye(size(RBM.W,1));
end

nVis = min(64,floor(sqrt(size(RBM.W/3,2))).^2);

subplot(331);
visPatches(RBM.X(RBM.batchIdx{jB},:)',invXForm);
title('Batch Data');

subplot(332);
visPatches(RBM.pVis',invXForm);
title('Fantasies');

subplot(333);
visPatches(RBM.dW(:,1:nVis),invXForm);
title ('Weight Gradients');

subplot(334);
hist(RBM.b);
title('Visible Bias');

subplot(335);
hist(RBM.aHid(:));
title(sprintf('E[hid]=%1.2f\nTarget Sparsity =%0.4f',mean(RBM.aHid(:)),RBM.sparsity))

subplot(336);
visPatches(RBM.W(:,1:nVis),invXForm);
title('Basis/Connection Weights');

subplot(337);
plot(RBM.log.err(1:iE));
title('Reconstruction errors');

subplot(338)
hist(RBM.W(:));
title('Connection Weights');

subplot(339);
semilogy(RBM.log.eta(1:iE));
title('Learning Rate');

drawnow

author: dustinstansbury

repo: medal


text:
function [A, b, C, d, results] = inflate_region(obstacles, A_bounds, b_bounds, start, varargin)
import iris.*;

p = inputParser();
p.addOptional('require_containment', false, @isscalar);
p.addOptional('error_on_infeasible_start', false, @isscalar);
p.addOptional('termination_threshold', 2e-2, @(x) x > 0);
p.addOptional('iter_limit', 100, @isnumeric);
p.parse(varargin{:});
options = p.Results;

if exist('+iris/inflate_regionmex', 'file')
  if ~iscell(obstacles)
    obs_cell = cell(1, size(obstacles, 3));
    for j = 1:size(obstacles, 3)
      obs_cell{j} = obstacles(:,:,j);
    end
    obstacles = obs_cell;

    % Note: could replace this with the following, but mat2cell is 100 times slower than the for loop
    % obstacles = mat2cell(obstacles, size(obstacles, 1), size(obstacles, 2), ones(1, size(obstacles, 3)));
  end

  if nargout > 4
    [A, b, C, d, p_history, e_history] = inflate_regionmex(obstacles, A_bounds, b_bounds, start, options);
    results = inflation_results();
    results.start = start;
    results.obstacles = obstacles;
    results.n_obs = numel(obstacles);
    results.e_history = e_history;
    results.p_history = p_history;
  else
    [A, b, C, d] = inflate_regionmex(obstacles, A_bounds, b_bounds, start, options);
  end
else
  if iscell(obstacles)
    padded = pad_obstacle_points(obstacles);
    obstacle_pts = cell2mat(reshape(padded, size(padded, 1), [], length(obstacles)));
  else
    obstacle_pts = obstacles;
  end
  [A, b, C, d, results] = inflate_region_fallback(obstacle_pts, A_bounds, b_bounds, start, options);
  results.obstacles = obstacles;
  if ~iscell(results.obstacles)
    results.obstacles = mat2cell(results.obstacles, size(results.obstacles, 1), size(results.obstacles, 2), ones(1, size(results.obstacles, 3)));
  end
end

author: rdeits

repo: iris-distro

save_path: github-repos/MATLAB/rdeits-iris-distro


text:
function h = plot3d(sF,varargin)
% plot spherical Function
%
% Syntax
%   plot3d(sF)
%
% See also
%   S2Fun/contour S2Fun/contourf S2Fun/pcolor S2Fun/plot
%

% plot the function values
h = plot(sF,'3d',varargin{:});

% remove output if not required
if nargout == 0, clear h; end

end

author: mtex-toolbox

repo: mtex

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


text:
function a = fourier_cosine ( n )

%*****************************************************************************80
%
%% FOURIER_COSINE returns the FOURIER_COSINE matrix.
%
%  Discussion:
%
%    FOURIER_COSINE is the discrete Fourier Cosine Transform matrix.
%
%  Example:
%
%    N = 5
%
%    0.447214      0.447214      0.447214      0.447214      0.447214
%    0.601501      0.371748      0.000000     -0.371748     -0.601501
%    0.511667     -0.195440     -0.632456     -0.195439      0.511667
%    0.371748     -0.601501      0.000000      0.601501     -0.371748
%    0.195439     -0.511667      0.632456     -0.511668      0.195439
%
%  Properties:
%
%    A * A' = I.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  a(1,1:n) = 1.0 / sqrt ( n );

  for i = 2 : n

    for j = 1 : n

      angle = ( i - 1 ) * ( 2 * j - 1 ) * pi / ( 2 * n );
      a(i,j) = sqrt ( 2.0 ) * cos ( angle ) / sqrt ( n );

    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function atlas_obj = atlas_add_L_R_to_labels(atlas_obj)
% Removes some strings indicating lateralization from atlas labels and adds new _L and _R suffixes for lateralized regions.
% - strings replaced are L_ R_ Left_ Right_ _L _R _Left _Right

labels = atlas_obj.labels;

mypatterns = {'L_' 'R_' 'Left_' 'Right_' '_L' '_R' '_Left' '_Right'};

for i = 1:length(mypatterns)
    labels = regexprep(labels, mypatterns{i}, '');
end

mysuffix = get_suffix(atlas_obj);

for i = 1:length(labels)

    labels{i} = [labels{i} mysuffix{i}];

end

atlas_obj.labels = labels;

end % function

function mysuffix = get_suffix(atlas_obj)

r = atlas2region(atlas_obj);

k = length(r);

for j = 1:k

    voxsign = sign(r(j).XYZmm(1, :)); % sign of each x coord. - = left
    modal_x = mode(voxsign);

    switch modal_x
        case -1
            labelstr = '_L';
        case 1
            labelstr = '_R';
        case 0
            labelstr = '';
    end

    % proportional asymmetry. 1 is vary lateralized. 0 is balanced
    % across L/R.  If symmetrical, then label _M
    prop_asym = abs(sum(voxsign == -1) - sum(voxsign == 1)) ./ length(voxsign);

    if prop_asym < .5, labelstr = ''; end

    mysuffix{j} = labelstr;

end % j

end % subfunction

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% Script demo_dim_info
% Exemplifies creation and usage of MrDimInfo class for retrieving and
% manipulating indices of multi-dimensional array
%
%  demo_dim_info
%
%
%   See also MrDimInfo

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2016-01-23
% Copyright (C) 2016 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

clear;
close all;
clc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1. Construct common dimInfo objects:
%   a) 4D EPI-fMRI array, with fixed TR
%   b) 5D multi-coil time series
%   c) 5D multi-echo time series
%   d) Create 5D multi-coil time series via nSamples and ranges
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  a) creates standard 4D dim info array from arraySize
% presets: units, dimLabels, resolutions
arraySize   = [64 50 33 100];
dimInfo     = MrDimInfo('nSamples', arraySize);
disp(dimInfo);

%% b) creates standard 5D dim info array from arraySize and resolutions
% presets of dimLabels, startingPoint = [1 1 1 1 1];
arraySize   = [64 50 33 100 8];
resolutions = [3 3 3 2.5 1];
units       = {'mm', 'mm', 'mm', 's', 'nil'};
dimInfo2    = MrDimInfo('nSamples', arraySize, 'resolutions', resolutions, ...
    'units', units);
disp(dimInfo2);
% no empty units allowed, will be overwritten by default!

%% c) creates standard 5D dim info array from arraySize, resolutions,
% startingPoints
% no presets
arraySize   = [64 50 33 8 3];
resolutions = [3 3 3 1 25];
units       = {'mm', 'mm', 'mm', 'nil', 'ms'};
dimLabels   = {'x', 'y', 'z', 'coil', 'echo_time'};
firstSamplingPoint = [-110, -110, -80, 0, 15];
dimInfo3    = MrDimInfo('nSamples', arraySize, 'resolutions', resolutions, ...
    'units', units, 'dimLabels', dimLabels, ...
    'firstSamplingPoint', firstSamplingPoint);
disp(dimInfo3.ranges);

%% d) Create 5D multi-coil time series via nSamples and ranges
% no presets, resolutions computed automatically
dimInfo4 = MrDimInfo(...
    'nSamples', [128 96 35 8 1000], ...
    'dimLabels', {'x', 'y', 'z', 'coil', 't'}, ...
    'units', {'mm', 'mm', 'mm', '', 's'}, ...
    'ranges', {[2 256], [2 192], [3 105], [1 8], [0 2497.5]});
disp(dimInfo4);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Get parameters of dimInfo via get_dims and dimInfo.'dimLabel'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% check out how neatly you can retrieve information about certain
% dimensions:
% dimInfo2.get_dims('x').resolutions
dimInfo2.x.resolutions
dimInfo2.get_dims('x')
dimInfo2.get_dims({'x' 'y' 'z'})

%% even cooler: use dimLabels directly for referencing!
dimInfo2.z.samplingPoints
dimInfo2.nSamples('z')
dimInfo2.nSamples({'z','y'})

dimInfo2('z')
dimInfo2({'z', 'y'})
dimInfo2([3 2])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 2. Modify dimInfo-dimensions via set_dims/add_dims-command
% a) Specify non-consecutive sampling-points (e.g. coil channels)
% b) Shift start sample of dimensions (e.g. centre FOV in x/y)
% c) Add a 6th dimension (e.g. additinal echoes)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  a) Specify non-consecutive sampling-points (e.g. coil channels)
disp(dimInfo3);
dimInfo3.set_dims('coil', 'samplingPoints', [2 3 4 7 8 10 11 12])
% Note that there is no concept of resolution here anymore, since there is
% equidistant spacing!
dimInfo3.resolutions
% However, samplingWidths is retained, since the data is still coming from
% one coil (a better way to think of it would be to select a number of
% slices; the resolution would be nan but each slice would still have a
% thickness = samplingWidth [which is the true slice thickness, resolution
% would be composed of samplingWidth + sliceGap])
dimInfo3.samplingWidths

% b) Shift start sample of dimensions (e.g. centre FOV in x/y)
disp(dimInfo4);
dimInfo4.set_dims([1 2], 'arrayIndex', [65 49], 'samplingPoint', [0 0]);
dimInfo4.ranges(:,1:2)

disp(dimInfo3);
dimInfo3.samplingPoints('coil') = {13:15};
dimInfo3.z.samplingPoints = {1:20};
dimInfo3.echo_time.nSamples = 5;

% c) Add a 6th dimension (e.g. additinal echoes)
disp(dimInfo2);
dimInfo2.add_dims(6, 'samplingPoints', [17 32], 'dimLabels', 'echo', 'units', 'ms');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 3. Display sampling points (=absolute indices with units) of
% selected first/center/last voxel
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

arrayIndexFirst     = [1,1,1,1,1];
arrayIndexLast      = [128 96 35 8 1000];
arrayIndexCenter    = [64, 48, 18, 4, 500];

arrayIndices = [
    arrayIndexFirst
    arrayIndexCenter
    arrayIndexLast
    ];

nVoxels = size(arrayIndices, 1);

samplingPointArray = dimInfo4.index2sample(...
    arrayIndices);

fprintf('===\ndimInfo.sample2index(arrayIndices): \n');
for iVoxel = 1:nVoxels
    fprintf('array Index, voxel %d:', iVoxel);
    fprintf('%5.1f  ', samplingPointArray(iVoxel,:));
    fprintf('\n');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 2. Display strings with units: sampling point of selected
%   first/center/last voxel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

indexLabelArray = dimInfo4.index2label(...
    arrayIndices);

fprintf('===\ndimInfo.index2label(arrayIndices): \n');
for iVoxel = 1:nVoxels
    fprintf('Voxel %d: ', iVoxel);
    fprintf('%s ',indexLabelArray{iVoxel}{:});
    fprintf('\n\n');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 3. Back-transform: Retrieve voxel index from absolute index
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

retrievedVoxelIndexArray = dimInfo4.sample2index(samplingPointArray);

fprintf('===\ndimInfo.sample2index(samplingPointArray): \n');
for iVoxel = 1:nVoxels
    fprintf('retrieved array index, voxel %d:', iVoxel);
    disp(retrievedVoxelIndexArray(iVoxel,:));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 4. dimInfo.select() - extract subset of dimension info from
%                        PropName/Value-pairs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  a) Select only data from a specific subset of one dimension (e.g. all
% data from some coils)
[selectionDimInfo, selectionIndexArray] = dimInfo4.select('type', 'index', ...
    'coil', [2 3 5 6]);

%% b) Select all data excluding subsets of one dimension (e.g. first volumes
% of time series data)
[selectionDimInfo2, selectionIndexArray2] = dimInfo4.select('type', 'index', ...
    'invert', true, 't', [1:5]);

%% c) Select subset of data array from multiple dimensions (e.g. cubic volume ROI,
% but see MrRoi for more sophisticated region definitions)
[selectionDimInfo3, selectionIndexArray3] = dimInfo4.select('type', 'index', ...
    'x', [20:40 80:100],  'y', [1:16 81:96], 'z', [10:15]);

%% d) Select subset of data array by specifying sampling-points
[selectionDimInfo4, selectionIndexArray4] = dimInfo4.select('type', 'sample', ...
    'x', [-128:2:0],  'y', [0:2:80]);

%% e) Combine selection into a nice structure, instead of
% ParameterName/Value-pairs (not all have to be given!)
selection.x = 1:10;
selection.t = 200:300;
[selectionDimInfo5, selectionIndexArray5] = dimInfo4.select(selection);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 5. dimInfo = MrDimInfo(fileName) - extract dimInfo directly from file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  3D Nifti
dataPath = tapas_uniqc_get_path('data');
niftiFile3D = fullfile(dataPath, 'nifti', 'rest', 'meanfmri.nii');
dimInfo3DFile = MrDimInfo(niftiFile3D);

% 4D Nifti
niftiFile4D = fullfile(dataPath, 'nifti', 'rest', 'fmri_short.nii');
dimInfo4DFile = MrDimInfo(niftiFile4D);

% 5D nifti
niftiFile5D = fullfile(dataPath, 'nifti', '5D', 'y_5d_deformation_field.nii');
dimInfo5DFile = MrDimInfo(niftiFile5D);

% several files in folder
niftiFolder = fullfile(dataPath, 'nifti', 'split', 'full');
dimInfoFolder = MrDimInfo(niftiFolder);

% par/rec
parRecFile = fullfile(dataPath, 'parrec', 'rest_feedback_7T', 'fmri1.par');
dimInfoParRec = MrDimInfo(parRecFile);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 6. Create MrDimInfo object from dimInfo struct (MrDimInfo(dimInfoStruct))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% creat MrDimInfo object
dimInfoObject = MrDimInfo('nSamples', [64 50 33 100], ...
    'units', {'u1', 'u2', 'u3', 'u4'}, ...
    'dimLabels', {'l1', 'l2', 'l3', 'l4'}, ...
    'resolutions', [0.3 2 0.85 5]);
disp(dimInfoObject);
% convert to struct (gives a warning)
dimInfoStruct = struct(dimInfoObject);
disp(dimInfoStruct);
% now create a new dimInfo object using the struct as input
newDimInfoObject = MrDimInfo(dimInfoStruct);

author: translationalneuromodeling

repo: tapas


text:
function b = backT(engine, f, t)

b.t = t;
b.obslik = f.obslik;
Q = length(f.alpha);
b.beta = ones(Q,1);
b.gamma = f.alpha;
if t > 1
  bb_t = b.obslik;
  b.xi = normalise((engine.transprob .* (f.past_alpha * bb_t'))); % T-1,T
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function dataset = ch_data_import(path)

tbl = readtable(path);

if sum(ismember(tbl.Properties.VariableNames,'AccX'))
    dataset.imu.acc(1,:) = tbl.AccX';
end

if sum(ismember(tbl.Properties.VariableNames,'AccY'))
    dataset.imu.acc(2,:) = tbl.AccY';
end

if sum(ismember(tbl.Properties.VariableNames,'AccZ'))
    dataset.imu.acc(3,:) = tbl.AccZ';
end

if sum(ismember(tbl.Properties.VariableNames,'GyrX'))
    dataset.imu.gyr(1,:) = tbl.GyrX';
end

if sum(ismember(tbl.Properties.VariableNames,'GyrY'))
    dataset.imu.gyr(2,:) = tbl.GyrY';
end

if sum(ismember(tbl.Properties.VariableNames,'GyrZ'))
    dataset.imu.gyr(3,:) = tbl.GyrZ';
end

if sum(ismember(tbl.Properties.VariableNames,'MagX'))
    dataset.imu.mag(1,:) = tbl.MagX';
end

if sum(ismember(tbl.Properties.VariableNames,'MagY'))
    dataset.imu.mag(2,:) = tbl.MagY';
end

if sum(ismember(tbl.Properties.VariableNames,'MagZ'))
    dataset.imu.mag(3,:) = tbl.MagZ';
end

if sum(ismember(tbl.Properties.VariableNames,'Roll'))
    dataset.eul.roll = tbl.Roll';
end

if sum(ismember(tbl.Properties.VariableNames,'Pitch'))
    dataset.eul.pitch = tbl.Pitch';
end

if sum(ismember(tbl.Properties.VariableNames,'Yaw'))
    dataset.eul.yaw = tbl.Yaw';
end

author: yandld


text:
% TRIAL2EEGPLOT - convert eeglab format to eeplot format of rejection window
%
% Usage:
%   >> eegplotarray = trial2eegplot(rej, rejE, points, color);
%
% Inputs:
%   rej   - rejection vector (0 and 1) with one value per trial
%   rejE  - electrode rejection array (size nb_elec x trials) also
%           made of 0 and 1.
%   points - number of points per trials
%   color  - color of the windows for EEGPLOT
%
% Outputs:
%   eegplotarray - array defining windows which is compatible with
%                  the function EEGPLOT
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: EEGTRESH, EEGLAB, EEGPLOT, POP_REJEPOCH

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

% ----------------------------------------------------------
function rejeegplot = trial2eegplot( rej, rejE, pnts, color)
	rej  = find(rej>0);
	rejE = rejE(:, rej)';
   	rejeegplot = zeros(length(rej), size(rejE,2)+5);
   	rejeegplot(:, 6:end) = rejE;
   	rejeegplot(:, 1) = (rej(:)-1)*pnts;
   	rejeegplot(:, 2) = rej(:)*pnts-1;
   	rejeegplot(:, 3:5) = ones(size(rejeegplot,1),1)*color;
return

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function c = rdivide(a,b)
% implements e ./ c and c ./ e
%
% Description
% If |e| is a matrix of embeddings and |c| is a matrix of coefficients
% then |e ./ c| is again a matrix of embeddings defined by
%
% $$ [\mathrm{e ./ c}]_{j\ell} = mathrm{e}_{j\ell} / \mathrm{c}_{j \ell}$$
%
% Syntax
%   out = e ./ c
%   out = c ./ e
%   out = e1 ./ e2
%
% Input
%  e, e1, e2 - @embedding
%  c - double
%
% Output
%  out- @embedding
%

if isa(a,'embedding')

  if isa(b,'embedding')

    for i = 1:length(a.u), a.u{i} = a.u{i} ./ b.u{i}; end

  else
    for i = 1:length(a.u), a.u{i} = a.u{i} ./ b; end
  end
  c = a;

else

  for i = 1:length(b.u), b.u{i} = a ./ b.u{i}; end
  c = b;

end

end


text:
%> @brief Convert a general cell array matrix to a string
%> @author Eric Cousineau <eacousineau@gmail.com>, AMBER Lab, under Dr.
%> Aaron Ames
%> @param X cell array
%> @param varargin Arguments to propagate to general2str()
function [str] = cell2str(X, varargin)
s = size(X);
func = @(x) general2str(x, varargin{:});
raw = cellfun(func, X, 'UniformOutput', false);
rows = cell(s(1), 1);
for i = 1:s(1)
    rows{i} = implode(raw(i, :), ', ');
end
str = ['{\n', indent(strjoin(rows', ';\n'), '\t'), ';\n}'];
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


%==============================================================================
% This code is part of the Matlab-based toolbox
% LagLDDDM - A Lagrangian Gauss--Newton--Krylov Solver for Mass- and
%                        Intensity-Preserving Diffeomorphic Image Registration
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM
%
% function [yc,dy,para] = getTrafoFromInstationaryVelocityRK4(vc,yc,varargin)
%
% compute transformation yc by integrating velocity field in time using a
% Runge-Kutta 4 method with fixed time steps. Velocity is time dependent
% here and assumed to be nodal in time. For stationary velocities, use
% getTrafoFromVelocityRK4.m
%
% For more details see Sec. 3 of the paper:
%
% @article{MangRuthotto2017,
%   Title = {A {L}agrangian {G}auss--{N}ewton--{K}rylov solver for mass- and intensity-preserving diffeomorphic image registration},
%   Year = {2017},
%   Journal = {SIAM Journal on Scientific Computing},
%   Author = {A. Mang, L. Ruthotto},
% }
%
% Input:
%  vc   - discrete velocity field (nodal, cell-centered, or staggered in space / nodal in time)
%  yc   - particle positions
%
% Additional REQUIRED Input (provided through varargin)
%
%  omega        - spatial domain (required!)
%  m            - number of cells in each direction (required!)
%
% Optional Input (provided through varargin)
%
%  doDerivative - compute derivative w.r.t. velocity
%  tspan        - time interval (default: [0 1])
%  N            - number of time steps (default: 5)
%  nt           - number of time points for velocity (default: compute from input)
%  storeInter   - store intermediate transformation (e.g., for visualization)
%
% Output:
%
%  yc           - end point of characteristics
%  dy           - derivative w.r.t. vc
%  para         - struct containing info
%
% =========================================================================
function [yc,dy,para] = getTrafoFromInstationaryVelocityRK4(vc,yc,varargin)

if nargin==0,
    runMinimalExample
    return;
end
doDerivative = (nargout>1);
omega = [];
m     = [];
tspan = [0,1];
N     = 5;
nt    = [];
storeInter = false;
for k=1:2:length(varargin)     % overwrites default parameter
    eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;
if isempty(omega) || isempty(m)
    error('%s - omega and m must be provided through varargin or trafo(''set'',''omega'',omega,''m'',m')
end
dim = numel(omega)/2;
if isempty(nt) % roughly estimate nt
    nt = round(numel(vc)/(prod(m)*dim))-1;
end
ht    = (tspan(2)-tspan(1))/nt;
tspan = (linspace(tspan(1),tspan(2),N)-min(tspan))/abs(ht)+1; % map to nodal indices
dt    = diff(tspan);
dt    = dt(1);
dy    = [];
if doDerivative
    dy = cell(nt+1,1);
    for j=1:nt+1, dy{j} = sparse(numel(yc),prod(m)*dim); end;
    dyi = dy;
end

para = struct('dt',dt,'N',N,'nt',nt,'omega',omega,'m',m);
if storeInter
    para.YC = zeros(numel(yc),N);
    para.YC(:,1) = yc;
end

vc = reshape(vc,[],nt+1);
for k=1:N-1
    % get v1 = v(yc,t(k));
    tk = tspan(k);
    p  = floor(tk); x = tk-p;  % split x into integer/remainder
    dvi      = zeros(nt+1,1);
    dvi(p)   = 1-x;
    if p<nt+1, dvi(p+1) = x; end
    vt = vc*dvi;
    [vi,dydvi] = linearInterGrid(vt,omega,m,yc,'doDerivative',doDerivative);
    if doDerivative
        Ty = getLinearInterGridMatrix(omega,m,yc);
        for j=1:nt+1
           dyi{j} = dvi(j)*Ty + dydvi*dy{j};
        end
        dytemp  = dyi;
    end
    yi = yc + .5*dt*vi;
	ytemp = vi;

    % get v2 = v(y1,.5*(t(k+1)+t(k)));
    tk1 = (tspan(k+1)+tspan(k))/2;
    p  = floor(tk1); x = tk1-p;  % split x into integer/remainder
    dvi      = zeros(nt+1,1);
    dvi(p)   = 1-x;
    if p<nt+1, dvi(p+1) = x; end
    vt = vc*dvi;
    [vi,dydvi] = linearInterGrid(vt,omega,m,yi,'doDerivative',doDerivative);
    if doDerivative
        Ty = getLinearInterGridMatrix(omega,m,yi);
        for j=1:nt+1
            dyi{j}    = dvi(j)*Ty + dydvi*(dy{j}+.5*dt*dyi{j});
            dytemp{j} = dytemp{j} + 2*dyi{j};
        end
    end
    yi = yc + .5*dt*vi;
	ytemp = ytemp + 2*vi;

    % get v3 = v(y2,.5*(t(k+1)+t(k)));
    [vi,dydvi] = linearInterGrid(vt,omega,m,yi,'doDerivative',doDerivative);
    if doDerivative
        Ty = getLinearInterGridMatrix(omega,m,yi);
        for j=1:nt+1
            dyi{j}    = dvi(j)*Ty + dydvi*(dy{j}+.5*dt*dyi{j});
            dytemp{j} = dytemp{j} + 2*dyi{j};
        end
    end
    yi  = yc + dt*vi;
	ytemp = ytemp + 2*vi;

    % get v4 = v(y3,t(k+1)));
    tk1 = tspan(k+1);
    p  = floor(tk1); x = tk1-p;  % split x into integer/remainder
    dvi      = zeros(nt+1,1);
    dvi(p)   = 1-x;
    if p<nt+1, dvi(p+1) = x; end
    vt = vc*dvi;
    [vi,dydvi] = linearInterGrid(vt,omega,m,yi,'doDerivative',doDerivative);
    if doDerivative
        Ty = getLinearInterGridMatrix(omega,m,yi);
        for j=1:nt+1
            dyi{j} = dvi(j)*Ty + dydvi*(dy{j}+dt*dyi{j});
            dytemp{j} = dytemp{j} + dyi{j};
            dy{j}  = dy{j} + (dt/6)*dytemp{j};
        end
    end
    ytemp = ytemp + vi;
    yc    = yc + (dt/6)*ytemp;
    if storeInter, para.YC(:,k+1) = yc; end
end
if doDerivative, dy = horzcat(dy{:}); end;

function runMinimalExample

omegaV = [-1 1 -1 1];
omega = .1*omegaV;
m     = [32 32];
tspan     = [8 0];
N     = 20;
yc    = getNodalGrid(omega,m)+.0002;

%
% % test cell-centered grid
regularizer('set','regularizer','mbCurvature','alpha',1)
xc      = reshape(getCellCenteredGrid(omegaV,m),[],2);
v0      = [.3*sign(xc(:,1)).*xc(:,1).^2; sin(pi*xc(:,2))];
vc      = v0*[1,1.1,1.5,2];
fctn    = @(vc) getTrafoFromInstationaryVelocityRK4(vc(:),yc,'omega',omegaV,'m',m,'T',tspan,'N',N);
checkDerivative(fctn,vc(:))


text:
function [TB] = GB2TB(GB)
% Convert computery things from gigabytes to terabytes.
% Chad A. Greene 2012
TB = GB*2^-10 ;


text:
function y = quad_form( x, Q, v, w )

%QUAD_FORM quadratic form.
%   QUAD_FORM(x,Q) is real(x'*Q*x) = x'*((Q+Q')/2)*x.
%   QUAD_FORM(x,Q,v,w) is real(x'*(Q*x+v)+w).
%
%   x must be a row or column vector, and Q must either be a scalar or
%   a square matrix with the same number of rows as x. If supplied, v must
%   be a scalar or a vector of the same size as x, and w must be a scalar.
%
%   NOTE: The use of QUAD_FORM can often be replaced by a call to NORM. For
%   example, if Q is positive definite, then the constraint
%       quad_form(x,Q) <= 1
%   is equivalent to
%       norm(sqrtm(Q)*x) <= 1
%   Generally speaking, the NORM version will be more reliable and more
%   accurate, so we encourage you to make similar conversions whenever
%   possible. We *strongly* discourage the QP-era practice of converting
%   NORM expressions into quadratic forms.
%
%   Disciplined convex programming information:
%       QUAD_FORM(x,Q,v,w) is neither convex nor concave in x and (Q,v)
%       jointly, so at least one of the two must be constant.
%
%       If (Q,v) is constant, then QUAD_FORM is convex if Q is positive
%       semidefinite, and concave if Q is negative semidefinite. An error
%       is generated if Q is indefinite (unless x is also constant).
%       QUAD_FORM is nonmonotonic in x, so x must be affine.
%
%       If x is constant, then QUAD_FORM is affine in Q, v, and w. The
%       signs of x will govern whether the elements of Q, v, and w may
%       be convex, concave, or affine.

narginchk(2,4);
sx = size( x );
if length( sx ) ~= 2 || all( sx ~= 1 ),
    error( 'The first argument must be a vector.' );
else
    sx = prod( sx );
end
if ndims( Q ) > 2 || size( Q, 1 ) ~= size( Q, 2 ), %#ok
    error( 'The second argument must be a scalar or a square matrix.' );
elseif all( size( Q, 1 ) ~= [ 1, sx ] ),
    error( 'The size of Q is incompatible with the size of x.' );
end
if nargin < 3,
    v = 0;
elseif ndims( v ) > 2 || all( size( v ) ~= 1 ), %#ok
    error( 'The third argument must be a vector.' );
elseif all( numel( v ) ~= [ 1, sx ] ),
    error( 'The size of v is incompatible with the size of x.' );
end
if nargin < 4,
    w = 0;
elseif numel( w ) > 1,
    error( 'The fourth argument must be a real scalar.' );
end
if sx == 0,
    y = real( w );
else
    x = x( : );
    v = v( : );
    y = real( x' * ( Q * x + v ) + w );
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function A = slinitarray(type, size)
%SLINITARRAY Initialize an array of specified type and size
%
% $ Syntax $
%   - A = slinitarray(type, size)
%
% $ Arguments $
%   - type          the string representing the element type of the array
%   - size          the vector representing the size of the array
%   - A             the initialized array
%
% $ Description $
%   - A = slinitarray(type, size) creates an array of specified type and
%     given size, with all elements being zeros.
%
%   - The typenames supported are listed below:
%       - 'double'
%       - 'single'
%       - 'float'
%       - 'uint8'
%       - 'uint16'
%       - 'uint32'
%       - 'uint64'
%       - 'int8'
%       - 'int16'
%       - 'int32'
%       - 'int64'
%       - 'logical'
%       - 'char'
%
% $ History $
%   - Created by Dahua Lin on Dec 7th, 2005
%

switch type
    case 'double'
        A = zeros(size);
    case {'single', 'float'}
        A = zeros(size, 'single');
    case {'uint8', 'uint16', 'uint32', 'uint64', ...
          'int8', 'int16', 'int32', 'int64'}
        A = zeros(size, type);
    case 'logical'
        A = false(size);
    case 'char'
        A = char(zeros(size, 'uint8'));
    otherwise
        error('sltoolbox:invalid_type', ...
            'Unsupported typename %s', typename);
end

author: lmthang


text:
% KSVD running file
% in this file a synthetic test of the K-SVD algorithm is performed. First,
% a random dictionary with normalized columns is being generated, and then
% a set of data signals, each as a linear combination of 3 dictionary
% element is created, with noise level of 20SNR. this set is given as input
% to the K-SVD algorithm.

% a different mode for activating the K-SVD algorithm is until a fixed
% error is reached in the Sparse coding stage, instead until a fixed number of coefficients is found
% (it was used by us for the
% denoising experiments). in order to switch between those two modes just
% change the param.errorFlag (0 - for fixed number of coefficients, 1 -
% until a certain error is reached).

param.L = 3;   % number of elements in each linear combination.
param.K = 50; % number of dictionary elements
param.numIteration = 50; % number of iteration to execute the K-SVD algorithm.

param.errorFlag = 0; % decompose signals until a certain error is reached. do not use fix number of coefficients.
%param.errorGoal = sigma;
param.preserveDCAtom = 0;

%%%%%%% creating the data to train on %%%%%%%%
N = 1500; % number of signals to generate
n = 20;   % dimension of each data
SNRdB = 20; % level of noise to be added
[param.TrueDictionary, D, x] = gererateSyntheticDictionaryAndData(N, param.L, n, param.K, SNRdB);
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%% initial dictionary: Dictionary elements %%%%%%%%
param.InitializationMethod =  'DataElements';

param.displayProgress = 1;
disp('Starting to  train the dictionary');

[Dictionary,output]  = KSVD(D,param);

disp(['The KSVD algorithm retrived ',num2str(output.ratio(end)),' atoms from the original dictionary']);

[Dictionary,output]  = MOD(D,param);

disp(['The MOD algorithm retrived ',num2str(output.ratio(end)),' atoms from the original dictionary']);

author: lbasek

repo: image-denoising-benchmark


text:
% mff_fileio_write - export Fieldtrip structures to MFF file
%
% Usage:
%   mff_fileio_write(mffFile, header, data, event);
%
% Inputs:
%   mffFile - [string] MFF file name
%   header   - fieldtrip data header
%   data     - fieldtrip raw data
%   event    - fieldtrip event structure (optional)

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

function mff_fileio_write(filename, hdr, data, evt)

if exist('pop_fileio')
    a = dir(which('pop_fileio.m'));
    b = dir(which('pop_fileio2.m'));
    if isempty(a) || isempty(b) || a.bytes ~= b.bytes
        fprintf([ '-----------------------------------------------------------------------------\n' ...
                  'WARNING: pop_fileio.m and pop_fileio2.m differ; likely you are using\n' ...
                  '         an old version of EEGLAB; The importer will try to import data using\n' ...
                  '         the pop_fileio.m function, and revert to pop_fileio2.m if it crashes\n' ...
                  '-----------------------------------------------------------------------------\n' ]);
    end
    try
        EEG = pop_fileio(hdr, data, evt); % maybe the wrong version?
    catch
        EEG = pop_fileio2(hdr, data, evt); % backup function included here
    end
else
    EEG = pop_fileio2(hdr, data, evt); % backup function included here
end

if isfield(hdr, 'orig')
    if isfield(hdr.orig, 'etc')
        EEG.etc = hdr.orig.etc;
    end
end

mff_export(EEG, filename);

% add code for eeg_emptyset
% add code for pop_fileio
% add code for finputcheck

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function indx = r8vec_sort_insert_index_a ( n, a )

%*****************************************************************************80
%
%% R8VEC_SORT_INSERT_INDEX_A ascending index sorts an R8VEC using insertion.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Algorithm 1.1,
%    Donald Kreher and Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998, page 11.
%
%  Parameters:
%
%    Input, integer N, the number of items in the vector.
%    N must be positive.
%
%    Input, real A(N), the array to be sorted.
%
%    Output, integer INDX(N), the sorted indices.  The array is sorted
%    when listed from A(INDX(1)) through A(INDX(N)).
%
  if ( n < 1 )
    return
  end

  indx = i4vec_indicator1 ( n );

  for i = 2 : n

    x = a(i);

    j = i - 1;

    while ( 1 <= j )

      if ( a(indx(j)) <= x )
        break
      end

      indx(j+1) = indx(j);
      j = j - 1;

    end

    indx(j+1) = i;

  end

  return
end


text:
function [q, pval] = ljungbox(data, lags)
% Ljung-Box tests for the presence of serial correlation in up to q lags.  Returns LAGS Ljung-Box
% statistics tests, one for tests for each lag between 1 and LAGS. Under the null of no serial
% correlation and assuming homoskedasticity, the Ljung-Box test statistic is asymptotically
% distributed X2(q)
%
% USAGE:
%  [Q,PVAL] = ljungbox(DATA,LAGS)
%
% INPUTS:
%  DATA      - A T by 1 vector of data
%  LAGS      - The maximum number of lags to compute the LB.  The statistic and pval will be
%                returned for all sets of lags up to and including LAGS
%
% OUTPUTS:
%  Q         - A LAGS by 1 vector of Q statistics
%  PVAL      - A LAGS by 1 set of appropriate pvals
%
% COMMENTS:
%  This test statistic is common but often inappropriate since it assumes homoskedasticity.  For a
%  heteroskedasticity consistent serial correlation test, see lmtest1
%
% SEE ALSO:
%  LMTEST1, SACF, SPACF

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 1/1/2007

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin~=2
    error('2 inputs required.')
end
T=length(data);
if T<=lags
    error('At least LAGS observations requires')
end
if size(data,1)~=T,
    data=data';
end
if size(data,2)~=1
    error('DATA must be a column vector')
end
if ~isscalar(lags) && lags>0 && floor(lags)==lags
    error('LAGS must be a positive integer')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ac = sacf(data,lags,0,0);

q = zeros(lags,1);
T = length(data);
for L=1:lags
    q(L)=T*(T+2)*sum(ac(1:L).^2./(T-(1:L)'));
end
pval = 1 - chi2cdf(q,(1:lags)');

author: bashtage

repo: mfe-toolbox


text:
function [u,flag,relres,iter,resvec] = minres(L,f,tol,maxit,R1,R2,u0,varargin)
%MINRES    A Preconditioned Minimal Residual method for ODEs.
%   U = MINRES(L,F) attempts to solve the linear ODE L(U) = F on DOMAIN(L) with
%   boundary conditions. The CHEBOP L must be a self-adjoint, uniformly elliptic
%   second-order differential operator of the form
%          L(U) = (a(x)*U')' + c(x)U,
%   with dirichlet boundary conditions. The righthand side F should be a chebfun
%   on DOMAIN(L). By default the indefinite integral operator is employed as a
%   preconditioner.
%
%   U = MINRES(L,F,TOL) specifies the tolerance of the method. If TOL is [] then
%   MINRES uses the default in cheboppref.
%
%   U = MINRES(L,F,TOL,MAXIT) specifies the maximum number of iterations. If
%   MAXIT is [] then MINRES uses the default in cheboppref.
%
%   U = MINRES(L,F,TOL,MAXIT,R1,R2) solves the preconditioned linear ODE of
%   (R2*L*R1)(V) = R2*f, where R2 must be the adjoint of R1. R1 and R2 must be
%   function handles. If R1 = [], then the default preconditioner is employed.
%   Only the default preconditioner is currently supported.
%
%   U = MINRES(L,F,TOL,MAXIT,R1,R2,U0) specifies the initial guess. If U0 is []
%   then MINRES MINRES the default, the zero function on DOMAIN(L).
%
%   [U,FLAG] = MINRES(L,F,...) also returns a convergence FLAG:
%    0 MINRES converged to the desired tolerance TOL within MAXIT iterations
%    1 MINRES iterated MAXIT times but did not converge.
%    2 preconditioner R1 was an unbounded operator.
%    3 MINRES stagnated (two consecutive iterates were the same).
%    4 one of the scalar quantities calculated during MINRES became too
%      small or too large to continue computing.
%
%   [U,FLAG,RELRES] = MINRES(L,F,...) also returns the relative residual
%    NORM(R2(F)-R2(L(U)),2)/NORM(F,2). If FLAG is 0, then RELRES <= TOL.
%
%   [U,FLAG,RELRES,ITER] = MINRES(L,F,...) also returns the iteration number at
%   which U was computed: 0 <= ITER <= MAXIT.
%
%   [U,FLAG,RELRES,ITER,RESVEC] = MINRES(L,F,...) returns a vector of estimated
%   residual norms at each iteration including NORM(R2(F)-R2(L(Uk)),2).
%
% See also CHEBOP/PCG and CHEBOP/GMRES.

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

% Only continue if L is a linear chebop:
if ( ~all(islinear(L)) )
    error('CHEBFUN:CHEBOP:pcg:nonlinear', ...
        'MINRES supports only linear CHEBOP instances.');
end

% At the moment, we need a second-order differential equation:
LinearOp = linop( L );
if ( LinearOp.diffOrder ~=2 )
    error('CHEBFUN:CHEBOP:pcg:DiffOrder', ...
        'MINRES supports only second-order ODEs.');
end

if ( nargin < 2 )
    error(message('chebop:minres:NotEnoughInputs'));
end

% Norm of righthand side function:
n2f = norm(f, 2);
n = length( f );
dom = domain( f );

% Grab boundary conditions from chebop:
if ( isa( L, 'chebop' ) )
    if isempty(L.lbc)
        left_bc = 0;
    else
        left_bc = L.lbcShow;

        % Current implementation requires Dirichlet boundary conditions:
        if ( ~isa(left_bc, 'double' ) )
            error('CHEBFUN:CHEBOP:pcg:leftbc', ...
                   'Currently, we require Dirichlet boundary conditions. Please supply N.lbc = double.');
        end
    end

    if ( isempty(L.rbc) )
        right_bc = 0;
    else
        right_bc = L.rbcShow;
        % Current implementation requires Dirichlet boundary conditions:
        if ( ~isa(right_bc, 'double' ) )
            error('CHEBFUN:CHEBOP:pcg:rightbc', ...
                   'Currently, we require Dirichlet boundary conditions. Please supply N.rbc = double.');
        end
    end

    L = L.op;
elseif ( isa( L, 'function_handle' ) )
    left_bc = 0;
    right_bc = 0;
else
    error(message('chebop:pcg:DiffOperatorIllDefined'));
end

% Ensure that L = @(x,u) ...
if ( nargin(L) == 1 )
    L = @(x,u) L(u);
elseif ( nargin(L) ~=2 )
    error(message('chebop:pcg:DiffOpNargin'));
end

% Data mine L for faster "mat-vec" computations in the CG iteration:
x = chebfun( @(x) x );
c = L(x, 1+0*x);
a = -L(x, x.^2/2) - (-L(x, x) + c.*x).*x + c.*(x.^2/2);
L = @(v) -diff( a.*diff( v ) ) + c.*v;

% Assign default values to unspecified parameters
if ( (nargin < 3) || isempty(tol) )
    tol = cheboppref().bvpTol;
end
warned = 0;
if ( tol <= eps )
    warning(message('chebop:pcg:tooSmallTolerance'));
    warned = 1;
    tol = eps;
elseif ( tol >= 1 )
    warning(message('chebop:pcg:tooBigTolerance'));
    warned = 1;
    tol = 1-eps;
end
if ( (nargin < 4) || isempty(maxit) )
    maxit = cheboppref().maxIter;
end

if ( (nargin < 5) || isempty(R1) )
    % Preconditioner:
    R1 = @(u) cumsum(u);
else
    % The method only works with the default preconditioner:
    error(message('chebop:pcg:OnlyDefaultPreconditionerAllowed'))
end

if ( (nargin < 6) || isempty(R2) )
    % Adjoint operator to R1:
    R2 = @(u) sum(u) - cumsum(u);
else
    % The method only works with the default preconditioner:
    error(message('chebop:pcg:OnlyDefaultPreconditionerAllowed'))
end

% Projection operator
Pi = @(g) g - mean(g);

% Build modified differential operator:
x = chebfun( @(x) x, dom );
T = @(v) Pi( R2( L( R1( v ) ) ) );

if ( (nargin >= 7) && ~isempty(u0) )
    if ( ~domainCheck(f, u0) )
        error(message('chebop:pcg:WrongInitGuessDomain'));
    else
        u = u0;
    end
    Tu = T(u);
else
    % Otherwise, initial guess is the zero function:
    u = 0*f;
    Tu = u;
end

if ( nargin > 7 )
    error(message('chebop:minres:TooManyInputs'));
end

% Ensure that rhs is within the correct space, and if not then solve a modified
% problem.
R2f = R2( f );
PiR2f = Pi( R2f );
if ( norm( R2f - PiR2f ) > tol || norm(left_bc)> tol || norm(right_bc) > tol )
    % f is NOT in the space Wn = {v in L_2: R1(v) in V_{n,0} }:
    basis = x.^(0:4);
    A = zeros(4, size(basis,2));
    for jj = 1:size(basis,2)
        A(1:2,jj) = feval( R1( R2( L(basis(:,jj)) ) ), dom([1,length(dom)])');
        A(3:4,jj) = feval(basis(:,jj),dom([1,length(dom)])');
    end
    b = [ feval( R1( R2f ), dom([1,length(dom)])') ; left_bc ; right_bc];
    warning('off','all')
    z = basis*(A\b);
    warning('on','all');
    % g is now in the space Wn = {v in L_2: R1(v) in V_{n,0} }:
    g = Pi( R2f - R2( L(z) ) );
else
    % f is in the space Wn = {v in L_2: R1(v) in V_{n,0} }:
    g = PiR2f;
    z = 0*f;
end

% Set up for the method
flag = 1;
iter = 0;
umin = u;                          % Iterate which has minimal residual so far
imin = 0;                          % Iteration at which xmin was computed
tolg = tol * normest(g);           % Relative tolerance
r = g - Tu;
normr = norm(r,2);                 % Norm of residual
normr_act = normr;

if ( normr <= tolg )               % Initial guess is a good enough solution
    flag = 0;
    relres = normr / n2f;
    resvec = normr;
    % Undo preconditioner:
    u = R1( Pi( u ) );
    return
end

resvec = zeros(maxit+1,1);         % Preallocate vector for MINRES residuals
resvec(1) = normr;                 % resvec(1) = norm(b-A*x0)
normrmin = normr;                  % Norm of minimum residual

vold = r;
v = vold;
beta1 = innerProduct(vold, v);
if ( beta1 <= 0 )
    flag = 5;
    relres = normr / n2b;
    resvec = resvec(1);
    u = R1( Pi( u ) );
    return
end
beta1 = sqrt(beta1);
snprod = beta1;
vv = v / beta1;
v = T(vv);
Amvv = v;
alpha = innerProduct(vv,v);
v = v - (alpha/beta1) * vold;

% Local reorthogonalization
numer = innerProduct(vv, v);
denom = innerProduct(vv, vv);
v = v - (numer/denom) * vv;
volder = vold;
vold = v;
betaold = beta1;
beta = innerProduct(v,v);
if ( beta < 0 )
    flag = 5;
    relres = normr / n2b;
    resvec = resvec(1);
    u = R1( Pi( u ) );
    return
end
iter = 1;
beta = sqrt(beta);
gammabar = alpha;
epsilon = 0;
deltabar = beta;
gamma = sqrt(gammabar^2 + beta^2);
mold = chebfun('0');
Amold = mold;
m = vv / gamma;
Am = Amvv / gamma;
cs = gammabar / gamma;
sn = beta / gamma;
u = u + snprod * cs * m;
snprod = snprod * sn;

normr = abs(snprod);
resvec(2,1) = normr;

% Check for convergence after first step.
if ( normr <= tolg )
    flag = 0;
    relres = normr / n2f;
    resvec = resvec(1:2);
    u = R1( Pi( u ) );
    return
end

stag = 0;                          % stagnation of the method
moresteps = 0;
maxmsteps = min([floor(n/50),5,n-maxit]);
maxstagsteps = 3;

% loop over maxit iterations (unless convergence or failure)
for ii = 2:maxit

    vv = v * (1/beta);
    v = T(vv);
    Amolder = Amold;
    Amold = Am;
    Am = v;
    v = v - (beta / betaold) * volder;
    alpha = innerProduct(vv, v);
    v = v - (alpha / beta) * vold;
    volder = vold;
    vold = v;
    betaold = beta;
    beta = innerProduct(v, v);
    if (beta < 0)
        flag = 5;
        break
    end
    beta = sqrt(beta);
    delta = cs * deltabar + sn * alpha;
    molder = mold;
    mold = m;
    m = vv - delta * mold - epsilon * molder;
    Am = Am - delta * Amold - epsilon * Amolder;
    gammabar = sn * deltabar - cs * alpha;
    epsilon = sn * beta;
    deltabar = - cs * beta;
    gamma = sqrt(gammabar^2 + beta^2);
    m = m / gamma;
    Am = Am / gamma;
    cs = gammabar / gamma;
    sn = beta / gamma;
    % Check for stagnation of the method
    if (snprod*cs == 0) || (abs(snprod*cs)*norm(m) < eps*norm(u))
        % increment the number of consecutive iterates which are the same
        stag = stag + 1;
    else
        stag = 0;
    end
    u = u + (snprod * cs) * m;
    snprod = snprod * sn;
    normr = abs(snprod);

    resvec(ii+1,1) = normr;

    % Check for convergence
    if ( ( normr <= tolg ) || ( stag >= maxstagsteps ) || moresteps)
        % double check residual norm is less than tolerance
        r = g - T( u );
        normr_act = norm(r);
        resvec(ii+1, 1) = normr_act;
        if (normr_act <= tolg)
            flag = 0;
            iter = ii;
            break
        else
            if stag >= maxstagsteps && moresteps == 0
                stag = 0;
            end
            moresteps = moresteps + 1;
            if moresteps >= maxmsteps
                if ~warned
                    warning('tooSmallTolerance');
                end
                flag = 3;
                iter = ii;
                break;
            end
        end
    end

    if ( normr < normrmin )        % update minimal norm quantities
        normrmin = normr;
        umin = u;
        imin = ii;
    end

    if ( stag >= maxstagsteps )    % 3 iterates are the same
        flag = 3;
        break
    end
end                                % for ii = 1 : maxit
if ( isempty(ii) )
    ii = 1;
end

% Returned solution is first with minimal residual:
if ( flag == 0 )
    relres = normr_act / n2f;
else
    r_comp = g - T(u);
    if ( norm(r_comp) <= normr_act )
        u = umin;
        iter = imin;
        relres = norm(r_comp) / n2f;
    else
        iter = ii;
        relres = normr_act / n2f;
    end
end

% Undo preconditioner:
u = R1( u ) + z;

% truncate the zeros from resvec
if ( (flag <= 1) || (flag == 3) )
    resvec = resvec(1:ii+1);
else
    resvec = resvec(1:ii);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function x0 = opti_guessX0(prob,ndec)
%Try choose a better x0 than zeros() if we have all finite bounds

if(~isempty(prob.lb) && ~isempty(prob.ub) && all(~isinf(prob.lb)) && all(~isinf(prob.ub)))
    x0 = (prob.ub-prob.lb)./2;
else
    x0 = zeros(ndec,1);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [w, infos] = subsamp_tr(problem, in_options)

% Trust Region Method
% Authors: Jonas Kohler and Aurelien Lucchi, 2017

% Minimize a continous, unconstrained function using the Trust Region method.
%
% References
% ----------
% Conn, A. R., Gould, N. I., & Toint, P. L. (2000). Trust region methods. Society for Industrial and Applied Mathematics.
%
%
% Original Python code was created by J. M. Kohler and A. Lucchi (https://github.com/dalab/subsampled_cubic_regularization)
%
%
% This file is part of SGDLibrary.
%
% Ported to MATLAB code by K.Yoshikawa and H.Kasai on March, 2018.
% Modified by H.Kasai on Apr. 17, 2018

    % set dimensions and samples
    d = problem.dim();
    n = problem.samples();

    %% set local options
    local_options.sampling_scheme = 'exponential'; % 'linear', 'exponential'
    local_options.penalty_increase_multiplier = 2;  % multiply by..
    local_options.penalty_derease_multiplier = 2;  % divide by..
    local_options.initial_tr_radius = 1;
    local_options.successful_treshold = 0.1;
    local_options.very_successful_treshold = 0.9;
    local_options.grad_tol = 1e-9;
    local_options.max_trust_radius = 1e4;

    % sampling
    local_options.Hessian_sampling = 1;
    local_options.gradient_sampling = 0;
    local_options.initial_sample_size_Hessian = 0.025;
    local_options.initial_sample_size_gradient = 0.25;

    % subproblem
    local_options.subproblem_solver = 'GLTR';

    local_options.krylov_tol = 1e-1;
    local_options.exact_tol = 1e-2;

    % merge options
    options = mergeOptions(get_default_options(d), local_options);
    options = mergeOptions(options, in_options);

    if options.verbose > 1
        fprintf('- Subproblem_solver:%s\n', options.subproblem_solver)
        fprintf('- Hessian_sampling:%d\n', options.Hessian_sampling)
        fprintf('- Gradient_sampling:%d\n', options.gradient_sampling)
        fprintf('- Sampling_scheme:%s\n\n', options.sampling_scheme)
    end

    % decide mode
    if options.gradient_sampling == 0 && options.Hessian_sampling == 0
        mode = 'TR';
    else
        mode = 'Subsamp TR';
    end

    % initialize
    iter = 0;
    grad_calc_count = 0;
    w = options.w_init;
    lambda_k = 0;
    successful_flag = false;
    tr_radius = options.initial_tr_radius;  % intial tr radius

    grad = problem.full_grad(w);

    % compute exponential growth constant such that full sample size is reached in n_iterations
    if strcmp(options.sampling_scheme, 'exponential')
        exp_growth_constant = ((1-options.initial_sample_size_Hessian)*n)^(1/options.max_iter);
    end

    % store first infos
    clear infos;
    [infos, f_val, optgap] = store_infos(problem, w, options, [], iter, grad_calc_count, 0);

    % display infos
    if options.verbose > 0
        fprintf('%s:        : Iter = %03d, cost = %.16e, optgap = %.4e\n', mode, iter, f_val, optgap);
    end

    % set start time
    start_time = tic();

    for i = 0 : options.max_iter-1

        %% I: Subsampling
        % a) determine batchsize
        if strcmp(options.sampling_scheme, 'exponential')

            sample_size_Hessian = options.Hessian_sampling*(fix(min([n, n*options.initial_sample_size_Hessian + exp_growth_constant^(i+1)]))+1) + (1-options.Hessian_sampling)*n;
            sample_size_gradient = options.gradient_sampling*(fix(min([n, n*options.initial_sample_size_gradient + exp_growth_constant^(i+1)]))+1) + (1-options.gradient_sampling)*n;

        elseif strcmp(options.sampling_scheme, 'linear')

            sample_size_Hessian = options.Hessian_sampling*fix(min([n, max([n*options.initial_sample_size_Hessian, n/options.max_iter*(i+1)])]))+(1-options.Hessian_sampling)*n;
            sample_size_gradient = options.gradient_sampling*fix(min([n, max([n*options.initial_sample_size_gradient, n/options.max_iter*(i+1)])]))+(1-options.gradient_sampling)*n;

        elseif strcmp(options.sampling_scheme, 'fix') % Added by HK

            sample_size_Hessian = options.Hessian_sampling*fix(min(n, n*options.initial_sample_size_Hessian))+(1-options.Hessian_sampling)*n;
            sample_size_gradient = options.gradient_sampling*fix(min(n, n*options.initial_sample_size_gradient))+(1-options.gradient_sampling)*n;

        else

            sample_size_Hessian = n;
            sample_size_gradient = n;

        end

        % b) draw batches
        if sample_size_Hessian < n
            int_idx_Hessian = randi([1, n], sample_size_Hessian,1); % KY
            bool_idx_Hessian = false(n,1);
            bool_idx_Hessian(int_idx_Hessian) = true;
            sub_hess_indices = find(bool_idx_Hessian);
        else
            sub_hess_indices = 1:n;
        end

        if sample_size_gradient < n
            int_idx_gradient = randi([1, n], sample_size_gradient,1);
            bool_idx_gradient = false(n,1);
            bool_idx_gradient(int_idx_gradient) = true;
            sub_grad_indices = find(bool_idx_gradient);
        else
            sub_grad_indices = 1:n;
        end

        n_samples_per_step = sample_size_Hessian+sample_size_gradient;

        % recompute gradient either because of accepted step or because of re-sampling
        if options.gradient_sampling == 1 || successful_flag == 1
            %grad = gradient_f(w, new_X2, new_Y2, alpha);
            grad = problem.grad(w, sub_grad_indices');
            grad_norm = norm(grad);
            if grad_norm < options.grad_tol
                fprintf('Norm of gradient (%e) reached: grad_tol = %g\n', grad_norm, options.grad_tol);
                break;
            end
        end

        %% II: Step computation
        % b) call subproblem solver
        [s, lambda_k] = tr_subsolver(problem, w, grad, tr_radius, sub_hess_indices, successful_flag, lambda_k, options.subproblem_solver,...
                                           options.exact_tol, options.krylov_tol);

        %sn = norm(s);

        %% III: Regularization Update
        f_prev = problem.cost(w);
        f_curr = problem.cost(w+s);
        function_decrease = f_prev - f_curr; % f(w) - f(w+w)

        %Hs = Hv_f(w, new_X, new_Y, s, alpha);
        Hs = problem.hess_vec(w, s, sub_hess_indices);

        model_decrease = -((grad'*s) + 0.5 * (s'*Hs));

        rho = function_decrease / model_decrease;

        if model_decrease >=0
            if options.verbose > 2
                fprintf('\tNegative model decrease (%e). This should not have happened\n', model_decrease);
            end
        end

        % update w if step s is successful
        if rho >= options.successful_treshold
            w = w + s;
            loss = f_curr;
            successful_flag = true;
            accstr = 'ACC';
        else
            loss = f_prev;
            accstr = 'REJ';
        end

        % Update trust region radius
        tr_radius_ = tr_radius;
        if rho < options.successful_treshold
            tr_radius = tr_radius * (1/options.penalty_increase_multiplier);
            %fprintf('unscuccesful iteration\n');
            successful_flag = false;
            trstr = 'TR-';
        elseif rho > options.very_successful_treshold && ((norm(s) - tr_radius) < 1e-10)
            tr_radius = min([options.penalty_derease_multiplier * tr_radius, options.max_trust_radius]);
            trstr = 'TR+';
        else
            trstr = '   ';
        end

%         % recompute gradient either because of accepted step or because of re-sampling
%         if options.gradient_sampling == 1 || successful_flag == 1
%             %grad = gradient_f(w, new_X2, new_Y2, alpha);
%             grad = problem.grad(w, sub_grad_indices');
%             grad_norm = norm(grad);
%             if grad_norm < options.grad_tol
%                 break;
%             end
%         end

        % measure elapsed time
        elapsed_time = toc(start_time);

        % count gradient evaluations
        grad_calc_count = grad_calc_count + n_samples_per_step;
        iter = iter + 1;

        % store infos
        [infos, f_val, optgap] = store_infos(problem, w, options, infos, iter, grad_calc_count, elapsed_time);

        % display infos
        if options.verbose > 0
            fprintf('%s: %s %s: Iter = %03d, cost = %.16e, optgap = %.4e', mode, accstr, trstr, iter, f_val, optgap);
            if options.verbose > 1
                fprintf(', sample (H,G)= (%d, %d)\n', sample_size_Hessian, sample_size_gradient)
            else
                fprintf('\n');
            end
        end

        if optgap < options.tol_optgap
            fprintf('Optimality gap tolerance reached: tol_optgap = %g\n', options.tol_optgap);
            break;
        end
    end

    if iter == options.max_iter
        fprintf('Max epoch reached: max_iter = %g\n', options.max_iter);
    end

end

author: hiroyuki-kasai

repo: SGDLibrary


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: compact version of matrix-free regularization
%
%  S(y) = alpha/2 * norm(B*(y-yRef))^2,
%
% where
%  alpha regularization parameter, weights regularization versus
%        distance in the joint objective function, alpha = 1 here
%  yRef  is a reference configuration, e.g. yRef = x or a
%        pre-registration result
%  B     a discrete partial differential operator either in explicit
%        matrix form or as a structure containing the necessary
%        parameters to compute B*y
% see also regularizer E8_regularization_MB
%==============================================================================

clear, close all, help(mfilename);

% initialize the regularization and create a starting  point
regularizer('reset','regularizer','mfElastic','alpha',1,'mu',1,'lambda',0);
y0 = @(omega,m) randn(size(getStaggeredGrid(omega,m)));

% 2D example, initialize physical domain and number of discretization points
omega = [0,1,0,1]; m = [16,12];   %

% test derivative of 2D implementation
fctn = @(yc) regularizer(yc,omega,m);  checkDerivative(fctn,y0(omega,m));

% 3D example, initialize physical domain and number of discretization points
omega = [0,1,0,1,0,1]; m  = [16,12,8];

% test derivative of 3D implementation
fctn = @(yc) regularizer(yc,omega,m);
checkDerivative(fctn,y0(omega,m));

%==============================================================================


text:
function [costs, grad, charInfo] = lstmCostGrad(model, trainData, params, isTest)
%%%
%
% Compute cost/grad for LSTM.
% When params.predictPos>0, returns costs.pos and costs.word
% If isTest==1, this method only computes cost (for testing purposes).
%
% Thang Luong @ 2014, 2015, <lmthang@stanford.edu>
%
%%%

  %%%%%%%%%%%%
  %%% INIT %%%
  %%%%%%%%%%%%
  curBatchSize = size(trainData.tgtInput, 1);
  if params.isBi
    srcMaxLen = trainData.srcMaxLen;
  else
    srcMaxLen = 1;
  end

  params.curBatchSize = curBatchSize;
  params.srcMaxLen = srcMaxLen;

  [params] = setAttnParams(params);
  [grad, params] = initGrad(model, params);
  zeroBatch = zeroMatrix([params.lstmSize, params.curBatchSize], params.isGPU, params.dataType);

  % init states
  [zeroState] = createZeroState(params);

  % init costs
  costs = initCosts(params);

  tgtOutput = trainData.tgtOutput;
  tgtInput = trainData.tgtInput;
  % char
  if params.charOpt
    % tgt
    if params.charTgtGen
      origTgtOutput = trainData.tgtOutput;
    end

    tgtOutput(tgtOutput > params.tgtCharShortList) = params.tgtUnk;
    tgtInput(tgtInput > params.tgtCharShortList) = params.tgtUnk;
  end

  %%%%%%%%%%%%%%%%%%%%
  %%% FORWARD PASS %%%
  %%%%%%%%%%%%%%%%%%%%
  %% encoder
  lastEncState = zeroState;
  if params.isBi
    [encStates, lastEncState, encRnnFlags, trainData, srcCharData] = encoderLayerForward(model, zeroState, trainData, params, isTest);

    charInfo.trainNumSrcUnkTypes = 0;
    charInfo.trainNumSrcUnkTokens = 0;
    if params.charSrcRep
      charInfo.trainNumSrcUnkTypes = srcCharData.numRareWords;
      charInfo.trainNumSrcUnkTokens = sum(srcCharData.rareFlags(:));
    end
  end

  %% decoder
  decRnnFlags = struct('decode', 1, 'test', isTest, 'attn', params.attnFunc, 'feedInput', params.feedInput, 'charSrcRep', params.charSrcRep, ...
      'charTgtGen', params.charTgtGen, 'initEmb', []);
  [decStates, trainData, attnInfos] = rnnLayerForward(model.W_tgt, model.W_emb_tgt, lastEncState, tgtInput, trainData.tgtMask, ...
    params, decRnnFlags, trainData, model, []);

  %% softmax
  [costs.word, grad.W_soft, dec_top_grads, costs.indLosses] = softmaxCostGrad(decStates, model.W_soft, tgtOutput, trainData.tgtMask, ...
    params, isTest);
  costs.total = costs.word;

  %% tgt char foward / backprop %%
  charInfo.numChars = 0;
  charInfo.trainNumTgtUnkTokens = 0;
  if params.charTgtGen
    [costs.char, tgtCharGrad, charInfo.numChars] = tgtCharCostGrad(decStates, attnInfos, model, origTgtOutput, trainData.tgtMask, params.tgtCharMap, params, isTest);
    costs.total = costs.total + costs.char;
    charInfo.trainNumTgtUnkTokens = tgtCharGrad.numRareWords;
    if isTest==0
      if tgtCharGrad.numRareWords>0
        grad.W_soft_char = tgtCharGrad.W_soft;
        grad.W_tgt_char = tgtCharGrad.W_tgt;
        grad.W_emb_tgt_char = tgtCharGrad.W_emb_tgt_char;
        grad.indices_tgt_char = tgtCharGrad.indices_tgt_char;

        if params.charFeedOpt
          assert(params.attnFunc>0);
          grad.W_h_char = tgtCharGrad.W_h_char;
        else
          % add top grads from tgt char
          rareCount = 0;
          decLen = length(dec_top_grads);
          for tt=1:decLen
            rareIndices = find(tgtCharGrad.rareFlags(:, tt));
            dec_top_grads{tt}(:, rareIndices) = dec_top_grads{tt}(:, rareIndices) + tgtCharGrad.initEmb(:, rareCount+1:rareCount+length(rareIndices));
            rareCount = rareCount + length(rareIndices);
          end
          assert(rareCount == tgtCharGrad.numRareWords);
        end
      else
        grad.indices_tgt_char = [];
      end

      % clear tgtCharGrad;
      if params.debug
        fprintf(2, '# after clearing tgtCharGrad, %s\n', gpuInfo(params.gpu));
      end
    end
  else
    tgtCharGrad = [];
  end

  if isTest==1 % don't compute grad
    return;
  end

  %%%%%%%%%%%%%%%%%%%%%
  %%% BACKWARD PASS %%%
  %%%%%%%%%%%%%%%%%%%%%
  % h_t and c_t gradients accumulate over time per layer
  zeroGrad = cell(params.numLayers, 1);
  for ll=1:params.numLayers
    zeroGrad{ll} = zeroBatch;
  end

  %% decoder
  [dc, dh, grad.W_tgt, grad.W_emb_tgt, grad.indices_tgt, attnGrad, grad.srcHidVecs] = rnnLayerBackprop(model.W_tgt, ...
    decStates, lastEncState, dec_top_grads, zeroGrad, zeroGrad, tgtInput, trainData.tgtMask, params, decRnnFlags, ...
    attnInfos, trainData, model, tgtCharGrad);
  if params.attnFunc % copy attention grads
    [grad] = copyStruct(attnGrad, grad);
  end

  %% encoder
  if params.isBi
    enc_top_grads = cell(srcMaxLen - 1, 1);
    for tt=1:params.numSrcHidVecs % attention
      enc_top_grads{tt} = grad.srcHidVecs(:,:,tt);
    end

    [~, ~, grad.W_src, grad.W_emb_src, grad.indices_src, ~, ~] = rnnLayerBackprop(model.W_src, encStates, zeroState, ...
    enc_top_grads, dc, dh, trainData.srcInput, trainData.srcMask, params, encRnnFlags, attnInfos, trainData, model, []);

    % char backprop
    if params.charSrcRep
      if srcCharData.numRareWords > 0
        if isTest==0 && params.charSrcSample > 0 % also include frequent words
          rareFlags = ismember(grad.indices_src, srcCharData.rareWords);
        else
          rareFlags = grad.indices_src > params.srcCharShortList;
        end

        % rare
        srcCharGrad.embs = grad.W_emb_src(:, rareFlags);
        srcCharGrad.indices = grad.indices_src(rareFlags);

        % frequent
        grad.W_emb_src = grad.W_emb_src(:, ~rareFlags);
        grad.indices_src = grad.indices_src(~rareFlags);

        [grad.W_src_char, grad.W_emb_src_char, grad.indices_src_char] = srcCharLayerBackprop(model.W_src_char, srcCharData, srcCharGrad);
      else
        grad.indices_src_char = [];
      end
    end
  end

  % remove unused variables
  if params.attnFunc>0
    grad = rmfield(grad, 'srcHidVecs');
  end
end

function [grad, params] = initGrad(model, params)
  %% grad
  for ii=1:length(params.varsDenseUpdate)
    field = params.varsDenseUpdate{ii};
    if iscell(model.(field))
      for jj=1:length(model.(field)) % cell, like W_src, W_tgt
        grad.(field){jj} = zeroMatrix(size(model.(field){jj}), params.isGPU, params.dataType);
      end
    else
      grad.(field) = zeroMatrix(size(model.(field)), params.isGPU, params.dataType);
    end
  end

  %% backprop to src hidden states for attention and positional models
  if params.attnFunc>0
    % we extract trainData.srcHidVecs later, which contains all src hidden states, lstmSize * curBatchSize * numSrcHidVecs
    grad.srcHidVecs = zeroMatrix([params.lstmSize, params.curBatchSize, params.numSrcHidVecs], params.isGPU, params.dataType);
  end
end