text:
%Summary simulation logger. Implements simulation logger interface.
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 1/10/2011
classdef SummaryLogger < edu.stanford.covert.cell.sim.util.Logger
    %options
    properties (Access = protected)
        stepSizeSec
        verbosity
        outputDirectory
    end

    %stored data
    properties (SetAccess = protected)
        log
        figureHandle
    end

    properties (Access = protected)
        timer
    end

    %indices into simulation states, processes
    properties (Access = protected)
        stateIndex_time
        stateIndex_mass
        stateIndex_rna
        stateIndex_metabolite
        stateIndex_monomer
        stateIndex_complex
        stateIndex_metabolicReaction
        stateIndex_chromosome
        stateIndex_geometry
        stateIndex_ftsZRing

		processIndex_transcription
		processIndex_translation
        processIndex_replicationInitiation
        processIndex_replication
        processIndex_ftsZPolymerization
    end

    properties(Constant = true)
        SIM_STATUS_DIDNT_START = -3;
        SIM_STATUS_STILL_RUNNING = -2;
        SIM_STATUS_TERMINATED_WITH_ERROR = -1;
        SIM_STATUS_COMPLETED_WITH_DIVISION = 0;
        SIM_STATUS_COMPLETED_WITHOUT_DIVISION = 1;

        SIM_STATUS_INDEX_STATUS = 1;
        SIM_STATUS_INDEX_REP_INIT_TIME = 2;
        SIM_STATUS_INDEX_INIT_REP_TIME = 3;
        SIM_STATUS_INDEX_REP_TIME = 4;
        SIM_STATUS_INDEX_CYTOKINESIS_TIME = 5;
        SIM_STATUS_INDEX_INIT_GROWTH = 6;
        SIM_STATUS_INDEX_FIN_GROWTH = 7;
        SIM_STATUS_INDEX_CUM_GROWTH = 8;
        SIM_STATUS_INDEX_FIN_MASS = 9;
        SIM_STATUS_INDEX_INIT_MASS = 10;
        SIM_STATUS_INDEX_MASS_DOUBLING_TIME = 11;
    end

    %ids
    properties
        ntpWholeCellModelIDs
        aaWholeCellModelIDs
    end

    methods
        function this = SummaryLogger(stepSizeSec, verbosity, outputDirectory)
            if nargin >= 1
                this.stepSizeSec = stepSizeSec;
            else
                this.stepSizeSec = 1;
            end

            if nargin >= 2
                this.verbosity = verbosity;
            else
                this.verbosity = 2;
            end

            if nargin >= 3
                this.outputDirectory = outputDirectory;
            else
                this.outputDirectory = [];
            end
        end
    end

    methods
        function setOptions(this, options)
            metaClass = metaclass(this);
            fields = intersect(fieldnames(options), cellfun(@(x) x.Name, metaClass.Properties, 'UniformOutput',false));
            for i = 1:numel(fields)
                this.(fields{i}) = options.(fields{i});
            end
        end

        %Allocates space for small summary stats gathered after each segment.
        function this = initialize(this, sim)
            %indices
            this.stateIndex_time = sim.stateIndex('Time');
            this.stateIndex_mass = sim.stateIndex('Mass');
            this.stateIndex_rna = sim.stateIndex('Rna');
            this.stateIndex_metabolite = sim.stateIndex('Metabolite');
            this.stateIndex_monomer = sim.stateIndex('ProteinMonomer');
            this.stateIndex_complex = sim.stateIndex('ProteinComplex');
            this.stateIndex_metabolicReaction = sim.stateIndex('MetabolicReaction');
            this.stateIndex_chromosome = sim.stateIndex('Chromosome');
            this.stateIndex_geometry = sim.stateIndex('Geometry');
            this.stateIndex_ftsZRing = sim.stateIndex('FtsZRing');

			this.processIndex_transcription = sim.processIndex('Transcription');
			this.processIndex_translation = sim.processIndex('Translation');
            this.processIndex_replicationInitiation = sim.processIndex('ReplicationInitiation');
            this.processIndex_replication = sim.processIndex('Replication');
            this.processIndex_ftsZPolymerization = sim.processIndex('FtsZPolymerization');

            %ids
            this.ntpWholeCellModelIDs = sim.state('Metabolite').wholeCellModelIDs(sim.state('Metabolite').ntpIndexs);
            this.aaWholeCellModelIDs = sim.state('Metabolite').wholeCellModelIDs(sim.state('Metabolite').aminoAcidIndexs);

            %initialize timer
            this.timer = tic;

            %allocate memory
            nSeg = sim.lengthSec / this.stepSizeSec + 1;
            validateattributes(nSeg, {'numeric'}, {'real', 'positive', 'integer'});

            this.log = struct();
            this.log.runTime               = zeros(1, nSeg, 'double');         %run time of simulation
            this.log.growth                = zeros(1, nSeg, 'double');         %growth rate
            this.log.mass                  = zeros(1, nSeg, 'double');         %cell mass
            this.log.time                  = 0:this.stepSizeSec:sim.lengthSec; %time (seconds)
            this.log.metabolites           = zeros(13, nSeg, 'int32');         %total numbers of ATP, ADP, AMP
            this.log.metaboliteUsages      = zeros(2, nSeg, 'int32');         %total numbers of ATP, ADP, AMP
            this.log.metaboliteProductions = zeros(2, nSeg, 'int32');         %total numbers of ATP, ADP, AMP
            this.log.rnas                  = zeros(7, nSeg, 'int32');          %total numbers of RNA; mature mRNA, rRNA, sRNA, tRNA
            this.log.proteins              = zeros(7, nSeg, 'int32');          %total number of protein monomers and complexes
            this.log.rnaPolymerases        = zeros(1, nSeg, 'int16');          %total number of RNA polymerases
            this.log.ribosomes             = zeros(1, nSeg, 'int16');          %total number of ribosomes
            this.log.supercoiled           = zeros(2, nSeg, 'double');         %number of supercoiled regions
            this.log.dnaA                  = zeros(2, nSeg, 'int16');          %numbers of free and total DnaA monomers
            this.log.dnaAFunctionalBoxes   = zeros(5, nSeg, 'int8');           %polymerization status of R1-5 functional DnaA boxes
            this.log.replisome             = zeros(6, nSeg, 'int32');          %positions of helicase and leading and lagging polymerases
            this.log.ftsZ                  = zeros(2, nSeg, 'int16');          %numbers of free and total FtsZ monomers
            this.log.ftsZRing              = zeros(4, nSeg, 'int32');          %size of FtsZ ring
            this.log.pinchedDiameter       = zeros(1, nSeg, 'double');         %pinched diameter
            this.log.ploidy                = zeros(1, nSeg, 'double');         %ploidy (genetic content--should range from 1 to 2)

            %initialize print, plot
            this.initializePrint();
            this.initializePlot();

            %append log
            this.append(sim);
        end

        %Appends a summary of the current state at a specific index.
        function this = append(this, sim)
            t = sim.states{this.stateIndex_time};
            i = t.values / this.stepSizeSec + 1;
            if mod(i, 1) ~= 0
                return;
            end

            %append log
            this.appendHelper(i, sim);

            %print and plot
            this.print(i);
            this.plot(i);
        end

        function appendHelper(this, i, sim)
            %references
            t = sim.states{this.stateIndex_time};
            mass = sim.states{this.stateIndex_mass};
            r = sim.states{this.stateIndex_rna};
            m = sim.states{this.stateIndex_metabolite};
            pm = sim.states{this.stateIndex_monomer};
            pc = sim.states{this.stateIndex_complex};
            mr = sim.states{this.stateIndex_metabolicReaction};
            c = sim.states{this.stateIndex_chromosome};
            geometry = sim.states{this.stateIndex_geometry};
            ftsZRing = sim.states{this.stateIndex_ftsZRing};

            %indices
            iCytosol = sim.compartment.cytosolIndexs;
            iMembrane = sim.compartment.membraneIndexs;

            %append
            this.log.runTime(1, i) = toc(this.timer); this.timer = tic;
            this.log.time(1, i) = t.values;

            this.log.mass(1, i) = sum(mass.cellDry) / mass.cellInitialDryWeight;
            this.log.growth(1, i) = mr.growth;

            this.log.metabolites(1, i) = sum(m.counts(m.atpIndexs, iCytosol));
            this.log.metabolites(2, i) = sum(m.counts(m.adpIndexs, iCytosol));
            this.log.metabolites(3, i) = sum(m.counts(m.ampIndexs, iCytosol));
            this.log.metabolites(4, i) = sum(m.counts(m.ntpIndexs, iCytosol));
            [this.log.metabolites(5, i), this.log.metabolites(6, i)] = min(m.counts(m.ntpIndexs, iCytosol));
            this.log.metabolites(7, i) = sum(m.counts(m.aminoAcidIndexs(1:20), iCytosol));
            [this.log.metabolites(8, i), this.log.metabolites(9, i)] = min(m.counts(m.aminoAcidIndexs(1:20), iCytosol));
            this.log.metabolites(10,i) = sum(m.counts(m.lipidIndexs, iMembrane));
            this.log.metabolites(11,i) = sum(m.counts(m.carbohydrateIndexs, iCytosol));
            this.log.metabolites(12,i) = sum(m.counts(m.dntpIndexs, iCytosol));
            this.log.metabolites(13,i) = sum(m.counts(m.antibioticIndexs, iCytosol));

            if ~isempty(m.processUsages)
                tmp1 = m.processUsages(sub2ind(size(m.counts), m.ntpIndexs(1), iCytosol), :);
                tmp2 = m.processUsages(sub2ind(size(m.counts), m.ntpIndexs(3), iCytosol), :);
                this.log.metaboliteUsages(1, i) = sum(max(0, -tmp1), 2);
                this.log.metaboliteUsages(2, i) = sum(max(0, -tmp2), 2);
                this.log.metaboliteProductions(1, i) = sum(max(0, tmp1), 2);
                this.log.metaboliteProductions(2, i) = sum(max(0, tmp2), 2);
            end

            this.log.rnas(1, i) = sum(r.counts(:));
            this.log.rnas(2, i) = sum(...
                + r.counts(r.matureIndexs(r.matureMRNAIndexs), iCytosol)...
                + r.counts(r.misfoldedIndexs(r.matureMRNAIndexs), iCytosol)...
                + r.counts(r.damagedIndexs(r.matureMRNAIndexs), iCytosol)...
                );
            this.log.rnas(3, i) = sum(...
                + r.counts(r.matureIndexs(r.matureRRNAIndexs), iCytosol)...
                + r.counts(r.misfoldedIndexs(r.matureRRNAIndexs), iCytosol)...
                + r.counts(r.damagedIndexs(r.matureRRNAIndexs), iCytosol)...
                );
            this.log.rnas(4, i) = sum(...
                + r.counts(r.matureIndexs(r.matureSRNAIndexs), iCytosol) ...
                + r.counts(r.boundIndexs(r.matureSRNAIndexs), iCytosol) ...
                + r.counts(r.aminoacylatedIndexs(r.matureSRNAIndexs), iCytosol) ...
                + r.counts(r.misfoldedIndexs(r.matureSRNAIndexs), iCytosol) ...
                + r.counts(r.damagedIndexs(r.matureSRNAIndexs), iCytosol) ...
                );
            this.log.rnas(5, i) = sum(...
                + r.counts(r.matureIndexs(r.matureTRNAIndexs), iCytosol) ...
                + r.counts(r.aminoacylatedIndexs(r.matureTRNAIndexs), iCytosol) ...
                + r.counts(r.misfoldedIndexs(r.matureTRNAIndexs), iCytosol) ...
                + r.counts(r.damagedIndexs(r.matureTRNAIndexs), iCytosol) ...
                );
            this.log.rnas(6, i) = this.log.rnas(1, i) - sum(this.log.rnas(2:5, i));
            this.log.rnas(7, i) = sum(r.counts(r.damagedIndexs, iCytosol));

            this.log.proteins(1, i) = sum(pm.counts(:));
            this.log.proteins(2, i) = sum(sum(...
                + pm.counts(pm.matureIndexs, :) ...
                + pm.counts(pm.boundIndexs, :) ...
                + pm.counts(pm.inactivatedIndexs, :) ...
                + pm.counts(pm.damagedIndexs, :) ...
                + pm.counts(pm.misfoldedIndexs, :) ...
                ));
            this.log.proteins(3, i) = this.log.proteins(1, i) - this.log.proteins(2, i);
            this.log.proteins(4, i) = sum(pc.counts(:));
            this.log.proteins(5, i) = sum(sum(...
                + pc.counts(pc.matureIndexs, :) ...
                + pc.counts(pc.boundIndexs, :) ...
                + pc.counts(pc.inactivatedIndexs, :) ...
                + pc.counts(pc.damagedIndexs, :) ...
                + pc.counts(pc.misfoldedIndexs, :) ...
                ));
            this.log.proteins(6, i) = this.log.proteins(4, i) - this.log.proteins(5, i);
            this.log.proteins(7, i) = ...
                + sum(sum(pm.counts(pm.damagedIndexs, :))) ...
                + sum(sum(pc.counts(pc.damagedIndexs, :)));

            if this.processIndex_transcription
                transcription = sim.processes{this.processIndex_transcription};
                this.log.rnaPolymerases(1, i) = ...
                    + transcription.enzymes(transcription.enzymeIndexs_rnaPolymerase) ...
                    + transcription.enzymes(transcription.enzymeIndexs_rnaPolymeraseHoloenzyme) ...
                    + transcription.boundEnzymes(transcription.enzymeIndexs_rnaPolymerase) ...
                    + transcription.boundEnzymes(transcription.enzymeIndexs_rnaPolymeraseHoloenzyme);
            end

            if this.processIndex_translation
                translation = sim.processes{this.processIndex_translation};
                this.log.ribosomes(1, i) = ...
                    + translation.enzymes(translation.enzymeIndexs_ribosome30S) ...
                    + translation.enzymes(translation.enzymeIndexs_ribosome30SIF3) ...
                    + translation.enzymes(translation.enzymeIndexs_ribosome50S) ...
                    + translation.enzymes(translation.enzymeIndexs_ribosome70S) ...
                    + translation.boundEnzymes(translation.enzymeIndexs_ribosome30S) ...
                    + translation.boundEnzymes(translation.enzymeIndexs_ribosome30SIF3) ...
                    + translation.boundEnzymes(translation.enzymeIndexs_ribosome50S) ...
                    + translation.boundEnzymes(translation.enzymeIndexs_ribosome70S);
            end

            this.log.supercoiled(1, i) = collapse(c.supercoiled) / 2;
            this.log.supercoiled(2, i) = (collapse(c.linkingNumbers) - collapse(c.polymerizedRegions)/c.relaxedBasesPerTurn) / (collapse(c.polymerizedRegions)/c.relaxedBasesPerTurn);

            if this.processIndex_replicationInitiation
                repInit = sim.processes{this.processIndex_replicationInitiation};

                this.log.dnaA(:, i) = [
                    repInit.enzymes(repInit.enzymeIndexs_DnaA_1mer_ATP)
                    [1 1:7 1:7] * (...
                    + repInit.enzymes([repInit.enzymeIndexs_DnaA; repInit.enzymeIndexs_DnaA_Nmer_ATP; repInit.enzymeIndexs_DnaA_Nmer_ADP]) ...
                    + repInit.boundEnzymes([repInit.enzymeIndexs_DnaA; repInit.enzymeIndexs_DnaA_Nmer_ATP; repInit.enzymeIndexs_DnaA_Nmer_ADP]) ...
                    )];
                dnaAPol = [
                    repInit.calculateDnaAR1234Polymerization()
                    repInit.calcuateIsDnaAR5Occupied()
                    ];
                this.log.dnaAFunctionalBoxes(:, i) = dnaAPol(:, 1);
            end
            if this.processIndex_replication
                rep = sim.processes{this.processIndex_replication};

                this.log.replisome(:, i) = [
                    rep.helicasePosition rep.leadingPosition rep.laggingPosition
                    ];
            end

            this.log.ploidy(1,i) = c.ploidy;

            if this.processIndex_ftsZPolymerization
                ftsZPol = sim.processes{this.processIndex_ftsZPolymerization};

                freeFtsZMonomers = ftsZPol.enzymes(ftsZPol.enzymeIndexs_FtsZ);
                totFtsZMonomers = ...
                    + ftsZPol.enzymes(ftsZPol.enzymeIndexs_FtsZ) ...
                    + ftsZPol.enzymes(ftsZPol.enzymeIndexs_FtsZ_GDP) ...
                    + ftsZPol.enzymes(ftsZPol.enzymeIndexs_FtsZ_activated)' * (1:9)' ...
                    + ftsZPol.boundEnzymes(ftsZPol.enzymeIndexs_FtsZ) ...
                    + ftsZPol.boundEnzymes(ftsZPol.enzymeIndexs_FtsZ_GDP) ...
                    + ftsZPol.boundEnzymes(ftsZPol.enzymeIndexs_FtsZ_activated)' * (1:9)';
                this.log.ftsZ(:, i) = [
                    freeFtsZMonomers
                    totFtsZMonomers
                    ];
            end
            this.log.ftsZRing(:, i) = [
                ftsZRing.numEdgesOneStraight
                ftsZRing.numEdgesTwoStraight
                ftsZRing.numEdgesTwoBent
                ftsZRing.numResidualBent
                ];
            this.log.pinchedDiameter(1, i) = geometry.pinchedDiameter;
        end

        function this = finalize(this, sim)
            t = sim.states{this.stateIndex_time};
            tmpi = t.values / this.stepSizeSec + 1;
            i = ceil(tmpi);

            %trim log
            fields = fieldnames(this.log);
            for j = 1:numel(fields)
                this.log.(fields{j}) = this.log.(fields{j})(:, 1:i);
            end

            %append
            this.appendHelper(i, sim);

            %print and plot
            if rem(tmpi, 1) ~= 0
                this.print(i);
            end
            this.finalizePrint();
            this.plot(i);

            %save
            this.save();
        end

        function this = save(this)
            if isempty(this.outputDirectory)
                return;
            end

            log = this.log; %#ok<PROP,NASGU>
            save([this.outputDirectory filesep 'summary.mat'], '-struct', 'log');

            if this.verbosity > 1
                try
                    saveas(this.figureHandle, [this.outputDirectory filesep 'summary.pdf']);
                catch exception
                    warning('WholeCell:warning', 'Summary figure could not be saved: %s.', exception.message);
                end
            end
        end
    end

    methods (Access = protected)
        function initializePrint(this)
            if this.verbosity <= 0
                return;
            end

            fprintf('%-5s %-3s %-4s %-10s %-89s %-29s %-11s %-11s %-3s %-7s %-9s %-41s %-6s %-7s %-11s %-10s\n', ...
                '', ...
                '', ...
                '', ...
                '', ...
                'Metabolites',...
                'RNAs', ...
                'Monomers', 'Complexes', ...
                '', '', '', ...
                'Replisome', ...
                '', ...
                '', '', '' ...
                );
            fprintf('%-5s %-3s %-4s %-10s %-89s %-29s %-11s %-11s %-3s %-7s %-9s %-41s %-6s %-7s %-11s %-10s\n', ...
                '', ...
                '', ...
                '', ...
                '', ...
                repmat('=', 1, 89),...
                repmat('=', 1, 29), ...
                repmat('=', 1, 11), repmat('=', 1, 11), ...
                '', '', '', ...
                repmat('=', 1, 41), ...
                '', ...
                '', '', '' ...
                );
            fprintf('%-5s %-3s %-4s %-10s %-8s %-8s %-8s %-8s %-9s %-8s %-9s %-8s %-15s %-4s %-4s %-4s %-4s %-4s %-4s %-5s %-5s %-5s %-5s %-3s %-7s %-9s %-13s %-13s %-13s %-6s %-7s %-11s %-10s\n',...
                'Time', ...
                'RT', ...
                'Mass', ...
                'Growth', ...
                'ATP', 'ADP', 'AMP', 'NTPs', 'Min NTP', 'AAs', 'Min AA', 'Lipids', 'Polysaccharides', ...
                'Tot', 'Nasc', 'mRNA', 'rRNA', 'sRNA', 'tRNA', ...
                'Nasct', 'Matur', 'Nasct', 'Matur', ...
                'Sup', 'DnaA', 'R1-5', ...
                'Helicase Pos', 'Leading Pos', 'Lagging Pos', ...
                'Ploidy', ...
                'FtsZ', 'FtsZ Ring', 'Pinch Diam' ...
                );
            fprintf('%-5s %-3s %-4s %-10s %-8s %-8s %-8s %-8s %-9s %-8s %-9s %-8s %-15s %-4s %-4s %-4s %-4s %-4s %-4s %-5s %-5s %-5s %-5s %-3s %-7s %-9s %-13s %-13s %-13s %-6s %-7s %-11s %-10s\n',...
                repmat('=', 1, 5), ...
                repmat('=', 1, 3), ...
                repmat('=', 1, 4), ...
                repmat('=', 1, 10), ...
                repmat('=', 1, 8), repmat('=', 1, 8), repmat('=', 1, 8), repmat('=', 1, 8), repmat('=', 1, 9), repmat('=', 1, 8), repmat('=', 1, 9), repmat('=', 1, 8), repmat('=', 1, 15), ...
                repmat('=', 1, 4), repmat('=', 1, 4), repmat('=', 1, 4), repmat('=', 1, 4), repmat('=', 1, 4), repmat('=', 1, 4), ...
                repmat('=', 1, 5), repmat('=', 1, 5), repmat('=', 1, 5), repmat('=', 1, 5), ...
                repmat('=', 1, 3), repmat('=', 1, 7), repmat('=', 1, 9), ...
                repmat('=', 1, 13), repmat('=', 1, 13), repmat('=', 1, 13), ...
                repmat('=', 1,6), ...
                repmat('=', 1, 7), repmat('=', 1, 11), repmat('=', 1, 10) ...
                );
        end

        function print(this, i)
            if this.verbosity <= 0
                return;
            end

            if ispc
                format = '%5d %3.0f %4.2f %5.3e %8d %8d %8d %8d %5d %3s %8d %5d %3s %8d %15d %4d %4d %4d %4d %4d %4d %5d %5d %5d %5d %3d %3d %3d %d %d %d %d %d %6d %6d %6d %6d %6d %6d %6.3f %3d %3d %2d %2d %2d %2d %5.3e\n';
            else
                format = '%5d %3.0f %4.2f %6.3e %8d %8d %8d %8d %5d %3s %8d %5d %3s %8d %15d %4d %4d %4d %4d %4d %4d %5d %5d %5d %5d %3d %3d %3d %d %d %d %d %d %6d %6d %6d %6d %6d %6d %6.3f %3d %3d %2d %2d %2d %2d %6.3e\n';
            end
            fprintf(format,...
                this.log.time(1, i),...
                this.log.runTime(1, i), ...
                this.log.mass(1, i), ...
                this.log.growth(1, i), ...
                this.log.metabolites(1, i),...
                this.log.metabolites(2, i),...
                this.log.metabolites(3, i),...
                this.log.metabolites(4, i),...
                this.log.metabolites(5, i),...
                this.ntpWholeCellModelIDs{this.log.metabolites(6, i)},...
                this.log.metabolites(7, i),...
                this.log.metabolites(8, i),...
                this.aaWholeCellModelIDs{this.log.metabolites(9, i)},...
                this.log.metabolites(10,i),...
                this.log.metabolites(11,i),...
                this.log.rnas(1, i),...
                this.log.rnas(6, i),...
                this.log.rnas(2, i),...
                this.log.rnas(3, i),...
                this.log.rnas(4, i),...
                this.log.rnas(5, i),...
                this.log.proteins(3, i),...
                this.log.proteins(2, i),...
                this.log.proteins(6, i),...
                this.log.proteins(5, i),...
                this.log.supercoiled(1, i), ...
                this.log.dnaA(1, i), ...
                this.log.dnaA(2, i), ...
                this.log.dnaAFunctionalBoxes(1, i), ...
                this.log.dnaAFunctionalBoxes(2, i), ...
                this.log.dnaAFunctionalBoxes(3, i), ...
                this.log.dnaAFunctionalBoxes(4, i), ...
                this.log.dnaAFunctionalBoxes(5, i), ...
                this.log.replisome(1, i), ...
                this.log.replisome(2, i), ...
                this.log.replisome(3, i), ...
                this.log.replisome(4, i), ...
                this.log.replisome(5, i), ...
                this.log.replisome(6, i), ...
                this.log.ploidy(1,i), ...
                this.log.ftsZ(1, i), ...
                this.log.ftsZ(2, i), ...
                this.log.ftsZRing(1, i), ...
                this.log.ftsZRing(2, i), ...
                this.log.ftsZRing(3, i), ...
                this.log.ftsZRing(4, i), ...
                this.log.pinchedDiameter(1, i) ...
                );
        end

        function finalizePrint(this)
            if this.verbosity <= 0
                return;
            end

            fprintf('Total simulation time: %d s\n', this.log.time(1, end));
            fprintf('Total runtime: %.3f s\n', sum(this.log.runTime));
            fprintf('\n');
        end

        function initializePlot(this)
            if this.verbosity <= 1
                return;
            end

            %create new figure
            this.figureHandle = figure();
        end

        %Plot progress of evolution of state of organism.
        function plot(this, ~)
            if this.verbosity <= 1
                return;
            end

            %run time
            axesHandle = subplot(2, 4, 1, 'Parent', this.figureHandle);
            plot(axesHandle, 0:numel(this.log.time)-1, this.log.runTime/60);
            title(axesHandle, 'Runtime')
            xlabel(axesHandle, 'Segment');
            ylabel(axesHandle, 'Run Time (m)');
            xlim(axesHandle, [0 numel(this.log.time)-1]);

            %metabolites
            axesHandle = subplot(2, 4, 2, 'Parent', this.figureHandle);
            h = plot(axesHandle, this.log.time / 3600, this.log.metabolites([1 2 3 4 5 7 8], :)');
            title(axesHandle, 'Metabolites')
            xlabel(axesHandle, 'Time (h)');
            ylabel(axesHandle, 'Metabolites');
            xlim(axesHandle, [0 this.log.time(end)] / 3600);
            legend(h, {'ATP', 'ADP', 'AMP', 'Tot NTPs', 'Min NTP', 'Tot AAs', 'Min AA'});

            %RNAs
            axesHandle = subplot(2, 4, 3, 'Parent', this.figureHandle);
            h = plot(axesHandle, this.log.time / 3600, this.log.rnas([1 6 2 3 4 5], :)');
            title(axesHandle, 'RNA')
            xlabel(axesHandle, 'Time (h)');
            ylabel(axesHandle, 'RNA');
            xlim(axesHandle, [0 this.log.time(end)] / 3600);
            legend(h, {'Total', 'Nascent', 'mRNA', 'rRNA', 'sRNA', 'tRNA'});

            %proteins
            axesHandle = subplot(2, 4, 4, 'Parent', this.figureHandle);
            h = plot(axesHandle, this.log.time/3600, this.log.proteins([1 3 2 4 6 5], :)');
            title(axesHandle, 'Protein')
            xlabel(axesHandle, 'Time (h)');
            ylabel(axesHandle, 'Protein');
            xlim(axesHandle, [0 this.log.time(end)]/3600);
            legend(h, {'Total Monomers', 'Nascent Monomers', 'Mature Monomers', 'Total Complexes', 'Nascent Complexes', 'Mature Complexes'});

            %Supercoiled
            axesHandle = subplot(2, 4, 5, 'Parent', this.figureHandle);
            plot(axesHandle, this.log.time/3600, this.log.supercoiled(1, :)');
            title(axesHandle, 'Supercoiled')
            xlabel(axesHandle, 'Time (h)');
            ylabel(axesHandle, 'No. Supercoiled Regions');
            xlim(axesHandle, [0 this.log.time(end)]/3600);

            %DnaA
            axesHandle = subplot(2, 4, 6, 'Parent', this.figureHandle);
            h = plot(axesHandle, this.log.time/3600, [this.log.dnaA/max(1, this.log.dnaA(1,1)); int16(this.log.dnaAFunctionalBoxes)]');
            title(axesHandle, 'DnaA')
            xlabel(axesHandle, 'Time (h)');
            ylabel(axesHandle, 'DnaA Count');
            xlim(axesHandle, [0 this.log.time(end)]/3600);
            legend(h, {'Free DnaA-ATP', 'Total DnaA', 'R1', 'R2', 'R3', 'R4', 'R5'});

            %Replisome
            axesHandle = subplot(2, 4, 7, 'Parent', this.figureHandle);
            h = plot(axesHandle, this.log.time/3600, this.log.replisome');
            title(axesHandle, 'Replisome')
            xlabel(axesHandle, 'Time (h)');
            ylabel(axesHandle, 'Position');
            xlim(axesHandle, [0 this.log.time(end)]/3600);
            legend(h, {'Helicase-1', 'Helicase-2', 'Leading-1', 'Leading-2', 'Lagging-1', 'Lagging-2'});

            %FtsZ
            axesHandle = subplot(2, 4, 8, 'Parent', this.figureHandle);
            h = plot(axesHandle, this.log.time/3600, [double(this.log.ftsZ/this.log.ftsZ(1, 1)); double(this.log.ftsZRing); this.log.pinchedDiameter/this.log.pinchedDiameter(1)]');
            title(axesHandle, 'FtsZ')
            xlabel(axesHandle, 'Time (h)');
            ylabel(axesHandle, 'FtsZ');
            xlim(axesHandle, [0 this.log.time(end)]/3600);
            legend(h, {'Free FtsZ-GTP', 'Total FtsZ', 'Ring-1strt', 'Ring-1bent', 'Ring-2strt', 'Ring-2 bent', 'Pinched Diameter'});

            %force display of figures
            drawnow;
        end
    end

    % Analyzes summary output of a batch of simulations. Stores analysis in mat file.
    methods (Static = true)
        function [simStats] = summarizeSimulations(baseDir)
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;
            import edu.stanford.covert.cell.sim.util.SummaryLogger;

            simStats = SummaryLogger.getSimulationStatistics(baseDir);

            nSimulations = size(simStats, 1);

            summaryFileName = [baseDir filesep 'summary.html'];
            stateFileName = [baseDir filesep 'initialAndFinalStates.html'];
            errorFileName = [baseDir filesep 'errors.html'];

            %% calc statistics
            stateFid = fopen(stateFileName, 'w');
            errorFid = fopen(errorFileName, 'w');
            fprintf(stateFid, '<html>\n<head>\n');
            fprintf(stateFid, '<style type="text/css">\n');
            fprintf(stateFid, '  table.states {width:12000px;}\n');
            fprintf(stateFid, '  table.states thead tr th, table.states tbody tr td{width:auto;padding-left:10px;}\n');
            fprintf(stateFid, '  table.states thead tr th:first-child, table.states tbody tr td:first-child{width:auto;padding-left:0px;}\n');
            fprintf(stateFid, '  table.states thead tr:first-child th {border-top:1px solid black;}\n');
            fprintf(stateFid, '  table.states thead tr:last-child th {border-bottom:1px solid black;}\n');
            fprintf(stateFid, '  table.states tbody tr:last-child td {border-bottom:1px solid black;}\n');
            fprintf(stateFid, '  table.states tr th:first-child, table.states tr td:first-child {border-right:1px solid black}\n');
            fprintf(stateFid, '  table.states thead tr:last-child th:nth-child(%d), table.states tbody tr td:nth-child(%d) {border-right:1px solid black}\n', 1+6, 1+6);
            fprintf(stateFid, '  table.states thead tr:last-child th:nth-child(%d), table.states tbody tr td:nth-child(%d) {border-right:1px solid black}\n', 1+6+48, 1+6+48);
            fprintf(stateFid, '  table.states thead tr th:first-child, table.states tbody tr td:first-child {border-left:1px solid black}\n');
            fprintf(stateFid, '  table.states thead tr th:last-child, table.states tbody tr td:last-child {border-right:1px solid black}\n');
            fprintf(stateFid, '  table.states thead tr th {background-color: #CCCCCC;}\n');
            fprintf(stateFid, '  table.states tbody tr:nth-child(even) td {background-color: #DDDDDD;}\n');
            fprintf(stateFid, '  table.states tbody tr:nth-child(odd) td {background-color: #EEEEEE;}\n');
            fprintf(stateFid, '</style>\n');
            fprintf(stateFid, '</head>\n');
            fprintf(stateFid, '<body>\n<table class="states" cellpadding="0" cellspacing="0">\n<thead>\n%s\n</thead>', SummaryLogger.formatColNames());
            fprintf(errorFid, '<html>\n<body>\n');
            for i = 1:nSimulations
                %calc simulation status
                if ~exist(sprintf('%s/%d/state-0.mat', baseDir, i), 'file')
                elseif ~exist(sprintf('%s/%d/err.mat', baseDir, i), 'file') && ~exist(sprintf('%s/%d/summary.mat', baseDir, i), 'file')
                elseif exist(sprintf('%s/%d/err.mat', baseDir, i), 'file')
                    tmpFid = fopen(sprintf('%s/%d/err.log', baseDir, i), 'r');
                    error = reshape(fread(tmpFid, '*char'), 1, []);
                    error = strrep(error, sprintf('\n\n\n\n'), sprintf('\n'));
                    error = strrep(error, sprintf('\n\n'), sprintf('\n'));
                    fclose(tmpFid);
                    fprintf(errorFid, '<p>\n<b>Simulation #%d</b><br/>\n<pre>\n%s\n</pre>\n</p>\n', i, error);
                else
                    summary = load(sprintf('%s/%d/summary.mat', baseDir, i));
                end

                if simStats(i, SummaryLogger.SIM_STATUS_INDEX_STATUS) >= 0
                    fprintf(stateFid, '<tr>\n<td>%d</td>\n%s\n%s\n%s\n</tr>\n', i, SummaryLogger.formatStats(simStats(i, :)), SummaryLogger.formatState(summary, 1), SummaryLogger.formatState(summary, numel(summary.time)));
                else
                    fprintf(stateFid, '<tr>\n<td>%d</td>\n%s\n%s\n%s\n</tr>\n', i, SummaryLogger.formatStats(simStats(i, :)), SummaryLogger.formatState(), SummaryLogger.formatState());
                end
            end
            fprintf(stateFid, '</tbody>\n</table>\n</body>\n</html>\n');
            fprintf(errorFid, '</body>\n</html>\n');
            fclose(stateFid);
            fclose(errorFid);

            if ~any(simStats(:, SummaryLogger.SIM_STATUS_INDEX_STATUS) == -1)
                errorFid = fopen(errorFileName, 'w');
                fprintf(errorFid, '<html>\n<body>\n');
                fprintf(errorFid, 'No run-time errors were encountered.\n');
                fprintf(errorFid, '</body>\n</html>\n');
                fclose(errorFid);
            end

            fid = fopen(summaryFileName, 'w');
            fwrite(fid, SummaryLogger.getSummaryReport(simStats));
            fclose(fid);

            %% plot
            figHandle = figure();

            %initial vs. cumulative growth
            plot(simStats(:, SummaryLogger.SIM_STATUS_INDEX_INIT_GROWTH)'*3600, ...
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_CUM_GROWTH)', '.');
            xlabel('Initial growth rate (cell/h)');
            ylabel('Cumulative growth rate (cell)');
            saveas(gcf, [baseDir filesep 'initialVsCumulativeGrowth.pdf']);

            %replication initiation time vs replication, cytokinesis times
            h = plot(simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_INIT_TIME)' / 3600, [...
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_TIME) ...
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_CYTOKINESIS_TIME)]'/ 3600, '.');
            xlabel('Replication Initiation Time (h)');
            ylabel('Time (h)');
            legend(h, {'Replication', 'Cytokinesis'});
            saveas(gcf, [baseDir filesep 'cellCyclePhaseLengths.pdf']);

            %replication initiation duration vs replication, cytokinesis
            %durations
            h = plot(simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_INIT_TIME)' / 3600, [...
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_TIME)-simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_INIT_TIME) ...
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_CYTOKINESIS_TIME)-simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_TIME)]'/ 3600, '.');
            xlabel('Replication Initiation Time (h)');
            ylabel('Duration (h)');
            legend(h, {'Replication', 'Cytokinesis'});
            saveas(gcf, [baseDir filesep 'cellCyclePhaseDurations.pdf']);

            %initial growth rate vs replication initiation, replication,
            %cytokinesis times
            h = plot(simStats(:, SummaryLogger.SIM_STATUS_INDEX_INIT_GROWTH)'*3600, [...
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_INIT_TIME) ...
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_TIME) ...
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_CYTOKINESIS_TIME)]'/ 3600, '.');
            xlabel('Initial growth rate (cell/h)');
            ylabel('Time (h)');
            legend(h, {'Replication initiation', 'Replication', 'Cytokinesis'});
            saveas(gcf, [baseDir filesep 'initialGrowthVsCellCyclePhaseTimes.pdf']);

            %cumulative growth vs replication initiation, replication,
            %cytokinesis times
            h = plot(simStats(:, SummaryLogger.SIM_STATUS_INDEX_CUM_GROWTH), [...
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_INIT_TIME) ...
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_TIME) simStats(:, SummaryLogger.SIM_STATUS_INDEX_CYTOKINESIS_TIME)]'/ 3600, '.');
            xlabel('Cumulative growth rate (cell)');
            ylabel('Time (h)');
            legend(h, {'Replication initiation', 'Replication', 'Cytokinesis'});
            saveas(gcf, [baseDir filesep 'cumulativeGrowthVsCellCyclePhaseTimes.pdf']);

            %initial growth rate vs replication initiation, replication,
            %cytokinesis durations
            h = plot(simStats(:, SummaryLogger.SIM_STATUS_INDEX_INIT_GROWTH)'*3600, [...
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_INIT_TIME) ....
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_TIME)-simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_INIT_TIME) ...
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_CYTOKINESIS_TIME)-simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_TIME)]'/ 3600, '.');
            xlabel('Initial growth rate (cell/h)');
            ylabel('Duration (h)');
            legend(h, {'Replication initiation', 'Replication', 'Cytokinesis'});
            saveas(gcf, [baseDir filesep 'initialGrowthVsCellCyclePhaseDurations.pdf']);

            %cumulative growth vs replication initiation, replication,
            %cytokinesis durations
            h = plot(simStats(:, SummaryLogger.SIM_STATUS_INDEX_CUM_GROWTH), [...
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_INIT_TIME) ...
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_TIME)-simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_INIT_TIME) ...
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_CYTOKINESIS_TIME)-simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_TIME)]'/ 3600, '.');
            xlabel('Cumulative growth rate (cell)');
            ylabel('Duration (h)');
            legend(h, {'Replication initiation', 'Replication', 'Cytokinesis'});
            saveas(gcf, [baseDir filesep 'cumulativeGrowthVsCellCyclePhaseDurations.pdf']);

            close(figHandle);
        end

        function simStats = getSimulationStatistics(baseDir, selectedSimulations)
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;
            import edu.stanford.covert.cell.sim.util.SummaryLogger;

            if nargin < 2
                selectedSimulations = 1:SimulationDiskUtil.getNumSimulations(baseDir);
            end

            nSimulations = numel(selectedSimulations);

            %% calc statistics
            simStats = NaN(nSimulations, 11);
            for j = 1:nSimulations
                i = selectedSimulations(j);

                %calc simulation status
                if ~exist(sprintf('%s/%d/state-0.mat', baseDir, i), 'file')
                    simStats(j, SummaryLogger.SIM_STATUS_INDEX_STATUS) = SummaryLogger.SIM_STATUS_DIDNT_START;
                elseif ~exist(sprintf('%s/%d/err.mat', baseDir, i), 'file') && ~exist(sprintf('%s/%d/summary.mat', baseDir, i), 'file')
                    simStats(j, SummaryLogger.SIM_STATUS_INDEX_STATUS) = SummaryLogger.SIM_STATUS_STILL_RUNNING;
                elseif exist(sprintf('%s/%d/err.mat', baseDir, i), 'file')
                    simStats(j, SummaryLogger.SIM_STATUS_INDEX_STATUS) = SummaryLogger.SIM_STATUS_TERMINATED_WITH_ERROR;
                else
                    summary = load(sprintf('%s/%d/summary.mat', baseDir, i), 'time', 'mass', 'growth', 'dnaAFunctionalBoxes', 'replisome', 'pinchedDiameter');
                    if summary.pinchedDiameter(end) == 0
                        simStats(j, SummaryLogger.SIM_STATUS_INDEX_STATUS) = SummaryLogger.SIM_STATUS_COMPLETED_WITH_DIVISION;
                        simStats(j, SummaryLogger.SIM_STATUS_INDEX_CYTOKINESIS_TIME) = summary.time(end);
                    else
                        simStats(j, SummaryLogger.SIM_STATUS_INDEX_STATUS) = SummaryLogger.SIM_STATUS_COMPLETED_WITHOUT_DIVISION;
                    end

                    %calc simulation statistics
                    idx = find(all(summary.dnaAFunctionalBoxes == repmat([7; 7; 7; 7; 1], 1, numel(summary.time)), 1), 1, 'first');
                    if ~isempty(idx)
                        simStats(j, SummaryLogger.SIM_STATUS_INDEX_REP_INIT_TIME) = summary.time(idx);
                    end

                    idx = find(any(summary.replisome, 1), 1, 'first');
                    if ~isempty(idx)
                        simStats(j, SummaryLogger.SIM_STATUS_INDEX_INIT_REP_TIME) = summary.time(idx);
                    end
                    simStats(j, SummaryLogger.SIM_STATUS_INDEX_REP_INIT_TIME) = min(...
                        simStats(j, SummaryLogger.SIM_STATUS_INDEX_REP_INIT_TIME), ...
                        simStats(j, SummaryLogger.SIM_STATUS_INDEX_INIT_REP_TIME));

                    idx = find(any(summary.replisome, 1), 1, 'last')+1;
                    if ~isempty(idx) && idx <= size(summary.time, 2)
                        simStats(j, SummaryLogger.SIM_STATUS_INDEX_REP_TIME) = summary.time(idx);
                    end

                    simStats(j, SummaryLogger.SIM_STATUS_INDEX_INIT_GROWTH) = summary.growth(1);
                    cellCycleLength = 9.0 * 3600;
                    idx = find(summary.time > cellCycleLength, 1, 'first');
                    if ~isempty(idx)
                        simStats(j, SummaryLogger.SIM_STATUS_INDEX_FIN_GROWTH) = summary.growth(idx);
                        simStats(j, SummaryLogger.SIM_STATUS_INDEX_CUM_GROWTH) = diff(summary.time(1:2)) * sum(summary.growth(1:idx));
                    else
                        f = (2 - 1) / (exp(log(2) * summary.time(end) / cellCycleLength) - 1); %#ok<CPROP,PROP>
                        simStats(j, SummaryLogger.SIM_STATUS_INDEX_FIN_GROWTH) = f * summary.growth(end);
                        simStats(j, SummaryLogger.SIM_STATUS_INDEX_CUM_GROWTH) = f * diff(summary.time(1:2)) * sum(summary.growth);
                    end

                    idx = find(summary.mass > 2 * summary.mass(1), 1, 'first');
                    if ~isempty(idx)
                        simStats(j, SummaryLogger.SIM_STATUS_INDEX_MASS_DOUBLING_TIME) = summary.time(idx-1);
                    else
                        simStats(j, SummaryLogger.SIM_STATUS_INDEX_MASS_DOUBLING_TIME) = ...
                            cellCycleLength / log(2) * log(2 / summary.mass(end)) + summary.time(end); %#ok<CPROP,PROP>
                    end

                    simStats(j, SummaryLogger.SIM_STATUS_INDEX_FIN_MASS) = summary.mass(end);
                    simStats(j, SummaryLogger.SIM_STATUS_INDEX_INIT_MASS) = summary.mass(1);
                end
            end
        end

        function value = formatColNames()
            stateColNames = {
                'Time'
                'Growth'
                'Mass'
                'ATP'
                'ADP'
                'AMP'
                'Total NTPs'
                'Min NTP'
                'Min NTP'
                'Total AAs'
                'Min AA'
                'Min AA'
                'Total Lipids'
                'Total Carbohydrates'
                'Total RNA'
                'Total mature mRNA'
                'Total mature rRNA'
                'Total mature sRNA'
                'Total mature tRNA'
                'Total immature RNA'
                'Total Monomers'
                'Total Mature Monomers'
                'Total Immature Monomers'
                'Total Complexs'
                'Total Mature Complexs'
                'Total Immature Complexs'
                'Supercoiled'
                'Free'
                'Total'
                'R1'
                'R2'
                'R3'
                'R4'
                'R5'
                'Helicase Position - 1'
                'Helicase Position - 2'
                'Leading Polymerase Position - 1'
                'Leading Polymerase Position - 2'
                'Lagging Polymerase Position - 2'
                'Lagging Polymerase Position - 2'
                'Ploidy'
                'Free'
                'Total'
                'One straight'
                'One bent'
                'Two straight'
                'Residual'
                'Pinched Diameter'
                };
            colNames = [{
                'Simulation'
                'Status'
                'Replication initiation time'
                'Initiation of replication time'
                'Replication time'
                'Cytokinesis time'
                'Cumulative growth'
                };
                stateColNames; stateColNames];

            value = [];

            value = [value sprintf('<tr>\n')];
            value = [value sprintf('<th> </th>\n')];
            value = [value sprintf('<th colspan="6">%s</th>\n', 'Summary')];
            value = [value sprintf('<th colspan="%d">%s</th>\n', numel(stateColNames), 'Initial State')];
            value = [value sprintf('<th colspan="%d">%s</th>\n', numel(stateColNames), 'Final State')];
            value = [value sprintf('</tr>\n')];

            value = [value sprintf('<tr>\n')];
            value = [value sprintf('<th> </th>\n')];
            value = [value sprintf('<th colspan="6">%s</th>\n', ' ')];
            value = [value sprintf('<th colspan="%d">%s</th>\n', 1, ' ')]; %Time
            value = [value sprintf('<th colspan="%d">%s</th>\n', 1, ' ')]; %Growth
            value = [value sprintf('<th colspan="%d">%s</th>\n', 1, ' ')]; %Mass
            value = [value sprintf('<th colspan="%d">%s</th>\n', 11, 'Metabolite')];
            value = [value sprintf('<th colspan="%d">%s</th>\n', 6, 'RNA')];
            value = [value sprintf('<th colspan="%d">%s</th>\n', 6, 'Protein')];
            value = [value sprintf('<th colspan="%d">%s</th>\n', 1, ' ')]; %Supercoiled
            value = [value sprintf('<th colspan="%d">%s</th>\n', 2, 'DnaA')];
            value = [value sprintf('<th colspan="%d">%s</th>\n', 5, 'Functional DnaA Boxes')];
            value = [value sprintf('<th colspan="%d">%s</th>\n', 6, 'Replisome')];
            value = [value sprintf('<th colspan="%d">%s</th>\n', 1, ' ')]; %ploidy
            value = [value sprintf('<th colspan="%d">%s</th>\n', 2, 'FtsZ')];
            value = [value sprintf('<th colspan="%d">%s</th>\n', 4, 'FtsZ Ring')];
            value = [value sprintf('<th colspan="%d">%s</th>\n', 1, ' ')]; %Pinched Diameter
            value = [value sprintf('<th colspan="%d">%s</th>\n', 1, ' ')]; %Time
            value = [value sprintf('<th colspan="%d">%s</th>\n', 1, ' ')]; %Growth
            value = [value sprintf('<th colspan="%d">%s</th>\n', 1, ' ')]; %Mass
            value = [value sprintf('<th colspan="%d">%s</th>\n', 11, 'Metabolite')];
            value = [value sprintf('<th colspan="%d">%s</th>\n', 6, 'RNA')];
            value = [value sprintf('<th colspan="%d">%s</th>\n', 6, 'Protein')];
            value = [value sprintf('<th colspan="%d">%s</th>\n', 1, ' ')]; %Supercoiled
            value = [value sprintf('<th colspan="%d">%s</th>\n', 2, 'DnaA')];
            value = [value sprintf('<th colspan="%d">%s</th>\n', 5, 'Functional DnaA Boxes')];
            value = [value sprintf('<th colspan="%d">%s</th>\n', 6, 'Replisome')];
            value = [value sprintf('<th colspan="%d">%s</th>\n', 1, ' ')]; %ploidy
            value = [value sprintf('<th colspan="%d">%s</th>\n', 2, 'FtsZ')];
            value = [value sprintf('<th colspan="%d">%s</th>\n', 4, 'FtsZ Ring')];
            value = [value sprintf('<th colspan="%d">%s</th>\n', 1, ' ')]; %Pinched Diameter
            value = [value sprintf('</tr>\n')];

            value = [value sprintf('<tr>\n')];
            for i = 1:numel(colNames)
                value = [value sprintf('<th>%s</th>\n', colNames{i})]; %#ok<AGROW>
            end
            value = [value sprintf('</tr>\n')];
        end

        function value = formatStats(stats)
            import edu.stanford.covert.cell.sim.util.SummaryLogger;

            switch stats(SummaryLogger.SIM_STATUS_INDEX_STATUS)
                case -3, status = 'Not yet started';
                case -2, status = 'Not yet terminated';
                case -1, status = 'Terminated with error';
                case 0, status = 'Terminated with division';
                case 1, status = 'Terminated without division';
            end
            value = sprintf('<td>%s</td>\n', status);

            if ~isnan(stats(SummaryLogger.SIM_STATUS_INDEX_REP_INIT_TIME))
                value = [value sprintf('<td>%d</td>\n', stats(SummaryLogger.SIM_STATUS_INDEX_REP_INIT_TIME))];
            else
                value = [value sprintf('<td> </td>\n')];
            end

            if ~isnan(stats(SummaryLogger.SIM_STATUS_INDEX_INIT_REP_TIME))
                value = [value sprintf('<td>%d</td>\n', stats(SummaryLogger.SIM_STATUS_INDEX_INIT_REP_TIME))];
            else
                value = [value sprintf('<td> </td>\n')];
            end

            if ~isnan(stats(SummaryLogger.SIM_STATUS_INDEX_REP_TIME))
                value = [value sprintf('<td>%d</td>\n', stats(SummaryLogger.SIM_STATUS_INDEX_REP_TIME))];
            else
                value = [value sprintf('<td> </td>\n')];
            end

            if ~isnan(stats(SummaryLogger.SIM_STATUS_INDEX_CYTOKINESIS_TIME))
                value = [value sprintf('<td>%d</td>\n', stats(SummaryLogger.SIM_STATUS_INDEX_CYTOKINESIS_TIME))];
            else
                value = [value sprintf('<td> </td>\n')];
            end

            if ~isnan(stats(SummaryLogger.SIM_STATUS_INDEX_CUM_GROWTH))
                value = [value sprintf('<td>%e</td>\n', stats(SummaryLogger.SIM_STATUS_INDEX_CUM_GROWTH))];
            else
                value = [value sprintf('<td> </td>\n')];
            end
        end

        function value = formatState(summary, idx)
            fields = {
                'time' 1
                'growth' 1
                'mass' 1
                'metabolites' 11
                'rnas' 6
                'proteins' 6
                'supercoiled' 1
                'dnaA' 2
                'dnaAFunctionalBoxes' 5
                'replisome' 6
                'ploidy' 1
                'ftsZ' 2
                'ftsZRing' 4
                'pinchedDiameter' 1
                };
            value = [];
            for i = 1:size(fields, 1)
                for j = 1:fields{i, 2}
                    if nargin == 0 || ~isfield(summary, fields{i, 1}) || j > size(summary.(fields{i, 1}), 1)
                        value = [value sprintf('<td> </td>\n')]; %#ok<AGROW>
                    elseif summary.(fields{i, 1})(j, idx) == ceil(summary.(fields{i, 1})(j, idx))
                        value = [value sprintf('<td>%d</td>\n', summary.(fields{i, 1})(j, idx))]; %#ok<AGROW>
                    else
                        value = [value sprintf('<td>%e</td>\n', summary.(fields{i, 1})(j, idx))]; %#ok<AGROW>
                    end
                end
            end
        end

        %% build summary report
        function summary = getSummaryReport(simStats)
            import edu.stanford.covert.cell.sim.util.SummaryLogger;

            nSimulations = size(simStats, 1);
            nTerminatedWithDivision = sum(simStats(:, SummaryLogger.SIM_STATUS_INDEX_STATUS) == 0);
            nTerminatedWithoutDivision = sum(simStats(:, SummaryLogger.SIM_STATUS_INDEX_STATUS) == 1);
            nTerminatedWithError = sum(simStats(:, SummaryLogger.SIM_STATUS_INDEX_STATUS) == -1);
            nNotYetTerminated = sum(simStats(:, SummaryLogger.SIM_STATUS_INDEX_STATUS) == -2);
            nNotYetStarted = sum(simStats(:, SummaryLogger.SIM_STATUS_INDEX_STATUS) == -3);
            repInitTimes = simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_INIT_TIME) / 3600;
            initRepTimes = simStats(:, SummaryLogger.SIM_STATUS_INDEX_INIT_REP_TIME) / 3600;
            repTimes = simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_TIME) / 3600;
            cytTimes = simStats(:, SummaryLogger.SIM_STATUS_INDEX_CYTOKINESIS_TIME) / 3600;
            initGrowths = simStats(:, SummaryLogger.SIM_STATUS_INDEX_INIT_GROWTH) * 3600;
            finGrowths = simStats(:, SummaryLogger.SIM_STATUS_INDEX_FIN_GROWTH) * 3600;
            cumGrowths = simStats(:, SummaryLogger.SIM_STATUS_INDEX_CUM_GROWTH);
            initMass = simStats(:, SummaryLogger.SIM_STATUS_INDEX_INIT_MASS);
            finMass = simStats(:, SummaryLogger.SIM_STATUS_INDEX_FIN_MASS);

            %styles
            summary = [];
            summary = [summary sprintf('<html>\n')];
            summary = [summary sprintf('<head>\n')];
            summary = [summary sprintf('<style type="text/css">\n')];
            summary = [summary sprintf('  table.summary tr th {text-align:left; width:400px;}\n')];
            summary = [summary sprintf('  table.summary tr td {text-align:right; padding-left:10px}\n')];
            summary = [summary sprintf('</style>\n')];
            summary = [summary sprintf('</head>\n')];
            summary = [summary sprintf('<body>\n')];

            %status
            summary = [summary sprintf('<table class="summary" cellpadding="0" cellspacing="0">\n')];
            summary = [summary sprintf('<tr><th>%s</th><td>%d</td><td>%5.1f%%</td></tr>\n', 'No. simulations', nSimulations, 100)];
            summary = [summary sprintf('<tr><th>%s</th><td>%d</td><td>%5.1f%%</td></tr>\n', 'Terminated with division', nTerminatedWithDivision, nTerminatedWithDivision/nSimulations*100)];
            summary = [summary sprintf('<tr><th>%s</th><td>%d</td><td>%5.1f%%</td></tr>\n', 'Terminated without division',nTerminatedWithoutDivision, nTerminatedWithoutDivision/nSimulations*100)];
            summary = [summary sprintf('<tr><th>%s</th><td>%d</td><td>%5.1f%%</td></tr>\n', 'Terminated with error', nTerminatedWithError, nTerminatedWithError/nSimulations*100)];
            summary = [summary sprintf('<tr><th>%s</th><td>%d</td><td>%5.1f%%</td></tr>\n', 'Not yet terminated', nNotYetTerminated, nNotYetTerminated/nSimulations*100)];
            summary = [summary sprintf('<tr><th>%s</th><td>%d</td><td>%5.1f%%</td></tr>\n', 'Not yet started', nNotYetStarted, nNotYetStarted/nSimulations*100)];
            summary = [summary sprintf('</table>\n<br/>\n')];

            %times
            meanMassDoublingTime = 9.0 / nanmean(cumGrowths);
            stdMassDoublingTime = sqrt(9.0 / nanmean(cumGrowths)) * nanstd(cumGrowths);
            summary = [summary sprintf('<table class="summary" cellpadding="0" cellspacing="0">\n')];
            summary = [summary sprintf('<tr><th>%s</th><td>%5.1f ± %5.1f</td></tr>\n', 'Mass doubling time (h)', meanMassDoublingTime, stdMassDoublingTime)];
            summary = [summary sprintf('<tr><th>%s</th><td>%5.1f ± %5.1f</td></tr>\n', 'Cell cycle length (h)', nanmean(cytTimes), nanstd(cytTimes))];
            summary = [summary sprintf('<tr><th>%s</th><td>%5.1f ± %5.1f</td></tr>\n', 'Replication initiation time (h)', nanmean(repInitTimes), nanstd(repInitTimes))];
            summary = [summary sprintf('<tr><th>%s</th><td>%5.1f ± %5.1f</td></tr>\n', 'Initiation of replication time (h)', nanmean(initRepTimes), nanstd(initRepTimes))];
            summary = [summary sprintf('<tr><th>%s</th><td>%5.1f ± %5.1f</td></tr>\n', 'Replication time (h)', nanmean(repTimes), nanstd(repTimes))];
            summary = [summary sprintf('</table>\n<br/>\n')];

            %durations
            summary = [summary sprintf('<table class="summary" cellpadding="0" cellspacing="0">\n')];
            summary = [summary sprintf('<tr><th>%s</th><td>%5.1f ± %5.1f</td></tr>\n', 'Replication initiation duration (h)', ...
                nanmean(repInitTimes), nanstd(repInitTimes))];
            summary = [summary sprintf('<tr><th>%s</th><td>%5.1f ± %5.1f</td></tr>\n', 'Replication initiation-replisome formation duration (h)', ...
                nanmean(initRepTimes - repInitTimes), nanstd(initRepTimes - repInitTimes))];
            summary = [summary sprintf('<tr><th>%s</th><td>%5.1f ± %5.1f</td></tr>\n', 'Replication duration (h)', ...
                nanmean(repTimes - initRepTimes), nanstd(repTimes - initRepTimes))];
            summary = [summary sprintf('<tr><th>%s</th><td>%5.1f ± %5.1f</td></tr>\n', 'Cytokinesis duration (h)', ...
                nanmean(cytTimes - repTimes), nanstd(cytTimes - repTimes))];
            summary = [summary sprintf('</table>\n<br/>\n')];

            %growth
            summary = [summary sprintf('<table class="summary" cellpadding="0" cellspacing="0">\n')];
            summary = [summary sprintf('<tr><th>%s</th><td>%5.3f ± %5.3f</td></tr>\n', 'Growth at t==0 (1/cells/h)', ...
                nanmean(initGrowths), nanstd(initGrowths))];
            summary = [summary sprintf('<tr><th>%s</th><td>%5.3f ± %5.3f</td></tr>\n', 'Growth at t==τ (1/cells/h)', ...
                nanmean(finGrowths), nanstd(finGrowths))];
            summary = [summary sprintf('<tr><th>%s</th><td>%5.3f ± %5.3f</td></tr>\n', 'Cumulative Growth over τ (1/cells)', ...
                nanmean(cumGrowths), nanstd(cumGrowths))];
            summary = [summary sprintf('</table>\n<br/>\n')];

            %mass
            summary = [summary sprintf('<table class="summary" cellpadding="0" cellspacing="0">\n')];
            summary = [summary sprintf('<tr><th>%s</th><td>%5.3f ± %5.3f</td></tr>\n', 'Mean Mass at t==0 (cells)', ...
                nanmean(initMass), nanstd(initMass))];
            summary = [summary sprintf('<tr><th>%s</th><td>%5.3f</td></tr>\n', 'Median Mass at t==0; (cells)', ...
                nanmedian(initMass))];
            summary = [summary sprintf('<tr><th>%s</th><td>%5.3f ± %5.3f</td></tr>\n', 'Mass Mass at t==τ (cells)', ...
                nanmean(finMass), nanstd(finMass))];
            summary = [summary sprintf('<tr><th>%s</th><td>%5.3f</td></tr>\n', 'Median at t==τ (cells)', ...
                nanmedian(finMass))];
            summary = [summary sprintf('</table>\n')];

            summary = [summary sprintf('</body>\n')];
            summary = [summary sprintf('</html>\n')];
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function Prod_cell = mult_vec(F_cell, G_cell)

chkarg(istypesizeof(F_cell, 'complexcell', [1 0]) || istypesizeof(F_cell, 'complexcell', [0 1]), ...
	'"F_cell" should be row or column cell array whose each element is array with complex elements.');
chkarg(istypesizeof(G_cell, 'complexcell', [1 0]) || istypesizeof(G_cell, 'complexcell', [0 1]), ...
	'"F_cell" should be row or column cell array whose each element is array with complex elements.', Axis.count);
chkarg(length(F_cell)==length(G_cell), '"F_cell" and "G_cell" should be the same length.');

n = length(F_cell);

Prod_cell = cell(size(F_cell));
for w = 1:n
	chkarg(all(size(F_cell{w})==size(G_cell{w})), '%d-component of "F_cell" and "G_cell" should have the same size.');
	Prod_cell{w} = F_cell{w} .* G_cell{w};
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function T = ctranspose(T,varargin)
% conjugate of a tensor
%
% Input
%  T - @tensor
%
% Output
%  T - @tensor
%

switch T.rank

  case 1

  case 2

    T.M = conj(permute(T.M,[2 1 3:ndims(T.M)]));

  case 3

  case 4

    % check for symmetry

    % convert to a matrix
    M = tensor42(T.M,T.doubleConvention);

    % invert the matrix
    M = conj(permute(M,[2 1 3:ndims(T.M)]));

    % make some corrections
    % TODO: why this is needed?
    if ~check_option(varargin,'skipCorrection')
      w = 1./(1+((1:6)>3));
      w = w.' * w;
      M = M .* w;
    end

    % convert back to a 4 rank tensor
    T.M = tensor24(M,T.doubleConvention);
end

author: mtex-toolbox

repo: mtex

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


text:
% evaluate simulated data - ICA and NMF
function separ(sep0, offset0, path_data, path_res, prename, niter, savethis, hint, sep_how)
% sep_how: string i - ica, n - nmf -> method for separation
if ~exist('hint', 'var')
    hint = 0;
end

if ~exist('sep_how', 'var')
    sep_how = 'in';
end

fprintf('Separation of components... \n')
for rr = 1: length(offset0)
    fprintf('\n%g: ',rr)
    for ll=1 : length(sep0)
        fprintf('.')
        p.namedir = [prename num2str(100*sep0(ll)) 'offset_' num2str(offset0(rr))];
        cd ([path_data p.namedir])
        for mm=1:niter

            namefile = [p.namedir '-iter_' num2str(mm)];
            load ([namefile '.mat'])
            %         ims(psf);
            %         SaveImageFULL('psf', 'pf');

            if sum(sep_how == 'i')>0 %ICA
                [icasig{mm}, A{mm}, W{mm}] = fastica (dveccr, 'numOfIC', 2, 'g', 'tanh');
                icapixICA{mm} = reshape(A{mm},32, 32, 2);
            end
            if sum(sep_how == 'n')>0 %NMF
                ncomp = 2; %number of components to be separated
                if hint
                    dvec_ind = double(squeeze(reshape(array2im(dpixc_ind), p.nx*p.ny, 1, 2))); % vectors of resized images
                    %                     [out, bg(mm), bg_im]=backgroundoffset(dpixc);
                    [out, bg(mm), bg_im]=backgroundoffset(dpixc, 'no', 5, 20, 8); %empirical values...
                    dvec_bg = bg(mm)*ones(1, p.nx*p.ny);
                    %                     dvec_bg = p.offset*ones(1, p.nx*p.ny); %changed for offset 10...

                    % % %                     blinkmatrand = rand(p.Nt, ncomp);
                    blinkmatrand = blinkmat';
                    winit = [blinkmatrand,ones(p.Nt,1)];             %random weights will be assigned to firts two and bg fixed

                    hinit = [dvec_ind; dvec_bg];       %original 'true' points + background
% % %                     hinit = [rand(ncomp, p.nx*p.ny); dvec_bg];
                    ncomp = ncomp+1; %background added
                    [w{mm},h{mm}]=nmf(double(dveccr'),ncomp+1,1,winit,hinit, [3], [1 2 3]);

                else
                    [w{mm},h{mm}]=nmf(double(dveccr'),ncomp,1);
                end
                icapixNMF{mm} = shiftdim(reshape(h{mm},ncomp,32,32), 1);
            end
            %             imstiled(icapixICA{mm});
            %             SaveImageFULL([p.namedir 'ICA_' num2str(mm)], 'p');
            %             imstiled(icapixNMF{mm});
            %             SaveImageFULL([p.namedir 'NMF_' num2str(mm)], 'p');
            close all

        end

        p.path_data = path_data;
        p.path = path_res;

        if savethis == 1
            fprintf('saving data \n');
            if ~(strcmp('p.path', 'p.path_data')) %not identical
                mkdir ([p.path p.namedir]);
                cd ([p.path p.namedir]);
            end
            save ([p.namedir '_separ'])
            writedata([],[],p,[p.namedir '_param'])
        end
    end
end

fprintf('\n')

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%% HEAT2D_POOL runs the HEAT2D example interactively with MATLABPOOL.
%

%
%  Read a black and white image which will be our "initial condition".
%
  x = imread ( 'surfsup.tif' );
%
%  Display the image.
%
  figure ( 1 )
  imshow ( x )
%
%  Get 4 "labs" to work on the job.
%
  matlabpool open local 4
  y = heat2d_fun ( x );
  matlabpool close

author: johannesgerer

repo: jburkardt-m


text:
function [mWalkout, fR95, fProbR95, PHI, R] = calc_Schusterwalk(mCatalog)
% function  [mWalkout, fR95, fProbR95, PHI, R] = calc_Schusterwalk(mCatalog)
% ------------------------------------------------------------------------
% Calculate random walkout (Schuster's method). See Rydelek & Sacks, Nature
% 1989. Shows phazor plot for one catalog.
%
% Incoming variables:
% mCatalog : EQ catalog
%
% Outgoing variables:
% mWalkout : Random walkout phases
% fR95     : Critical radius at 95% level
% fProbR95 : Probability to obtain vector of length >= R
% PHI      : Sum of phase angles
% R        : Possible sum of vector length
%
% J. Woessner, woessner@seismo.ifg.ethz.ch
% last update: 12.02.03

% Get hours and minutes, calculate minutes
vH = mCatalog(:,8);
vMin = mCatalog(:,9);
vTime = vH*60+vMin;

% Calculate degrees and radians
vThetaDegree = vTime*180/720;
vThetaRad = deg2rad(vThetaDegree);

mWalkout(1,1)=0;
mWalkout(1,2)=0;
for i=1:length(mCatalog(:,8))
    mWalkout(i+1,1)=mWalkout(i,1)+sin(vThetaRad(i));
    mWalkout(i+1,2)=mWalkout(i,2)+cos(vThetaRad(i));
end

A=sum(sin(vThetaRad));
B=sum(cos(vThetaRad));
R=sqrt(A^2+B^2);
PHI=atan(B/A);
if mWalkout(length(mCatalog(:,8)),2)<0
    PHI=PHI+pi;
end

N=length(mCatalog(:,8));
% Critical radius at 95% level
fR95=1.37*sqrt(N); %Nature 1989
% Probability to obtain vector of length >= R
fProbR95=exp(-R^2/N);

figure_w_normalized_uicontrolunits('tag','schuster')
% Plot Schuster walkout
plot(mWalkout(:,1),mWalkout(:,2));
hold on;
polar([0:360]*pi/180,ones(1,361)*fR95,'r');
[x, y] =pol2cart(PHI,R);
plot([0,x],[0 y],'g');
plot ([0 0],[-1/5*fR95,1/5*fR95],'k');
plot ([-1/5*fR95,1/5*fR95],[0 0],'k');
hold off
text (0,2/5*fR95,'0:00','HorizontalAlignment','center');
text (2/5*fR95,0,'6:00','HorizontalAlignment','left');
text (0,-2/5*fR95,'12:00','HorizontalAlignment','center');
text (-2/5*fR95,0,'18:00','HorizontalAlignment','right');
text (0,fR95 ,'95KI','VerticalAlignment','bottom','HorizontalAlignment','center');
axis equal;
%xlim([-15 15])
%ylim([-15 15])
text (-fR95,-fR95,num2str(length(mCatalog(:,8))));

figure_w_normalized_uicontrolunits('tag','schuster2')
% Plot hourly histogram
subplot(2,1,1);
vTime = (mCatalog(:,8)*60+mCatalog(:,9))/60; % Calculate decimal hour
histogram(vTime,0:1:24)
% Plot FMD
subplot(2,1,2);
[mFMDC, mFMD] = calc_FMD(mCatalog);
plot(mFMDC(1,:),mFMDC(2,:),'ks',mFMD(1,:),mFMD(2,:),'k^');

author: CelsoReyes


text:
function [root] = bsdsRoot()
% function [root] = bsdsRoot()
%
% Return the root directory for the BSDS data.
%
% David R. Martin <dmartin@eecs.berkeley.edu>
% March 2003

root = '/home/cs/dmartin/public_html/segbench/BSDS300';

author: Cloud-CV

repo: object-proposals


text:
% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
function s = char(robot)
%
% S = R.char() is a string representation of the robot parameters.

	s = '';
	for j=1:length(robot)
		r = robot(j);

		% informational line
		if r.mdh
			convention = 'modDH';
		else
			convention = 'stdDH';
		end
		s = sprintf('%s (%d axis, %s, %s)', r.name, r.n, config(r), convention);

		% comment and other info
		line = '';
		if ~isempty(r.manufacturer)
			line = strcat(line, sprintf(' %s;', r.manufacturer));
		end
		if ~isempty(r.comment)
			line = strcat(line, sprintf(' %s;', r.comment));
		end
		s = strvcat(s, line);

		% link parameters
		s = strvcat(s, '+---+-----------+-----------+-----------+-----------+');
		s = strvcat(s, '| j |     theta |         d |         a |     alpha |');
		s = strvcat(s, '+---+-----------+-----------+-----------+-----------+');
		s = strvcat(s, char(r.links, true));
		s = strvcat(s, '+---+-----------+-----------+-----------+-----------+');

		% gravity, base, tool

		s_grav = horzcat(strvcat('grav = ', ' ', ' '), num2str(r.gravity));
		s_grav = strvcat(s_grav, " ");
		s_base = horzcat(strvcat("  base = ",' ',' ', ' '),num2str(r.base));

		s_tool = horzcat(strvcat('   tool =  ',' ',' ', ' '), num2str(r.tool));

		line = horzcat(s_grav, s_base, s_tool);

		s = strvcat(s, ' ', line);
		if j ~= length(robot)
			s = strvcat(s, ' ');
		end
	end
endfunction

function v = config(r)
            v = '';
            for i=1:r.n
                v(i) = r.links(i).RP;
            end
        end

author: petercorke

repo: robotics-toolbox-matlab


text:
function r = subsref(this, selector)
%SUBSREF   Access parameters of HS optical flow algorithm
%   R = SUBSREF(THIS, SEL) extracts and returns parameter of a HS
%   optical algorithm THIS using given selector SEL.
%   Available selectors:
%     - ".optical_flow":   Flow base object
%     - ".pyramid_levels": Levels of the image pyramid
%     - ".warping_mode":   'forward' for forward warping and 'backward'
%                          for backward warping
%     - ".solver":         Selector for the linear equation solver
%                          ('backslash' for the builtin solver and 'sor'
%                          for a MEX-based SOR solver)
%     - ".max_iters":      Maximum number of iterations at each pyramid
%                          level
%     - ".limit_update":   Flag whether to restrict the incremental flow
%                          at each level to [-1, 1]
%
%   This is a member function of the class 'hs_optical_flow'.
%
%   Author: Deqing Sun, Department of Computer Science, Brown University
%   Contact: dqsun@cs.brown.edu
%   $Date: 2007-11-30 $
%   $Revision: $
%
% Copyright 2007-2010, Brown University, Providence, RI. USA
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

  % Allow full access to all fields of the structure
  r = builtin('subsref', this, selector);

author: kristinbranson

repo: JAABA


text:
function DEMr = resample(DEM,target,method,swapzone,varargin)

%RESAMPLE change spatial resolution of a GRIDobj
%
% Syntax
%
%     DEMr = resample(DEM,cellsize)
%     DEMr = resample(DEM,GRID)
%     DEMr = resample(...,method)
%     DEMr = resample(DEM,GRID,method,swapzone)
%
% Description
%
%     resample changes the cellsize of a grid. The function uses the MATLAB
%     function imtransform. If an instance of GRIDobj is supplied as
%     second argument, resample interpolates values in DEM to match the
%     spatial reference of GRID.
%
% Input arguments
%
%     DEM       grid object (GRIDobj)
%     cellsize  cellsize of resampled grid
%     GRID      other grid object
%     method    'bicubic', 'bilinear', or 'nearest'
%     swapzone  true or false. If true and if DEM and GRID have different
%               projected coordinate systems, the function will attempt to
%               reproject and resample the DEM in one step. Note that this
%               requires the mapping toolbox. In case, the DEM is in a
%               geographic coordinate system, please use the function
%               reproject2utm(DEM,GRID).
%
% Output arguments
%
%     DEMr    grid object (GRIDobj)
%
% Example
%
%     DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
%     DEMr = resample(DEM,100);
%     imagesc(DEMr)
%
%
% See also: griddedInterpolant, imtransform
%
% Author:  Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 8. August, 2015

% check input arguments
% narginchk(2,6)
validateattributes(target,{'double' 'GRIDobj'},{'scalar'})
if nargin == 2;
    method = 'bilinear';
    swapzone = false;
elseif nargin == 3;
    method = validatestring(method,{'bicubic', 'bilinear', 'nearest' });
    swapzone = false;
else
    method = validatestring(method,{'bicubic', 'bilinear', 'nearest' });
end

% check underlying class
if islogical(DEM.Z)
    method = 'nearest';
end

if swapzone && isa(target,'GRIDobj');
    if ~isequal(DEM.georef.GeoKeyDirectoryTag.ProjectedCSTypeGeoKey,...
            target.georef.GeoKeyDirectoryTag.ProjectedCSTypeGeoKey)
        warning('UTM zones differ. Will attempt to match.');
        swapzone = true;
    else
        swapzone = false;
    end
end

% get coordinate vectors
[u,v] = getcoordinates(DEM);

% tform
T = maketform('affine',[1 0 0; 0 1 0; 0 0 1]);

% Fillvalues

if isinteger(DEM.Z)
    fillval = 0;
elseif islogical(DEM.Z)
    fillval = 0;
else
    fillval = nan;
end

p = inputParser;
p.FunctionName = 'Paras';
% optional
% default values.
addParameter(p,'fillval',fillval);
parse(p,varargin{:});
fillval = p.Results.fillval;

if isa(target,'GRIDobj')
    % the target is another GRIDobj

    % check whether both grids have the same projection
    if swapzone
        mstructsource = DEM.georef.mstruct;
        mstructtarget = target.georef.mstruct;
        T = maketform('custom', 2, 2, ...
                @FWDTRANS, ...
                @INVTRANS, ...
                []);
    end

    DEMr    = target;
    [xn,yn] = getcoordinates(DEMr);

    DEMr.Z = imtransform(DEM.Z,T,method,...
        'Udata',[u(1) u(end)],'Vdata',[v(1) v(end)],...
        'Xdata',[xn(1) xn(end)],'Ydata',[yn(1) yn(end)],...
        'Size',DEMr.size,...
        'FillValues',fillval);
    DEMr.name = [DEM.name ' (resampled)'];

else
    csnew   = target;
    DEMr    = GRIDobj([]);
    [DEMr.Z,xn,yn] = imtransform(DEM.Z,T,method,...
        'Udata',[u(1) u(end)],'Vdata',[v(1) v(end)],...
        'Xdata',[u(1) u(end)],'Ydata',[v(1) v(end)],...
        'XYscale',[csnew csnew],...
        'FillValues',fillval);

    % new referencing matrix
    DEMr.refmat = [0 -csnew;...
                   csnew 0; ...
                   xn(1)-csnew ...
                   yn(1)+csnew];
    % size of the resampled grid
    DEMr.size    = size(DEMr.Z);
    DEMr.cellsize = csnew;
    DEMr.georef = DEM.georef;
end

DEMr.name    = [DEM.name ' (resampled)'];

if ~isempty(DEMr.georef) && ~isa(target,'GRIDobj');
    DEMr.georef.RefMatrix = DEMr.refmat;
    DEMr.georef.Height = DEMr.size(1);
    DEMr.georef.Width  = DEMr.size(2);

    DEMr.georef.SpatialRef = refmatToMapRasterReference(DEMr.refmat, DEMr.size);

end

%%

    function x = FWDTRANS(u,~)
        % invtrans first
        [lati,long] = minvtran(mstructsource,u(:,1),u(:,2));
        [x,y] = mfwdtran(mstructtarget,lati,long);
        x = [x y];
    end

    function u = INVTRANS(x,~)
        [lati,long] = minvtran(mstructtarget,x(:,1),x(:,2));
        [x,y] = mfwdtran(mstructsource,lati,long);
        u = [x y];

    end
end

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
% Demo for aggregate channel features object detector on KAIST dataset.
%
% See also acfReadme.m
%
% Piotr's Computer Vision Matlab Toolbox      Version 3.40
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]
%
% 2015.06.02. Modified by Soonmin Hwang [smhwang-at-rcv.kaist.ac.kr]
% 2017.10.11. Update

%% extract training and testing images and ground truth
dataDir='data/kaist-rgbt/';
addpath( genpath( 'pdollar-toolbox' ) );
addpath( genpath( 'libs' ) );

%% set up opts for training detector (see acfTrain)
opts=acfTrain2(); opts.modelDs=[50 20.5]; opts.modelDsPad=[64 32];
opts.pPyramid.smooth=.5;
opts.pPyramid.pChns.pColor.smooth=0; opts.nWeak=[32 128 256 512 1024 2048];
opts.pBoost.pTree.maxDepth=4;
opts.pBoost.pTree.fracFtrs=1/16;
opts.nNeg=50000;
opts.nAccNeg = 100000;
opts.nPerNeg = 20;
opts.pPyramid.pChns.pGradHist.softBin=1; opts.pJitter=struct('flip',1);

% opts.trainSet = [dataDir 'imageSets/train20.txt'];
opts.trainSet = [dataDir 'imageSets/train04.txt'];
opts.posGtDir=[dataDir 'annotations'];
opts.posImgDir=[dataDir 'images'];

pLoad={'lbls',{'person'},'ilbls',{'people','person?','cyclist'},'squarify',{3,.41}};
opts.pLoad = [pLoad 'hRng',[45 inf], 'vType', {'none'} ];

%% Set thermal features
opts = acfSetTchFeatures(opts, [1 0 1 0]);  % T, TM, THOG, TM+TO

%% train detector (see acfTrain)
detector = acfTrain2( opts );

%% modify detector (see acfModify)
pModify=struct('cascThr',0,'cascCal',.025);
detector=acfModify(detector,pModify);

%% run detector on a sample image (see acfDetect)
imgDir = [dataDir 'images']; gtDir = [dataDir 'annotations'];
testSet = [dataDir 'imageSets/test20.txt'];

imgNms =bbGt2('getSubsetFiles',{imgDir, imgDir, gtDir}, testSet);

idx = 102;
I = imreadHistEq(imgNms{1,idx}, imgNms{2,idx});
tic, bbs=acfDetect2(I,detector); toc
figure(1); imshow(I(:,:,1:3)); bbApply('draw',bbs); pause(.1);

%% test detector and plot roc (see acfTest)
% Reasonable
[~,~,gt,dt]=acfTest2('name',opts.name,'imgDir',imgDir,...
  'gtDir',gtDir,'pLoad',[pLoad, 'hRng',[45 inf],...
  'vType',{{'none','partial'}},'xRng',[5 635],'yRng',[5 475]],...
  'pModify',pModify,'reapply',0,'show',2,...
  'lims', [3.1e-3 3.1e1 .2 1], 'clr', 'g', 'subset', testSet, ...
  'figName', 'Reasonable-all');

author: SoonminHwang

repo: rgbt-ped-detection


text:
function f_x = ParFor14(in1)
%PARFOR14
%    F_X = PARFOR14(IN1)

%    This function was generated by the Symbolic Math Toolbox version 8.2.
%    20-Sep-2019 09:35:38

u = in1(:,1);
ux = in1(:,4);
uxxx = in1(:,6);
f_x = ((ux.*-2.895410174450939e19+uxxx.*2.0816405633874e18+u.*ux.*2.276963278976451e19).*-1.402326054557642e-16)./u;

author: dynamicslab

repo: SINDy-PI


text:
% This file is part of TREEQSM.
%
% TREEQSM 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.
%
% TREEQSM 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 TREEQSM.  If not, see <http://www.gnu.org/licenses/>.

function [R,dR1,dR2] = form_rotation_matrices(theta)

% --------------------------------------------------------------------------
% FORM_ROTATION_MATRICES.M      Forms rotation matrices R = R2*R1 and its
%                                   derivatives
%
% Input
% theta    Plane rotation angles (t1, t2)
%
% Output
% R        Rotation matrix
% R1       Plane rotation [1 0 0; 0 c1 -s1; 0 s1 c1]
% R2       Plane rotation [c2 0 s2; 0 1 0; -s2 0 c2]

c = cos(theta);
s = sin(theta);

R1 = [1 0 0; 0 c(1) -s(1); 0 s(1) c(1)];
R = R1;

R2 = [c(2) 0 s(2); 0 1 0; -s(2) 0 c(2)];
R = R2*R;

if nargout > 1
    dR1 = [0 0 0; 0 -R1(3,2) -R1(2,2); 0 R1(2,2) -R1(3,2)];
end

if nargout > 2
    dR2 = [-R2(1,3) 0 R2(1,1); 0 0 0; -R2(1,1) 0 -R2(1,3)];
end

author: InverseTampere

repo: TreeQSM


text:
classdef FullInnerMeshCreator < handle

    properties (Access = protected) % Inputs
        type
        unfittedMesh
    end

    methods (Static, Access = public)

        function ime = create(cParams)
            switch cParams.type
                case 'Matlab'
                    ime = FullInnerMeshCreator_Matlab(cParams);
                case 'GiD'
                    ime = FullInnerMeshCreator_GiD(cParams);
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [calscores,w0] = parallel_cal(w,scores,wfuse)
%

if nargin==0
    test_this();
    return;
end

if ~exist('scores','var') || isempty(scores)
    calscores = sprintf(['parallel calibration:',repmat(' %g',1,length(w))],w);
    return;
end

[m,n] = size(scores);

if nargout>1, w0 = init_w0(wfuse); end

calscores = linTrans(w,@(w)map_this(w),@(w)transmap_this(w));

    function w0 = init_w0(wfuse)
        assert(length(wfuse)-1==m);
        scal = wfuse(1:end-1);
        offs = wfuse(end);
        W = [scal*(m+1);((m+1)/m)*offs*ones(m,1)];
        w0 = W(:);
    end

    function y = map_this(w)
        w = reshape(w,m,2);
        y = bsxfun(@times,scores,w(:,1));
        y = bsxfun(@plus,y,w(:,2));
    end

    function w = transmap_this(y)
        y = reshape(y,m,n);
        w = [sum(y.*scores,2),sum(y,2)];
    end

end

function test_this()

    scores = randn(4,10);
    [sys,w0] = parallel_cal([],scores,(1:5)');
    test_MV2DF(sys,w0);

end

author: nesl


text:
classdef substrTest < TestCase
    methods
        function this = substrTest(methodName)
            this = this@TestCase(methodName);
        end

        function testEmptyString(~)
            s = '';
            assertEqual(s(1:0), substr('', 1, 0));
            assertExceptionThrown(@() substr('', 1, 1), 'MATLAB:badsubscript');
            assertExceptionThrown(@() substr('', 0, 0), 'MATLAB:badsubscript');
        end

        function testWholeString(~)
            assertEqual('abc', substr('abc', 1, 0));
        end

        function testPartOfString(~)
            assertEqual('bc', substr('abcd', 2, 3));
        end

        function testPartOfStringFromEnd(~)
            assertEqual('bc', substr('abcd', 2, -1));
        end
    end
end

function result = substr(s, i, j)
    result = edu.stanford.covert.util.substr(s, i, j);
end


text:
classdef GeographicElementSet < AbstractElementSet
    %GeographicElementSet Summary of this class goes here
    %   Detailed explanation goes here

    properties
        lat%(1,1) double %rad
        long%(1,1) double %rad
        alt%(1,1) double %km
        velAz%(1,1) double %rad NEZ frame
        velEl%(1,1) double %rad NEZ frame
        velMag%(1,1) double %km/s NEZ frame

        optVar GeographicElementSetVariable
    end

    properties(Constant)
        typeEnum = ElementSetEnum.GeographicElements;
    end

    methods
        function obj = GeographicElementSet(time, lat, long, alt, velAz, velEl, velMag, frame)
            if(nargin > 0)
                obj.time = time;
                obj.lat = lat;
                obj.long = long;
                obj.alt = alt;
                obj.velAz = velAz;
                obj.velEl = velEl;
                obj.velMag = velMag;
                obj.frame = frame;
            end
        end

        %vectorized
        function cartElemSet = convertToCartesianElementSet(obj)
            radius = NaN(1, length(obj));
            for(i=1:length(obj))
                radius(i) = obj(i).frame.getOriginBody().radius;
            end
            r = radius + [obj.alt];

            x = r.*cos([obj.lat]).*cos([obj.long]);
            y = r.*cos([obj.lat]).*sin([obj.long]);
            z = r.*sin([obj.lat]);

            rVect = [x;y;z];

            sezVVectAz = pi - [obj.velAz];
            sezVVectEl = [obj.velEl];
            sezVVectMag = [obj.velMag];

            [x,y,z] = sph2cart(sezVVectAz, sezVVectEl, sezVVectMag);
            vectorSez = [x;y;z];
            vVect = rotSEZVectToECEFCoords(rVect, vectorSez); %not necessarily ECEF coords here but the transform holds

%             cartElemSet = CartesianElementSet(obj.time, rVect, vVect, obj.frame);
%             cartElemSet = repmat(CartesianElementSet.getDefaultElements(), size(obj));
%             for(i=1:length(obj))
%                 cartElemSet(i) = CartesianElementSet(obj(i).time, rVect(:,i), vVect(:,i), obj(i).frame);
%             end

            obj = obj(:)';
            cartElemSet = CartesianElementSet([obj.time], rVect, vVect, [obj.frame]);
        end

        %vectorized
        function kepElemSet = convertToKeplerianElementSet(obj)
            kepElemSet = convertToKeplerianElementSet(convertToCartesianElementSet(obj));
        end

        %vectorized
        function geoElemSet = convertToGeographicElementSet(obj)
            geoElemSet = obj;
        end

        %vectorized
        function univElemSet = convertToUniversalElementSet(obj)
            univElemSet = convertToUniversalElementSet(convertToKeplerianElementSet(obj));
        end

        function elemVect = getElementVector(obj)
            elemVect = [rad2deg(obj.lat),rad2deg(obj.long),obj.alt,rad2deg(obj.velAz),rad2deg(obj.velEl),obj.velMag];
        end

        function newElemSet = copyWithoutOptVar(obj)
            newElemSet = GeographicElementSet(obj.time, obj.lat, obj.long, obj.alt, obj.velAz, obj.velEl, obj.velMag, obj.frame);
        end
    end

    methods(Access=protected)
        function displayScalarObject(obj)
            fprintf('Geographic State \n\tTime: %0.3f sec UT \n\tLatitude: %0.3f deg \n\tLongitude: %0.3f deg \n\tAltitude: %0.3f km \n\tVel Az (NEZ): %0.3f deg \n\tVel El (NEZ): %0.3f deg \n\tVel Mag (NEZ): %0.3f km/s \n\tFrame: %s\n', ...
                    obj.time, ...
                    rad2deg(obj.lat), ...
                    rad2deg(obj.long), ...
                    obj.alt, ...
                    rad2deg(obj.velAz), ...
                    rad2deg(obj.velEl), ...
                    obj.velMag, ...
                    obj.frame.getNameStr());
        end
    end

    methods(Static)
        function elemSet = getDefaultElements()
            elemSet = GeographicElementSet();
        end

        function errMsg = validateInputOrbit(errMsg, hLat, hLong, hAlt, hVelAz, hVelEl, hVelMag, bodyInfo, bndStr, checkElement)
            if(isempty(bndStr))
                bndStr = '';
            else
                bndStr = sprintf(' (%s Bound)', bndStr);
            end

            if(checkElement(1))
                lat = str2double(get(hLat,'String'));
                enteredStr = get(hLat,'String');
                numberName = ['Latitude', bndStr];
                lb = -90;
                ub = 90;
                isInt = false;
                errMsg = validateNumber(lat, numberName, lb, ub, isInt, errMsg, enteredStr);
            end

            if(checkElement(2))
                long = str2double(get(hLong,'String'));
                enteredStr = get(hLong,'String');
                numberName = ['Longitude', bndStr];
                lb = -360;
                ub = 360;
                isInt = false;
                errMsg = validateNumber(long, numberName, lb, ub, isInt, errMsg, enteredStr);
            end

            if(checkElement(3))
                alt = str2double(get(hAlt,'String'));
                enteredStr = get(hAlt,'String');
                numberName = ['Altitude', bndStr];
                lb = -bodyInfo.radius;
                ub = Inf;
                isInt = false;
                errMsg = validateNumber(alt, numberName, lb, ub, isInt, errMsg, enteredStr);
            end

            if(checkElement(4))
                bfVx = str2double(get(hVelAz,'String'));
                enteredStr = get(hVelAz,'String');
                numberName = ['Body-Fixed Velocity Azimuth', bndStr];
                lb = -360;
                ub = 360;
                isInt = false;
                errMsg = validateNumber(bfVx, numberName, lb, ub, isInt, errMsg, enteredStr);
            end

            if(checkElement(5))
                bfVy = str2double(get(hVelEl,'String'));
                enteredStr = get(hVelEl,'String');
                numberName = ['Body-Fixed Velocity Elevation', bndStr];
                lb = -90;
                ub = 90;
                isInt = false;
                errMsg = validateNumber(bfVy, numberName, lb, ub, isInt, errMsg, enteredStr);
            end

            if(checkElement(6))
                bfVz = str2double(get(hVelMag,'String'));
                enteredStr = get(hVelMag,'String');
                numberName = ['Body-Fixed Velocity Magnitude', bndStr];
                lb = 0;
                ub = Inf;
                isInt = false;
                errMsg = validateNumber(bfVz, numberName, lb, ub, isInt, errMsg, enteredStr);
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function u = u_function ( time )

%*****************************************************************************80
%
%% U_FUNCTION evaluates the time-dependent boundary values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 December 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real TIME, the current time.
%
%    Output, real U(2), the boundary values at the left and right endpoints.
%
  ua =   sin ( pi * time );
  ub = - sin ( pi * time );

  u = [ ua; ub ];

  return
end


text:
function normqqplot(y,class)

%NORMQQPLOT produces a Quantile-Quantile plot in which the vector y is plotted against
% the quantiles of a standard normal distribution.
%
% Required input arguments:
%       y : row or column vector
%
% Optional input arguments:
%    class : a string used for the y-label and the title(default: ' ')
%
% This function is part of LIBRA: the Matlab library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html
%
%Written by Nele Smets
% Last update: 12/03/2004

set(gcf,'Name', 'Normal QQ-plot', 'NumberTitle', 'off');
if nargin==1
    class=' ';
end
n=length(y);
for i=1:n
    normalquantile(i)=norminv((i-1/3)/(n+1/3),0,1);
end
y=sort(y);
plot(normalquantile,y,'bo')
xlabel('Quantiles of the standard normal distribution');
ylabel(['Standardized ',class,' residual']);
title(class);
box on

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
%% h = HEXSCATTER( x, y, ... )
% Gordon Bean, February 2014
% A scatter-plot substitute - generate a density plot using hexagonal
% patches.
%
% Syntax
% hexscatter(xdata, ydata)
% hexscatter(xdata, ydata, 'Name', Value, ...)
% h = hexscatter(...)
%
% Description
% hexscatter(xdata, ydata) creates a density plot of the ydata versus the
% xdata using hexagonal tiles. xdata and ydata should be vectors. NaN
% values (and their corresponding values in the other vector) are ignored.
%
% hexscatter(xdata, ydata, 'Name', Value, ...) accepts name-value pairs of
% arguments from the following list (defaults in {}):
%  'xlim' { [min(xdadta(:) max(xdata(:))] } - a 2-element vector containing
%  the lower and upper bounds of the 2nd dimension of the grid.
%  'ylim' { [min(ydadta(:) max(ydata(:))] } - a 2-element vector containing
%  the lower and upper bounds of the 1st dimension of the grid.
%  'res' { 50 } - the resolution, or number of bins in each dimension. The
%  total number of bins will be the resolution squared.
%  'drawEdges' { false } - if true, edges are drawn around each hexagonal
%  patch.
%  'showZeros' { false } - if true, bins with 0 counts are shaded; if
%  false, only bins with non-zero counts are colored.
%
% h = hexscatter( ... ) returns the object handle to the patch object
% created.
%
% Examples
% hexscatter(rand(2000,1), rand(2000,1))
%
% hexscatter(rand(2000,1), rand(2000,1), 'res', 90)
%
% Also available in the Bean Matlab Toolkit:
% https://github.com/brazilbean/bean-matlab-toolkit

function h = hexscatter( xdata, ydata, varargin )
    params = default_param( varargin, ...
        'xlim', [min(xdata(:)) max(xdata(:))], ...
        'ylim', [min(ydata(:)) max(ydata(:))], ...
        'res', 50, ...
        'drawEdges', false, ...
        'showZeros', false);

    if params.drawedges
        ec = 'flat';
    else
        ec = 'none';
    end

    %% Determine grid
    xl = params.xlim;
    yl = params.ylim;

    xbins = linspace(xl(1), xl(2), params.res);
    ybins = linspace(yl(1), yl(2), params.res);
    dy = diff(ybins([1 2]))*0.5;

    [X, Y] = meshgrid(xbins, ybins);
    n = size(X,1);
    Y(:,1:fix(end/2)*2) = ...
        Y(:,1:fix(end/2)*2) + repmat([0 dy],[n,fix(n/2)]);

    %% Map points to boxes
    nix = isnan(xdata) | isnan(ydata);
    xdata = xdata(~nix);
    ydata = ydata(~nix);

    % Which pair of columns?
    dx = diff(xbins([1 2]));
    foox = floor((xdata - xbins(1)) ./ dx)+1;
    foox(foox > length(xbins)) = length(xbins);

    % Which pair of rows?
    % Use the first row, which starts without an offset, as the standard
    fooy = floor((ydata - ybins(1)) ./ diff(ybins([1 2])))+1;
    fooy(fooy > length(ybins)) = length(ybins);

    % Which orientation
    orientation = mod(foox,2) == 1;

    % Map points to boxes
    foo = [xdata - xbins(foox)', ydata - ybins(fooy)'];

    % Which layer
    layer = foo(:,2) > dy;

    % Convert to block B format
    toflip = layer == orientation;
    foo(toflip,1) = dx - foo(toflip,1);

    foo(layer==1,2) = foo(layer==1,2) - dy;

    % Find closest corner
    dist = sqrt(sum(foo.^2,2));
    dist2 = sqrt(sum(bsxfun(@minus, [dx dy], foo).^2, 2));

    topright = dist > dist2;

    %% Map corners back to bins
    % Which x bin?
    x = foox + ~(orientation == (layer == topright));
    x(x > length(xbins)) = length(xbins);

    % Which y bin?
    y = fooy + (layer & topright);
    y(y > length(ybins)) = length(ybins);

    ii = sub2ind(size(X), y, x);

    %% Determine counts
    counts = sum(bsxfun(@eq, ii, 1:numel(X)),1);

    newplot;
    xscale = diff(xbins([1 2]))*2/3;
    yscale = diff(ybins([1 2]))*2/3;
    theta = 0 : 60 : 360;
    x = bsxfun(@plus, X(:), cosd(theta)*xscale)';
    y = bsxfun(@plus, Y(:), sind(theta)*yscale)';

    if params.showzeros
        h = patch(x, y, counts, 'edgeColor', ec);
    else
        jj = counts > 0;
        h = patch(x(:,jj), y(:,jj), counts(jj), 'edgeColor', ec);
    end

    if nargout == 0
        clear h;
    end

    %% Function: default_param
    % Gordon Bean, March 2012
    % Copied from https://github.com/brazilbean/bean-matlab-toolkit
    function params = default_param( params, varargin )
        if (iscell(params))
            params = get_params(params{:});
        end
        defaults = get_params(varargin{:});

        for f = fieldnames(defaults)'
            field = f{:};
            if (~isfield( params, lower(field) ))
                params.(lower(field)) = defaults.(field);
            end
        end
    end

    %% Function: get_params - return a struct of key-value pairs
    % Gordon Bean, January 2012
    %
    % Usage
    % params = get_params( ... )
    %
    % Used to parse key-value pairs in varargin - returns a struct.
    % Converts all keys to lower case.
    %
    % Copied from https://github.com/brazilbean/bean-matlab-toolkit
    function params = get_params( varargin )
        params = struct;

        nn = length(varargin);
        if (mod(nn,2) ~= 0)
            error('Uneven number of parameters and values in list.');
        end

        tmp = reshape(varargin, [2 nn/2]);
        for kk = 1 : size(tmp,2)
            params.(lower(tmp{1,kk})) = tmp{2,kk};
        end
    end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function meow(n, paws)
%
% This Matlab function plays n successive cat meows, where each of the n
% meows is randomly chosen from one of two available sounds.  The function
% pauses for paws seconds after each meow.  The default values of n and
% paws are 5 and 3, respectively.
%
% Dr. Phillip M. Feldman   18 June, 2009

if nargin<1, n= 5; end
if nargin<2, paws= 3; end

for i= 1 : n
   if (rand(1) <= 0.5)
      infile= 'meow1.wav';
   else
      infile= 'meow2.wav';
   end

   [Y,Fs,Nbits]= wavread(infile);
   sound(Y, Fs);
   pause(paws);
end

author: Sable

repo: mcbench-benchmarks


text:
report_this_filefun(mfilename('fullpath'));

l = isnan(tmap);
tmap(l) = 1;

[lat,lon] = meshgrat(tmap,tmapleg);
%[lat,lon] = meshgrat(vlat,vlon);

%[smap,smapleg] = country2mtx('switzerland',100);
%[lat0, lon0] = meshgrat(smap,smapleg);

% tmap = km2deg(tmap/1);
[X , Y]  = meshgrid(gx,gy);

%sw = interp2(lon0,lat0,smap,lon,lat);

ren = interp2(X,Y,re4,lon,lat);

mi = min(min(ren));
l =  isnan(ren);
ren(l) = mi-0.05;

figure_w_normalized_uicontrolunits('pos',[150 500 1000 700])

hold on; axis off
axesm('MapProjection','eqaconic','MapParallels',[],...
    'MapLatLimit',[s4_south s3_north],'MapLonLimit',[s2_west s1_east])

meshm(ren,tmapleg,size(tmap),tmap);

daspectm('m',50);
tightmap
view([0 90])
camlight; lighting phong
set(gca,'projection','perspective');

% load usalo
% h = displaym(usalo('state'));

%set(h(1),'color',[0.9 0.9 0.9],'Linewidth',2)
% h = displaym(gtlakevec); set(h(1),'color',[0.9 0.9 0.9],'Linewidth',2)

% h2 = displaym(PPpoint);
%h = displaym(PPtext); trimcart(h);

%pl = plotm(lima(:,2),lima(:,1),'hw');
%set(pl,'LineWidth',1.5,'MarkerSize',12,...
% 'MarkerFaceColor','y','MarkerEdgeColor','k')

pl = plotm(dam(:,2),dam(:,1),'^w');
set(pl,'LineWidth',1.5,'MarkerSize',12,...
    'MarkerFaceColor','w','MarkerEdgeColor','k')
pl = plotm(coastline(:,2), coastline(:,1),'w','Linewidth',2);
set(pl,'LineWidth',2);
zdatam(handlem('allline'),10000) % keep line on surface
%zdatam(handlem('alltext'),10000) % keep line on surface

j = jet;
%j = j(64:-1:1,:);
j = [ [ 0.85 0.9 0.9] ; j];

colormap(j); brighten(0.1);

axis off; set(gcf,'color','k')

setm(gca,'ffacecolor','k')
setm(gca,'fedgecolor','y','flinewidth',3);

setm(gca,'mlabellocation',2)
setm(gca,'meridianlabel','on')
setm(gca,'plabellocation',2)
setm(gca,'parallellabel','on')
setm(gca,'Fontcolor','w','Fontweight','bold','FontSize',12,'Labelunits','dm')

h5 = colorbar;
set(h5,'position',[0.82 0.35 0.01 0.3],'TickDir','out','Ycolor','w','Xcolor','w',...
    'Fontweight','FontSize',12);
set(gcf,'Inverthardcopy','off');