text:
%CellState
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last Updated: 9/12/2010
classdef CellState
    %printing
    methods (Static)
        function run(sim, fileName)
            import edu.stanford.covert.cell.sim.analysis.CellState;
            import edu.stanford.covert.cell.sim.util.PlotUtil;
            import edu.stanford.covert.cell.sim.util.PrintUtil;

            %excel file
            [content, colLabels, indentation] = CellState.state(sim);
            if nargin == 1
                PrintUtil.printToStdIO(content, colLabels, struct('indentation', indentation));
            else
                PrintUtil.printToFile(content, colLabels, [fileName '.xls'], 'State', struct('indentation', indentation));
            end

            [content, colLabels, indentation] = CellState.processes(sim);
            if nargin == 1
                PrintUtil.printToStdIO(content, colLabels, struct('indentation', indentation));
            else
                PrintUtil.printToFile(content, colLabels, [fileName '.xls'], 'Process', struct('indentation', indentation));
            end

            % plots
            if nargin == 1
                CellState.plotStateOnChromosome(sim, PlotUtil.newAxesHandle());
            else
                [axesHandle, figHandle] = PlotUtil.newAxesHandle();

                cla(axesHandle);
                CellState.plotStateOnChromosome(sim, axesHandle);
                saveas(figHandle, [fileName '-OnChromosome.pdf']);

                close(figHandle);
            end
        end

        function [content, colLabels, indentation] = state(sim)
            %import classes
            import edu.stanford.covert.util.ConstantUtil;
            import edu.stanford.covert.cell.sim.Simulation;

            content = cell(0, 4);
            colLabels = {'Quantity', 'Value', 'Units'};

            %% time
            time = sim.state('Time');
            content = [content;{
                0 'Time' time.values(1)/ConstantUtil.secondsPerHour 'h';
                }];

            %% mass
            mass = sim.state('Mass');
            cellDryMass = sum(mass.cellDry);
            content = [content;{
                0 'Mass' sum(mass.cell) 'g';
                    1 'Dry' sum(mass.cellDry) 'g';
                    1 'Metabolite' sum(mass.metaboliteWt) 'g';
                        2 'Percent Dry' sum(mass.metaboliteWt) / cellDryMass * 100 '%';
                    1 'DNA' sum(mass.dnaWt) 'g';
                        2 'Percent Dry' sum(mass.dnaWt) / cellDryMass * 100 '%';
                    1 'RNA' sum(mass.rnaWt) 'g';
                        2 'Percent Dry' sum(mass.rnaWt) / cellDryMass * 100 '%';
                    1 'Protein' sum(mass.proteinWt) 'g';
                        2 'Percent Dry' sum(mass.proteinWt) / cellDryMass * 100 '%';
                    1 'No. cells' sum(mass.cellDry)/mass.cellInitialDryWeight 'cells';
                }];

            %% growth
            s = sim.state('MetabolicReaction');
            time = sim.state('Time');
            content = [content;{
                0 'Growth' s.growth 'cells/s';
                    1 'Relative to expected rate' s.growth / (1 / time.cellCycleLength * log(2) * exp(time.values / time.cellCycleLength * log(2))) [];
                    1 'Doubling Time', s.doublingTime / ConstantUtil.secondsPerHour 'h';
                }];

            %% shape
            s = sim.state('Geometry');
            content = [content;{
                0 'Shape' 'Cyclinder with hemispherical caps' []
                    1 'Width' s.width 'm'
                    1 'Pinched Diameter' s.pinchedDiameter 'm'
                    1 'Cylindrical length' s.cylindricalLength 'm'
                    1 'Total length' s.totalLength 'm'
                    1 'Pinched circumference' s.pinchedCircumference 'm'
                    1 'Volume' s.volume 'L'
                    1 'Surface area' s.surfaceArea' 'm^2'
                    1 'Pinched' s.pinched 'true/false'
                }];

            %% FtsZ ring
            s = sim.state('FtsZRing');
            content = [content;{
                0 'FtsZ ring' [] []
                    1 'No. edges' s.numEdges []
                    1 'No. edges bound by 1 straight FtsZ 9-mer' s.numEdgesOneStraight []
                    1 'No. edges bound by 2 straight FtsZ 9-mers' s.numEdgesTwoStraight []
                    1 'No. edges bound by 2 bent FtsZ 9-mers' s.numEdgesTwoBent []
                    1 'No. edges bound by 1 bent FtsZ 9-mer' s.numResidualBent []
                }];

            %% stimuli
            content = [content;{
                0 'Stimuli' '' []
                }];
            stim = sim.state('Stimulus');
            for i = 1:numel(stim.wholeCellModelIDs)
                value = stim.values(i, stim.values(i, :) ~= 0);
                if numel(value) == 0
                    value = 0;
                elseif numel(value) > 1
                    value = mean(value);
                end
                content = [content;{
                    1 stim.wholeCellModelIDs{i} value []
                    }]; %#ok<AGROW>
            end

            %% metabolites
            c = sim.state('Chromosome');
            s = sim.state('Metabolite');

            freedNTPs = s.counts(s.dntpIndexs, sim.compartment.cytosolIndexs);
            freeNTPs = s.counts(s.ntpIndexs, sim.compartment.cytosolIndexs);
            freeNDPs = s.counts(s.ndpIndexs, sim.compartment.cytosolIndexs);
            freeNMPs = s.counts(s.nmpIndexs, sim.compartment.cytosolIndexs);
            freeAAs = s.counts(s.aminoAcidIndexs, sim.compartment.cytosolIndexs);

            [mindNTPVal, mindNTPIdx] = min(freedNTPs);
            [minNTPVal, minNTPIdx] = min(freeNTPs);
            [minNDPVal, minNDPIdx] = min(freeNDPs);
            [minNMPVal, minNMPIdx] = min(freeNMPs);
            [minAAVal, minAAIdx] = min(freeAAs);

            freeNTPMass = freeNTPs' * s.molecularWeights(s.ntpIndexs) / ConstantUtil.nAvogadro;
            freeAAMass = freeAAs' * s.molecularWeights(s.aminoAcidIndexs) / ConstantUtil.nAvogadro;

            content = [content;{
                0 'Metabolites' [] []
                    1 'Free' sum(s.counts(:)) 'molecules'
                    1 'dNTPs' [] []
                        2 'Free' sum(freedNTPs) 'molecules'
                        2 'Min free' s.wholeCellModelIDs{s.dntpIndexs(mindNTPIdx)} []
                        2 'Min free' mindNTPVal 'molecules'
                        2 'Incorporated into DNA' (1 - sum(freedNTPs) / (sum(freedNTPs) + collapse(c.polymerizedRegions)))*100 '%'
                    1 'NTPs' [] []
                        2 'Free' sum(freeNTPs) 'molecules'
                        2 'Min free' s.wholeCellModelIDs{s.ntpIndexs(minNTPIdx)} []
                        2 'Min free' minNTPVal 'molecules'
                        2 'Incorporated into RNA' (1 - freeNTPMass / (freeNTPMass + sum(sim.state('Mass').rnaWt)))*100 '%'
                    1 'NDPs' [] []
                        2 'Free' sum(freeNDPs) 'molecules'
                        2 'Min free' s.wholeCellModelIDs{s.ndpIndexs(minNDPIdx)} []
                        2 'Min free' minNDPVal 'molecules'
                    1 'NMPs' [] []
                        2 'Free' sum(freeNMPs) 'molecules'
                        2 'Min free' s.wholeCellModelIDs{s.nmpIndexs(minNMPIdx)} []
                        2 'Min free' minNMPVal 'molecules'
                    1 'Amino acids' [] []
                        2 'Free' sum(freeAAs) 'molecules'
                        2 'Min free' s.wholeCellModelIDs{s.aminoAcidIndexs(minAAIdx)} []
                        2 'Min free' minAAVal 'molecules'
                        2 'Incorporated into Protein' (1 - freeAAMass / (freeAAMass + sum(sim.state('Mass').proteinWt)))*100 '%'
                    1 'Energy' [] []
                        2 'ATP' sum(s.counts(s.getIndexs({'ATP'}), :)) 'molecules'
                        2 'ADP' sum(s.counts(s.getIndexs({'ADP'}), :)) 'molecules'
                        2 'AMP' sum(s.counts(s.getIndexs({'AMP'}), :)) 'molecules'
                        2 'PPi' sum(s.counts(s.getIndexs({'PPI'}), :)) 'molecules'
                        2 'Pi'  sum(s.counts(s.getIndexs({'PI'}),  :)) 'molecules'
                        2 'H2O' sum(s.counts(s.getIndexs({'H2O'}), :)) 'molecules'
                        2 'H'   sum(s.counts(s.getIndexs({'H'}),   :)) 'molecules'
                }];

            %% DNA
            c = sim.state('Chromosome');
            pm = sim.state('ProteinMonomer');
            pc = sim.state('ProteinComplex');
            ri = sim.process('ReplicationInitiation');
            dr = sim.process('DNARepair');
            dnaABoxBound = ...
                ri.calculateDnaABoxStatus() == ri.dnaABoxStatus_DnaAATPBound | ...
                ri.calculateDnaABoxStatus() == ri.dnaABoxStatus_DnaAADPBound;
            [polATPs, polADPs] = ri.calculateDnaAR1234Polymerization();
            oriCComplexSize = ri.calculateDnaAR1234ComplexSize(polATPs, polADPs);
            [~, boundMonomers] = find(c.monomerBoundSites);
            [~, boundComplexs] = find(c.complexBoundSites);
            occupiedBases = sum(c.monomerDNAFootprints(boundMonomers) + c.complexDNAFootprints(boundComplexs));

            methylatedRMSites = nnz(c.damagedBases([
                dr.RM_MunI_RecognitionSites(:, dr.RM_MunI_MethylatedPositions(1))   ones(size(dr.RM_MunI_RecognitionSites, 1), 1);
                dr.RM_MunI_RecognitionSites(:, dr.RM_MunI_MethylatedPositions(2)) 2*ones(size(dr.RM_MunI_RecognitionSites, 1), 1);
                dr.RM_MunI_RecognitionSites(:, dr.RM_MunI_MethylatedPositions(1)) 3*ones(size(dr.RM_MunI_RecognitionSites, 1), 1);
                dr.RM_MunI_RecognitionSites(:, dr.RM_MunI_MethylatedPositions(2)) 4*ones(size(dr.RM_MunI_RecognitionSites, 1), 1);
                ])) / 2;
            polymerizedRMSites = sum(c.isRegionPolymerized([
                dr.RM_MunI_RecognitionSites(:, 1)   ones(size(dr.RM_MunI_RecognitionSites, 1), 1);
                dr.RM_MunI_RecognitionSites(:, 1) 2*ones(size(dr.RM_MunI_RecognitionSites, 1), 1);
                dr.RM_MunI_RecognitionSites(:, 1) 3*ones(size(dr.RM_MunI_RecognitionSites, 1), 1);
                dr.RM_MunI_RecognitionSites(:, 1) 4*ones(size(dr.RM_MunI_RecognitionSites, 1), 1);
                ], 6, false)) / 2;

            content = [content;{
                0 'No. chromosomes' collapse(c.polymerizedRegions)/c.sequenceLen [];
                    1 'Bases' collapse(c.polymerizedRegions) 'nt';
                    1 'Linking number' collapse(c.linkingNumbers)/2 'Links';
                    1 'Superhelical density', (collapse(c.linkingNumbers)  - collapse(c.linkingNumbers_minFreeEnergy)) / collapse(c.linkingNumbers_minFreeEnergy) [];
                    1 'Bound proteins' nnz(c.monomerBoundSites) + nnz(c.complexBoundSites) 'molecules'
                        2 'Monomers' nnz(c.monomerBoundSites) 'molecules'
                            3 'GntR' nnz(c.monomerBoundSites == pm.getIndexs({'MG_101_MONOMER'})) 'molecules'
                        2 'Complexes' nnz(c.monomerBoundSites) 'molecules'
                            3 'DnaA' nnz(dnaABoxBound) 'molecules'
                                4 '7-mer' nnz(dnaABoxBound(ri.dnaABoxIndexs_7mer, :)) 'molecules'
                                4 '8-mer' nnz(dnaABoxBound(ri.dnaABoxIndexs_8mer, :)) 'molecules'
                                4 '9-mer' nnz(dnaABoxBound(ri.dnaABoxIndexs_9mer, :)) 'molecules'
                                4 'OriC Complex' oriCComplexSize(1) 'molecules'
                                4 'R1' nnz(dnaABoxBound(ri.dnaABoxIndexs_R12345(1), :)) 'molecules'
                                4 'R2' nnz(dnaABoxBound(ri.dnaABoxIndexs_R12345(2), :)) 'molecules'
                                4 'R3' nnz(dnaABoxBound(ri.dnaABoxIndexs_R12345(3), :)) 'molecules'
                                4 'R4' nnz(dnaABoxBound(ri.dnaABoxIndexs_R12345(3), :)) 'molecules'
                                4 'R5' nnz(dnaABoxBound(ri.dnaABoxIndexs_R12345(5), :)) 'molecules'
                            3 'Gyrase' nnz(c.complexBoundSites == pc.getIndexs({'DNA_GYRASE'})) 'molecules'
                            3 'LuxR' nnz(c.complexBoundSites == pc.getIndexs({'MG_428_DIMER'})) 'molecules'
                            3 'SMC' nnz(c.complexBoundSites == pc.getIndexs({'MG_213_214_298_6MER_ADP'})) 'molecules'
                                4 'Average spacing' collapse(c.polymerizedRegions)/2 / nnz(c.complexBoundSites == pc.getIndexs({'MG_213_214_298_6MER_ADP'})) 'nt'
                        2 'Occupied' occupiedBases 'nt'
                        2 'Occupied' occupiedBases / collapse(c.polymerizedRegions) / 2 * 100 '%'
                    1 'Damage' [] []
                        2 'Gap sites' nnz(c.gapSites) 'nt'
                        2 'Abasic sites' nnz(c.abasicSites) 'nt'
                        2 'Damaged bases' nnz(c.damagedBases) 'nt'
                            3 'm6AD' nnz(c.damagedBases == c.metabolite.m6ADIndexs) 'nt'
                            3 'methylated R/M sites' methylatedRMSites 'nt'
                            3 'methylated R/M sites' methylatedRMSites/polymerizedRMSites*100 '%'
                        2 'Damaged sugar-phosphate' nnz(c.damagedSugarPhosphates) 'nt'
                        2 'Strand breaks' nnz(c.strandBreaks) 'nt'
                        2 'Intrastrand crosslinks' nnz(c.intrastrandCrossLinks) 'nt'
                        2 'Holliday junctions' nnz(c.hollidayJunctions) 'nt'
                }];

            %% RNA
            r = sim.state('Rna');

            matureRNAs = sum(r.counts(r.matureIndexs, :), 2);
            aminoacylatedRNAs = sum(r.counts(r.aminoacylatedIndexs, :), 2);
            content = [content;{
                0 'RNA' sum(r.counts(:)) 'molecules'
                    1 'mRNA' sum(matureRNAs(r.matureMRNAIndexs)) 'molecules'
                        2 'No. species expressed' nnz(r.counts(r.matureIndexs(r.matureMRNAIndexs), :)) []
                    1 'rRNA' sum(matureRNAs(r.matureRRNAIndexs)) 'molecules'
                        2 '5S'  matureRNAs(r.matureRibosomalRRNAIndexs(1)) 'molecules'
                        2 '15S' matureRNAs(r.matureRibosomalRRNAIndexs(2)) 'molecules'
                        2 '23S' matureRNAs(r.matureRibosomalRRNAIndexs(3)) 'molecules'
                    1 'sRNA' sum(matureRNAs(r.matureSRNAIndexs) + aminoacylatedRNAs(r.matureSRNAIndexs)) 'molecules'
                        2 'scRNA'   matureRNAs(r.matureSRNAIndexs(1))+aminoacylatedRNAs(r.matureSRNAIndexs(1)) 'molecules'
                        2 'MG170'   matureRNAs(r.matureSRNAIndexs(2))+aminoacylatedRNAs(r.matureSRNAIndexs(2)) 'molecules'
                        2 'RNase P' matureRNAs(r.matureSRNAIndexs(3))+aminoacylatedRNAs(r.matureSRNAIndexs(3)) 'molecules'
                        2 'tmRNA'   matureRNAs(r.matureSRNAIndexs(4))+aminoacylatedRNAs(r.matureSRNAIndexs(4)) 'molecules'
                            3 'Free' matureRNAs(r.matureSRNAIndexs(4)) 'molecules'
                            4 'Aminoacylated' aminoacylatedRNAs(r.matureSRNAIndexs(4)) 'molecules'
                    1 'tRNA' sum(matureRNAs(r.matureTRNAIndexs) + aminoacylatedRNAs(r.matureTRNAIndexs)) 'molecules'
                        2 'Free' sum(matureRNAs(r.matureTRNAIndexs)) 'molecules'
                        2 'Aminoacylated' sum(aminoacylatedRNAs(r.matureTRNAIndexs)) 'molecules'
                        2 'No. species expressed' nnz(matureRNAs(r.matureTRNAIndexs) + aminoacylatedRNAs(r.matureTRNAIndexs)) 'molecules'
                    1 'Immature' nnz(sim.process('Transcription').transcripts.boundTranscriptionUnits) 'molecules'
                    1 'Aborted' size(sim.state('Transcript').abortedTranscripts, 1) 'molecules'
                }];

            %% Protein
            m = sim.state('ProteinMonomer');
            c = sim.state('ProteinComplex');

            ribosomalMonomerIdxs = sum(c.proteinComplexComposition(sim.gene.mRNAIndexs, c.getIndexs({'RIBOSOME_70S'}), :),3);
            ribosomalMonomers = sum(m.counts(m.matureIndexs(ribosomalMonomerIdxs > 0), :), 2);
            [minRibosomalMonomers, minRibosomalMonomerIdx] = min(ribosomalMonomers);
            minRibosomalMonomerID = m.wholeCellModelIDs{m.matureIndexs(minRibosomalMonomerIdx)};

            [idxs, ~, counts] = find(sum(m.counts(m.inactivatedIndexs, :), 2));
            inactiveMonomers = cell(numel(idxs), 4);
            inactiveMonomers(:, 1) = {2};
            inactiveMonomers(:, 2) = m.wholeCellModelIDs(m.inactivatedIndexs(idxs));
            inactiveMonomers(:, 3) = num2cell(counts);

            [idxs, ~, counts] = find(sum(c.counts(c.inactivatedIndexs, :), 2));
            inactiveComplexs = cell(numel(idxs), 4);
            inactiveComplexs(:, 1) = {2};
            inactiveComplexs(:, 2) = c.wholeCellModelIDs(c.inactivatedIndexs(idxs));
            inactiveComplexs(:, 3) = reshape(num2cell(counts), [], 1);

            content = [content;{
                0 'Protein' sum(m.counts(:)) + sum(c.counts(:)) 'molecules'
                    1 'Monomers' sum(m.counts(:)) 'molecules'
                        2 'Ribosomal proteins' sum(ribosomalMonomers) 'molecules'
                            3 'Min expressed species' minRibosomalMonomers 'molecules'
                            3 'Min expressed species' minRibosomalMonomerID []
                    1 'Complexs' sum(c.counts(:)) 'molecules'
                    1 'Bound' sum(sum(m.counts(m.boundIndexs, :)))+sum(sum(c.counts(c.boundIndexs, :))) 'molecules'
                        2 'Chromosome' nnz(sim.state('Chromosome').monomerBoundSites)+nnz(sim.state('Chromosome').complexBoundSites) 'molecules'
                        2 'RNA polymerase' nnz(sim.process('Transcription').transcripts.boundTranscriptionUnits) 'molecules'
                        2 'Ribosome' nnz(sim.state('Ribosome').boundMRNAs) 'molecules'
                    1 'Immature' nnz(sim.state('Ribosome').boundMRNAs) 'molecules'
                    1 'Inactive' sum(sum(m.counts(m.inactivatedIndexs, :)))+sum(sum(c.counts(c.inactivatedIndexs, :))) 'molecules'
                        };
                        inactiveMonomers;
                        inactiveComplexs;
                        {
                    1 'Damaged' sum(sum(m.counts(m.damagedIndexs, :)))+sum(sum(c.counts(c.damagedIndexs, :))) 'molecules'
                    1 'Aborted' size(sim.state('Polypeptide').abortedPolypeptides, 1) 'molecules'
                    1 'Cytosol' sum(m.counts(:, sim.compartment.cytosolIndexs))+ sum(c.counts(:, sim.compartment.cytosolIndexs)) 'molecules'
                    1 'Membrane' sum(m.counts(:, sim.compartment.membraneIndexs))+ sum(c.counts(:, sim.compartment.membraneIndexs)) 'molecules'
                    1 'Term. Org. Cytosol' sum(m.counts(:, sim.compartment.terminalOrganelleCytosolIndexs))+ sum(c.counts(:, sim.compartment.terminalOrganelleCytosolIndexs)) 'molecules'
                    1 'Term. Org. Membrane' sum(m.counts(:, sim.compartment.terminalOrganelleMembraneIndexs))+ sum(c.counts(:, sim.compartment.terminalOrganelleMembraneIndexs)) 'molecules'
                    1 'Extracellular' sum(m.counts(:, sim.compartment.extracellularIndexs))+ sum(c.counts(:, sim.compartment.extracellularIndexs)) 'molecules'
                }];

            %% DNA Polymerase
            r = sim.process('Replication');
            free = r.enzymes(r.enzymeIndexs_core) + r.enzymeComposition(r.enzymeIndexs_core, :) * r.enzymes;
            bound = r.boundEnzymes(r.enzymeIndexs_core) + r.enzymeComposition(r.enzymeIndexs_core, :) * r.boundEnzymes;
            content = [content;{
                0 'DNA polymerase' free+bound 'molecules'
                    1 'Free' free 'molecules'
                    1 'Bound' bound 'molecules'
                }];

            %% Transcription enzymes
            t = sim.process('Transcription');
            s = sim.state('RNAPolymerase');
            rna = t.rnaPolymerases;
            rnaPolymeraseStateOccupancies = s.stateOccupancies*100;
            content = [content;{
                0 'Transcription enzymes' [] []
                    1 'RNA polymerase' t.enzymes(t.enzymeIndexs_rnaPolymerase)+t.boundEnzymes(t.enzymeIndexs_rnaPolymerase) 'molecules'
                        2 'Active' rnaPolymeraseStateOccupancies(rna.activelyTranscribingIndex) ' %'
                        2 'Specifically bound' rnaPolymeraseStateOccupancies(rna.specificallyBoundIndex) ' %'
                        2 'Non-specifically bound' rnaPolymeraseStateOccupancies(rna.nonSpecificallyBoundIndex) ' %'
                        2 'Free' rnaPolymeraseStateOccupancies(rna.freeIndex) ' %'
                    1 'Initiation/Elongation factors' sum(t.enzymes(t.enzymeIndexs_transcriptionFactors)+t.boundEnzymes(t.enzymeIndexs_transcriptionFactors)) 'molecules'
                        2 'Free' sum(t.enzymes(t.enzymeIndexs_transcriptionFactors)) 'molecules'
                        2 'Bound' sum(t.boundEnzymes(t.enzymeIndexs_transcriptionFactors)) 'molecules'
                }];

            %% tRNA aminoacylation enzymes
            t = sim.process('tRNAAminoacylation');
            totEnzymes = t.enzymes + t.boundEnzymes;
            content = [content;{
                0 'tRNA aminoacylation enzymes' sum(totEnzymes) 'molecules'
                    1 'tRNA synthetases' sum(totEnzymes(t.enzymeIndexs_tRNASynthetases)) 'molecules'
                    1 'tRNA transferases' sum(totEnzymes(t.enzymeIndexs_tRNATransferases)) 'molecules'
                }];

            %% Translation enzymes
            t = sim.process('Translation');
            r = sim.state('Ribosome');
            totEnzymes = t.enzymes + t.boundEnzymes;
            content = [content;{
                0 'Translation enzymes' [] []
                    1 'Ribosome' [] []
                        2 'Ribosome 30S' totEnzymes(t.enzymeIndexs_ribosome30S) 'molecules'
                        2 'Ribosome 30S-IF3' totEnzymes(t.enzymeIndexs_ribosome30SIF3) 'molecules'
                        2 'Ribosome 50S' totEnzymes(t.enzymeIndexs_ribosome50S) 'molecules'
                        2 'Ribosome 70S' totEnzymes(t.enzymeIndexs_ribosome70S) 'molecules'
                            3 'Free' t.enzymes(t.enzymeIndexs_ribosome70S) 'molecules'
                            3 'Bound' t.boundEnzymes(t.enzymeIndexs_ribosome70S) 'molecules'
                            3 'Active' nnz(r.boundMRNAs) 'molecules'
                    1 'Translation factors' sum(totEnzymes(t.enzymeIndexs_translationFactors)) 'molecules'
                        1 'Free' sum(t.enzymes(t.enzymeIndexs_translationFactors)) 'molecules'
                        1 'Bound' sum(t.boundEnzymes(t.enzymeIndexs_translationFactors)) 'molecules'
                }];

            %% FtsZ
            f = sim.process('FtsZPolymerization');
            totEnzymes = f.enzymes + f.boundEnzymes;
            content = [content;{
                0 'FtsZ' [1 1 1 2:9] * totEnzymes 'molecules'
                    1 'Free monomers' totEnzymes(f.enzymeIndexs_FtsZ) 'molecules'
                    1 'GDP 1-mer' totEnzymes(f.enzymeIndexs_FtsZ_GDP) 'molecules'
                    1 'GTP 1-mer' totEnzymes(f.enzymeIndexs_FtsZ_GTP) 'molecules'
                    1 'GTP 2-mer' totEnzymes(f.enzymeIndexs_FtsZ_activated(2)) 'molecules'
                    1 'GTP 3-mer' totEnzymes(f.enzymeIndexs_FtsZ_activated(3)) 'molecules'
                    1 'GTP 4-mer' totEnzymes(f.enzymeIndexs_FtsZ_activated(4)) 'molecules'
                    1 'GTP 5-mer' totEnzymes(f.enzymeIndexs_FtsZ_activated(5)) 'molecules'
                    1 'GTP 6-mer' totEnzymes(f.enzymeIndexs_FtsZ_activated(6)) 'molecules'
                    1 'GTP 7-mer' totEnzymes(f.enzymeIndexs_FtsZ_activated(7)) 'molecules'
                    1 'GTP 8-mer' totEnzymes(f.enzymeIndexs_FtsZ_activated(8)) 'molecules'
                    1 'GTP 9-mer' totEnzymes(f.enzymeIndexs_FtsZ_activated(9)) 'molecules'
                    1 'Bound 9-mer' f.boundEnzymes(f.enzymeIndexs_FtsZ_activated(9)) 'molecules'
                }];

            %% format output
            indentation = cell2mat(content(:, 1));
            content = content(:, 2:end);
        end

        %diagnostics
        function [content, colLabels, indentation] = processes(sim)
            content = cell(0, 4);
            colLabels = {'Quantity', 'Value', 'Units'};

            %% Transcription
            t = sim.process('Transcription');
            s = sim.state('RNAPolymerase');

            %building blocks
            availableNTPs = t.substrates(t.substrateIndexs_ntp);
            rnaPolymeraseMaxNTPIncorporation = s.nActive * t.rnaPolymeraseElongationRate;
            transcriptionFactors = t.enzymes(t.enzymeIndexs_transcriptionFactors);

            %limits
            transcriptionNTPLimit = min(availableNTPs, ...
                availableNTPs / sum(availableNTPs) * rnaPolymeraseMaxNTPIncorporation);
            if abs(rnaPolymeraseMaxNTPIncorporation - sum(transcriptionNTPLimit)) < 1e-6
                transcriptionLimitingFractor = 'Active RNA Polymerases';
            else
                transcriptionLimitingFractor = 'Available NMPs, NTPs';
            end

            %limiting blocks
            [minNTPs, minNTP] = min(availableNTPs);
            minNTP = t.substrateWholeCellModelIDs{t.substrateIndexs_ntp(minNTP)};

            [minTranscriptionFactors, minTranscriptionFactor] = min(transcriptionFactors);
            minTranscriptionFactor = t.enzymeWholeCellModelIDs{t.enzymeIndexs_transcriptionFactors(minTranscriptionFactor)};

            %print
            content = [content;{
                0 'Transcription' [] []
                    1 'Allocated ATP' t.substrates(t.substrateIndexs_atp) 'molecules'
                    1 'Available NTPs' sum(availableNTPs) 'molecules'
                    1 'Min NTP' minNTPs 'molecules'
                    1 'Min NTP' minNTP []
                    1 'Active RNA Polymerases' s.nActive 'molecules'
                    1 'RNA Polymerase Elongation Rate' t.rnaPolymeraseElongationRate 'nt/s'
                    1 'Max RNA Polymerase NTP Incorporation'  rnaPolymeraseMaxNTPIncorporation 'nt/s'
                    1 'Transcription factors' sum(transcriptionFactors) 'molecules'
                    1 'Min transcription factors' minTranscriptionFactors 'molecules'
                    1 'Min transcription factors' minTranscriptionFactor []
                    1 'Max NTP Incorporation' sum(transcriptionNTPLimit) 'nt/s'
                    1 'Limiting Factor' transcriptionLimitingFractor []
             }];

            %% Translation
            t = sim.process('Translation');
            r = sim.state('Ribosome');

            %building blocks
            energy = t.substrates(t.substrateIndexs_gtp);

            aminoAcids = t.substrates(t.substrateIndexs_aminoAcids);
            ribosomeMaxAAIncorporation = r.nActive * t.ribosomeElongationRate;

            tRNAs = t.freeTRNAs + t.aminoacylatedTRNAs;
            translationFactors = t.enzymes(t.enzymeIndexs_translationFactors) + t.boundEnzymes(t.enzymeIndexs_translationFactors);

            %limiting blocks
            [minAminoAcids, minAminoAcid] = min(aminoAcids);
            minAminoAcid = t.substrateWholeCellModelIDs{t.substrateIndexs_aminoAcids(minAminoAcid)};

            [minTRNAs, minTRNA] = min(tRNAs);
            minTRNA = t.freeTRNAWholeCellModelIDs{minTRNA};

            [minTranslationFactors, minTranslationFactor]=min(translationFactors);
            minTranslationFactor = t.enzymeWholeCellModelIDs{t.enzymeIndexs_translationFactors(minTranslationFactor)};

            %limits
            content = [content;{
                0 'Translation' [] []
                    1 'Allocated ATP' energy 'molecules'
                    1 'Available AAs'  sum(aminoAcids) 'molecules'
                    1 'Min Amino Acid' minAminoAcids 'molecules'
                    1 'Min Amino Acid' minAminoAcid []
                    1 'tRNAs' sum(tRNAs) 'molecules'
                    1 'Min tRNAs' minTRNAs 'molecules'
                    1 'Min tRNAs' minTRNA []
                    1 'Active Ribosomes' r.nActive 'molecules'
                    1 'Ribosome Elongation Rate' t.ribosomeElongationRate 'aa/s'
                    1 'Max Ribosome AA Incorporation' ribosomeMaxAAIncorporation 'aa/s'
                    1 'Translation factors' sum(translationFactors) 'molecules'
                    1 'Min translation factors' minTranslationFactors 'molecules'
                    1 'Min translation factors' minTranslationFactor []
                }];

            %% format output
            indentation = cell2mat(content(:, 1));
            content = content(:, 2:end);
        end
    end

    %plotting
    methods (Static)
        function plotStateOnChromosome(sim, axesHandle)
            c = sim.state('Chromosome');

            cla(axesHandle);
            hold on;
            legendHandles = [];

            %base chromosome
            for i = 1:4
                [r, xoffset] = calcRadiusXOffset(i);
                h = plot(axesHandle, xoffset + r*cos(pi/180 * (1:360)), r*sin(pi/180 * (1:360)), 'Color', [0.7 0.7 0.7], 'LineStyle', ':');
            end
            legendHandles = [legendHandles; h];

            %polymerized regions
            for i = 1:4
                [pos, lens] = find(c.polymerizedRegions(:, i));
                [r, xoffset] = calcRadiusXOffset(i);
                for j = 1:size(pos, 1)
                    theta = 2*pi * (pos(j,1) + (0:lens(j)/360:lens(j)))/c.sequenceLen;
                    h = plot(axesHandle, xoffset + r*cos(theta), r*sin(theta), 'Color', [0.7 0.7 0.7], 'LineWidth', 2);
                end
            end
            legendHandles = [legendHandles; h];

            %polymerized transcription units
            for i = 1:2
                for j = 1:numel(c.transcriptionUnitWholeCellModelIDs)
                    if ~c.polymerizedTranscriptionUnits(j, i); continue; end;
                    [r, xoffset] = calcRadiusXOffset(2*(i-1) + c.transcriptionUnitStrands(j));
                    theta = 2*pi * (c.transcriptionUnitStartCoordinates(j) + (0:c.transcriptionUnitLengths(j)/360:c.transcriptionUnitLengths(j))) / c.sequenceLen;
                    h = plot(axesHandle, xoffset + r*cos(theta), r*sin(theta), 'Color', [0 1 0], 'LineWidth', 2);
                end
            end
            legendHandles = [legendHandles; h];

            %polymerized genes
            for i = 1:2
                for j = 1:numel(sim.gene.wholeCellModelIDs)
                    if ~c.polymerizedGenes(j, i); continue; end;
                    [r, xoffset] = calcRadiusXOffset(2*(i-1) + sim.gene.strands(j));
                    theta = 2*pi * (sim.gene.startCoordinates(j) + (0:sim.gene.lengths(j)/360:sim.gene.lengths(j))) / c.sequenceLen;
                    h = plot(axesHandle, xoffset + r*cos(theta), r*sin(theta), 'Color', [0 0 1], 'LineWidth', 2);
                end
            end
            legendHandles = [legendHandles; h];

            %bound proteins
            [monomerPosStrnds, monomers] = find(c.monomerBoundSites);
            [complexPosStrnds, complexs] = find(c.complexBoundSites);
            posStrnds = [monomerPosStrnds; complexPosStrnds];
            lens = [c.monomerDNAFootprints(monomers); c.complexDNAFootprints(complexs)];
            ssStranded = [c.monomerDNAFootprintRegionStrandedness(monomers); c.complexDNAFootprintRegionStrandedness(complexs)] == c.dnaStrandedness_ssDNA;
            for i = 1:size(posStrnds, 1)
                [r, xoffset] = calcRadiusXOffset(posStrnds(i,2));
                theta = 2*pi * (posStrnds(i,1) + (0:lens(i)/360:lens(i))) / c.sequenceLen;
                h = plot(axesHandle, xoffset + r*cos(theta), r*sin(theta), 'Color', [1 1 0], 'LineWidth', 4);
                if ~ssStranded(i)
                    [r, xoffset] = calcRadiusXOffset(posStrnds(i, 2) -iseven(posStrnds(i, 2)) + isodd(posStrnds(i, 2)));
                    plot(axesHandle, xoffset + r*cos(theta), r*sin(theta), 'Color', [1 1 0], 'LineWidth', 4);
                end
            end
            legendHandles = [legendHandles; h];

            %damage
            posStrnds = find(c.getDamagedSites(1,1,1,1,1,1,1));
            for i = 1:size(posStrnds, 1)
                [r, xoffset] = calcRadiusXOffset(posStrnds(i,2));
                theta = 2*pi * posStrnds(i,1) / c.sequenceLen;
                h = plot(axesHandle, xoffset + r*cos(theta), r*sin(theta), 'Color', [1 0 0], 'LineWidth', 4);
            end
            legendHandles = [legendHandles; h];

            %expressed RNAs
            r = sim.state('Rna');
            genes = r.matureRNAGeneComposition * sum(r.counts(r.matureIndexs, :) + r.counts(r.aminoacylatedIndexs, :), 2);
            for i = 1:numel(sim.gene.wholeCellModelIDs)
                if genes(i) == 0; continue; end;

                r = 1 + [0.05;0.25] * iseven(sim.gene.strands(i)) - [0.05;0.25] * isodd(sim.gene.strands(i));
                xoffset = 0;

                theta = 2*pi * (sim.gene.startCoordinates(i) + sim.gene.lengths(i)/2) / c.sequenceLen;
                h = plot(axesHandle, xoffset + r*cos(theta), r*sin(theta), 'Color', [1 0 1], 'LineWidth', 2);
            end
            legendHandles = [legendHandles; h];

            %expressed proteins
            pm = sim.state('ProteinMonomer');
            pc = sim.state('ProteinComplex');

            monomers = ...
                sum(pm.counts(pm.matureIndexs, :) + pm.counts(pm.boundIndexs, :) + pm.counts(pm.inactivatedIndexs, :), 2) + ...
                sum(pc.proteinComplexComposition(sim.gene.mRNAIndexs, :, :), 3) * sum(pc.counts(pc.matureIndexs, :) + pc.counts(pc.boundIndexs, :) + pc.counts(pc.inactivatedIndexs, :), 2);
            for i = 1:numel(monomers)
                if monomers(i) == 0; continue; end;

                r = 1 + [0.05;0.25] * iseven(sim.gene.strands(sim.gene.mRNAIndexs(i))) - [0.05;0.25] * isodd(sim.gene.strands(sim.gene.mRNAIndexs(i)));
                xoffset = 0;

                theta = 2*pi * (sim.gene.startCoordinates(sim.gene.mRNAIndexs(i)) + sim.gene.lengths(sim.gene.mRNAIndexs(i))/2) / c.sequenceLen;
                h = plot(axesHandle, xoffset + r*cos(theta), r*sin(theta), 'Color', [0 1 1], 'LineWidth', 2);
            end
            legendHandles = [legendHandles; h];

            %ribosomes
            ribosome = sim.state('Ribosome');
            for i = 1:numel(ribosome.boundMRNAs)
                if ribosome.boundMRNAs(i) == 0; continue; end;

                r = 1 + [0.05;0.25] * iseven(sim.gene.strands(sim.gene.mRNAIndexs(ribosome.boundMRNAs(i)))) - [0.05;0.25] * isodd(sim.gene.strands(sim.gene.mRNAIndexs(ribosome.boundMRNAs(i))));
                xoffset = 0;

                theta = 2*pi * (sim.gene.startCoordinates(sim.gene.mRNAIndexs(ribosome.boundMRNAs(i))) + ribosome.states(i)) / c.sequenceLen;
                plot(axesHandle, xoffset + r*cos(theta), r*sin(theta), 'Color', [1 1 0], 'LineWidth', 2);
            end

            %legend
            legend(legendHandles, {
                'Chromosome'
                'Polymerized'
                '  Gene'
                '  TU'
                'Bnd Prot'
                'Damage'
                'Exp RNA'
                'Exp Prot'}, ...
                'Orientation', 'Horizontal', ...
                'Location', 'SouthOutside');

            %format
            xlim([-1.4 2.5+1.4]);
            ylim([-1.4 1.4]);

            set(axesHandle, 'DataAspectRatio', [1 1 1])
            set(axesHandle, 'Visible', 'off')
            set(get(axesHandle, 'Parent'), 'Color', [1 1 1]);
            position = get(gca, 'Position');
            scale = 1.35;
            set(gca, 'Position', [position(1:2)-position(3:4)*(scale-1)/2 position(3:4)*scale]);

            function [r, xoffset] = calcRadiusXOffset(strnd)
                r = 1 + 0.05 * iseven(strnd) - 0.05 * isodd(strnd);
                xoffset = 2.5 * floor((strnd-1)/2);
            end
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [data] = ft_dipolesimulation(cfg)

% FT_DIPOLESIMULATION simulates channel-level time-series data that consists of the
% the spatial distribution of the the field or potential of one or multiple dipoles.
%
% Use as
%   data = ft_dipolesimulation(cfg)
% which will return a raw data structure that resembles the output of
% FT_PREPROCESSING.
%
% The dipoles position and orientation have to be specified with
%   cfg.sourcemodel.pos        = [Rx Ry Rz] (size Nx3)
%   cfg.sourcemodel.mom        = [Qx Qy Qz] (size 3xN)
%   cfg.sourcemodel.unit       = string, can be 'mm', 'cm', 'm' (default is automatic)
%
% The timecourse of the dipole activity is given as a cell-array with one
% dipole signal per trial
%   cfg.sourcemodel.signal     = cell-array with one dipole signal per trial
% or by specifying the parameters of a sine-wave signal
%   cfg.sourcemodel.frequency  =   in Hz
%   cfg.sourcemodel.phase      =   in radians
%   cfg.sourcemodel.amplitude  =   per dipole
%
% The number of trials and the time axes of the trials can be specified by
%   cfg.fsample     = simulated sample frequency (default = 1000)
%   cfg.trllen      = length of simulated trials in seconds (default = 1)
%   cfg.numtrl      = number of simulated trials (default = 10)
%   cfg.baseline    = number (default = 0.3)
% or by
%   cfg.time        = cell-array with one time axis per trial, for example obtained from an existing dataset
%
% Random white noise can be added to the data in each trial, either by
% specifying an absolute or a relative noise level
%   cfg.relnoise    = add noise with level relative to data signal
%   cfg.absnoise    = add noise with absolute level
%   cfg.randomseed  = 'yes' or a number or vector with the seed value (default = 'yes')
%
% Optional input arguments are
%   cfg.channel    = Nx1 cell-array with selection of channels (default = 'all'),
%                    see FT_CHANNELSELECTION for details
%   cfg.dipoleunit = units for dipole amplitude (default nA*m)
%   cfg.chanunit   = Nx1 cell-array with units for the channel data
%
% Optionally, you can modify the leadfields by reducing the rank, i.e. remove the weakest orientation
%   cfg.reducerank    = 'no', or number (default = 3 for EEG, 2 for MEG)
%   cfg.backproject   = 'yes' or 'no',  determines when reducerank is applied whether the
%                       lower rank leadfield is projected back onto the original linear
%                       subspace, or not (default = 'yes')
%
% The volume conduction model of the head should be specified as
%   cfg.headmodel     = structure with volume conduction model, see FT_PREPARE_HEADMODEL
%
% The EEG or MEG sensor positions should be specified as
%   cfg.elec          = structure with electrode positions or filename, see FT_READ_SENS
%   cfg.grad          = structure with gradiometer definition or filename, see FT_READ_SENS
%
% See also FT_SOURCEANALYSIS, FT_DIPOLEFITTING, FT_TIMELOCKSIMULATION,
% FT_FREQSIMULATION, FT_CONNECTIVITYSIMULATION

% Undocumented local options
% cfg.feedback
% cfg.previous
% cfg.version

% Copyright (C) 2004, Robert Oostenveld
%
% 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$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble provenance
ft_preamble randomseed

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'channels'}); % prevent accidental typos, see issue 1729
cfg = ft_checkconfig(cfg, 'renamed',    {'elecfile', 'elec'});
cfg = ft_checkconfig(cfg, 'renamed',    {'gradfile', 'grad'});
cfg = ft_checkconfig(cfg, 'renamed',    {'optofile', 'opto'});
cfg = ft_checkconfig(cfg, 'renamed',    {'hdmfile', 'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'vol',     'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'dip',    'sourcemodel'});

% for consistency with FT_TIMELOCKSIMULUATION and FT_FREQSIMULATION
cfg = ft_checkconfig(cfg, 'createsubcfg', 'sourcemodel');
cfg = ft_checkconfig(cfg, 'renamed', {'ntrials', 'numtrl'});
cfg = ft_checkconfig(cfg, 'renamed', {'triallength', 'trllen'});

% set the defaults
cfg.sourcemodel         = ft_getopt(cfg, 'sourcemodel', []);
cfg.sourcemodel.pos     = ft_getopt(cfg.sourcemodel, 'pos', [-5 0 15]);
cfg.sourcemodel.mom     = ft_getopt(cfg.sourcemodel, 'mom', [1 0 0]');
cfg.sourcemodel.time    = ft_getopt(cfg.sourcemodel, 'time', {});
cfg.sourcemodel.signal  = ft_getopt(cfg.sourcemodel, 'signal', {});
cfg.fsample             = ft_getopt(cfg, 'fsample', 250);
cfg.relnoise            = ft_getopt(cfg, 'relnoise', 0);
cfg.absnoise            = ft_getopt(cfg, 'absnoise', 0);
cfg.feedback            = ft_getopt(cfg, 'feedback', 'text');
cfg.channel             = ft_getopt(cfg, 'channel',  'all');
cfg.dipoleunit          = ft_getopt(cfg, 'dipoleunit', 'nA*m');
cfg.chanunit            = ft_getopt(cfg, 'chanunit', {});

% collect and preprocess the electrodes/gradiometer and head model
% this will also update cfg.channel to match the electrodes/gradiometers
[headmodel, sens, cfg] = prepare_headmodel(cfg, []);

% construct the low-level options for the leadfield computation as key-value pairs, these are passed to FT_COMPUTE_LEADFIELD
leadfieldopt = {};
leadfieldopt = ft_setopt(leadfieldopt, 'reducerank',     ft_getopt(cfg, 'reducerank'));
leadfieldopt = ft_setopt(leadfieldopt, 'backproject',    ft_getopt(cfg, 'backproject'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalize',      ft_getopt(cfg, 'normalize'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalizeparam', ft_getopt(cfg, 'normalizeparam'));
leadfieldopt = ft_setopt(leadfieldopt, 'weight',         ft_getopt(cfg, 'weight'));

cfg.sourcemodel = fixdipole(cfg.sourcemodel);
Ndipoles = size(cfg.sourcemodel.pos,1);

% in case no time or signal was given, set some additional defaults
if ~isempty(cfg.sourcemodel.time) && ~isempty(cfg.sourcemodel.signal)
  assert(length(cfg.sourcemodel.signal)==length(cfg.sourcemodel.time)); % these must match
  cfg.numtrl    = length(cfg.sourcemodel.time);
  cfg.fsample   = 1/mean(diff(cfg.sourcemodel.time{1}));  % determine from time-axis
  cfg.trllen    = length(cfg.sourcemodel.time{1})/cfg.fsample;
  cfg.baseline  = -cfg.sourcemodel.time{1}(1);
elseif ~isempty(cfg.sourcemodel.time)
  cfg.numtrl    = length(cfg.sourcemodel.time);
  cfg.fsample   = 1/mean(diff(cfg.sourcemodel.time{1}));  % determine from time-axis
  cfg.trllen    = length(cfg.sourcemodel.time{1})/cfg.fsample;
  cfg.baseline  = -cfg.sourcemodel.time{1}(1);
elseif ~isempty(cfg.sourcemodel.signal)
  cfg.numtrl    = length(cfg.sourcemodel.signal);
  cfg.fsample   = ft_getopt(cfg, 'fsample', 1000);
  cfg.trllen    = length(cfg.sourcemodel.signal{1})/cfg.fsample;
  cfg.baseline  = ft_getopt(cfg, 'baseline', 0);
else
  cfg.numtrl    = ft_getopt(cfg, 'numtrl', 10);
  cfg.fsample   = ft_getopt(cfg, 'fsample', 1000);
  cfg.trllen    = ft_getopt(cfg, 'trllen', 1);
  cfg.baseline  = ft_getopt(cfg, 'baseline', 0);
end

% no signal was given, set some additional defaults
if isempty(cfg.sourcemodel.signal)
  cfg.sourcemodel.frequency   = ft_getopt(cfg, 'frequency', ones(Ndipoles,1)*10);
  cfg.sourcemodel.phase       = ft_getopt(cfg, 'phase', zeros(Ndipoles,1));
  cfg.sourcemodel.amplitude   = ft_getopt(cfg, 'amplitude', ones(Ndipoles,1));
end

if isfield(cfg.sourcemodel, 'frequency')
  % this should be a column vector
  cfg.sourcemodel.frequency = cfg.sourcemodel.frequency(:);
end

if isfield(cfg.sourcemodel, 'phase')
  % this should be a column vector
  cfg.sourcemodel.phase = cfg.sourcemodel.phase(:);
end

if ~isempty(cfg.sourcemodel.time)
  % use the user-supplied time vectors
  diptime = cfg.sourcemodel.time;
else
  % construct a time axis for every trial
  nsample = round(cfg.trllen*cfg.fsample);
  diptime = cell(1, cfg.numtrl);
  for iTr = 1:cfg.numtrl
    diptime{iTr} = (((1:nsample)-1)/cfg.fsample) - cfg.baseline;
  end
end

if ~isempty(cfg.sourcemodel.signal)
  % use the user-supplied signal for the dipoles
  dipsignal = cfg.sourcemodel.signal;
else
  dipsignal = cell(1, cfg.numtrl);
  for iTr = 1:cfg.numtrl
    % compute a cosine signal with the desired frequency, phase and amplitude for each dipole
    for i=1:Ndipoles
      dipsignal{iTr}(i,:) = cos(cfg.sourcemodel.frequency(i)*diptime{iTr}*2*pi + cfg.sourcemodel.phase(i)) * cfg.sourcemodel.amplitude(i);
    end
  end
end

dippos = cfg.sourcemodel.pos;
dipmom = cfg.sourcemodel.mom;

if ~iscell(dipmom)
  dipmom = {dipmom};
end

if ~iscell(dippos)
  dippos = {dippos};
end

if length(dippos)==1
  dippos = repmat(dippos, 1, cfg.numtrl);
elseif length(dippos)~=cfg.numtrl
  ft_error('incorrect number of trials specified in the dipole position');
end

if length(dipmom)==1
  dipmom = repmat(dipmom, 1, cfg.numtrl);
elseif length(dipmom)~=cfg.numtrl
  ft_error('incorrect number of trials specified in the dipole moment');
end

data.time  = diptime;
data.trial = {};

ft_progress('init', cfg.feedback, 'computing data data');
for trial=1:cfg.numtrl
  ft_progress(trial/cfg.numtrl, 'computing data data for trial %d\n', trial);
  if numel(cfg.chanunit) == numel(cfg.channel)
    lf = ft_compute_leadfield(dippos{trial}, sens, headmodel, 'dipoleunit', cfg.dipoleunit, 'chanunit', cfg.chanunit, leadfieldopt{:});
  else
    lf = ft_compute_leadfield(dippos{trial}, sens, headmodel, leadfieldopt{:});
  end
  nsamples = size(dipsignal{trial},2);
  nchannels = size(lf,1);
  data.trial{trial} = zeros(nchannels,nsamples);
  for i = 1:3
    data.trial{trial} = data.trial{trial} + ...
      lf(:,i:3:end) * (repmat(dipmom{trial}(i:3:end),1,nsamples) .* dipsignal{trial});
  end
end
ft_progress('close');

if ft_senstype(sens, 'meg')
  data.grad = sens;
elseif ft_senstype(sens, 'eeg')
  data.elec = sens;
end

% determine RMS value of data data
ss = 0;
sc = 0;
for trial=1:cfg.numtrl
  ss = ss + sum(data.trial{trial}(:).^2);
  sc = sc + length(data.trial{trial}(:));
end
rms = sqrt(ss/sc);
ft_info('RMS value of data data is %g\n', rms);

% add noise to the data data
for trial=1:cfg.numtrl
  relnoise = randn(size(data.trial{trial})) * cfg.relnoise * rms;
  absnoise = randn(size(data.trial{trial})) * cfg.absnoise;
  data.trial{trial} = data.trial{trial} + relnoise + absnoise;
end

data.fsample = cfg.fsample;
data.label   = sens.label;

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble randomseed
ft_postamble provenance data
ft_postamble history    data
ft_postamble savevar    data

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function ival = c8mat_is_hermitian ( m, n, a )

%*****************************************************************************80
%
%% C8MAT_IS_HERMITIAN checks if a complex matrix is hermitian.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of A.
%
%    Input, complex A(M,N), the matrix.
%
%    Output, integer IVAL:
%    -1, the matrix is not symmetric because M /= N.
%    -2, the matrix is not symmetric because A(I,J) = conjg ( A(J,I) ) fails.
%    1, the matrix is symmetric.
%
  if ( m ~= n )
    ival = -1;
    return
  end

  for i = 1 : n
    ffor j = 1 : i - 1
      if ( a(i,j) ~= conj ( a(j,i) ) )
        ival = -2;
      end
    end
  end

  ival = 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [CellResp,data,dimCR] = LoadFileFromParts(data_dir,filename,namestring)
disp('loading file from parts...');
% tic
fishdir = fullfile(data_dir,filename);
load(fishdir,'data');
load(fishdir,'dimCR');
CellResp = zeros(dimCR);
num = 0;
nParts = round(dimCR(1)*dimCR(2)/(10^8));
disp(['in ' num2str(nParts) ' parts']);
for i = 1:nParts,
    disp(num2str(i));
    load(fishdir,[namestring,'_' num2str(i)]);
    eval(['len = size(',namestring,'_' num2str(i) ',1);']);
    eval(['CellResp(num+1:num+len,:) = ',namestring,'_' num2str(i) ';']);
    num = num+len;
end
% toc

end

author: xiuyechen

repo: FishExplorer


text:
function res = reconWUpyr(pyr, pind, daub_order);

% RES = reconWUpyr(PYR, INDICES, DAUB_ORDER)

% Reconstruct image from its separable undecimated orthonormal QMF/wavelet pyramid
% representation, as created by buildWUpyr.
%
% PYR is a vector containing the N pyramid subbands, ordered from fine
% to coarse.  INDICES is an Nx2 matrix containing the sizes of
% each subband.
%
% DAUB_ORDER: specifies the order of the daubechies wavelet filter used

% JPM, Univ. de Granada, 03/2003, based on Rice Wavelet Toolbox
% functions "mrdwt" and  "mirdwt", and on Matlab Pyrtools from Eero Simoncelli.

Nor = 3;
Nsc = (size(pind,1)-2)/Nor-1;
h = daubcqf(daub_order);

yh = [];

nband = 1;
last = prod(pind(1,:)); % empty "high pass residual band" for compatibility with full steerpyr 2
for nsc = 1:Nsc+1,  % The number of scales corresponds to the number of pyramid levels (also for compatibility)
    for nor = 1:Nor,
        nband = nband +1;
        first = last + 1;
        last = first + prod(pind(nband,:)) - 1;
        band = pyrBand(pyr,pind,nband);
        sh = (daub_order/2 - 1)*2^nsc;  % approximate phase compensation
        if nsc > 2,
            band = expand(band, 2^(nsc-2));
        end
        if nor == 1,        % horizontal
            band = shift(band, [-sh -2^(nsc-1)]);
        elseif nor == 2,    % vertical
            band = shift(band, [-2^(nsc-1) -sh]);
        else
            band = shift(band, [-sh -sh]);    % diagonal
        end
        yh = [yh band];
    end
end

nband = nband + 1;
band = pyrBand(pyr,pind,nband);
lpr = expand(band,2^Nsc);

res= mirdwt(lpr,yh,h,Nsc+1);

author: ricedsp


text:
function scale_change_factor = scale_filter_track(im, pos, base_target_sz, currentScaleFactor, scale_filter, params)

% Track the scale using the scale filter.

% Get scale filter features
scales = currentScaleFactor*scale_filter.scaleSizeFactors;
xs = extract_scale_sample(im, pos, base_target_sz, scales, params.scale_model_sz, params.use_mexResize);

% Project
xs = feature_projection_scale(xs, scale_filter.basis, scale_filter.window);

% Get scores
xsf = fft(xs, [], 2);
scale_responsef = sum(scale_filter.sf_num .* xsf, 1) ./ (scale_filter.sf_den + params.lambda);
interp_scale_response = ifft(resizeDFT(scale_responsef, params.number_of_interp_scales), 'symmetric');

recovered_scale_index = find(interp_scale_response == max(interp_scale_response(:)), 1);

if params.do_poly_interp
    % Fit a quadratic polynomial to get a refined scale
    % estimate.
    id1 = mod(recovered_scale_index -1 -1,params.number_of_interp_scales)+1;
    id2 = mod(recovered_scale_index +1 -1,params.number_of_interp_scales)+1;

    poly_x = [scale_filter.interpScaleFactors(id1), scale_filter.interpScaleFactors(recovered_scale_index), scale_filter.interpScaleFactors(id2)];
    poly_y = [interp_scale_response(id1), interp_scale_response(recovered_scale_index), interp_scale_response(id2)];

    poly_A_mat = [poly_x(1)^2, poly_x(1), 1;...
        poly_x(2)^2, poly_x(2), 1;...
        poly_x(3)^2, poly_x(3), 1 ];

    poly = poly_A_mat\poly_y';

    scale_change_factor = -poly(2)/(2*poly(1));
else
    scale_change_factor = scale_filter.interpScaleFactors(recovered_scale_index);
end

author: martin-danelljan

repo: ECO


text:
classdef Net < handle
    %NET  Create and manipulate comprehensive artificial neural networks
    %
    % This module contains:
    %
    % - API for new layers creation, layers are building bricks of neural
    %   networks;
    % - set of built-in most-useful Layers;
    % - API to constuct and modify comprehensive neural networks from layers;
    % - functionality for loading serialized networks models from different
    %   frameworks.
    %
    % Functionality of this module is designed only for forward pass
    % computations (i. e. network testing). A network training is in principle
    % not supported.
    %
    % [Wiki](https://github.com/opencv/opencv/wiki/Deep-Learning-in-OpenCV)
    %
    % ## Net class
    % Neural network is presented as directed acyclic graph (DAG), where
    % vertices are Layer instances, and edges specify relationships between
    % layers inputs and outputs.
    %
    % Each network layer has unique integer id and unique string name inside
    % its network. LayerId can store either layer name or layer id.
    %
    % See also: cv.Net.Net, nnet.cnn.layer.Layer, trainNetwork,
    %  SeriesNetwork, importCaffeNetwork, importCaffeLayers, alexnet, vgg16,
    %  vgg19
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    %% Constructor/destructor
    methods
        function this = Net(varargin)
            %NET  Constructor and importer of trained serialized models from different dnn-frameworks
            %
            %     net = cv.Net()
            %
            %     net = cv.Net('Caffe', prototxt)
            %     net = cv.Net('Caffe', prototxt, caffeModel)
            %
            %     net = cv.Net('Tensorflow', modelmodel)
            %     net = cv.Net('Tensorflow', model, config)
            %
            %     net = cv.Net('Torch', filename)
            %     net = cv.Net('Torch', filename, isBinary)
            %
            %     net = cv.Net('Darknet', cfgFile)
            %     net = cv.Net('Darknet', cfgFile, darknetModel)
            %
            % ## Input
            % * __prototxt__ path to the `.prototxt` file with text
            %   description of the network architecture.
            % * __caffeModel__ (optional) path to the `.caffemodel` file with
            %   learned network. Empty by default.
            % * __model__ path to the `.pb` file with binary protobuf
            %   description of the network architecture. Binary serialized
            %   TensorFlow graph includes weights.
            % * __config__ Optional path to the `.pbtxt` file that contains
            %   text graph definition in protobuf format. Resulting net is
            %   built by text graph using weights from a binary one. This is
            %   more flexible than binary format and may be used to build the
            %   network using binary format only as a weights storage. This
            %   approach is similar to Caffe's `.prorotxt` and `.caffemodel`.
            % * __filename__ path to the file, dumped from Torch by using
            %   `torch.save()` function.
            % * __isBinary__ specifies whether the network was serialized in
            %   ascii mode or binary. default true.
            % * __cfgFile__ path to the `.cfg` file with text description of
            %   the network architecture.
            % * __darknetModel__ (optional) path to the `.weights` file with
            %   learned network.
            %
            % The first variant creates an empty network.
            %
            % The second variant reads a network model stored in
            % [Caffe](http://caffe.berkeleyvision.org) framework's format.
            %
            % The third variant reads a network model stored in
            % [TensorFlow](https://www.tensorflow.org/) framework's format.
            %
            % The fourth variant reads a network model stored in
            % [Torch7](http://torch.ch) framework's format.
            %
            % The fifth variant reads a network model stored in
            % [Darknet](https://pjreddie.com/darknet/) model files.
            %
            % The importers first create a net, add loaded layers into it, and
            % set connections between them.
            %
            % ### Notes for Torch
            %
            % NOTE: ASCII mode of Torch serializer is more preferable, because
            % binary mode extensively use `long` type of C language, which has
            % various bit-length on different systems.
            %
            % The loading file must contain serialized
            % [`nn.Module`](https://github.com/torch/nn/blob/master/doc/module.md)
            % object with importing network. Try to eliminate a custom objects
            % from serialazing data to avoid importing errors.
            %
            % List of supported layers (i.e. object instances derived from
            % Torch `nn.Module` class):
            % - `nn.Sequential`
            % - `nn.Parallel`
            % - `nn.Concat`
            % - `nn.Linear`
            % - `nn.SpatialConvolution`
            % - `nn.SpatialMaxPooling`, `nn.SpatialAveragePooling`
            % - `nn.ReLU`, `nn.TanH`, `nn.Sigmoid`
            % - `nn.Reshape`
            % - `nn.SoftMax`, `nn.LogSoftMax`
            %
            % Also some equivalents of these classes from cunn, cudnn, and
            % fbcunn may be successfully imported.
            %
            % See also: cv.Net, cv.Net.forward
            %
            this.id = Net_(0, 'new', varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     net.delete()
            %
            % See also: cv.Net
            %
            if isempty(this.id), return; end
            Net_(this.id, 'delete');
        end
    end

    %% Net (set/get blobs and params, forward pass)
    methods
        function setInput(this, blob, name)
            %SETINPUT  Sets the new value for the layer output blob
            %
            %     net.setInput(blob)
            %     net.setInput(blob, name)
            %
            % ## Input
            % * __blob__ new blob, constructed from an image or an array of
            %   images.
            % * __name__ descriptor of the updating layer output blob. See
            %   cv.Net.connect to know format of the descriptor.
            %
            % NOTE: If updating blob is not empty then `blob` must have the
            % same shape, because network reshaping is not implemented yet.
            %
            % The blob (4-dimensional blob, so-called batch) is constructed
            % from image or array of images. Image is a 2-dimensional
            % multi-channel or 3-dimensional single-channel image (or array of
            % such images).
            %
            % See also: cv.Net.forward, cv.Net.blobFromImages
            %
            if nargin > 2
                Net_(this.id, 'setInput', blob, name);
            else
                Net_(this.id, 'setInput', blob);
            end
        end

        function setParam(this, layerId, numParam, blob)
            %SETPARAM  Sets the new value for the learned param of the layer
            %
            %     net.setParam(layerId, numParam, blob)
            %
            % ## Input
            % * __layerId__ name or id of the layer.
            % * __numParam__ index of the layer parameter in the blobs array.
            % * __blob__ the new value.
            %
            % NOTE: If shape of the new blob differs from the previous shape,
            % then the following forward pass may fail.
            %
            % See also: cv.Net.getParam
            %
            Net_(this.id, 'setParam', layerId, numParam, blob);
        end

        function blob = getParam(this, layerId, numParam)
            %GETPARAM  Returns parameter blob of the layer
            %
            %     blob = net.getParam(layerId)
            %     blob = net.getParam(layerId, numParam)
            %
            % ## Input
            % * __layerId__ name or id of the layer.
            % * __numParam__ index of the layer parameter in the blobs array.
            %   default 0.
            %
            % ## Output
            % * __blob__ returned parameter blob.
            %
            % Parameters are the weights and biases.
            %
            % See also: cv.Net.setParam
            %
            if nargin > 2
                blob = Net_(this.id, 'getParam', layerId, numParam);
            else
                blob = Net_(this.id, 'getParam', layerId);
            end
        end

        function blob = forward(this, varargin)
            %FORWARD  Runs forward pass
            %
            %     blob = net.forward()
            %     blob = net.forward(outputName)
            %
            %     blobs = net.forward(outBlobNames)
            %
            % ## Input
            % * __outputName__ name for layer which output is needed to get.
            % * __outBlobNames__ names for layers which outputs are needed to
            %   get.
            %
            % ## Output
            % * __blob__ blob for first output of specified layer.
            % * __blobs__ blobs for first outputs of specified layers
            %   (cell array).
            %
            % The first form runs forward pass to compute output of layer
            % with name `outputName`. By default (`outputName` not specified)
            % runs forward pass for the whole network.
            % (i.e `names = net.getLayerNames(); outputName = names(end);`).
            % It returns blob for first output of specified layer.
            %
            % The second form runs forward pass to compute outputs of layers
            % listed in `outBlobNames`. It returns blobs for first outputs of
            % specified layers.
            %
            % See also: cv.Net.forwardAndRetrieve, cv.Net.Net
            %
            blob = Net_(this.id, 'forward', varargin{:});
        end

        function blobs = forwardAndRetrieve(this, varargin)
            %FORWARDANDRETRIEVE  Runs forward pass
            %
            %     blobs = net.forwardAndRetrieve()
            %     blobs = net.forwardAndRetrieve(outputName)
            %
            %     blobsArr = net.forwardAndRetrieve(outBlobNames)
            %
            % ## Input
            % * __outputName__ name for layer which output is needed to get.
            % * __outBlobNames__ names for layers which outputs are needed to
            %   get.
            %
            % ## Output
            % * __blobs__ contains all output blobs for specified layer
            %   (cell array)
            % * __blobsArr__ contains all output blobs for each layer
            %   specified in `outBlobNames` (cell array of cell arrays).
            %
            % The first form runs forward pass to compute output of layer
            % with name `outputName`. By default (`outputName` not specified)
            % runs forward pass for the whole network
            % (i.e `names = net.getLayerNames(); outputName = names(end);`).
            % It returns all output blobs for specified layer.
            %
            % The second form runs forward pass to compute outputs of layers
            % listed in `outBlobNames`. It returns all output blobs for each
            % layer specified in `outBlobNames`.
            %
            % See also: cv.Net.forward, cv.Net.Net
            %
            blobs = Net_(this.id, 'forwardAndRetrieve', varargin{:});
        end

        function [timings, total] = getPerfProfile(this)
            %GETPERFPROFILE  Returns overall time for inference and timings (in ticks) for layers
            %
            %     [timings, total] = net.getPerfProfile()
            %
            % ## Output
            % * __timings__ vector for tick timings for all layers.
            % * __total__ overall ticks for model inference.
            %
            % Indexes in returned vector correspond to layers ids. Some layers
            % can be fused with others, in this case zero ticks count will be
            % return for that skipped layers.
            %
            % See also: cv.Net.forward, cv.TickMeter
            %
            [timings, total] = Net_(this.id, 'getPerfProfile');
        end
    end

    %% Net (network architecture)
    methods
        function b = empty(this)
            %EMPTY  Returns true if there are no layers in the network.
            %
            %     b = net.empty()
            %
            % ## Output
            % * __b__ Boolean.
            %
            % See also: cv.Net.Net
            %
            b = Net_(this.id, 'empty');
        end

        function id = addLayer(this, name, layerType, params)
            %ADDLAYER  Adds new layer to the net
            %
            %     id = net.addLayer(name, layerType, params)
            %
            % ## Input
            % * __name__ unique name of the adding layer.
            % * __layerType__ typename of the adding layer (type must be
            %   registered).
            % * __params__ parameters which will be used to initialize the
            %   creating layer. Scalar structure with the following fields:
            %   * __dict__ name-value dictionary as struct, values are scalar
            %     values (or arrays) of one of the following type: double,
            %     integer, or string.
            %   * __blobs__ List of learned parameters stored as blobs.
            %   * __name__ Name of the layer instance (optional, can be used
            %     internal purposes).
            %   * __type__ Type name which was used for creating layer by
            %     layer factory (optional).
            %
            % ## Output
            % * __id__ unique identifier of created layer, or -1 if a failure
            %   will happen.
            %
            % A LayerParams provides all data needed to initialize layer. It
            % includes dictionary with scalar params (`params.dict` struct),
            % blob params `params.blobs` and optional meta information
            % `params.name` and `params.type` of layer instance.
            %
            % Built-in layers listed below partially reproduce functionality
            % of corresponding Caffe and Torch7 layers. In partuclar, the
            % following layers and Caffe importer were tested to reproduce
            % [Caffe](http://caffe.berkeleyvision.org/tutorial/layers.html)
            % functionality:
            % - Convolution
            % - Deconvolution
            % - Pooling
            % - InnerProduct
            % - TanH, ReLU, Sigmoid, BNLL, Power, AbsVal
            % - Softmax
            % - Reshape, Flatten, Slice, Split
            % - LRN
            % - MVN
            % - Dropout (since it does nothing on forward pass)
            %
            % See also: cv.Net.addLayerToPrev, cv.Net.deleteLayer, cv.Net.connect
            %
            id = Net_(this.id, 'addLayer', name, layerType, params);
            id = int32(id);
        end

        function id = addLayerToPrev(this, name, layerType, params)
            %ADDLAYERTOPREV  Adds new layer and connects its first input to the first output of previously added layer
            %
            %     id = net.addLayerToPrev(name, layerType, params)
            %
            % ## Input
            % * __name__ unique name of the adding layer.
            % * __layerType__ typename of the adding layer (type must be
            %   registered).
            % * __params__ parameters which will be used to initialize the
            %   creating layer.
            %
            % ## Output
            % * __id__ unique identifier of created layer, or -1 if a failure
            %   will happen.
            %
            % See also: cv.Net.addLayer, cv.Net.deleteLayer, cv.Net.connect
            %
            id = Net_(this.id, 'addLayerToPrev', name, layerType, params);
            id = int32(id);
        end

        function id = getLayerId(this, name)
            %GETLAYERID  Converts string name of the layer to the integer identifier
            %
            %     id = net.getLayerId(name)
            %
            % ## Input
            % * __name__ string name of the layer.
            %
            % ## Output
            % * __id__ id of the layer, or -1 if the layer wasn't found.
            %
            % See also: cv.Net.getLayer, cv.Net.getLayerNames
            %
            id = Net_(this.id, 'getLayerId', name);
            id = int32(id);
        end

        function names = getLayerNames(this)
            %GETLAYERNAMES  Get layer names
            %
            %     names = net.getLayerNames()
            %
            % ## Output
            % * __names__ names of layers.
            %
            % See also: cv.Net.getLayerId, cv.Net.getLayer
            %
            names = Net_(this.id, 'getLayerNames');
        end

        function layer = getLayer(this, layerId)
            %GETLAYER  Returns layer with specified id or name which the network use
            %
            %     layer = net.getLayer(layerId)
            %
            % ## Input
            % * __layerId__ layer name or layer id.
            %
            % ## Output
            % * __layer__ returned layer. Scalar structure with the following
            %   fields:
            %   * __blobs__ List of stored learned parameters as returned by
            %     cv.Net.getParam.
            %   * __name__ name of the layer instance, can be used for logging
            %     or other internal purposes.
            %   * __type__ Type name which was used for creating layer by
            %     layer factory.
            %   * __preferableTarget__ preferred target for layer forwarding
            %     (see cv.Net.setPreferableTarget).
            %
            % Layers are the building blocks of networks.
            %
            % See also: cv.Net.getLayerId
            %
            layer = Net_(this.id, 'getLayer', layerId);
        end

        function layers = getLayerInputs(this, layerId)
            %GETLAYERINPUTS  Returns input layers of specific layer
            %
            %     layers = net.getLayerInputs(layerId)
            %
            % ## Input
            % * __layerId__ layer name or layer id.
            %
            % ## Output
            % * __layers__ returned layers, struct array.
            %
            % See also: cv.Net.getLayerId, cv.Net.getLayer
            %
            layers = Net_(this.id, 'getLayerInputs', layerId);
        end

        function deleteLayer(this, layerId)
            %DELETELAYER  Delete layer for the network
            %
            %     net.deleteLayer(layerId)
            %
            % ## Input
            % * __layerId__ layer name or layer id.
            %
            % Warning: Not yet implemented.
            %
            % See also: cv.Net.addLayer
            %
            Net_(this.id, 'deleteLayer', layerId);
        end

        function connect(this, varargin)
            %CONNECT  Connects output of the first layer to input of the second layer
            %
            %     net.connect(outPin, inpPin)
            %     net.connect(outLayerId, outNum, inpLayerId, inpNum)
            %
            % ## Input
            % * __outPin__ descriptor of the first layer output. See below.
            % * __inpPin__ descriptor of the second layer input. See below.
            %
            % ## Input
            % * __outLayerId__ identifier of the first layer.
            % * __outNum__ number of the first layer output.
            % * __inpLayerId__ identifier of the second layer.
            % * __inpNum__ number of the second layer input.
            %
            % Descriptors have the following template
            % `<layer_name>[.input_number]`:
            % - the first part of the template `layer_name` is sting name of
            %   the added layer. If this part is empty then the network input
            %   pseudo layer will be used;
            % - the second optional part of the template `input_number` is
            %   either number of the layer input, either label one. If this
            %   part is omitted then the first layer input will be used.
            %
            % See also: cv.Net.setInputsNames, cv.Net.addLayer
            %
            Net_(this.id, 'connect', varargin{:});
        end

        function setInputsNames(this, inputBlobNames)
            %SETINPUTSNAMES  Sets outputs names of the network input pseudo layer
            %
            %     net.setInputsNames(inputBlobNames)
            %
            % ## Input
            % * __inputBlobNames__ blob names.
            %
            % Each net always has special own the network input pseudo layer
            % with `id=0`. This layer stores the user blobs only and don't
            % make any computations. In fact, this layer provides the only way
            % to pass user data into the network. As any other layer, this
            % layer can label its outputs and this function provides an easy
            % way to do this.
            %
            % See also: cv.Net.connect, cv.Net.setInput
            %
            Net_(this.id, 'setInputsNames', inputBlobNames);
        end

        function indices = getUnconnectedOutLayers(this)
            %GETUNCONNECTEDOUTLAYERS  Returns indexes of layers with unconnected outputs
            %
            %     indices = net.getUnconnectedOutLayers()
            %
            % ## Output
            % * __indices__ vector of indices.
            %
            % See also: cv.Net.getLayer
            %
            indices = Net_(this.id, 'getUnconnectedOutLayers');
        end

        function layersTypes = getLayerTypes(this)
            %GETLAYERTYPES  Returns list of types for layer used in model
            %
            %     layersTypes = net.getLayerTypes()
            %
            % ## Output
            % * __layersTypes__ layer types.
            %
            % See also: cv.Net.getLayersCount
            %
            layersTypes = Net_(this.id, 'getLayerTypes');
        end

        function count = getLayersCount(this, layerType)
            %GETLAYERSCOUNT  Returns count of layers of specified type
            %
            %     count = net.getLayersCount(layerType)
            %
            % ## Input
            % * __layerType__ type.
            %
            % ## Output
            % * __count__ count of layers.
            %
            % See also: cv.Net.getLayerTypes
            %
            count = Net_(this.id, 'getLayersCount', layerType);
        end

        function enableFusion(this, fusion)
            %ENABLEFUSION  Enables or disables layer fusion in the network
            %
            %     net.enableFusion(fusion)
            %
            % ## Input
            % * __fusion__ true to enable the fusion, false to disable. The
            %   fusion is enabled by default.
            %
            % See also: cv.Net.connect
            %
            Net_(this.id, 'enableFusion', fusion);
        end

        function setHalideScheduler(this, scheduler)
            %SETHALIDESCHEDULER  Compile Halide layers
            %
            %     net.setHalideScheduler(scheduler)
            %
            % ## Input
            % * __scheduler__ scheduler Path to YAML file with scheduling
            %   directives.
            %
            % Schedule layers that support Halide backend. Then compile them
            % for specific target. For layers that not represented in
            % scheduling file or if no manual scheduling used at all,
            % automatic scheduling will be applied.
            %
            % See also: cv.Net.setPreferableBackend
            %
            Net_(this.id, 'setHalideScheduler', scheduler);
        end

        function setPreferableBackend(this, backend)
            %SETPREFERABLEBACKEND  Ask network to use specific computation backend where it supported
            %
            %     net.setPreferableBackend(backend)
            %
            % ## Input
            % * __backend__ computation backend supported by layers, one of:
            %   * __Default__
            %   * __Halide__ Halide language backend.
            %   * __InferenceEngine__ Intel's Deep Learning Inference Engine.
            %
            % See also: cv.Net.setPreferableTarget, cv.Net.setHalideScheduler
            %
            Net_(this.id, 'setPreferableBackend', backend);
        end

        function setPreferableTarget(this, target)
            %SETPREFERABLETARGET  Ask network to make computations on specific target device
            %
            %     net.setPreferableTarget(target)
            %
            % ## Input
            % * __target__ target device for computations, one of:
            %   * __CPU__
            %   * __OpenCL__
            %
            % See also: cv.Net.setPreferableBackend
            %
            Net_(this.id, 'setPreferableTarget', target);
        end
    end

    %% Auxiliary functions
    methods (Static)
        function blob = readTorchBlob(filename, varargin)
            %READTORCHBLOB  Loads blob which was serialized as torch.Tensor object of Torch7 framework
            %
            %     blob = cv.Net.readTorchBlob(filename)
            %     blob = cv.Net.readTorchBlob(filename, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __filename__ path to the blob file.
            %
            % ## Output
            % * __blob__ output blob.
            %
            % ## Options
            % * __IsBinary__ specifies whether blob file was serialized in
            %   ascii mode or binary. default true.
            %
            % This function has the same limitations as cv.Net.Net with
            % regards to the Torch importer.
            %
            % See also: cv.Net.setInput, cv.Net.blobFromImages
            %
            blob = Net_(0, 'readTorchBlob', filename, varargin{:});
        end

        function blob = blobFromImages(img, varargin)
            %BLOBFROMIMAGES  Creates 4-dimensional blob from image or series of images
            %
            %     blob = cv.Net.blobFromImages(img)
            %     blob = cv.Net.blobFromImages(imgs)
            %     blob = cv.Net.blobFromImages(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __img__ input image (with 1-, 3- or 4-channels).
            % * __imgs__ input images (all with 1-, 3- or 4-channels).
            %
            % ## Output
            % * __blob__ 4-dimansional array with NCHW dimensions order.
            %
            % ## Options
            % * __Size__ spatial size for output image `[w,h]`. default [0,0]
            %   (in which case input image size is used)
            % * __Mean__ scalar with mean values which are subtracted from
            %   channels. Values are intended to be in
            %   (mean-R, mean-G, mean-B) order if image has BGR ordering and
            %   `SwapRB` is true. default [0,0,0]
            % * __ScaleFactor__ multiplier for images values. default 1.0
            % * __SwapRB__ flag which indicates that swap first and last
            %   channels in 3-channel image is necessary. For instance, Caffe
            %   models are usually trained on BGR images, while TensorFlow
            %   models expect RGB images as input. default true
            % * __Crop__ flag which indicates whether image will be cropped
            %   after resize or not. default true
            %
            % Creates blob and optionally resizes and crops the images from
            % center, subtracts mean values, scales values, and swaps blue and
            % red channels.
            %
            % If `Crop` is true, input image is resized so one side after
            % resize is equal to corresponding dimension in `Size` and another
            % one is equal or larger. Then, crop from the center is performed.
            % If `Crop` is false, direct resize without cropping and
            % preserving aspect ratio is performed.
            %
            % A blob is a 4-dimensional matrix (so-called batch) with the
            % following shape: `[num, cn, rows, cols]`.
            %
            % See also: cv.Net.setInput
            %
            blob = Net_(0, 'blobFromImages', img, varargin{:});
        end

        function imgs = imagesFromBlob(blob)
            %IMAGESFROMBLOB  Parse a 4D blob and output the images it contains
            %
            %     imgs = cv.Net.imagesFromBlob(blob)
            %
            % ## Input
            % * __blob__ 4-dimensional array `(images, channels, height, width)`
            %   in floating-point precision (`single`) from which you would
            %   like to extract the images.
            %
            % ## Output
            % * __imgs__ cell-array of matrices containing the images
            %   extracted from the blob in floating-point precision (`single`).
            %   They are non-normalized neither mean-added. The number of
            %   returned images equals the first dimension of the blob
            %   (batch size). Every image has a number of channels equals to
            %   the second dimension of the blob (depth).
            %
            % See also: cv.Net.blobFromImages
            %
            imgs = Net_(0, 'imagesFromBlob', blob);
        end

        function shrinkCaffeModel(src, dst, varargin)
            %SHRINKCAFFEMODEL  Convert all weights of Caffe network to half precision floating point
            %
            %     cv.Net.shrinkCaffeModel(src, dst)
            %     cv.Net.shrinkCaffeModel(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __src__ Path to origin model from Caffe framework contains
            %   single precision floating point weights (usually has
            %   `.caffemodel` extension).
            % * __dst__ Path to destination model with updated weights.
            %
            % ## Options
            % * __LayersTypes__ Set of layers types which parameters will be
            %   converted. By default (not set), converts only Convolutional
            %   and Fully-Connected layers' weights,
            %   i.e `{'Convolution', 'InnerProduct'}`.
            %
            % Note: Shrinked model has no origin `float32` weights so it can't
            % be used in origin Caffe framework anymore. However the structure
            % of data is taken from NVidia's
            % <https://github.com/NVIDIA/caffe Caffe fork>. So the resulting
            % model may be used there.
            %
            Net_(0, 'shrinkCaffeModel', src, dst, varargin{:});
        end

        function indices = NMSBoxes(bboxes, scores, score_threshold, nms_threshold, varargin)
            %NMSBOXES  Performs non-maximum suppression given boxes and corresponding scores
            %
            %     indices = cv.Net.NMSBoxes(bboxes, scores, score_threshold, nms_threshold)
            %     indices = cv.Net.NMSBoxes(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __bboxes__ a set of bounding boxes to apply NMS.
            % * __scores__ a set of corresponding confidences.
            % * **score_threshold** a threshold used to filter boxes by score.
            % * **nms_threshold** a threshold used in non maximum suppression.
            %
            % ## Output
            % * __indices__ the kept indices of bboxes after NMS.
            %
            % ## Options
            % * __Eta__ a coefficient in adaptive threshold formula:
            %   `nms_threshold_{i+1} = eta * nms_threshold_{i}`. default 1.0
            % * __TopK__ if `> 0`, keep at most `TopK` picked indices.
            %   default 0
            %
            % See also: cv.groupRectangles
            %
            indices = Net_(0, 'NMSBoxes', bboxes, scores, score_threshold, nms_threshold, varargin{:});
        end
    end
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%compute ang between vector tailsminfl and vector inflheadsm

function [data,units]=compute_angihsmvsangtsmi(trx,n)

larvae = trx.exp2flies{n};
numlarvae = numel(larvae);
angihsmvsangtsmi=cell(1,numlarvae);
absangihsmvsangtsmi=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    absangihsmvsangtsmi{1,i}=real(acos(cos(trx(larva).tailsminflang).*cos(trx(larva).inflheadsmang)+sin(trx(larva).tailsminflang).*sin(trx(larva).inflheadsmang)));
    temp=trx(larva).tailsminflang-pi/2;
    cosperp=sign(cos(temp).*cos(trx(larva).inflheadsmang)+sin(temp).*sin(trx(larva).inflheadsmang));
    angihsmvsangtsmi{1,i}=bsxfun(@times,absangihsmvsangtsmi{1,i},cosperp);
end

units=parseunits('rad');
data=angihsmvsangtsmi;

author: kristinbranson

repo: JAABA


text:
function a_heap = r8vec_heap_a ( n, a )

%*****************************************************************************80
%
%% R8VEC_HEAP_A reorders an R8VEC into an ascending heap.
%
%  Discussion:
%
%    An ascending heap is an array A with the property that, for every index J,
%    A(J) <= A(2*J) and A(J) <= A(2*J+1), (as long as the indices
%    2*J and 2*J+1 are legal).
%
%  Diagram:
%
%                  A(1)
%                /      \
%            A(2)         A(3)
%          /     \        /  \
%      A(4)       A(5)  A(6) A(7)
%      /  \       /   \
%    A(8) A(9) A(10) A(11)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Albert Nijenhuis, Herbert Wilf,
%    Combinatorial Algorithms,
%    Academic Press, 1978, second edition,
%    ISBN: 0-12-519260-6.
%
%  Parameters:
%
%    Input, integer N, the size of the input array.
%
%    Input, real A(N), an unsorted array.
%
%    Output, real A_HEAP(N), the values of the array have been
%    reordered into a heap.
%
  a_heap(1:n) = a(1:n);
%
%  Only nodes N/2 down to 1 can be "parent" nodes.
%
  for i = floor ( n/2 ) : -1 : 1
%
%  Copy the value out of the parent node.
%  Position IFREE is now "open".
%
    key = a_heap(i);
    ifree = i;

    while ( 1 )
%
%  Positions 2*IFREE and 2*IFREE + 1 are the descendants of position
%  IFREE.  (One or both may not exist because they exceed N.)
%
      m = 2 * ifree;
%
%  Does the first position exist?
%
      if ( n < m )
        break;
      end
%
%  Does the second position exist?
%
      if ( m + 1 <= n )
%
%  If both positions exist, take the smaller of the two values,
%  and update M if necessary.
%
        if ( a_heap(m+1) < a_heap(m) )
          m = m + 1;
        end

      end
%
%  If the small descendant is smaller than KEY, move it up,
%  and update IFREE, the location of the free position, and
%  consider the descendants of THIS position.
%
      if ( key <= a_heap(m) )
        break;
      end

      a_heap(ifree) = a_heap(m);
      ifree = m;

    end
%
%  Once there is no more shifting to do, KEY moves into the free spot.
%
    a_heap(ifree) = key;

  end

  return
end


text:
function model=sortModel(model,sortReversible,sortMetName,sortReactionOrder)
% sortModel
%   Sorts a model based on metabolite names and compartments
%
%   model             a model structure
%   sortReversible    sorts the reversible reactions so the the metabolite
%                     that is first in lexiographical order is a reactant
%                     (opt, default true)
%   sortMetName       sort the metabolite names in the equation, also uses
%                     compartment abbreviation (opt, default false)
%   sortReactionOrder sorts the reaction order within each subsystem so that
%                     reactions consuming some metabolite comes efter
%                     reactions producing it. This overrides the
%                     sortReversible option and reactions are sorted so that
%                     the production direction matches the consumption
%                     direction (opt, default false)
%
%   model             an updated model structure
%
%   Usage: model=sortModel(model,sortReversible,sortMetName,sortReactionOrder)

if nargin<2
    sortReversible=true;
end
if nargin<3
    sortMetName=false;
end
if nargin<4
    sortReactionOrder=false;
end

if sortMetName==true
    %Assuming that metComps are the indexes. Should be changed at one point
    [~, metIndexes]=sort(strcat(model.metNames,'[',model.comps(model.metComps),']'));
    model=permuteModel(model,metIndexes,'mets');
end

if sortReversible==true && sortReactionOrder==false
    %Get all reversible reactions
    revIndexes=find(model.rev);

    %Loop through them
    for i=1:numel(revIndexes)
        %Create a cell array with all the metabolite names
        mets=find(model.S(:,revIndexes(i)));
        metNames=strcat(model.metNames(mets),model.comps(model.metComps(mets)));

        if iscellstr(metNames)
            [~, indexes]=sort(metNames);

            if model.S(mets(indexes(1)),revIndexes(i))>0
                model.S(:,revIndexes(i))=model.S(:,revIndexes(i))*-1;
            end
        end
    end
end

if sortReactionOrder==true
    %Check if the model has sub-systems, otherwise throw an error
    if ~isfield(model,'subSystems')
        EM='The model must contain a subSystems field in order to sort reaction order';
        dispEM(EM);
    end

    subsystemsUnique='';
    subsystemsConcatenated='';
    for i=1:numel(model.subSystems)
        subsystemsConcatenated{i,1}=strjoin(model.subSystems{i,1},';');
        if ~isempty(model.subSystems{i,1})
            for j=1:numel(model.subSystems{i,1})
                subsystemsUnique{numel(subsystemsUnique)+1,1}=model.subSystems{i,1}{1,j};
            end
        end
    end
    subsystemsUnique=unique(subsystemsUnique);
    for i=1:numel(subsystemsUnique)
        %Get all reactions for that subsystem
        rxns=find(~cellfun(@isempty,regexp(subsystemsConcatenated,subsystemsUnique(i))));

        %Temporarily ignore large subsystems because of inefficient
        %implementation
        if numel(rxns)<2 || numel(rxns)>250
            continue;
        end

        nRxns=numel(rxns);
        revRxns=rxns(model.rev(rxns)~=0);

        %This variable will hold the current reversibility directions of
        %the reversible reactions. 1 means the same direction as in the
        %original model and -1 means the opposite direction.
        oldRev=ones(numel(revRxns),1);

        %The problem could probably be solved analytically but a simple
        %random method is implemented here instead. Two reactions are
        %chosen randomly and their positions are switched. A score is
        %calculated based on the number of metabolites that are produced
        %before they are consumed. If the perturbed model has a better or
        %equal score than the original the reaction order is switched. If
        %no increase in score has been seen after 1000*rxnsInSubsystem then
        %the optimization is terminated

        rxnOrder=1:nRxns;
        oldScore=-inf;
        counter=0;
        firstIter=true;
        while 1==1
            counter=counter+1;
            if counter==100*nRxns
                break;
            end

            newRxnOrder=rxnOrder;
            rev=oldRev;

            if firstIter==false
                y=randperm(nRxns,2);

                %Switch the order
                newRxnOrder(y(1))=rxnOrder(y(2));
                newRxnOrder(y(2))=rxnOrder(y(1));

                %With a 50% chance, also switch the reversibility of one of
                %the reactions
                if rand()>0.5 && numel(rev)>1
                    n=randperm(numel(rev),1);
                    rev(n)=rev(n)*-1;
                end
            end
            firstIter=false;

            tempS=model.S;

            %Switch the directionalities
            for j=1:numel(rev)
                if rev(j)==-1
                    tempS(:,revRxns(j))=tempS(:,revRxns(j)).*-1;
                end
            end

            %Get the metabolites that are involved and when they are
            %produced/consumed
            s=tempS(:,newRxnOrder);

            %Remove mets that aren't used in both directions
            s=s(any(s,2),:);

            %Add so that all mets are produced and consumed in the end
            s=[s ones(size(s,1),1) ones(size(s,1),1)*-1];

            %For each metabolite, find the reaction where it's first
            %produced and the reaction where it's first consumed
            s1=s>0;
            r1=arrayfun(@(x) find(s1(x,:),1,'first'),1:size(s1,1));
            s2=s<0;
            r2=arrayfun(@(x) find(s2(x,:),1,'first'),1:size(s2,1));

            score=sum(r1<r2);

            if score>=oldScore
                if score>oldScore
                    counter=0;
                end
                oldScore=score;
                oldRev=rev;
                rxnOrder=newRxnOrder;
            end
        end

        %Update the model for this subsystem
        for j=1:numel(oldRev)
            if oldRev(j)==-1
                model.S(:,revRxns(j))=model.S(:,revRxns(j)).*-1;
            end
        end
        order=1:numel(model.rxns);
        order(rxns)=rxns(rxnOrder);
        model=permuteModel(model, order, 'rxns');
    end
end
end

author: SysBioChalmers

repo: RAVEN


text:
function sliderc(r,L1,rec,sp)
% SLIDERC Slider Crank Mechanism
% sliderc(r,L1,rec) plots the slider crank mechanism specified and rotates it,
% where:
%
% r   : Radius of crank
% L1  : Length of piston
% rec : Side lengths of slider
% sp  : Speed of rotation
%
% Example 1:
%
% sliderc(1,3,[1,.5],.5)
%
%
% Example 2:
%
% sliderc(2,5,[5,.5],2)
%
%
% Uses the function Draw a Circle by Zhenhai Wang on the MATLAB File Exchange.
%
% numandina@gmail.com

L2=rec(1); % get horizontal side of slider
L3=rec(2); % get vertical side of slider
g=circle([0 0],r,100); % draw the crank
grid on
hold on
t=0; % originally, the angle theta is zero

xL1=r*cos(t); % the first piston end x location
yL1=r*sin(t); % the first piston end y location

xL2=r*cos(t)+(L1^2-yL1^2)^(1/2); % the second piston end x location (Pythagoras theorem)
yL2=0;                           % the second piston end y location

ln=line([xL1 xL2],[yL1 yL2]); % draw the piston
set([ln,g],'linewidth',3,'color','k') % format the lines
plot(0,0,'k+') % centre of crank
set(gcf,'color','w','menubar','none') % more formatting

while true % rotate indefinitely

	% piston start and end x and y locations, just like before
	xL1=r*cos(t);
	yL1=r*sin(t);
	xL2=r*cos(t)+(L1^2-yL1^2)^(1/2);
	yL2=0;

	set(ln,'xdata',[xL1 xL2],'ydata',[yL1 yL2]); % draw the piston accordingly

	% put coloured points at piston start and end locations and update them
	delete(findobj('marker','.'))
	plot([xL1 xL2],[yL1 yL2],'r.')

	% slider location points

	xp1=xL2; % first end x location is the same as the piston second end x location, since they are connected
	xp2=xL2+L2; % second end x location depends on the horizontal side length

	yp1=yL2-L3/2; %	the vertical side lengths determine the rest of the slider shape.
	yp2=yL2+L3/2;

	gr=[.5 .5 .5]; % RGB colour vector, [1 1 1] is white, [0 0 0] is black, so [.5 .5 .5] is grey

 	pl=patch([xp1,xp1,xp2,xp2],[yp1,yp2,yp2,yp1],gr); % draw the slider and paint it grey

 	axis([-(r+1) r+L1+L2+1 min([-(r+1),-(L3/2+1)]) max([r+1,L3/2+1])]) % adjust the axis limits
	axis equal % so that the circle would appear proper
	drawnow()  % update axes
 	delete(pl) % delete slider so it would be updated in next iteration
	t=t+.1*sp; % adjust theta for next iteration
end

end

author: Sable

repo: mcbench-benchmarks