text:
%Simulation unit test cases
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Author: Jared Jacobs, jmjacobs@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 9/14/2010
classdef Simulation_Test < TestCase
    properties
        simulation
        outputDir
    end

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

        function setUp(this)
            this.simulation = this.newSimulation();
            this.outputDir = '';
        end

        function tearDown(this)
            if ~isempty(this.outputDir) && exist(this.outputDir, 'dir')
                rmdir(this.outputDir, 's');
            end
        end
    end

    %tests
    methods
        function testConstants(this)
            sim = this.simulation;

            assertAllEqual(false, isnan(sim.state('Metabolite').molecularWeights));
            assertEqual(1, sim.compartment.cytosolIndexs);
        end

        function testAllocateMemoryForState(this)
            this.simulation.allocateMemoryForState(1);  %no run-time errors
        end

        function testApplyPerturbations(this)
            sim = this.simulation;
            g =  sim.gene;

            % 1 mRNA genes
            this.helper_testApplyPerturbations(g.mRNAIndexs(1));

            % 2 mRNA genes
            this.helper_testApplyPerturbations(g.mRNAIndexs([1 101]));

            % mRNA genes
            this.helper_testApplyPerturbations(g.mRNAIndexs);

            % rRNA genes
            this.helper_testApplyPerturbations(g.rRNAIndexs);

            % sRNA genes
            this.helper_testApplyPerturbations(g.sRNAIndexs);

            % tRNA genes
            this.helper_testApplyPerturbations(g.tRNAIndexs);

            % all genes
            this.helper_testApplyPerturbations((1:numel(g.wholeCellModelIDs))');
        end

        function helper_testApplyPerturbations(this, geneIdxs)
            %% references
            sim = this.simulation;
            g =  sim.gene;
            r = sim.state('Rna');
            pm = sim.state('ProteinMonomer');
            pc = sim.state('ProteinComplex');
            ftsZRing = sim.state('FtsZRing');
            c = sim.state('Chromosome');
            rnaPol = sim.state('RNAPolymerase');
            rib = sim.state('Ribosome');
            pol = sim.state('Polypeptide');
            transcript = sim.state('Transcript');

            %% initialize state
            r.decayRates(:) = 0;
            pm.decayRates(:) = 0;
            pc.decayRates(:) = 0;
            r.counts(:) = 1;
            pm.counts(:) = 1;
            pc.counts(:) = 1;

            r.counts(r.boundIndexs(r.matureSRNAIndexs), :) = 0;

            c.initialize();
            vals = setdiff((1:numel(pm.boundIndexs))', pm.translationFactorIndexs);
            monPosStrnds = [100 * (1:numel(vals))'  ones(size(vals))];
            c.monomerBoundSites(monPosStrnds) = vals;
            vals = setdiff((1:numel(pc.boundIndexs))', [
                pc.ribosome70SIndexs
                pc.translationFactorIndexs
                pc.ftsZGTPIndexs
                pc.ftsZGDPIndexs
                pc.dnaAPolymerIndexs
                pc.replisomeIndexs
                ]);
            cpxPosStrnds = [1000 * (1:numel(vals))' + monPosStrnds(end, 1)  ones(size(vals))];
            c.complexBoundSites(cpxPosStrnds) = vals;
            assertIn(cpxPosStrnds(end, 1), [1 c.sequenceLen-1000]);
            pc.counts(pc.boundIndexs([pc.dnaAPolymerIndexs; pc.replisomeIndexs]), :) = 0;

            rnaPol.states = [
                rnaPol.nonSpecificallyBoundValue
                rnaPol.specificallyBoundValue
                ];
            rnaPol.positionStrands = [
                find(c.complexBoundSites == pc.rnaPolymeraseIndexs(1))
                find(c.complexBoundSites == pc.rnaPolymeraseIndexs(2))
                ];
            transcript.boundTranscriptionUnits = [0; 1];
            transcript.boundTranscriptProgress = [0; 0];
            transcript.boundTranscriptChromosome = [0; 0];

            rib.states = rib.activeValue;
            rib.boundMRNAs = 1;
            rib.mRNAPositions = 1;
            rib.tmRNAPositions = 0;
            pol.boundMRNAs = rib.boundMRNAs;
            pol.nascentMonomerLengths = rib.mRNAPositions;
            pol.proteolysisTagLengths = rib.tmRNAPositions;

            pc.counts(pc.boundIndexs([pc.ftsZGTPIndexs; pc.ftsZGDPIndexs]), :) = 0;
            ftsZRing.numEdgesOneStraight = 0;
            ftsZRing.numEdgesTwoStraight = 0;
            ftsZRing.numEdgesTwoBent = 0;
            ftsZRing.numResidualBent = 0;

            %% apply perturbations
            monIdxs = find(ismember(g.mRNAIndexs, geneIdxs));
            matRnaIdxs = find(any(r.matureRNAGeneComposition(setdiff(geneIdxs, g.mRNAIndexs), :), 1));
            cpxIdxs = find(any(any(pc.proteinComplexComposition(geneIdxs, :, :), 3), 1));

            otherMatRnaIdxs = setdiff(1:numel(r.matureIndexs), matRnaIdxs);
            otherMonIdxs = setdiff(1:numel(pm.matureIndexs), monIdxs);
            otherCpxIdxs = setdiff(1:numel(pc.matureIndexs), cpxIdxs);

            rnaDecayRates = r.decayRates;
            monDecayRates = pm.decayRates;
            cpxDecayRates = pc.decayRates;

            rnaCounts = r.counts;
            monCounts = pm.counts;
            cpxCounts = pc.counts;

            sim.applyOptions('geneticKnockouts', g.wholeCellModelIDs(geneIdxs));
            sim.applyPerturbations();

            %% assertions
            %decay rates
            decayRate = log(2)/(log(2)/1e6);

            assertAllEqual(0, r.decayRates(r.nascentIndexs));
            assertAllEqual(0, r.decayRates(r.intergenicIndexs));
            if ~isempty(matRnaIdxs)
                assertAllEqual(decayRate, r.decayRates(r.processedIndexs(matRnaIdxs)));
                assertAllEqual(decayRate, r.decayRates(r.matureIndexs(matRnaIdxs)));
                assertAllEqual(decayRate, r.decayRates(r.boundIndexs(matRnaIdxs)));
                assertAllEqual(decayRate, r.decayRates(r.misfoldedIndexs(matRnaIdxs)));
                assertAllEqual(decayRate, r.decayRates(r.damagedIndexs(matRnaIdxs)));
                assertAllEqual(decayRate, r.decayRates(r.aminoacylatedIndexs(matRnaIdxs)));
            end
            assertEqual(rnaDecayRates(r.processedIndexs(otherMatRnaIdxs)), r.decayRates(r.processedIndexs(otherMatRnaIdxs)));
            assertEqual(rnaDecayRates(r.matureIndexs(otherMatRnaIdxs)), r.decayRates(r.matureIndexs(otherMatRnaIdxs)));
            assertEqual(rnaDecayRates(r.boundIndexs(otherMatRnaIdxs)), r.decayRates(r.boundIndexs(otherMatRnaIdxs)));
            assertEqual(rnaDecayRates(r.misfoldedIndexs(otherMatRnaIdxs)), r.decayRates(r.misfoldedIndexs(otherMatRnaIdxs)));
            assertEqual(rnaDecayRates(r.damagedIndexs(otherMatRnaIdxs)), r.decayRates(r.damagedIndexs(otherMatRnaIdxs)));
            assertEqual(rnaDecayRates(r.aminoacylatedIndexs(otherMatRnaIdxs)), r.decayRates(r.aminoacylatedIndexs(otherMatRnaIdxs)));

            if ~isempty(monIdxs)
                assertAllEqual(decayRate, pm.decayRates(pm.nascentIndexs(monIdxs)));
                assertAllEqual(decayRate, pm.decayRates(pm.signalSequenceIndexs(monIdxs)));
                assertAllEqual(decayRate, pm.decayRates(pm.processedIIndexs(monIdxs)));
                assertAllEqual(decayRate, pm.decayRates(pm.processedIIIndexs(monIdxs)));
                assertAllEqual(decayRate, pm.decayRates(pm.foldedIndexs(monIdxs)));
                assertAllEqual(decayRate, pm.decayRates(pm.matureIndexs(monIdxs)));
                assertAllEqual(decayRate, pm.decayRates(pm.inactivatedIndexs(monIdxs)));
                assertAllEqual(decayRate, pm.decayRates(pm.boundIndexs(monIdxs)));
                assertAllEqual(decayRate, pm.decayRates(pm.misfoldedIndexs(monIdxs)));
                assertAllEqual(decayRate, pm.decayRates(pm.damagedIndexs(monIdxs)));
            end
            assertEqual(monDecayRates(pm.nascentIndexs(otherMonIdxs)), pm.decayRates(pm.nascentIndexs(otherMonIdxs)));
            assertEqual(monDecayRates(pm.signalSequenceIndexs(otherMonIdxs)), pm.decayRates(pm.signalSequenceIndexs(otherMonIdxs)));
            assertEqual(monDecayRates(pm.processedIIndexs(otherMonIdxs)), pm.decayRates(pm.processedIIndexs(otherMonIdxs)));
            assertEqual(monDecayRates(pm.processedIIIndexs(otherMonIdxs)), pm.decayRates(pm.processedIIIndexs(otherMonIdxs)));
            assertEqual(monDecayRates(pm.foldedIndexs(otherMonIdxs)), pm.decayRates(pm.foldedIndexs(otherMonIdxs)));
            assertEqual(monDecayRates(pm.matureIndexs(otherMonIdxs)), pm.decayRates(pm.matureIndexs(otherMonIdxs)));
            assertEqual(monDecayRates(pm.inactivatedIndexs(otherMonIdxs)), pm.decayRates(pm.inactivatedIndexs(otherMonIdxs)));
            assertEqual(monDecayRates(pm.boundIndexs(otherMonIdxs)), pm.decayRates(pm.boundIndexs(otherMonIdxs)));
            assertEqual(monDecayRates(pm.misfoldedIndexs(otherMonIdxs)), pm.decayRates(pm.misfoldedIndexs(otherMonIdxs)));
            assertEqual(monDecayRates(pm.damagedIndexs(otherMonIdxs)), pm.decayRates(pm.damagedIndexs(otherMonIdxs)));

            if ~isempty(cpxIdxs)
                assertAllEqual(decayRate, pc.decayRates(pc.nascentIndexs(cpxIdxs)));
                assertAllEqual(decayRate, pc.decayRates(pc.matureIndexs(cpxIdxs)));
                assertAllEqual(decayRate, pc.decayRates(pc.inactivatedIndexs(cpxIdxs)));
                assertAllEqual(decayRate, pc.decayRates(pc.boundIndexs(cpxIdxs)));
                assertAllEqual(decayRate, pc.decayRates(pc.misfoldedIndexs(cpxIdxs)));
                assertAllEqual(decayRate, pc.decayRates(pc.damagedIndexs(cpxIdxs)));
            end
            assertEqual(cpxDecayRates(pc.nascentIndexs(otherCpxIdxs)), pc.decayRates(pc.nascentIndexs(otherCpxIdxs)));
            assertEqual(cpxDecayRates(pc.matureIndexs(otherCpxIdxs)), pc.decayRates(pc.matureIndexs(otherCpxIdxs)));
            assertEqual(cpxDecayRates(pc.inactivatedIndexs(otherCpxIdxs)), pc.decayRates(pc.inactivatedIndexs(otherCpxIdxs)));
            assertEqual(cpxDecayRates(pc.boundIndexs(otherCpxIdxs)), pc.decayRates(pc.boundIndexs(otherCpxIdxs)));
            assertEqual(cpxDecayRates(pc.misfoldedIndexs(otherCpxIdxs)), pc.decayRates(pc.misfoldedIndexs(otherCpxIdxs)));
            assertEqual(cpxDecayRates(pc.damagedIndexs(otherCpxIdxs)), pc.decayRates(pc.damagedIndexs(otherCpxIdxs)));

            %counts
            assertAllEqual(1, r.counts(r.nascentIndexs, :));
            assertAllEqual(1, r.counts(r.intergenicIndexs, :));
            if ~isempty(matRnaIdxs)
                assertAllEqual(0, r.counts(r.processedIndexs(matRnaIdxs), :));
                assertAllEqual(0, r.counts(r.matureIndexs(matRnaIdxs), :));
                assertAllEqual(0, r.counts(r.boundIndexs(matRnaIdxs), :));
                assertAllEqual(0, r.counts(r.misfoldedIndexs(matRnaIdxs), :));
                assertAllEqual(0, r.counts(r.damagedIndexs(matRnaIdxs), :));
                assertAllEqual(0, r.counts(r.aminoacylatedIndexs(matRnaIdxs), :));
            end
            assertEqual(rnaCounts(r.processedIndexs(otherMatRnaIdxs, :)), r.counts(r.processedIndexs(otherMatRnaIdxs, :)));
            assertEqual(rnaCounts(r.matureIndexs(otherMatRnaIdxs, :)), r.counts(r.matureIndexs(otherMatRnaIdxs, :)));
            assertEqual(rnaCounts(r.boundIndexs(otherMatRnaIdxs, :)), r.counts(r.boundIndexs(otherMatRnaIdxs, :)));
            assertEqual(rnaCounts(r.misfoldedIndexs(otherMatRnaIdxs, :)), r.counts(r.misfoldedIndexs(otherMatRnaIdxs, :)));
            assertEqual(rnaCounts(r.damagedIndexs(otherMatRnaIdxs, :)), r.counts(r.damagedIndexs(otherMatRnaIdxs, :)));
            assertEqual(rnaCounts(r.aminoacylatedIndexs(otherMatRnaIdxs, :)), r.counts(r.aminoacylatedIndexs(otherMatRnaIdxs, :)));

            if ~isempty(monIdxs)
                assertAllEqual(0, pm.counts(pm.nascentIndexs(monIdxs), :));
                assertAllEqual(0, pm.counts(pm.signalSequenceIndexs(monIdxs), :));
                assertAllEqual(0, pm.counts(pm.processedIIndexs(monIdxs), :));
                assertAllEqual(0, pm.counts(pm.processedIIIndexs(monIdxs), :));
                assertAllEqual(0, pm.counts(pm.foldedIndexs(monIdxs), :));
                assertAllEqual(0, pm.counts(pm.matureIndexs(monIdxs), :));
                assertAllEqual(0, pm.counts(pm.inactivatedIndexs(monIdxs), :));
                assertAllEqual(0, pm.counts(pm.boundIndexs(monIdxs), :));
                assertAllEqual(0, pm.counts(pm.misfoldedIndexs(monIdxs), :));
                assertAllEqual(0, pm.counts(pm.damagedIndexs(monIdxs), :));
            end
            assertEqual(monCounts(pm.nascentIndexs(otherMonIdxs, :)), pm.counts(pm.nascentIndexs(otherMonIdxs, :)));
            assertEqual(monCounts(pm.signalSequenceIndexs(otherMonIdxs, :)), pm.counts(pm.signalSequenceIndexs(otherMonIdxs, :)));
            assertEqual(monCounts(pm.processedIIndexs(otherMonIdxs, :)), pm.counts(pm.processedIIndexs(otherMonIdxs, :)));
            assertEqual(monCounts(pm.processedIIIndexs(otherMonIdxs, :)), pm.counts(pm.processedIIIndexs(otherMonIdxs, :)));
            assertEqual(monCounts(pm.foldedIndexs(otherMonIdxs, :)), pm.counts(pm.foldedIndexs(otherMonIdxs, :)));
            assertEqual(monCounts(pm.matureIndexs(otherMonIdxs, :)), pm.counts(pm.matureIndexs(otherMonIdxs, :)));
            assertEqual(monCounts(pm.inactivatedIndexs(otherMonIdxs, :)), pm.counts(pm.inactivatedIndexs(otherMonIdxs, :)));
            assertEqual(monCounts(pm.boundIndexs(otherMonIdxs, :)), pm.counts(pm.boundIndexs(otherMonIdxs, :)));
            assertEqual(monCounts(pm.misfoldedIndexs(otherMonIdxs, :)), pm.counts(pm.misfoldedIndexs(otherMonIdxs, :)));
            assertEqual(monCounts(pm.damagedIndexs(otherMonIdxs, :)), pm.counts(pm.damagedIndexs(otherMonIdxs, :)));

            if ~isempty(cpxIdxs)
                assertAllEqual(0, pc.counts(pc.nascentIndexs(cpxIdxs), :));
                assertAllEqual(0, pc.counts(pc.matureIndexs(cpxIdxs), :));
                assertAllEqual(0, pc.counts(pc.inactivatedIndexs(cpxIdxs), :));
                assertAllEqual(0, pc.counts(pc.boundIndexs(cpxIdxs), :));
                assertAllEqual(0, pc.counts(pc.misfoldedIndexs(cpxIdxs), :));
                assertAllEqual(0, pc.counts(pc.damagedIndexs(cpxIdxs), :));
            end
            assertEqual(cpxCounts(pc.nascentIndexs(otherCpxIdxs, :)), pc.counts(pc.nascentIndexs(otherCpxIdxs, :)));
            assertEqual(cpxCounts(pc.matureIndexs(otherCpxIdxs, :)), pc.counts(pc.matureIndexs(otherCpxIdxs, :)));
            assertEqual(cpxCounts(pc.inactivatedIndexs(otherCpxIdxs, :)), pc.counts(pc.inactivatedIndexs(otherCpxIdxs, :)));
            assertEqual(cpxCounts(pc.boundIndexs(otherCpxIdxs, :)), pc.counts(pc.boundIndexs(otherCpxIdxs, :)));
            assertEqual(cpxCounts(pc.misfoldedIndexs(otherCpxIdxs, :)), pc.counts(pc.misfoldedIndexs(otherCpxIdxs, :)));
            assertEqual(cpxCounts(pc.damagedIndexs(otherCpxIdxs, :)), pc.counts(pc.damagedIndexs(otherCpxIdxs, :)));

            %% test membrane protein decay
            %monomer
            sim = edu.stanford.covert.cell.sim.SimulationFixture.load([], true);
            comp = sim.compartment;
            g = sim.gene;
            pm = sim.state('ProteinMonomer');
            protDcy = sim.process('ProteinDecay');

            idx = find(pm.compartments(pm.matureIndexs) == comp.membraneIndexs, 1, 'first');

            sim.applyOptions('geneticKnockouts', g.wholeCellModelIDs(g.mRNAIndexs(idx)));
            sim.applyPerturbations();

            monIdx = pm.matureIndexs(idx);
            pm.counts(monIdx, :) = 0;
            pm.counts(monIdx, comp.membraneIndexs) = 1;

            protDcy.copyFromState();
            protDcy.enzymes(:) = 1e6;
            protDcy.substrates(:) = 1e6;
            protDcy.evolveState();
            protDcy.copyToState();

            assertEqual(0, pm.counts(monIdx, comp.membraneIndexs));

            %complex
            sim = edu.stanford.covert.cell.sim.SimulationFixture.load([], true);
            comp = sim.compartment;
            g = sim.gene;
            pc = sim.state('ProteinComplex');
            protDcy = sim.process('ProteinDecay');

            idx = find(pc.compartments(pc.matureIndexs) == comp.membraneIndexs, 1, 'first');
            geneIdx = find(any(pc.proteinComplexComposition(:, idx, :), 3), 1, 'first');

            sim.applyOptions('geneticKnockouts', g.wholeCellModelIDs(geneIdx));
            sim.applyPerturbations();

            cpxIdx = pc.matureIndexs(idx);
            pc.counts(cpxIdx, :) = 0;
            pc.counts(cpxIdx, comp.membraneIndexs) = 1;

            protDcy.copyFromState();
            protDcy.enzymes(:) = 1e6;
            protDcy.substrates(:) = 1e6;
            protDcy.evolveState();
            protDcy.copyToState();

            assertEqual(0, pc.counts(cpxIdx, comp.membraneIndexs));
        end

        function testEvolveState(this)
            sim = this.simulation;
            met = sim.state('Metabolite');

            t0 = sim.state('Time').values;

            %mechanism for recording a call sequence
            calls = {};
            function r = recorder(s, value)
                if ~exist('value', 'var')
                    value = [];
                end
                function result = onCall(varargin)
                    calls{end+1} = s;
                    if isa(value, 'function_handle')
                        value(varargin{:}); %#ok<NOEFF>
                        result = [];
                    else
                        result = value;
                    end
                end
                r = @onCall;
            end

            %mechanism for using substrates
            function u = user(substrates)
                function use(process)
                    process.substrates = process.substrates - substrates;
                end
                u = @use;
            end

            %metabolite indices of interest (no special significance)
            M1 = 7; M2 = 123; M3 = 432;

            %three fake processes
            sim.setForTest('processes', {
                %simple case: all metabolites in same compartment, all
                %substrates are metabolites
                edu.stanford.covert.test.mock(...
                    'substrates', [0 0]',...
                    'substrateMetaboliteLocalIndexs', [1 2]',...
                    'substrateMetaboliteGlobalIndexs', [M1 M3]',...
                    'substrateMetaboliteCompartmentIndexs', [1 1]',...
                    'substrateMetaboliteGlobalCompartmentIndexs', sub2ind(size(met.counts), [M1 M3]', [1 1]'), ...
                    'calcResourceRequirements_Current', recorder('r1', [10 10]'),...
                    'copyFromState', recorder('f1'),...
                    'copyToState', recorder('t1'),...
                    'evolveState', recorder('e1', user([3 4]')))
                %more complex: metabolites in various compartments, some
                %substrates are not metabolites
                edu.stanford.covert.test.mock(...
                    'substrates', [0 0 0 0 0]',...
                    'substrateMetaboliteLocalIndexs', [1 2 5]',...
                    'substrateMetaboliteGlobalIndexs', [M1 M2 M3]',...
                    'substrateMetaboliteCompartmentIndexs', [3 1 1]',...
                    'substrateMetaboliteGlobalCompartmentIndexs', sub2ind(size(met.counts), [M1 M2 M3]', [3 1 1]'), ...
                    'calcResourceRequirements_Current', recorder('r2', [5 5 0 0 5]'),...
                    'copyFromState', recorder('f2'),...
                    'copyToState', recorder('t2'),...
                    'evolveState', recorder('e2', user([1 0 0 0 -1]')))
                %each metabolite in multiple compartments
                edu.stanford.covert.test.mock(...
                    'substrates', [0 0 0 0; 0 0 0 0]',...
                    'substrateMetaboliteLocalIndexs', [1 3 4]',...
                    'substrateMetaboliteGlobalIndexs', [M2 M3 M1; M2 M3 M1]',...
                    'substrateMetaboliteCompartmentIndexs', [5 4 5; 3 1 3]',...
                    'substrateMetaboliteGlobalCompartmentIndexs', sub2ind(size(met.counts), [M2 M3 M1; M2 M3 M1]', [5 4 5; 3 1 3]'), ...
                    'calcResourceRequirements_Current', recorder('r3', [5 0 5 0; 0 0 5 5]'),...
                    'copyFromState', recorder('f3'),...
                    'copyToState', recorder('t3'),...
                    'evolveState', recorder('e3'))
                });
            sim.setForTest('processEvalOrderIndexs', [3 2 1]');
            sim.setForTest('processesInEvalOrder', sim.processes(sim.processEvalOrderIndexs));
            met.counts(M1,:) = 20;
            met.counts(M2,:) = 0;
            met.counts(M3,:) = 10;

            sim.evolveState();
            assertEqual(t0 + sim.stepSizeSec, sim.state('Time').values);
            assertEqual(...
                'f1 r1 f2 r2 f3 r3 f3 e3 t3 f2 e2 t2 f1 e1 t1',...
                strjoin(' ', calls{:}));

            %processes receive proporational to their allocation
            assertEqual(...
                [17 1],...
                sim.processes{1}.substrates');
            assertEqual(...
                [9 0 0 0 3],...
                sim.processes{2}.substrates');
            assertEqual(...
                [0 0 10 0;
                 0 0 2 10],...
                sim.processes{3}.substrates');

            %verify that metabolite/substrate usage/production propagate back
            assertEqual([17 20 19 20 20 20], met.counts(M1,:));
            assertEqual([ 0  0  0  0  0  0], met.counts(M2,:));
            assertEqual([ 7 10 10 10 10 10], met.counts(M3,:));
        end

        function testRunLogging(this)
            sim = this.simulation;
            sim.applyOptions('lengthSec', 10);
            sim.state('MetabolicReaction').initialGrowthFilterWidth = Inf;

            %1 logger
            logger = edu.stanford.covert.cell.sim.util.SummaryLogger(2, 0);
            sim.run(logger);

            %2 loggers
            logger1 = edu.stanford.covert.cell.sim.util.SummaryLogger(2, 0);
            logger2 = edu.stanford.covert.cell.sim.util.SummaryLogger(2, 0);
            sim.run({logger1 logger2});
        end

        function testMassBalance(this)
            sim = this.simulation;
            sim.applyOptions('lengthSec', 100);
            sim.state('MetabolicReaction').initialGrowthFilterWidth = Inf;

            sim.allocateMemoryForState(1);
            sim.initializeState();
            sim.applyPerturbations();

            mass = sim.state('Mass');
            initialMass = sum(mass.cell);

            sim.run();

            assertElementsAlmostEqual(initialMass, sum(mass.cell));
        end

        function testPlots(this)
            sim = this.simulation;

            %create figure, add axis
            handle = figure();
            axesHandle = subplot(1,1,1,'Parent',handle);
            legend(axesHandle,'off');
            set(axesHandle,'Box','on');

            %inputs
            time = 1;
            compartments = 1;

            %call all plotting functions
            simulationMetaData = metaclass(sim);
            failedPlotMethodNames = {};
            for i = 1:length(simulationMetaData.Methods)
                if length(simulationMetaData.Methods{i}.Name) > 4 && ...
                   strcmp(simulationMetaData.Methods{i}.Access, 'public') && ...
                   strcmp(simulationMetaData.Methods{i}.Name(1:4), 'plot') && ...
                   numel(simulationMetaData.Methods{i}.InputNames) == 4
                    cla(axesHandle);
                    try
                        sim.(simulationMetaData.Methods{i}.Name)(axesHandle, time, compartments)
                    catch %#ok<CTCH>
                        failedPlotMethodNames{end+1}=simulationMetaData.Methods{i}.Name; %#ok<AGROW>
                    end
                end
            end

            %close figure
            close(handle);

            assert(isempty(failedPlotMethodNames), ...
                ['The following plot methods failed: ' strjoin(', ',failedPlotMethodNames{:})]);
        end

        function testOptions(this)
            sim = this.simulation;
            assertEqual(0, sim.getOptions.verbosity);
            assertEqual(1, sim.getOptions.stepSizeSec);

            sim.applyOptions();
            assertEqual(0, sim.getOptions.verbosity);
            assertEqual(1, sim.getOptions.stepSizeSec);

            sim.applyOptions('verbosity', 3, 'stepSizeSec', 0.5);
            assertEqual(3, sim.getOptions.verbosity);
            assertEqual(0.5, sim.getOptions.stepSizeSec);

            sim.applyOptions('verbosity', 4);
            assertEqual(4, sim.getOptions.verbosity);
            assertEqual(0.5, sim.getOptions.stepSizeSec);
        end
    end

    %helpers
    methods (Access = private)
        function sim = newSimulation(~)
            sim = edu.stanford.covert.cell.sim.SimulationFixture.load();
        end

        function prepareOutputDir(this)
            this.outputDir = ['tmp/Simulation_Integrated_Test_' this.MethodName];
            if ~exist(this.outputDir, 'dir')
                mkdir(this.outputDir);
            end
            delete([this.outputDir '/*']);
        end

        function result = createMetadata(this)
            result = struct;
            result.shortDescription = 'Wild-type simulation';
            result.longDescription  = 'Wild-type simulation';
            result.email            = 'jkarr@stanford.edu';
            result.firstName        = 'Jonathan';
            result.lastName         = 'Karr';
            result.affiliation      = 'Covert Lab, Department of Bioengineering, Stanford University';
            result.knowledgeBaseWID = this.simulation.knowledgeBaseWID;
            result.downsampleStep   = 5000;
            result.outputDirectory  = this.outputDir;
            result.startTime        = datestr(clock,'yyyy-mm-dd HH:MM:SS');
            result.endTime          = datestr(clock,'yyyy-mm-dd HH:MM:SS');
            [result.userName, result.hostName, result.ipAddress] = ...
                edu.stanford.covert.util.computerInfo();
            result.revision = 234;
            result.differencesFromRevision = '';
        end
    end

    %static helpers
    methods (Static, Access = private)
        function assertPropertyValuesEqual(...
                metadata, options, parameters, fittedConstants, timeCourses, sim)
            %meta data
            metadataNames = fieldnames(sim.metadata);
            for i = 1:length(metadataNames)
                name = metadataNames{i};
                assertEqual(metadata.(name), sim.metadata.(name), name);
            end

            %options
            optionNames = sim.optionNames;
            for i = 1:length(optionNames)
                name = optionNames{i};
                switch class(options.Core__.(name))
                    case {'cell','struct'}
                        assertEqual(options.Core__.(name), sim.(name));
                    otherwise
                        assertElementsAlmostEqual(...
                            options.Core__.(name), sim.(name),...
                            'relative', 1e-4, name);
                end

            end
            for i = 1:length(sim.processes)
                process = sim.processes{i};
                optionNames = process.optionNames;
                for j = 1:length(optionNames)
                    name = optionNames{j};
                    switch class(options.Core__.(name))
                        case {'cell','struct'}
                            assertEqual(options.(process.wholeCellModelID).(name), process.(name));
                        otherwise
                            assertElementsAlmostEqual(...
                                options.(process.wholeCellModelID).(name),...
                                process.(name),...
                                'relative', 1e-4, name);
                    end
                end
            end

            %parameters
            parameterNames = sim.parameterNames;
            simulationfieldnames = fieldnames(sim);
            for i = 1:length(parameterNames)
                name = parameterNames{i};
                if ismember(name, simulationfieldnames)
                    if isnumeric(parameters.Core__.(name))
                        assertElementsAlmostEqual(...
                            parameters.Core__.(name), sim.(name),...
                            'relative', 1e-4, name);
                    else
                        assertEqual(parameters.Core__.(name), sim.(name), name);
                    end
                end
            end
            for i = 1:length(sim.processes)
                process = sim.processes{i};
                parameterNames = process.parameterNames;
                processFieldnames = fieldnames(process);
                for j = 1:length(parameterNames)
                    name = parameterNames{j};
                    if ismember(name, processFieldnames)
                        if isnumeric(parameters.(process.wholeCellModelID).(name))
                            assertElementsAlmostEqual(...
                                parameters.(process.wholeCellModelID).(name), ...
                                process.(name),...
                                'relative', 1e-4, name);
                        else
                            assertEqual(...
                                parameters.(process.wholeCellModelID).(name), ...
                                process.(name), name);
                        end
                    end
                end
            end

            %fit constants
            fittedConstantNames = sim.fittedConstantNames;
            for i = 1:length(fittedConstantNames)
                name = fittedConstantNames{i};
                assertElementsAlmostEqual(...
                    fittedConstants.Core__.(name), sim.(name),...
                    'relative', 1e-4, name);
            end
            for i = 1:length(sim.processes)
                process = sim.processes{i};
                fittedConstantNames = process.fittedConstantNames;
                for j = 1:length(fittedConstantNames)
                    name = fittedConstantNames{j};
                    assertElementsAlmostEqual(...
                        fittedConstants.(process.wholeCellModelID).(name),...
                        process.(name),...
                        'relative', 1e-4, name);
                end
            end

            %time courses
            timeCourseNames = sim.timeCourseNames;

            for i = 1:length(timeCourseNames)
                name = timeCourseNames{i};
                assertElementsAlmostEqual(...
                    timeCourses.Core__.(name), sim.(name),...
                    'relative', 1e-4, name);
            end

            for i = 1:length(sim.processes)
                process = sim.processes{i};
                timeCourseNames = process.timeCourseNames;
                for j = 1:length(timeCourseNames)
                    name = timeCourseNames{j};
                    if isa(timeCourses.(sim.processWholeCellModelIDs{i}).(name), 'edu.stanford.covert.util.CircularSparseMat')
                        assertTrue(...
                            isalmostequal(timeCourses.Core__.(name), sim.(name)),...
                            name);
                    else
                        assertElementsAlmostEqual(...
                            timeCourses.(process.wholeCellModelID).(name), ...
                            process.(name),...
                            'relative', 1e-4,...
                            name);
                    end
                end
            end
        end
    end

    %test division
    methods
        function testDivideState(~)
            import edu.stanford.covert.cell.sim.constant.Condition;
            import edu.stanford.covert.cell.sim.util.CachedSimulationObjectUtil;
            import edu.stanford.covert.cell.sim.util.DiskLogger;
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil

            %load simulation
            %sim = CachedSimulationObjectUtil.load();

            %baseOutDir = [SimulationDiskUtil.getBaseDir() filesep '2012_10_24_00_49_53'];
            %simDir = [baseOutDir filesep '2'];
            %md = DiskLogger.loadMetadata(simDir);
            %sim = loadSimulation(simDir, md.lengthSec, md.lengthSec, 1, 'extract');

            tmp = load('src_test/+edu/+stanford/+covert/+cell/+sim/fixtures/Simulation_EndOfCellCycle.mat');
            sim = tmp.simulation;

            %references to constants, states, processes
            c = sim.compartment;
            chr = sim.state('Chromosome');
            met = sim.state('Metabolite');
            r = sim.state('Rna');
            pm = sim.state('ProteinMonomer');
            pc = sim.state('ProteinComplex');
            trl = sim.process('Translation');

            %divide cell into 2 daughters
            daughters = sim.divideState();
            daughter1 = daughters{1};
            daughter2 = daughters{2};

            %test options, parameters copied correctly
            assertEqual(sim.getOptions(), daughter1.getOptions());
            assertEqual(sim.getParameters(), daughter1.getParameters());

            %% test
            % mass
            assertElementsAlmostEqual(sim.state('Mass').cell, ...
                + daughter1.state('Mass').cell ...
                + daughter2.state('Mass').cell, 'relative', 1e-100);

            %bound proteins
            d1_chr = daughter1.state('Chromosome');
            d2_chr = daughter2.state('Chromosome');
            d1_pm = daughter1.state('ProteinMonomer');
            d2_pm = daughter2.state('ProteinMonomer');
            d1_pc = daughter1.state('ProteinComplex');
            d2_pc = daughter2.state('ProteinComplex');

            [~, vals] = find(chr.monomerBoundSites);
            [~, vals1] = find(d1_chr.monomerBoundSites);
            [~, vals2] = find(d2_chr.monomerBoundSites);
            cnts = histc(vals, (1:numel(pm.boundIndexs)));
            cnts1 = histc(vals1, (1:numel(pm.boundIndexs)));
            cnts2 = histc(vals2, (1:numel(pm.boundIndexs)));
            cnts = cnts(:);
            cnts1 = cnts1(:);
            cnts2 = cnts2(:);
            assertEqual(cnts, cnts1 + cnts2);
            trlFactorIdxs = setdiff(pm.getIndexs(trl.enzymeWholeCellModelIDs(trl.enzymeIndexs_elongationFactors)), 0);
            idxs = setdiff((1:numel(pm.boundIndexs))', trlFactorIdxs);
            assertEqual(cnts(idxs), pm.counts(pm.boundIndexs(idxs), c.cytosolIndexs));

            [~, vals] = find(chr.complexBoundSites);
            [~, vals1] = find(d1_chr.complexBoundSites);
            [~, vals2] = find(d2_chr.complexBoundSites);
            cnts = histc(vals, (1:numel(pc.boundIndexs)));
            cnts1 = histc(vals1, (1:numel(pc.boundIndexs)));
            cnts2 = histc(vals2, (1:numel(pc.boundIndexs)));
            cnts = cnts(:);
            cnts1 = cnts1(:);
            cnts2 = cnts2(:);
            assertEqual(cnts, cnts1 + cnts2);
            trlFactorIdxs = setdiff(pc.getIndexs(trl.enzymeWholeCellModelIDs(trl.enzymeIndexs_elongationFactors)), 0);
            idxs = setdiff((1:numel(pc.boundIndexs))', [
                trlFactorIdxs
                pc.getIndexs({'MG_224_9MER_GTP'; 'MG_224_9MER_GDP'})
                pc.ribosome70SIndexs
                ]);
            assertEqual(cnts(idxs), pc.counts(pc.boundIndexs(idxs), c.cytosolIndexs));

            %total metabolites, RNA, protein except, FtsZ
            d1_met = daughter1.state('Metabolite');
            d2_met = daughter2.state('Metabolite');
            d1_r = daughter1.state('Rna');
            d2_r = daughter2.state('Rna');

            assertEqual(met.counts(:, setdiff(1:end, c.extracellularIndexs)), ...
                + d1_met.counts(:, setdiff(1:end, c.extracellularIndexs)) ...
                + d2_met.counts(:, setdiff(1:end, c.extracellularIndexs)));

            mediaIdxs = met.setCounts(:, Condition.objectIndexs);
            nonMediaIdxs = setdiff((1:numel(met.wholeCellModelIDs))', mediaIdxs);
            assertEqual(met.counts(nonMediaIdxs, c.extracellularIndexs), ...
                + d1_met.counts(nonMediaIdxs, c.extracellularIndexs) ...
                + d2_met.counts(nonMediaIdxs, c.extracellularIndexs));
            assertEqual(...
                d1_met.counts(mediaIdxs, c.extracellularIndexs), ...
                d2_met.counts(mediaIdxs, c.extracellularIndexs));

            assertEqual(r.counts, ...
                + d1_r.counts ...
                + d2_r.counts);
            assertEqual(pm.counts, ...
                + d1_pm.counts ...
                + d2_pm.counts);

            tmp = pc.getIndexs({'MG_224_9MER_GTP'; 'MG_224_9MER_GDP'});
            ftsZIdxs = [pc.matureIndexs(tmp) pc.boundIndexs(tmp)];
            idxs = setdiff(1:numel(pc.wholeCellModelIDs), ftsZIdxs(:));
            assertEqual(pc.counts(idxs, :), ...
                + d1_pc.counts(idxs, :) ...
                + d2_pc.counts(idxs, :));
            assertEqual(sum(pc.counts(ftsZIdxs(1, :), :), 1), ...
                + sum(d1_pc.counts(ftsZIdxs(1, :), :), 1) ...
                + sum(d2_pc.counts(ftsZIdxs(1, :), :), 1));
            assertEqual(sum(pc.counts(ftsZIdxs(2, :), :), 1), ...
                + sum(d1_pc.counts(ftsZIdxs(2, :), :), 1) ...
                + sum(d2_pc.counts(ftsZIdxs(2, :), :), 1));
            assertAllEqual(0, d1_pc.counts(ftsZIdxs(:, 2), :));
            assertAllEqual(0, d2_pc.counts(ftsZIdxs(:, 2), :));

            %metabolite, RNA, protein counts non-negative
            validateattributes(daughter1.state('Metabolite').counts, {'numeric'}, {'nonnegative'});
            validateattributes(daughter2.state('Metabolite').counts, {'numeric'}, {'nonnegative'});
            validateattributes(d1_r.counts, {'numeric'}, {'nonnegative'});
            validateattributes(d2_r.counts, {'numeric'}, {'nonnegative'});
            validateattributes(d1_pm.counts, {'numeric'}, {'nonnegative'});
            validateattributes(d2_pm.counts, {'numeric'}, {'nonnegative'});
            validateattributes(d1_pc.counts, {'numeric'}, {'nonnegative'});
            validateattributes(d2_pc.counts, {'numeric'}, {'nonnegative'});

            %total transcripts, RNA polymerase
            rnaPol = sim.state('RNAPolymerase');
            d1_rnaPol = daughter1.state('RNAPolymerase');
            d2_rnaPol = daughter2.state('RNAPolymerase');
            transcript = sim.state('Transcript');
            d1_transcript = daughter1.state('Transcript');
            d2_transcript = daughter2.state('Transcript');

            tmp1 = [rnaPol.states rnaPol.positionStrands transcript.boundTranscriptionUnits transcript.boundTranscriptProgress transcript.boundTranscriptChromosome];
            tmp2 = [
                d1_rnaPol.states d1_rnaPol.positionStrands(:, 1) d1_rnaPol.positionStrands(:, 2)                                                 d1_transcript.boundTranscriptionUnits d1_transcript.boundTranscriptProgress d1_transcript.boundTranscriptChromosome
                d2_rnaPol.states d2_rnaPol.positionStrands(:, 1) d2_rnaPol.positionStrands(:, 2)+2*(d2_rnaPol.states ~= rnaPol.freeValue) d2_transcript.boundTranscriptionUnits d2_transcript.boundTranscriptProgress 2*d2_transcript.boundTranscriptChromosome
                ];
            tmp1 = tmp1(tmp1(:, 1) ~= rnaPol.notExistValue, :);
            tmp2 = tmp2(tmp2(:, 1) ~= rnaPol.notExistValue, :);
            tmp1 = sortrows(tmp1, 1:6);
            tmp2 = sortrows(tmp2, 1:6);
            assertEqual(tmp1, tmp2);

            for i = 1:numel(daughters)
                d_pc = daughters{i}.state('ProteinComplex');
                d_rnaPol = daughters{i}.state('RNAPolymerase');
                d_transcript = daughters{i}.state('Transcript');

                assertEqual(...
                    sum(d_pc.counts(pc.boundIndexs(pc.rnaPolymeraseIndexs), c.cytosolIndexs)), ...
                    d_rnaPol.nActive + d_rnaPol.nNonSpecificallyBound + d_rnaPol.nSpecificallyBound);
                assertEqual(...
                    d_pc.counts(pc.boundIndexs(pc.rnaPolymeraseIndexs(2), c.cytosolIndexs)), ...
                    sum(d_rnaPol.states == rnaPol.activelyTranscribingValue) + d_rnaPol.nSpecificallyBound);
                assertEqual(...
                    d_pc.counts(pc.boundIndexs(pc.rnaPolymeraseIndexs(1), c.cytosolIndexs)), ...
                    sum(d_rnaPol.states > rnaPol.activelyTranscribingValue) + d_rnaPol.nNonSpecificallyBound);

                assertEqual(...
                    sum(d_pc.counts(pc.matureIndexs(pc.rnaPolymeraseIndexs), c.cytosolIndexs)), ...
                    d_rnaPol.nFree);

                assertEqual(double(d_rnaPol.states >= rnaPol.activelyTranscribingValue | d_rnaPol.states == rnaPol.specificallyBoundValue)', ...
                    d_transcript.boundTranscriptChromosome');
            end

            %total polypeptides, ribosomes
            rib = sim.state('Ribosome');
            d1_rib = daughter1.state('Ribosome');
            d2_rib = daughter2.state('Ribosome');

            tmp1 = [rib.states rib.boundMRNAs rib.mRNAPositions rib.tmRNAPositions];
            tmp2 = [
                d1_rib.states d1_rib.boundMRNAs d1_rib.mRNAPositions d1_rib.tmRNAPositions
                d2_rib.states d2_rib.boundMRNAs d2_rib.mRNAPositions d2_rib.tmRNAPositions
                ];
            tmp1 = tmp1(tmp1(:, 1) ~= rib.notExistValue, :);
            tmp2 = tmp2(tmp2(:, 1) ~= rib.notExistValue, :);
            tmp1 = sortrows(tmp1, 1:4);
            tmp2 = sortrows(tmp2, 1:4);
            assertEqual(tmp1, tmp2)

            assertEqual(d1_pc.counts(pc.boundIndexs(pc.ribosome70SIndexs), c.cytosolIndexs), d1_rib.nActive + d1_rib.nStalled);
            assertEqual(d2_pc.counts(pc.boundIndexs(pc.ribosome70SIndexs), c.cytosolIndexs), d2_rib.nActive + d2_rib.nStalled);

            assertEqual(rib.nStalled, r.counts(r.boundIndexs(trl.tmRNAGlobalIndex), c.cytosolIndexs));
            assertEqual(d1_rib.nStalled, d1_r.counts(r.boundIndexs(trl.tmRNAGlobalIndex), c.cytosolIndexs));
            assertEqual(d2_rib.nStalled, d2_r.counts(r.boundIndexs(trl.tmRNAGlobalIndex), c.cytosolIndexs));

            %chromosome
            for i = 1:numel(chr.stateNames)
                sName = chr.stateNames{i};
                if isequal(sName, 'segregated')
                    continue;
                end
                assertEqual(chr.(sName), ...
                    [daughter1.state('Chromosome').(sName)(:, [1 2]) ...
                    daughter2.state('Chromosome').(sName)(:, [1 2])]);
            end
            assertEqual(false, daughter1.state('Chromosome').segregated)
            assertEqual(false, daughter2.state('Chromosome').segregated)

            %stimulus
            assertEqual(sim.state('Stimulus').values, daughter1.state('Stimulus').values)
            assertEqual(sim.state('Stimulus').values, daughter2.state('Stimulus').values)

            %geometry
            assertElementsAlmostEqual(sim.state('Geometry').volume, ...
                + daughter1.state('Geometry').volume ...
                + daughter2.state('Geometry').volume, ...
                'relative', 1e-24);

            %FtsZ ring
            assertEqual(0, daughter1.state('FtsZRing').numEdgesOneStraight);
            assertEqual(0, daughter1.state('FtsZRing').numEdgesTwoStraight);
            assertEqual(0, daughter1.state('FtsZRing').numEdgesTwoBent);
            assertEqual(0, daughter1.state('FtsZRing').numResidualBent);

            assertEqual(0, daughter2.state('FtsZRing').numEdgesOneStraight);
            assertEqual(0, daughter2.state('FtsZRing').numEdgesTwoStraight);
            assertEqual(0, daughter2.state('FtsZRing').numEdgesTwoBent);
            assertEqual(0, daughter2.state('FtsZRing').numResidualBent);

            %host
            sNames = sim.state('Host').stateNames;
            for i = 1:numel(sNames)
                assertEqual(sim.state('Host').(sNames{i}), daughter1.state('Host').(sNames{i}));
            end

            %time
            assertEqual(0, daughter1.state('Time').values)
            assertEqual(0, daughter2.state('Time').values)
        end

        function testInitializeFromDaughter(~)
            import edu.stanford.covert.cell.sim.util.CachedSimulationObjectUtil;

            %sim = CachedSimulationObjectUtil.load();
            tmp = load('src_test/+edu/+stanford/+covert/+cell/+sim/fixtures/Simulation_EndOfCellCycle.mat');
            sim = tmp.simulation;

            %divide cell
            [~, state_daughters] = sim.divideState();

            %setup daughter cell simulation
            tmpDirName = 'output/runMediumTests/initializeFromDaughter';
            if exist(tmpDirName, 'dir')
                rmdir(tmpDirName, 's')
            end

            mkdir(tmpDirName)

            state_daughter1 = state_daughters(1);
            save([tmpDirName filesep 'initialConditions.mat'], '-struct', 'state_daughter1');
            copyfile('src_test/+edu/+stanford/+covert/+cell/+sim/fixtures/initializeFromDaughter.xml', ...
                [tmpDirName filesep 'conditions.xml']);

            %run simulation
            daughter = runSimulation('outDir', tmpDirName, 'logToDisk', true);

            %test options, parameters copied correctly
            sim.setForTest('lengthSec', 100);
            sim.setForTest('seed', daughter.seed);
            sim.setForTest('media', reshape(sim.media, [], 6));
            sim.setForTest('stimulus', reshape(sim.stimulus, [], 6));
            sim.setForTest('geneticKnockouts', reshape(sim.geneticKnockouts, [], 1));

            for i = 1:numel(sim.states)
                s = sim.states{i};
                s.seed = daughter.states{i}.seed;
            end
            for i = 1:numel(sim.processes)
                p = sim.processes{i};
                p.seed = daughter.processes{i}.seed;
            end

            assertEqual(sim.getOptions(), daughter.getOptions());
            assertEqual(sim.getParameters(), daughter.getParameters());

            %test state copied correctly
            tmp = load([tmpDirName filesep 'state-0.mat']);
            assertEqual(state_daughter1, tmp)

            %clean up
            rmdir(tmpDirName, 's')
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function val=cosInt(x)
%%COSINT Evaluate the cosine integral of a real non-negative value x. The
%        cosine integral of x is defined to be the negative of the integral
%        from x to infinity of cos(t)/t dt. The derivative of the cosine
%        integral function is cos(x)/x.
%
%INPUTS:    x A matrix of non-negative real values, at which one wishes to
%             evaluate the cosine integral. x does not have to be finite.
%
%OUTPUTS: val A matrix of the real values of the cosine integral evaluated
%             at the points in x.
%
%The cosine integral can be related to the exponential integral as
%discussed in [1]. This function uses Matlab's built-in expint function to
%evaluate the cosine integral using the expint function.
%
%The cosine integral asymptotically approaches zero. If the results from
%the exponential integral functions are not finite, then it is assumed that
%the corresponding values of x were too large and the asymptotic result (0)
%is substituted.
%
%x should be non-negative. If negative values of x are (incorrectly)
%provided, then the value of cosInt(abs(x)) is returned. Strict definitions
%of the cosine integral are usually such that negative values return
%complex numbers.
%
%REFERENCES:
%[1] Weisstein, Eric W. "Cosine Integral." From MathWorld--A Wolfram Web
%    Resource. http://mathworld.wolfram.com/CosineIntegral.html
%
%June 2014 David F.Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    %The real function is in case of finite precision problems.
    val=-0.5*real(expint(1j*x)+expint(-1j*x));
    val(~isfinite(val))=0;

    %Make sure that NaN inputs remain NaN on the output.
    val(isnan(x))=NaN;
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function cnn = cnnupdateweights(cnn)
%CNNUPDATEWEIGHTS Updates the weights and biases of a cnn.
%   CNN = CNNUPDATEWEIGHTS(CNN) Updates the kernel weights and biases of
%   convolutional neural net CNN.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Number of layers.
Lc = numel([cnn.NumLayerKernels]);

% All the alphas are the same.
alpha = cnn(1).Alpha;

% Update cnn weights and biases.
for k = 1:Lc
   for j = 1:cnn(k).NumOutputMaps
      for i = 1:cnn(k).NumInputMaps
         % The following implement Eqs. (14-65) and (14-66) in DIPUM3E.
         cnn(k).Kernel{i}{j} = cnn(k).Kernel{i}{j} - alpha*cnn(k).WeightGradient{i}{j};
      end
         cnn(k).Bias{j} = cnn(k).Bias{j} - alpha*cnn(k).BiasGradient{j};
   end
end

author: dipum

repo: dipum-toolbox

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


text:
function dispsig(signalMatrix, range, titlestr);
%DISPSIG - deprecated!
%
% Please use icaplot instead.
%
%   See also ICAPLOT

% @(#)$Id$

fprintf('\nNote: DISPSIG is now deprecated! Please use ICAPLOT.\n');

if nargin < 3, titlestr = ''; end
if nargin < 2, range = 1:size(signalMatrix, 1); end

icaplot('dispsig',signalMatrix',0,range,range,titlestr);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function gl = spm_get_volumes(P)
% Compute total volumes from tissue segmentations
% FORMAT gl = spm_get_volumes(P)
% P  - a matrix of image filenames
% gl - a vector of volumes (in litres)
%__________________________________________________________________________
% Copyright (C) 2006-2011 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_get_volumes.m 5670 2013-10-04 16:48:38Z ged $

warning('spm:deprecated', ...
    ['spm_get_volumes will be removed in the future, please use ' ...
    'the Tissue Volumes Utility in the Batch interface, or:\n\t' ...
    'spm_summarise(P, ''all'', ''litres'')']);

if ~nargin
    [P,sts] = spm_select(Inf,'image','Select images');
    if ~sts, gl = []; return; end
end

V = spm_vol(P);
if spm_check_orientations(V, false)
    gl = spm_summarise(V, 'all', 'litres');
else
    N = numel(V);
    gl = nan(N, 1);
    for n = 1:N
        gl(n) = spm_summarise(V(n), 'all', 'litres');
    end
end

author: spm


text:
function [DF] = df_build_fast(V,C)
% DF_BUILD_FAST Build a volumetric distance field and spatial index that can be
% used to speed up approximate NN queries in 3D This code uses the kdtree
% toolbox
%
% [DF] = df_build_fast(V,C)
%
% Input:
%   V  Coordinates of the points
%   C  Column vector that defines the number of cells for every
%      dimension
% Output:
%   DF  structure to be passed to df_query to execute queries

if ~exist('C','var')
    C=[10;10;10];
end

MIN = min(V,[],1);
MAX = max(V,[],1);

S = (MAX-MIN)./C';

D = zeros(C'+1);
N = zeros(C'+1);

%% Build kd-tree
tree = kdtree_build( V );

%% Queries
for x=1:size(D,1)
    progressbar(x,size(D,1));
    for y=1:size(D,2)
        for z=1:size(D,3)
            i = [x,y,z];
            p = MIN + (i-1).*S;

            N(x,y,z) = kdtree_nearest_neighbor(tree,p);
            D(x,y,z) = normrow(p-V(N(x,y,z),:));
        end
    end
end

kdtree_delete(tree);

% Prepare DF struct
DF.MIN = MIN;
DF.MAX = MAX;
DF.D = D;
DF.N = N;
DF.C = C;
DF.S = S;
end

author: alecjacobson

repo: gptoolbox


%==============================================================================
% 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:  Landmark Based Registration, TPS
%
% - load data (see setup2DhandData)
% - setup  viewer (viewImage2D), interpolator (splineInter),
% - setup landmarks (LM)
% - run TPS
%==============================================================================

clear, close all, help(mfilename)

%% setup hand data
setup2DhandData

if FAIRinput('','set new landmarks ? ',0),
  [LM,fig] = getLandmarks(dataT,dataR,omega,m);
  close(fig);
end;

omegaT = omega(1,:);
omegaR = omega(end,:);
xT = getCellCenteredGrid(omegaT,m);
xR = getCellCenteredGrid(omegaR,m);
Tc = imgModel(dataT,omegaT,xT);
Rc = imgModel(dataR,omegaR,xR);

%% visualize data
FAIRfigure(1,'figname',mfilename); clf;
subplot(1,3,1); viewImage(Tc,omegaT,m); hold on;
ph = plotLM(LM(:,1:2),'numbering','on','color','r');
set(ph,'linewidth',2,'markersize',20);
title(sprintf('%s','T&LM'),'fontsize',20);

subplot(1,3,2); viewImage(Rc,omegaR,m); hold on;
ph = plotLM(LM(:,3:4),'numbering','on','color','g','marker','+');
set(ph,'linewidth',2,'markersize',20);
title(sprintf('%s','R&LM'),'fontsize',20);

%% run TPS registration
[yc,LM] = LMreg('TPS',LM(:,1:4),xR);
TLM = imgModel(dataT,omegaT,yc);

subplot(1,3,3); cla; viewImage(TLM,omegaR,m); hold on;
ph = plotLM(LM(:,3:4),'numbering','off','color','g','marker','+');
qh = plotLM(LM(:,7:8),'numbering','off','color','m','marker','x');
rh = plot(LM(:,[3,7])',LM(:,[4,8])','m-','linewidth',3);
set([ph;qh;rh],'linewidth',2,'markersize',20);
title(sprintf('%s','T(y^{TPS})&LM'),'fontsize',20);
%==============================================================================


text:
% FT_POSTAMBLE_HASTOOLBOX is executed at the end of each FieldTrip
% function to remove other toolboxes that have been added automatically
% by FT_HASTOOLBOX during execution of the specific function.
%
% Use as
%   ft_postamble hastoolbox
%
% See also FT_PREAMBLE, FT_POSTAMBLE, FT_HASTOOLBOX

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

global ft_default

if ~isempty(ft_default) && isfield(ft_default, 'toolbox') && isfield(ft_default.toolbox, 'cleanup')
  while ~isempty(ft_default.toolbox.cleanup)
    toolbox = ft_default.toolbox.cleanup{end};
    ft_warning('off','backtrace');
    ft_warning('removing %s toolbox from your MATLAB path', toolbox);
    ft_warning('on','backtrace');
    rmpath(genpath(toolbox));
    ft_default.toolbox.cleanup = ft_default.toolbox.cleanup(1:end-1);
  end
end


text:
function fg = mk_fgraph(G, node_sizes, factors, varargin)
% MK_FGRAPH Make a factor graph
% fg = mk_fgraph(G, node_sizes, factors, ...)
%
% A factor graph is a bipartite graph, with one side containing variables,
% and the other containing functions of (subsets of) these variables.
% For details, see "Factor Graphs and the Sum-Product Algorithm",
%  F. Kschischang and B. Frey and H-A. Loeliger,
%  IEEE Trans. Info. Theory, 2001
%
% G(i,j) = 1 if there is an arc from variable i to factor j
%
% node_sizes(i) is the number of values node i can take on,
%   or the length of node i if i is a continuous-valued vector.
%
% 'factors' is the list of factors (kernel functions)
%
% The list below gives optional arguments [default value in brackets].
%
% equiv_class - equiv_class(i)=j  means factor node i gets its params from factors{j} [1:F]
% discrete - the list of nodes which are discrete random variables [1:N]
%
% e.g., fg = mk_fgraph(G, [2 2], {bnet.CPD{1},bnet.CPD{2}}, 'discrete', [1 2])

fg.G = G;
fg.node_sizes = node_sizes;
fg.factors = factors;
[fg.nvars fg.nfactors] = size(G);

% default values for parameters
fg.equiv_class = 1:fg.nfactors;
fg.dnodes = 1:fg.nvars;

if nargin >= 4
  args = varargin;
  nargs = length(args);
  for i=1:2:nargs
    switch args{i},
     case 'equiv_class', fg.equiv_class = args{i+1};
     case 'discrete',    fg.dnodes = args{i+1};
     otherwise,
      error(['invalid argument name ' args{i}]);
    end
  end
end

% so that determine_pot_type will work...
fg.utility_nodes = [];
%fg.decision_nodes = [];
%fg.chance_nodes = fg.nvars;

fg.dom = cell(1, fg.nfactors);
for f=1:fg.nfactors
  fg.dom{f} = find(G(:,f));
end
fg.dep = cell(1, fg.nvars);
for x=1:fg.nvars
  fg.dep{x} = find(G(x,:));
end
fg.cnodes = mysetdiff(1:fg.nvars, fg.dnodes);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [x,esq,j] = kmeanlbg(d,k)
%KMEANLBG Vector quantisation using the Linde-Buzo-Gray algorithm [X,ESQ,J]=(D,K)
%
%Inputs:
% D contains data vectors (one per row)
% K is number of centres required
%
%Outputs:
% X is output row vectors (K rows)
% ESQ is mean square error
% J indicates which centre each data vector belongs to
%
%  Implements LBG K-means algorithm:
% Linde, Y., A. Buzo, and R. M. Gray,
% "An Algorithm for vector quantiser design,"
% IEEE Trans Communications, vol. 28, pp.84-95, Jan 1980.

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

nc=size(d,2);
[x,esq,j]=kmeans(d,1);
m=1;
while m<k
   n=min(m,k-m);
   m=m+n;
   e=1e-4*sqrt(esq)*rand(1,nc);
   [x,esq,j]=kmeans(d,m,[x(1:n,:)+e(ones(n,1),:); x(1:n,:)-e(ones(n,1),:); x(n+1:m-n,:)]);
end

author: decouples


text:
classdef TestBm3dDenoising
    %TestBm3dDenoising

    methods (Static)
        function test_1
            img = cv.imread(fullfile(mexopencv.root(),'test','cat.jpg'), ...
                'Grayscale',true, 'ReduceScale',2);
            out = cv.bm3dDenoising(img, 'H',10);
            validateattributes(out, {class(img)}, {'size',size(img)});
        end

        function test_error_argnum
            try
                cv.bm3dDenoising();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% Example using the k-means algorithm function on an image and clustering
% the letters as seperate clusters

test = imread('Example.bmp');
test = test>=100;
test = test(:,:,1);

imagesc(test); colormap('gray')

% initial guess
means = [100,100; 200, 200; 150, 300; 50, 300];

hold on

% before the first step
plot(means(:,2),means(:,1),'o','Color',[1,0,0],'MarkerSize',10,'LineWidth',4)

% one step
means = kmean(test,means);
plot(means(:,2),means(:,1),'o','Color',[0,1,0],'MarkerSize',10,'LineWidth',4)

% second step
means = kmean(test,means);
plot(means(:,2),means(:,1),'o','Color',[0,0,1],'MarkerSize',10,'LineWidth',4)

author: Sable

repo: mcbench-benchmarks


text:
function obj = t_port_inj_power_acp(quiet)
%T_PORT_INJ_POWER_ACP  Tests of port_inj_power() derivatives wrt polar V.

%   MATPOWER
%   Copyright (c) 2019-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

if nargin < 1
    quiet = 0;
end

tc = struct( ...        %% test cases
        'name', {'1', '2'}, ...
        'ec', { @mp.nme_gizmo_acp, ...
                {   {@mp.nme_gen_acp_nln, 'mp.nme_gen'}, ...
                    {@mp.nme_load_acp_nln, 'mp.nme_load'}, ...
                    {@mp.nme_branch_acp_nln, 'mp.nme_branch'}, ...
                    {@mp.nme_shunt_acp_nln, 'mp.nme_shunt'}, ...
                    @mp.nme_gizmo_acp_nln }    } ...
    );

t_begin(87*length(tc), quiet);

define_constants;
if quiet
    verbose = 0;
else
    verbose = 1;
end

casefile = 't_case9_gizmo';
mpopt = mpoption('out.all', 0, 'verbose', 0);
dmc = mp.dm_converter_mpc2().modify_element_classes(@mp.dmce_gizmo_mpc2).build();

for c = 1:length(tc)
    %% create network model object
    mpc = loadcase(casefile);
    dm = mp.data_model().modify_element_classes(@mp.dme_gizmo).build(mpc, dmc);
    ac = mp.net_model_acp().modify_element_classes(@mp.nme_gizmo_acp).build(dm);
    C = ac.C;
    D = ac.D;
    np = ac.np;
    nv = ac.nv/2;
    nz = ac.nz;
    A = [   C sparse(nv, nz);
            sparse(nz, np) D    ];
    A2 = [A sparse(nv+nz, np+nz); sparse(nv+nz, np+nz) A];

    %% other parameters
    dx = 1e-8;
    idx = randperm(np, fix(0.67*np))';
    lam = (1.5*rand(np, 1) + 0.5); k = randperm(np, fix(np/2)); lam(k) = -lam(k);
    e0 = zeros(np, 1);
    e1 = ones(np, 1);

    %% get ref bus index
    ref = find(mpc.bus(:, BUS_TYPE) == REF);

    %% construct initial system v1, v2, zr, zi, v_, z_, x_
    t = sprintf('%s : construct initial system v_, z_', tc(c).name);
    sv1 = ac.params_var('va');
    sv2 = ac.params_var('vm');
    szr = ac.params_var('zr');
    szi = ac.params_var('zi');

    %% randomize voltages a bit
    sv1 = sv1 + (0.6*rand(size(sv1)) - 0.3);    sv1(ref) = 0;
    sv2 = sv2 + (0.06*rand(size(sv1)) - 0.03);  sv2(ref) = 1;

    %% adjust values of z_
    szr(1) = 0.67;
    szi(1:3) = [0.1; 0.2; 0.3];

    %% initialize v_, z_, x_
    sv = sv2 .* exp(1j * sv1);
    sz = szr + 1j * szi;
    sx = [sv; sz];
    nx = length(sx);
    t_is(nx, nv+nz, 12, t);

    %%-----  tests using system voltages  -----
    t = sprintf('%s : ac.port_inj_power(x_) : ', tc(c).name);
    v10 = sv1; v20 = sv2; zr0 = szr; zi0 = szi; %% init w/ system v_, z_ components
    v0 = v20 .* exp(1j * v10);
    z0 = zr0 + 1j * zi0;
    x0 = [v0; z0];
    Nv = length(v0);
    Nz = length(z0);
    [S0, Sv1, Sv2, Szr, Szi] = ac.port_inj_power(x0);   %% analytical

    %% check matrix input/output
    SS0 = ac.port_inj_power(x0*ones(1,Nv));
    t_is(SS0, S0*ones(1,Nv), 12, [t 'matrix input']);

    %% Sv1
    v_ = (v20*ones(1,Nv)) .* exp(1j * (v10*ones(1,Nv) + dx*eye(Nv,Nv)));
    z_ = (zr0 + 1j * zi0) * ones(1,Nv);
    x_ = [v_; z_];
    SS = ac.port_inj_power(x_);
    num_Sv1b = (SS - SS0) / dx;
    t_is(full(Sv1), num_Sv1b, 6, [t 'Sv1']);

    %% Sv2
    v_ = (v20*ones(1,Nv) + dx*eye(Nv,Nv)) .* exp(1j * (v10*ones(1,Nv)));
    z_ = (zr0 + 1j * zi0) * ones(1,Nv);
    x_ = [v_; z_];
    SS = ac.port_inj_power(x_);
    num_Sv2b = (SS - SS0) / dx;
    t_is(full(Sv2), num_Sv2b, 6, [t 'Sv2']);

    SS0 = ac.port_inj_power(x0*ones(1,Nz));

    %% Szr
    v_ = (v20*ones(1,Nz)) .* exp(1j * (v10*ones(1,Nz)));
    z_ = (zr0 * ones(1,Nz) + dx*eye(Nz,Nz)) + 1j * (zi0 * ones(1,Nz));
    x_ = [v_; z_];
    SS = ac.port_inj_power(x_);
    num_Szrb = (SS - SS0) / dx;
    t_is(full(Szr), num_Szrb, 6, [t 'Szr']);

    %% Szi
    v_ = (v20*ones(1,Nz)) .* exp(1j * (v10*ones(1,Nz)));
    z_ = (zr0 * ones(1,Nz)) + 1j * (zi0 * ones(1,Nz) + dx*eye(Nz,Nz));
    x_ = [v_; z_];
    SS = ac.port_inj_power(x_);
    num_Szib = (SS - SS0) / dx;
    t_is(full(Szi), num_Szib, 6, [t 'Szi']);

    t = sprintf('%s : ac.port_inj_power(x_, 1, idx) : ', tc(c).name);
    [iS0, iSv1, iSv2, iSzr, iSzi] = ac.port_inj_power(x0, 1, idx);
    t_is(iS0, S0(idx), 12, [t 'S0']);
    t_is(iSv1, Sv1(idx, :), 12, [t 'Sv1']);
    t_is(iSv2, Sv2(idx, :), 12, [t 'Sv2']);
    t_is(iSzr, Szr(idx, :), 12, [t 'Szr']);
    t_is(iSzi, Szi(idx, :), 12, [t 'Szi']);

    t = sprintf('%s : ac.port_inj_power_hess(x0, ek) == ac.p_i_p_h(x0, 1, 1, k) : ', tc(c).name);
    for k = 1:length(lam)
        ek = e0; ek(k) = 1;
        H1 = ac.port_inj_power_hess(x0, ek);
        H2 = ac.port_inj_power_hess(x0, 1, 1, k);
        t_is(H1, H2, 12, sprintf('%s%d', t, k));
    end

    t = sprintf('%s : ac.port_inj_power_hess(x_, lam) : ', tc(c).name);
    H = ac.port_inj_power_hess(x0, lam);
    HH = sparse(size(H, 1), size(H, 2));
    for k = 1:length(lam)
        ek = e0; ek(k) = 1;
        HH = HH + lam(k) * ac.port_inj_power_hess(x0, ek);
    end
    t_is(H, HH, 12, [t 'weighted sum indiv Hessians']);

    t = sprintf('%s : ac.port_inj_power_hess(x_, lam, 1, idx) : ', tc(c).name);
    H = ac.port_inj_power_hess(x0, lam(idx), 1, idx);
    HH = sparse(size(H, 1), size(H, 2));
    for k = 1:length(idx)
        ek = e0; ek(idx(k)) = 1;
        HH = HH + lam(idx(k)) * ac.port_inj_power_hess(x0, ek);
    end
    t_is(H, HH, 12, [t 'weighted sum indiv Hessians']);

    t = sprintf('%s : ac.port_inj_power_hess(x_, lam) : ', tc(c).name);
    H = ac.port_inj_power_hess(x0, lam);
    [S0, Sv1, Sv2, Szr, Szi] = ac.port_inj_power(x0);
    numH = zeros(2*nx, 2*nx);
    for k = 1:Nv
        v1 = v10; v1(k) = v1(k) + dx;
        v_ = v20 .* exp(1j * v1);
        x_ = [v_; z0];
        [S0p, Sv1p, Sv2p, Szrp, Szip] = ac.port_inj_power(x_);
        numH(:, k) = ([Sv1p, Sv2p, Szrp, Szip]- [Sv1, Sv2, Szr, Szi]).' * lam / dx;

        v2 = v20; v2(k) = v2(k) + dx;
        v_ = v2 .* exp(1j * v10);
        x_ = [v_; z0];
        [S0p, Sv1p, Sv2p, Szrp, Szip] = ac.port_inj_power(x_);
        numH(:, Nv+k) = ([Sv1p, Sv2p, Szrp, Szip]- [Sv1, Sv2, Szr, Szi]).' * lam / dx;
    end
    for k = 1:Nz
        z_ = zr0 + 1j * zi0;
        z_(k) = z_(k) + dx;
        x_ = [v0; z_];
        [S0p, Sv1p, Sv2p, Szrp, Szip] = ac.port_inj_power(x_);
        numH(:, 2*Nv+k) = ([Sv1p, Sv2p, Szrp, Szip]- [Sv1, Sv2, Szr, Szi]).' * lam / dx;

        z_ = zr0 + 1j * zi0;
        z_(k) = z_(k) + 1j*dx;
        x_ = [v0; z_];
        [S0p, Sv1p, Sv2p, Szrp, Szip] = ac.port_inj_power(x_);
        numH(:, 2*Nv+Nz+k) = ([Sv1p, Sv2p, Szrp, Szip]- [Sv1, Sv2, Szr, Szi]).' * lam / dx;
    end
    t_is(full(H), numH, 5, [t 'numerical Hessian']);

    %%-----  tests using port voltages  -----
    t = sprintf('%s : ac.port_inj_power(x_, 0) : ', tc(c).name);
    v10 = C'*sv1; v20 = C'*sv2; zr0 = D'*szr; zi0 = D'*szi; %% init w/ port v_, z_ components
    v0 = v20 .* exp(1j * v10);
    z0 = zr0 + 1j * zi0;
    x0 = [v0; z0];
    Nv = length(v0);
    Nz = length(z0);
    [S0, Sv1, Sv2, Szr, Szi] = ac.port_inj_power(x0, 0);    %% analytical

    %% check matrix input/output
    SS0 = ac.port_inj_power(x0*ones(1,Nv), 0);
    t_is(SS0, S0*ones(1,Nv), 12, [t 'matrix input']);

    %% Sv1
    v_ = (v20*ones(1,Nv)) .* exp(1j * (v10*ones(1,Nv) + dx*eye(Nv,Nv)));
    z_ = (zr0 + 1j * zi0) * ones(1,Nv);
    x_ = [v_; z_];
    SS = ac.port_inj_power(x_, 0);
    num_Sv1b = (SS - SS0) / dx;
    t_is(full(Sv1), num_Sv1b, 6, [t 'Sv1']);

    %% Sv2
    v_ = (v20*ones(1,Nv) + dx*eye(Nv,Nv)) .* exp(1j * (v10*ones(1,Nv)));
    z_ = (zr0 + 1j * zi0) * ones(1,Nv);
    x_ = [v_; z_];
    SS = ac.port_inj_power(x_, 0);
    num_Sv2b = (SS - SS0) / dx;
    t_is(full(Sv2), num_Sv2b, 6, [t 'Sv2']);

    SS0 = ac.port_inj_power(x0*ones(1,Nz), 0);

    %% Szr
    v_ = (v20*ones(1,Nz)) .* exp(1j * (v10*ones(1,Nz)));
    z_ = (zr0 * ones(1,Nz) + dx*eye(Nz,Nz)) + 1j * (zi0 * ones(1,Nz));
    x_ = [v_; z_];
    SS = ac.port_inj_power(x_, 0);
    num_Szrb = (SS - SS0) / dx;
    t_is(full(Szr), num_Szrb, 6, [t 'Szr']);

    %% Szi
    v_ = (v20*ones(1,Nz)) .* exp(1j * (v10*ones(1,Nz)));
    z_ = (zr0 * ones(1,Nz)) + 1j * (zi0 * ones(1,Nz) + dx*eye(Nz,Nz));
    x_ = [v_; z_];
    SS = ac.port_inj_power(x_, 0);
    num_Szib = (SS - SS0) / dx;
    t_is(full(Szi), num_Szib, 6, [t 'Szi']);

    t = sprintf('%s : ac.port_inj_power(x_, 0, idx) : ', tc(c).name);
    [iS0, iSv1, iSv2, iSzr, iSzi] = ac.port_inj_power(x0, 0, idx);
    t_is(iS0, S0(idx), 12, [t 'S0']);
    t_is(iSv1, Sv1(idx, :), 12, [t 'Sv1']);
    t_is(iSv2, Sv2(idx, :), 12, [t 'Sv2']);
    t_is(iSzr, Szr(idx, :), 12, [t 'Szr']);
    t_is(iSzi, Szi(idx, :), 12, [t 'Szi']);

    t = sprintf('%s : ac.port_inj_power_hess(x0, ek, 0) == ac.p_i_p_h(x0, 1, 0, k) : ', tc(c).name);
    for k = 1:length(lam)
        ek = e0; ek(k) = 1;
        H1 = ac.port_inj_power_hess(x0, ek, 0);
        H2 = ac.port_inj_power_hess(x0, 1, 0, k);
        t_is(H1, H2, 12, sprintf('%s%d', t, k));
    end

    t = sprintf('%s : ac.port_inj_power_hess(x_, lam, 0) : ', tc(c).name);
    H = ac.port_inj_power_hess(x0, lam, 0);
    HH = sparse(size(H, 1), size(H, 2));
    for k = 1:length(lam)
        ek = e0; ek(k) = 1;
        HH = HH + lam(k) * ac.port_inj_power_hess(x0, ek, 0);
    end
    t_is(H, HH, 12, [t 'weighted sum indiv Hessians']);

    t = sprintf('%s : ac.port_inj_power_hess(x_, lam, 0, idx) : ', tc(c).name);
    H = ac.port_inj_power_hess(x0, lam(idx), 0, idx);
    HH = sparse(size(H, 1), size(H, 2));
    for k = 1:length(idx)
        ek = e0; ek(idx(k)) = 1;
        HH = HH + lam(idx(k)) * ac.port_inj_power_hess(x0, ek, 0);
    end
    t_is(H, HH, 12, [t 'weighted sum indiv Hessians']);

    t = sprintf('%s : ac.port_inj_power_hess(x_, lam, 0) : ', tc(c).name);
    H = ac.port_inj_power_hess(x0, lam, 0);
    [S0, Sv1, Sv2, Szr, Szi] = ac.port_inj_power(x0, 0);
    Nx = 2*Nv+2*Nz;
    numH = zeros(Nx, Nx);
    for k = 1:Nv
        v1 = v10; v1(k) = v1(k) + dx;
        v_ = v20 .* exp(1j * v1);
        x_ = [v_; z0];
        [S0p, Sv1p, Sv2p, Szrp, Szip] = ac.port_inj_power(x_, 0);
        numH(:, k) = ([Sv1p, Sv2p, Szrp, Szip]- [Sv1, Sv2, Szr, Szi]).' * lam / dx;

        v2 = v20; v2(k) = v2(k) + dx;
        v_ = v2 .* exp(1j * v10);
        x_ = [v_; z0];
        [S0p, Sv1p, Sv2p, Szrp, Szip] = ac.port_inj_power(x_, 0);
        numH(:, Nv+k) = ([Sv1p, Sv2p, Szrp, Szip]- [Sv1, Sv2, Szr, Szi]).' * lam / dx;
    end
    for k = 1:Nz
        z_ = zr0 + 1j * zi0;
        z_(k) = z_(k) + dx;
        x_ = [v0; z_];
        [S0p, Sv1p, Sv2p, Szrp, Szip] = ac.port_inj_power(x_, 0);
        numH(:, 2*Nv+k) = ([Sv1p, Sv2p, Szrp, Szip]- [Sv1, Sv2, Szr, Szi]).' * lam / dx;

        z_ = zr0 + 1j * zi0;
        z_(k) = z_(k) + 1j*dx;
        x_ = [v0; z_];
        [S0p, Sv1p, Sv2p, Szrp, Szip] = ac.port_inj_power(x_, 0);
        numH(:, 2*Nv+Nz+k) = ([Sv1p, Sv2p, Szrp, Szip]- [Sv1, Sv2, Szr, Szi]).' * lam / dx;
    end
    t_is(full(H), numH, 5, [t 'numerical Hessian']);
end

t_end;

if nargout
    obj = ac;
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function N = readNEIGH(filename)
  % READNEIGH  Read tetrahedral neighbor information from a .neigh file (as
  % produced by tetgen)
  %
  % N = readNEIGH(filename)
  %
  % Inputs:
  %   filename  path to .neigh file
  % Outputs:
  %   N  #simplices by #size-of-simplex neighborhood information (-1) indicates
  %     boundary. T(i,j) *should* indicate the neighbor to the jth face of the
  %     ith tet. *However* tetgen does not seem consistent. Consider
  %     post-processing with fixNEIGH.m
  %
  % See also: tt, tetgen, readNODE, readELE, fixNEIGH
  %

  fp = fopen(filename);
  line = fscanf(fp,' %[^\n]s');
  [header,count] = sscanf(line,'%d %d',2);
  if count~=2
    fclose(fp);
    error('Bad header');
  end

  % number of elements
  n = header(1);
  % size of an element
  size_e = header(2);

  parser = '%d';
  % append to parser enough to read all entries in element + 1 for index
  parser = [parser repmat(' %d',1,size_e+1)];
  N = fscanf(fp,parser,[size_e+1 n])';
  fclose(fp);

  % get rid of row indices and make one indexed
  N = N(:,2:end) + 1;
  N(~N) = -1;
end


text:
function Y = modelSamp(model, X)

% MODELSAMP Give a sample from a model for given X.

% MLTOOLS

fhandle = str2func([model.type 'Samp']);
Y = fhandle(model, X);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [day_weight]=rmr_day_weight(data_close,data,t1,day_weight, w, epsilon)
%to calculate the day_weight at the t+1 day with robust mean reversion method
%Input:
%    data
%    t1         ---the new day
%    day_weight ---the weight before the new day
%    w          ---the length of window
%    epsilon      ---the parameter to control the reversion threshold
%Output:
%    day_weight ---the weight at the new day
%w=min(w,t1-1);
%x_t1 = l1median_VaZh_z(data((t1-w):(t1-1),:));
if t1<w+2
    x_t1=data(t1-1,:);
else
   x_t1 = l1median_VaZh_z(data_close((t1-w):(t1-1),:))./data_close(t1-1,:);
end
%   if sum(x_t1)==size(data,2)
%       x_t1 = median(data_close((t1-w):(t1-1),:))./data_close(t1-1,:);
%   end

if (norm(x_t1-mean(x_t1)))^2==0
    tao = 0;
else
    tao = min(0,(x_t1*day_weight-epsilon)/(norm(x_t1-mean(x_t1)))^2);
end
day_weight = day_weight - tao*(x_t1-mean(x_t1)*ones(size(x_t1)))';
day_weight = simplex_projection(day_weight,1);

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function [a] = mci_logistic_act (P,M,U)
% Activations of logistic model
% FORMAT [a] = mci_logistic_act (P,M,U)
%
% P         parameters
% M         model structure
% U         contains rewards and times
%
% a         activations of logistic model
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

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

a = U.X*P;


function exportAndCropFigure(h,name,saveDir)
    if nargin < 3
        return;
    end
    [status,message,messageid] = mkdir(saveDir);
    [status,message,messageid] = mkdir(fullfile(saveDir,'png'));
    [status,message,messageid] = mkdir(fullfile(saveDir,'pdf'));
    [status,message,messageid] = mkdir(fullfile(saveDir,'fig'));

    saveas(h,fullfile(saveDir,name));
    saveas(h,fullfile(saveDir,'png',[name '.png']));
    saveas(h,fullfile(saveDir,'fig',[name '.fig']));
    saveas(h,fullfile(saveDir,'pdf',[name '.pdf']));

    %- auto crop pdf -%
    cmdStr = ['pdfcrop.pl' ' ' fullfile(saveDir,'pdf',name)  ...
    ' ' fullfile(saveDir,'pdf',name)];
    [status, result] = system(cmdStr);

end


author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% DEMYALESVARGPLVM4 Run the Shared Var. GP-LVM on a subset of the Yale
% faces.
% DESC Run the Shared Var. GP-LVM on a subset of the Yale faces. The code
% for creating this subset out of raw images exists in comments. Unlike
% demYaleSvargplvm1, this demo is not a wrapper, it can be used as a
% standalone demo.
%
% VARGPLVM

%---
% ca;clear; trainModel=true; indPoints = 30;  experimentNo=2; % demSitranArtificialLight1 % BAD
% ca;clear; trainModel=true; indPoints = 60;  experimentNo=3; demSitranArtificialLight1
% ca;clear; trainModel=true; indPoints = 40;  experimentNo=4; demSitranArtificialLight1 % GOOD % Liv. Machines
% ca;clear; trainModel=true; indPoints = 100; experimentNo=5; demSitranArtificialLight1 % so and so
% ca;clear; trainModel=true; indPoints = 35;  experimentNo=6;
        % initVardistIters = 1200; itNo=[500 8000 5000];demSitranArtificialLight1 % GOOD!!

%--

% Fix seeds
randn('seed', 1e4);
rand('seed', 1e4);

% Create the dataset out of the images.
baseDir=[localDatasetsDirectorySmall 'sharedVargplvm' filesep 'sitranFaces2'];
selDirs = {'arif','ricardo','zhenwen','andreas','fariba','max'};

for d=1:length(selDirs)
    dirFrom=[baseDir filesep selDirs{d}];
    a=dir(dirFrom);
    counter = 0;
    for i=1:length(a)
        if length(a(i).name)>4 & strcmp(a(i).name(1:7),'cropped')
            im = imread([dirFrom filesep a(i).name]);
            im = im(:,:,1);
            %imagesc(im), colormap('gray'); title(a(i).name), pause
            counter = counter+1;
            Yall{d}(counter,:)=im(:)';
        end
    end
    Yall{d} = double(Yall{d});
end
height = size(im,1);
width = size(im,2);
for d=1:length(Yall)
    Yall{d} = util_artificialShadow(reshape(Yall{d},height,width),9,30,[],3000);
end

%{
for i=1:size(Yall{1},1)
    for d=1:length(Yall)
        subplot(2,3,d)
        imagesc(reshape(Yall{d}(i,:), height, width)); colormap('gray');
    end
    if i==1
        pause
    else
        pause(0.1)
    end
end
%}

fpm = length(Yall)/2;
Y = Yall(1:fpm);
identities{1} = [];
for i=1:fpm
    Y{end+1}=[];
    identities{1} = [identities{1}; i*ones(size(Y{1},1),1)];
    idTmp{i+fpm}=[];
end
for i=1:size(Yall{1},1)
    perm = randperm(fpm)+fpm;%perm=[1 2 3 4]+4;
    for j=1:fpm
        Y{j+fpm} = [Y{j+fpm}; Yall{perm(j)}(i,:)]; idTmp{j+fpm}=[idTmp{j+fpm}; perm(j)];
    end
end
clear 'Yall'
Yall{1}=[]; Yall{2}=[]; identities{2}=[];
for i=1:fpm
    Yall{1} = [Yall{1}; Y{i}];
    Yall{2} = [Yall{2}; Y{i+fpm}];
    identities{2} = [identities{2}; idTmp{i+fpm}];
end

%{
Y{1} = Yall{1};
Y{2} = Yall{2};
Y{3} = Yall{3};
Y{4} = Yall{4};
Y{5} = [];
Y{6} = [];
Y{7} = [];
Y{8} = [];
identities{1}=[ones(size(Y{1},1),1); 2*ones(size(Y{2},1),1); 3*ones(size(Y{3},1),1); 4*ones(size(Y{4},1),1)];
idTmp{5}=[];idTmp{6}=[];idTmp{7}=[];idTmp{8}=[];
for i=1:size(Yall{4},1)
    perm = randperm(4)+4;%perm=[1 2 3 4]+4;
    Y{5} = [Y{5}; Yall{perm(1)}(i,:)]; idTmp{5}=[idTmp{5}; perm(1)];
    Y{6} = [Y{6}; Yall{perm(2)}(i,:)]; idTmp{6}=[idTmp{6}; perm(2)];
    Y{7} = [Y{7}; Yall{perm(3)}(i,:)]; idTmp{7}=[idTmp{7}; perm(3)];
    Y{8} = [Y{8}; Yall{perm(4)}(i,:)]; idTmp{8}=[idTmp{8}; perm(4)];
end
clear 'Yall'
Yall{1} = [Y{1}; Y{2}; Y{3}; Y{4}];
Yall{2} = [Y{5}; Y{6}; Y{7}; Y{8}];
identities{2} = [idTmp{5};idTmp{6};idTmp{7};idTmp{8};];
%}

clear 'Y' 'idTmp' 'd'

numberOfDatasets = length(Yall);

 %{
for i=1:size(Yall{1},1)
    for d=1:length(Yall)
        subplot(1,2,d)
        imagesc(reshape(Yall{d}(i,:), height, width)); colormap('gray'); title(num2str(identities{d}(i)))
    end
    if i==1
        pause
    else
        pause(0.05)
    end
end
clear 'd'
 %}
%%

if ~exist('trainModel','var'), trainModel = false; end
if ~exist('itNo')         ,  itNo = [500 1500 1500];              end     % Default: 2000
if ~exist('indPoints')    ,  indPoints = 55;          end     % Default: 49
if ~exist('initVardistIters'), initVardistIters = 900;      end
if ~exist('mappingKern')   ,  mappingKern = 'rbfardjit'; end

% Set to empty value {} to work with toy data
if ~exist('dataSetNames')    ,    dataSetNames = {};    end
if ~exist('invWidthMult'),       invWidthMult = 5;                     end
if ~exist('dataType'), dataType = 'artificialLight'; end
if ~exist('latentDimPerModel'), latentDimPerModel = 7; end
if ~exist('experimentNo'), experimentNo = 1; end
if ~exist('doPredictions'), doPredictions = false; end
% If this is true, then the model is in "D > N" mode.
if ~exist('DgtN'), DgtN = true; end
% Create initial X by doing e.g. ppca in the concatenated model.m's or by
% doing ppca in the model.m's separately and concatenate afterwards?
if ~exist('initial_X'), initial_X = 'separately'; end % Other options: 'concatenated'
% Which indices to use for training, rest for test
if ~exist('indTr'), indTr = -1; end

enableParallelism = 0;

if exist('pyramid','var') % extract pyramid representation of the images
    if pyramid
        for e=1:size(Y,2)
            Y{e} = im2pyramid(Y{e}, lbls(1), lbls(2), 4);
        end
    end
end

%-- Load datasets
for i=1:numberOfDatasets
    Y = Yall{i};
    dims{i} = size(Y,2);
    N{i} = size(Y,1);
    if indTr == -1
        indTr = 1:N{i};
    end
    indTs = setdiff(1:size(Y,1), indTr);
    Ytr{i} = Y(indTr,:);
    Yts{i} = Y(indTs,:);
    t{i} = linspace(0, 2*pi, size(Y, 1)+1)'; t{i} = t{i}(1:end-1, 1);
    timeStampsTraining{i} = t{i}(indTr,1); %timeStampsTest = t(indTs,1);
    d{i} = size(Ytr{i}, 2);
end

for i=2:numberOfDatasets
    if N{i} ~= N{i-1}
        error('The number of observations in each dataset must be the same!');
    end
end

%%
if trainModel

    %-- Options for the models
    for i=1:numberOfDatasets
        % Set up models
        options{i} = vargplvmOptions('dtcvar');
        options{i}.kern = mappingKern; %{'rbfard2', 'bias', 'white'};
        %indPoints = 80; %%%%%
        options{i}.numActive = indPoints;
        options{i}.optimiser = 'scg2';
        if ~DgtN
            options{i}.enableDgtN = false;
        end
        % !!!!! Be careful to use the same type of scaling and bias for all
        % models!!!

        % scale = std(Ytr);
        % scale(find(scale==0)) = 1;
        %options.scaleVal = mean(std(Ytr));
        options{i}.scaleVal = sqrt(var(Ytr{i}(:)));
    end

    %-------------- INIT LATENT SPACE ---%
    for i=1:length(Ytr)
        % Compute m, the normalised version of Ytr (to be used for
        % initialisation of X)
        bias = mean(Ytr{i});
        scale = ones(1, d{i});

        if(isfield(options{i},'scale2var1'))
            if(options{i}.scale2var1)
                scale = std(Ytr{i});
                scale(find(scale==0)) = 1;
                if(isfield(options{i}, 'scaleVal'))
                    warning('Both scale2var1 and scaleVal set for GP');
                end
            end
        end
        if(isfield(options{i}, 'scaleVal'))
            scale = repmat(options{i}.scaleVal, 1, d{i});
        end

        % Remove bias and apply scale.
        m{i} = Ytr{i};
        for j = 1:d{i}
            m{i}(:, j) = m{i}(:, j) - bias(j);
            if scale(j)
                m{i}(:, j) = m{i}(:, j)/scale(j);
            end
        end

        %    mAll = [mAll m{i}]; % Concatenation (doesn't work if different sizes)
    end
    if strcmp(initial_X, 'separately')
        fprintf('# Initialising X by performing ppca in each observed (scaled) dataset separately and then concatenating...\n');
        X_init{1} = ppcaEmbed(m{1},latentDimPerModel);
        X_init{2} = ppcaEmbed(m{2},latentDimPerModel);
        X_init = [X_init{1} X_init{2}];
    else
        fprintf('# Initialising X by performing ppca in concatenated observed (scaled) data...\n');
        X_init = ppcaEmbed([m{1} m{2}], latentDimPerModel*2);
    end
    %-----------------

    latentDim = size(X_init,2);

    % Free up some memory
    clear('Y')

    %-- Create the sub-models: Assume that for each dataset we have one model.
    % This can be changed later, as long as we find a reasonable way to
    % initialise the latent spaces.
    for i=1:numberOfDatasets
        %---- Here put some code to assign X to the global common X which must
        % be created by doing pca in the concatenation of Y's...After this
        % point, model{i}.X will be the same for all i's. TODO...
        fprintf(1,'# Creating the model...\n');
        options{i}.initX = X_init;
        model{i} = vargplvmCreate(latentDim, d{i}, Ytr{i}, options{i});

        model{i}.X = X_init; %%%%%%%
        model{i} = vargplvmParamInit(model{i}, m{i}, model{i}.X);
        model{i}.X = X_init; %%%%%%%

        inpScales = invWidthMult./(((max(model{i}.X)-min(model{i}.X))).^2); % Default 5
        %inpScales(:) = max(inpScales); % Optional!!!!!
        model{i}.kern.comp{1}.inputScales = inpScales;

        if strcmp(model{i}.kern.type, 'rbfardjit')
            model{i}.kern.inputScales = model{i}.kern.comp{1}.inputScales;
        end
        params = vargplvmExtractParam(model{i});
        model{i} = vargplvmExpandParam(model{i}, params);
        model{i}.vardist.covars = 0.5*ones(size(model{i}.vardist.covars)) + 0.001*randn(size(model{i}.vardist.covars));

        model{i}.beta=1/(0.01*var(m{i}(:)));
        prunedModelInit{i} = vargplvmPruneModel(model{i});
        %disp(model{i}.vardist.covars)
    end

    %modelInit = model;%%%TEMP

    %--  Unify models into a structure
    svargplvm_init
    model = svargplvmModelCreate(model);
    model.dataSetNames = dataSetNames;
    model.experimentNo = experimentNo;
    model.dataType = dataType;
    %%---
    capName = dataType;
    capName(1) = upper(capName(1));
    modelType = model.type;
    modelType(1) = upper(modelType(1));
    fileToSave = ['dem' capName modelType num2str(experimentNo) '.mat'];
    %%---

    %-- Define what level of parallelism to use (w.r.t submodels or/and w.r.t
    % datapoints).
    %{
fprintf('# Parallel computations w.r.t the submodels!\n');
model.parallel = 1;
model = svargplvmPropagateField(model,'parallel', 1);
%
fprintf('# Parallel computations w.r.t the datapoints!\n');
model.vardist.parallel = 1;
for i=1:model.numModels
    model.comp{i}.vardist.parallel = 1;
end
    %}

    % Force kernel computations
    params = svargplvmExtractParam(model);
    model = svargplvmExpandParam(model, params);

    %%

    display = 1;
    %%%% Optimisation
    % do not learn beta and sigma_f for few iterations for intitialization
    if initVardistIters ~=0
        model.initVardist = 1; model.learnSigmaf = 0;
        model = svargplvmPropagateField(model,'initVardist', model.initVardist);
        model = svargplvmPropagateField(model,'learnSigmaf', model.learnSigmaf);
        fprintf(1,'# Intitiliazing the variational distribution for %s iters...\n',num2str(initVardistIters));
        model = svargplvmOptimise(model, display, initVardistIters); % Default: 20
        %fprintf(1,'1/b = %.4d\n',1/model.beta);

        modelInitVardist = model;
        model.initVardistIters=initVardistIters;
    end

    model.initVardist = 0; model.learnSigmaf=1;
    model = svargplvmPropagateField(model,'initVardist', model.initVardist);
    model = svargplvmPropagateField(model,'learnSigmaf', model.learnSigmaf);

    % Optimise the model.
    model.iters = 0;
    for i=1:length(itNo)
        iters = itNo(i); % default: 2000
        fprintf(1,'\n# Optimising the model for %d iterations (session %d)...\n',iters,i);
        model = svargplvmOptimise(model, display, iters);
        model.iters = model.iters + iters;
        % Save model
        prunedModel = svargplvmPruneModel(model);
        fprintf(1,'# Saving %s\n',fileToSave);
        save(fileToSave, 'prunedModel', 'prunedModelInit');
    end
else
    capName = dataType;
    capName(1) = upper(capName(1));
    modelType = 'svargplvm';
    modelType(1) = upper(modelType(1));
    fileToSave = ['dem' capName modelType num2str(experimentNo) '.mat'];
    load(fileToSave);
    model = svargplvmRestorePrunedModel(prunedModel, Ytr);clear('prunedModel');clear('prunedModelInit');
end

%%

if ~exist('resultsDynamic')  resultsDynamic = 0; end

for i=1:length(model.comp)
    if model.comp{i}.DgtN
        model.comp{i}.m = model.comp{i}.mOrig;
    end
end
for i=1:length(model.comp)
    model.comp{i}.vis.startDim = {1,2};
    model.comp{i}.vis.startPos = model.vardist.means(1,:);
end

%%
% v = 1;
% modelVis = model.comp{v};
% if resultsDynamic
%     %bar(model.comp{v}.kern.comp{1}.inputScales);
%     %figure
%     % The following causes OUTOFMEMORY exception except from when we prune the
%     % video dimensions: (No2te: also, lvmVisualise was changed a bit so that no
%     % dynamic slides is presented, because otherwise a strange error occurs)..
%     modelVis.y = Ytr{v};
%     reduction = 1; % 4
%     opt.showVariance = 1;    opt.showInducing =1 ;
%     %opt.showVariance=0; opt.showInducing=0;
%     [modelP, newHeight, newWidth] = vargplvmReduceVidModel(modelVis, height, width, reduction, reduction);
%     lvmVisualiseGeneral(modelP, [], 'imageMRDVisualise', 'imageMRDModify', opt, [newHeight newWidth],0,0,1);
%     clear modelP
%     figure;svargplvmShowScales(model);
% end

%%
if resultsDynamic
    if exist('modelOrig','var'), model = modelOrig; end
    modelOrig = model;
    clear global visualiseInfo
    for v=1:model.numModels
        model.comp{v}.y = Ytr{v};
        reduction = 2; % 4
        opt.showVariance = 1;    opt.showInducing =1 ;
        %opt.showVariance=0; opt.showInducing=0;
        [model.comp{v}, newHeight, newWidth] = vargplvmReduceVidModel(model.comp{v}, height, width, reduction, reduction);
    end
    lvmVisualiseMRD(model,[],{'imageMRDVisualise','imageMRDVisualise'},{'imageMRDModify','imageMRDModify'},opt,{{[newHeight newWidth],0,0,1},{[newHeight newWidth],0,0,1}});
    clear modelP
    figure;svargplvmShowScales(model);
    set(gca, 'FontSize', 18);
end

%% ---
%%{
axFs = 18;
titleFs = 0;
legFs = 20;

figure
sc = svargplvmShowScales(model,0);
maxScales1 = max(sc{1});
maxScales2 = max(sc{2});
sc{1} = sc{1}./maxScales1;
sc{2} = sc{2}./maxScales2;
sc{1}=sigmoid(sc{1}*5)-0.5;
sc{2}=sigmoid(sc{2}*5)-0.5;

x=1:size(sc{1},2);
K=0.75;
bar1=bar(x, sc{1}, 'FaceColor', 'b', 'EdgeColor', 'b');
set(gca, 'YtickLabel',[])
set(bar1,'BarWidth',K);
hold on;
bar2=bar(x, sc{2}, 'FaceColor', 'r', 'EdgeColor', 'r');
set(gca, 'YtickLabel',[])
set(bar2,'BarWidth',K/2.5);
hold off;
%lg=legend(['scales1'],['scales2']);
%set(lg, 'FontSize',legFs)
set(gca, 'FontSize', axFs);
%%}

%% Eigenfaces
%{
m=1;
Kuu=model.comp{m}.K_uu;
beta = model.comp{m}.beta;
Psi2 = model.comp{m}.Psi2;
Psi1 = model.comp{m}.Psi1;
h=newHeight;
w=newWidth;
mu = zeros(size(Kuu,1),h*w);
for i=1:size(Y,2)
    mu(:,i) = Kuu*pdinv(1/beta*Kuu+Psi2)*Psi1'*model.comp{m}.m(:,i);
end
x = mu';
x=bsxfun(@minus, x', mean(x'))';

% calculate covariance
s = cov(x');
% obtain eigenvalue & eigenvector
[V,D] = eig(s);
eigval = diag(D);
% sort eigenvalues in descending order
eigval = eigval(end:-1:1);
V = fliplr(V);
% show 0th through 15th principal eigenvectors
eig0 = reshape(mean(x,2), [h,w]);
figure,subplot(4,4,1)
imagesc(eig0)
colormap gray
for i = 1:15
    subplot(4,4,i+1)
    imagesc(reshape(V(:,i),h,w))
end

%
figure;imagesc(reshape(V(:,3), h,w)); axis off; colormap('gray'); axis('image')
figure;imagesc(reshape(V(:,5), h,w)); axis off; colormap('gray'); axis('image')
figure;imagesc(reshape(V(:,6), h,w)); axis off; colormap('gray'); axis('image')
figure;imagesc(reshape(V(:,3), h,w)); axis off; colormap('gray'); axis('image');colormap(flipud(colormap));
figure;imagesc(reshape(V(:,5), h,w)); axis off; colormap('gray'); axis('image');colormap(flipud(colormap));
figure;imagesc(reshape(V(:,6), h,w)); axis off; colormap('gray'); axis('image');colormap(flipud(colormap));
%}

%%

%---------------------------- PREDICTIONS ---------------
if ~doPredictions
    return
end

obsMod = 1; % one of the involved sub-models (the one for which we have the data)
infMod = setdiff(1:2, obsMod);

% Find the dimensions that are shared for obsMod and infMod
if ~exist('sharedDims')
    s1 = model.comp{obsMod}.kern.comp{1}.inputScales;
    s2 = model.comp{infMod}.kern.comp{1}.inputScales;
    % Normalise values between 0 and 1
    s1 = s1 / max(s1);
    s2 = s2 / max(s2);

    %  thresh = max(model.comp{obsMod}.kern.comp{1}.inputScales) * 0.001;
    thresh = 0.005;

    retainedScales{obsMod} = find(s1 > thresh);
    %thresh = max(model.comp{infMod}.kern.comp{1}.inputScales) * 0.001;
    retainedScales{infMod} = find(s2  > thresh);
    sharedDims = intersect(retainedScales{obsMod}, retainedScales{infMod});
end

% Find X_* only for the shared dimensions (Xs*):
if ~exist('privateDims')
    privateDims = setdiff(1:model.comp{obsMod}.q, sharedDims);
end

if ~exist('testOnTraining')
    testOnTraining=1;
end

numberTestPoints = 10;
if testOnTraining
    perm = randperm(model.N);
    testInd = perm(1:numberTestPoints);
else
    perm = randperm(size(Yts{obsMod},1));
    testInd = perm(1:numberTestPoints);
end

scrsz = get(0,'ScreenSize');

for i=1:length(testInd)
    curInd = testInd(i);
    fprintf('# Testing indice number %d ', curInd);
    if testOnTraining
        fprintf('taken from the training set\n');
        y_star = model.comp{obsMod}.y(curInd,:);
        x_star = model.comp{obsMod}.vardist.means(curInd,:);
        varx_star = model.comp{obsMod}.vardist.covars(curInd,:);
    else
        fprintf('taken from the test set\n');
        y_star = Yts{obsMod}(curInd,:);
        z_star = Yts{infMod}(curInd,:);
        dst = dist2(y_star, model.comp{obsMod}.y);
        [mind, mini] = min(dst);

        Init(i,:) = model.vardist.means(mini,:);
        vardistx = vardistCreate(model.comp{obsMod}.vardist.means(mini,:), model.q, 'gaussian');
        vardistx.covars = model.comp{obsMod}.vardist.covars(mini,:);
        model.comp{obsMod}.vardistx = vardistx;
        display=1;
        iters = 250;
        % Find p(X_* | Y_*) which is approximated by q(X_*)
        [x_star, varx_star, modelUpdated] = vargplvmOptimisePoint(model.comp{obsMod}, vardistx, y_star, display, iters);%%%
    end
    numberOfNN = 9;
    % Now we selected a datapoint X_* by taking into account only the
    % private dimensions for Y. Now, based on the shared dimensions of
    % that, we select the closest (in a NN manner) X from the training data.
    fprintf('# Finding the %d NN of X_* with the training X based only on the shared dims.\n', numberOfNN);
    [ind, distInd] = nn_class(model.X(:,sharedDims), x_star(:,sharedDims), numberOfNN, 'euclidean');

    ZpredMu = zeros(length(ind), size(model.comp{infMod}.y,2));
    ZpredSigma = zeros(length(ind), size(model.comp{infMod}.y,2));

    % Find p(y_*|x_*) for every x_* found from the NN
    fprintf('# Predicting images from the NN of X_* ');
    for k=1:numberOfNN
        fprintf('.');
        x_cur = model.X(ind(k),:);
        %x_cur(sharedDims) = x_star(sharedDims); %%% OPTIONAL!!!
        %[ZpredMu(k,:), ZpredSigma(k,:)] = vargplvmPosteriorMeanVar(model.comp{infMod}, model.X(ind(k),:));
        ZpredMu(k,:) = vargplvmPosteriorMeanVar(model.comp{infMod}, x_cur);
    end
    fprintf('\n\n');

    %-- Plots
    % Open a big figure (first 2 args control the position, last 2 control
    % the size)
    figure('Position',[scrsz(3)/100.86 scrsz(4)/6.666 scrsz(3)/1.0457 scrsz(4)/1.0682],...
        'Name',['Fig: ' num2str(i) ' (Exp: ' num2str(experimentNo) ')'],'NumberTitle','off')
    numRows = 2;

    if testOnTraining
        numCols = ceil((numberOfNN+1)/numRows);
        plotCounter = 1;
    else
        % For the real test image!
        numCols = ceil((numberOfNN+2)/numRows);
        plotCounter = 2;
    end
    subplot(numRows, numCols, 1)
    imagesc(reshape(y_star,height,width)), title(['Original y (image #' num2str(curInd) ')']), colormap('gray')

    if ~testOnTraining
        subplot(numRows, numCols, 2)
        imagesc(reshape(z_star,height,width)), title(['Corresponding z (image #' num2str(curInd) ')']), colormap('gray')
    end

    for k=1:numberOfNN
        subplot(numRows, numCols, k+plotCounter)
        imagesc(reshape(ZpredMu(k,:), height, width)), title(['NN #' num2str(k)]), colormap('gray');
    end

    %{
    indYnn = [];
    % Do a NN on the DATA space, for every predicted output.
    for j=1:size(ZpredMu,1)
        [indYnn(j), distInd] = nn_class(model.comp{infMod}.y, ZpredMu(j,:),1,'euclidean');
    end
    for j=1:length(indYnn)
        figure, imagesc(reshape(model.comp{infMod}.y(indYnn(j),:),height, width)), title([num2str(j)]), colormap('gray')
    end
    %}
end

if ~testOnTraining
    errsumFull = sum((ZpredMu - Yts).^2);
    errorFull = mean(errsumFull);
end

repo: vargplvm