text:
% Analyzes families of whole-cell simulations
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Created: 11/1/2012
classdef MultiGenerations
    properties (Constant = true)
        ANCESTRY_COLIDX_CELL = 1
        ANCESTRY_COLIDX_GEN = 2
        ANCESTRY_COLIDX_FAMILY = 3
        ANCESTRY_COLIDX_PARENT = 4
        ANCESTRY_COLIDX_CHILD1 = 5
        ANCESTRY_COLIDX_CHILD2 = 6
    end

    methods (Static = true)
        function run(simBatchDir, nGen, nCellFirstGen)
            %% import
            import edu.stanford.covert.cell.sim.analysis.MultiGenerations;
            import edu.stanford.covert.cell.sim.util.CachedSimulationObjectUtil;
            import edu.stanford.covert.cell.sim.util.PlotUtil;
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;
            import edu.stanford.covert.cell.sim.util.SimulationEnsemble;

            %%
            if nargin < 3
                simBatchDir = '2012_11_15_18_48_23';
                nGen = 3;
                nCellFirstGen = 8;
            end

            outDir = [SimulationDiskUtil.getBaseDir() filesep simBatchDir];

            if ~exist(outDir, 'dir')
                mkdir(outDir)
            end

            %% load constants
            sim = CachedSimulationObjectUtil.load();
            massState = sim.state('Mass');
            pm = sim.state('ProteinMonomer');
            pc = sim.state('ProteinComplex');

            %% load data
            ancestry = MultiGenerations.calcAncestry(nGen, nCellFirstGen);
            relations = MultiGenerations.calcAncestryRelations(nGen, nCellFirstGen);

            [~, simData, simStartTimes, divSimTfs, finSimTfs, finSimIdxs, propNames, ...
                dnaBndMons, dnaBndCpxs, otherBndCpxs, dnaBndMonTfs, dnaBndCpxTfs, otherBndCpxTfs] = ...
                MultiGenerations.cacheSimData(simBatchDir);
            simData = simData{1};
            simStartTimes = simStartTimes{1};
            divSimTfs = divSimTfs{1};
            finSimTfs = finSimTfs{1};
            finSimIdxs = finSimIdxs{1};
            dnaBndMons = dnaBndMons{1};
            dnaBndCpxs = dnaBndCpxs{1};
            otherBndCpxs = otherBndCpxs{1};
            dnaBndMonTfs = dnaBndMonTfs{1};
            dnaBndCpxTfs = dnaBndCpxTfs{1};
            otherBndCpxTfs = otherBndCpxTfs{1};

            %% population growth
            %growth traces
            for iFamily = 1:nCellFirstGen
                famIdxs = find(finSimTfs & ancestry(:, MultiGenerations.ANCESTRY_COLIDX_FAMILY) == iFamily);
                colors = MultiGenerations.calcRedGreenColors(numel(famIdxs));
                [~, order] = sort(simStartTimes(famIdxs));
                colors(order, :) = colors;
                tStep = 10;

                [axesHandle, figHandle] = PlotUtil.newAxesHandle();
                hold(axesHandle, 'on');
                for i = 1:numel(famIdxs)
                    tmp = SimulationEnsemble.load(simBatchDir, {'Time' 'values'; 'MetabolicReaction', 'growth'}, [], [], tStep, 'extract', famIdxs(i));
                    tmpTime = permute(tmp.Time.values, [1 3 2]) / 3600;
                    tmpGrowth = permute(tmp.MetabolicReaction.growth, [1 3 2]) * ...
                        massState.cellInitialDryWeight / (1 - massState.fractionWetWeight) * 3600 * 1e15;
                    plot(axesHandle, simStartTimes(famIdxs(i)) + tmpTime, tmpGrowth, ...
                        'Color', colors(i, :));
                end
                xlabel(axesHandle, 'Time (h)');
                ylabel(axesHandle, 'Growth (fg h^{-1})');
                saveas(figHandle, [outDir filesep 'Growth-Family-' num2str(iFamily) '.pdf']);
                close(figHandle);
            end

            %number cells
            nCells = NaN(nGen, 1);
            nSurvive = NaN(nGen, 1);
            for iGen = 0:nGen-1
                nCells(iGen + 1) = sum(ancestry(:, MultiGenerations.ANCESTRY_COLIDX_GEN) == iGen & finSimTfs);
                nSurvive(iGen + 1) = sum(ancestry(:, MultiGenerations.ANCESTRY_COLIDX_GEN) == iGen & divSimTfs);
            end
            fracSurvive = nSurvive ./ nCells;

            [axesHandle, figHandle] = PlotUtil.newAxesHandle();
            hold(axesHandle, 'on');

            plot(axesHandle, 0:nGen-1, nCells, 'b.')
            xlabel(axesHandle, 'Generation');
            ylabel(axesHandle, 'Population');
            xlim(axesHandle, [-0.5 nGen-0.5]);
            ylim(axesHandle, [0 75]);

            expFunc = @(a, b, x) a * exp(x / b);
            [f, gof] = fit((0:nGen-1)', nCells, ...
                fittype(expFunc),  ...
                fitoptions('Method', 'NonlinearLeastSquares', 'lower', [0.5 * nCellFirstGen  0.5], 'upper', [2.0 * nCellFirstGen  3], 'StartPoint', [nCellFirstGen 1/log(2)]));
            plot(axesHandle, (0:0.1:nGen-1)', expFunc(f.a, f.b, (0:0.1:nGen-1)'), 'Color', 'r');

            saveas(figHandle, [outDir filesep 'PopulationSize.pdf']);
            close(figHandle);

            %% survival
            %all families
            [axesHandle, figHandle] = PlotUtil.newAxesHandle();
            hold(axesHandle, 'on');

            errorbar(0:nGen-1, fracSurvive * 100, sqrt(fracSurvive .* (1 - fracSurvive)) ./ sqrt(nCells) * 100, 'b.', 'Parent', axesHandle);
            xlabel(axesHandle, 'Generation');
            ylabel(axesHandle, 'Survival (%)');
            xlim(axesHandle, [-0.5 nGen-0.5]);
            ylim(axesHandle, [0 100])

            [f, gof] = fit((0:nGen-1)', fracSurvive * 100, ...
                fittype('poly1'), ...
                fitoptions('Method', 'LinearLeastSquares', 'Weights', sqrt(fracSurvive .* (1 - fracSurvive)) ./ sqrt(nCells) * 100));
            plot(axesHandle, (0:0.1:nGen-1)', (0:0.1:nGen-1)' * f.p1 + f.p2, 'Color', 'r');

            saveas(figHandle, [outDir filesep 'Survival.pdf']);
            close(figHandle);

            %left/right
            nCellsLR = NaN(nGen, 2);
            nSurviveLR = NaN(nGen, 2);
            for iGen = 0:nGen-1
                nCellsLR(iGen + 1, 1) = sum(ancestry(1:2:end, MultiGenerations.ANCESTRY_COLIDX_GEN) == iGen & finSimTfs(1:2:end, 1));
                nCellsLR(iGen + 1, 2) = sum(ancestry(2:2:end, MultiGenerations.ANCESTRY_COLIDX_GEN) == iGen & finSimTfs(2:2:end, 1));
                nSurviveLR(iGen + 1, 1) = sum(ancestry(1:2:end, MultiGenerations.ANCESTRY_COLIDX_GEN) == iGen & divSimTfs(1:2:end, 1));
                nSurviveLR(iGen + 1, 2) = sum(ancestry(2:2:end, MultiGenerations.ANCESTRY_COLIDX_GEN) == iGen & divSimTfs(2:2:end, 1));
            end
            fracSurviveLR = nSurviveLR ./ nCellsLR;

            [axesHandle, figHandle] = PlotUtil.newAxesHandle();

            h = plot(axesHandle, 1:nGen-1, fracSurviveLR(2:end, :) * 100);
            legend(h, {'Left', 'Right'}, 'Location', 'NorthEastOutside');
            xlabel(axesHandle, 'Generation');
            ylabel(axesHandle, 'Survival (%)');
            xlim(axesHandle, [0.5 nGen-0.5]);
            ylim(axesHandle, [0 100])

            saveas(figHandle, [outDir filesep 'Survival-LeftRight.pdf']);
            close(figHandle);

            %by family
            [axesHandle, figHandle] = PlotUtil.newAxesHandle();
            hold(axesHandle, 'on');

            colors = MultiGenerations.calcRedGreenColors(nCellFirstGen);
            labels = cell(nCellFirstGen, 1);
            h = zeros(nCellFirstGen, 1);
            for i = 1:nCellFirstGen
                tfs = finSimTfs & ancestry(:, MultiGenerations.ANCESTRY_COLIDX_FAMILY) == i;
                fracSurvive = zeros(nGen, 1);
                for iGen = 0:nGen-1
                    fracSurvive(iGen + 1) = ...
                        sum(ancestry(divSimTfs & tfs, MultiGenerations.ANCESTRY_COLIDX_GEN) == iGen) / ...
                        sum(ancestry(tfs, MultiGenerations.ANCESTRY_COLIDX_GEN) == iGen);
                end
                h(i) = plot(axesHandle, 0:nGen-1, fracSurvive * 100, 'Color', colors(i, :));
                labels{i} = num2str(i);
            end

            legend(h, labels, 'Location', 'NorthEastOutside');
            xlabel(axesHandle, 'Generation');
            ylabel(axesHandle, 'Survival (%)');
            xlim(axesHandle, [-0.5 nGen-0.5]);
            ylim(axesHandle, [0 100])

            saveas(figHandle, [outDir filesep 'Survival-ByFamily.pdf']);
            close(figHandle);

            %% left/right chromosome bias
            %- bound protein
            %- methylation
            %- superhelical density

            gen_dnaBndMons = NaN(nGen, 2, sum(dnaBndMonTfs));
            gen_dnaBndCpxs = NaN(nGen, 2, sum(dnaBndCpxTfs));
            gen_otherBndCpxs = NaN(nGen, 2, sum(otherBndCpxTfs));

            for iGen = 0:nGen - 1
                lefts  = find( isodd(finSimIdxs) & ancestry(finSimIdxs, MultiGenerations.ANCESTRY_COLIDX_GEN) == iGen);
                rights = find(~isodd(finSimIdxs) & ancestry(finSimIdxs, MultiGenerations.ANCESTRY_COLIDX_GEN) == iGen);

                gen_dnaBndMons(iGen+1, 1, :) = mean(dnaBndMons(lefts, :), 1);
                gen_dnaBndMons(iGen+1, 2, :) = mean(dnaBndMons(rights, :), 1);

                gen_dnaBndCpxs(iGen+1, 1, :) = mean(dnaBndCpxs(lefts, :), 1);
                gen_dnaBndCpxs(iGen+1, 2, :) = mean(dnaBndCpxs(rights, :), 1);

                gen_otherBndCpxs(iGen+1, 1, :) = mean(otherBndCpxs(lefts, :), 1);
                gen_otherBndCpxs(iGen+1, 2, :) = mean(otherBndCpxs(rights, :), 1);
            end

            [axesHandle, figHandle] = PlotUtil.newAxesHandle();
            hold(axesHandle, 'on');
            h = plot(axesHandle, (1:nGen-1)', permute(gen_dnaBndMons(2:end, 1, :), [1 3 2]));
            g = plot(axesHandle, (1:nGen-1)', permute(gen_dnaBndMons(2:end, 2, :), [1 3 2]));
            colors = MultiGenerations.calcRedGreenColors(size(gen_dnaBndMons, 3));
            for i = 1:size(gen_dnaBndMons, 3)
                set(h(i), 'LineStyle', '-', 'Color', colors(i, :))
                set(g(i), 'LineStyle', ':', 'Color', colors(i, :))
            end
            legend(h, pm.wholeCellModelIDs(pm.boundIndexs(dnaBndMonTfs)), 'Location', 'NorthEastOutside', 'Interpreter', 'none');
            xlim(axesHandle, [0.5 nGen-0.5])
            xlabel(axesHandle, 'Generation');
            ylabel(axesHandle, 'Protein monomers')
            saveas(figHandle, [outDir filesep 'DNABoundMonomerLeftRightBias.pdf']);
            close(figHandle);

            [axesHandle, figHandle] = PlotUtil.newAxesHandle();
            hold(axesHandle, 'on');
            h = plot(axesHandle, (1:nGen-1)', permute(gen_dnaBndCpxs(2:end, 1, :), [1 3 2]));
            g = plot(axesHandle, (1:nGen-1)', permute(gen_dnaBndCpxs(2:end, 2, :), [1 3 2]));
            colors = MultiGenerations.calcRedGreenColors(size(gen_dnaBndCpxs, 3));
            for i = 1:size(gen_dnaBndCpxs, 3)
                set(h(i), 'LineStyle', '-', 'Color', colors(i, :))
                set(g(i), 'LineStyle', ':', 'Color', colors(i, :))
            end
            legend(h, pc.wholeCellModelIDs(pc.boundIndexs(dnaBndCpxTfs)), 'Location', 'NorthEastOutside', 'Interpreter', 'none');
            xlim(axesHandle, [0.5 nGen-0.5])
            xlabel(axesHandle, 'Generation');
            ylabel(axesHandle, 'Protein complexes')
            saveas(figHandle, [outDir filesep 'DNABoundComplexLeftRightBias.pdf']);
            close(figHandle);

            [axesHandle, figHandle] = PlotUtil.newAxesHandle();
            hold(axesHandle, 'on');
            h = plot(axesHandle, (1:nGen-1)', permute(gen_otherBndCpxs(2:end, 1, :), [1 3 2]));
            g = plot(axesHandle, (1:nGen-1)', permute(gen_otherBndCpxs(2:end, 2, :), [1 3 2]));
            colors = MultiGenerations.calcRedGreenColors(size(gen_otherBndCpxs, 3));
            for i = 1:size(gen_otherBndCpxs, 3)
                set(h(i), 'LineStyle', '-', 'Color', colors(i, :))
                set(g(i), 'LineStyle', ':', 'Color', colors(i, :))
            end
            legend(h, pc.wholeCellModelIDs(pc.boundIndexs(otherBndCpxTfs)), 'Location', 'NorthEastOutside', 'Interpreter', 'none');
            xlim(axesHandle, [0.5 nGen-0.5])
            xlabel(axesHandle, 'Generation');
            ylabel(axesHandle, 'Protein complexes')
            saveas(figHandle, [outDir filesep 'OtherBoundComplexLeftRightBias.pdf']);
            close(figHandle);

            %% growth, mass, cell cycle phase durations
            for iProp = 1:size(simData, 2)
                %all families
                [axesHandle, figHandle] = PlotUtil.newAxesHandle();
                boxplot(axesHandle, simData(finSimTfs, iProp), ancestry(finSimTfs, MultiGenerations.ANCESTRY_COLIDX_GEN));
                title(axesHandle, [propNames{iProp, 1} ' - All Cells'])
                xlabel(axesHandle, 'Generation');
                ylabel(axesHandle, propNames{iProp, 2});
                saveas(figHandle, [outDir filesep propNames{iProp, 1} 'BoxPlot-AllCells.pdf']);
                close(figHandle);

                [axesHandle, figHandle] = PlotUtil.newAxesHandle();
                boxplot(axesHandle, simData(divSimTfs, iProp), ancestry(divSimTfs, MultiGenerations.ANCESTRY_COLIDX_GEN));
                title(axesHandle, [propNames{iProp, 1} ' - Successful Divisions'])
                xlabel(axesHandle, 'Generation');
                ylabel(axesHandle, propNames{iProp, 2});
                saveas(figHandle, [outDir filesep propNames{iProp, 1} 'BoxPlot-SuccessfulCells.pdf']);
                close(figHandle);

                [axesHandle, figHandle] = PlotUtil.newAxesHandle();
                hold(axesHandle, 'on');
                h = zeros(nGen, 1);
                x = simData(finSimTfs, iProp);
                edges = linspace(min(x), max(x), 20);
                labels = cell(nGen, 1);
                colors = MultiGenerations.calcRedGreenColors(nGen);
                for iGen = 0:nGen-1
                    tfs = finSimTfs & ancestry(:, MultiGenerations.ANCESTRY_COLIDX_GEN) == iGen;
                    cnts = histc(simData(tfs, iProp), edges);
                    cnts = cnts / sum(tfs);
                    h(iGen + 1) = plot(axesHandle, edges, cnts, 'Color', colors(iGen + 1, :));
                    labels{iGen + 1} = num2str(iGen);
                end
                title(axesHandle, [propNames{iProp, 1} ' - All Cells'])
                xlabel(axesHandle, propNames{iProp, 2});
                ylabel(axesHandle, 'Frequency');
                legend(h, labels);
                saveas(figHandle, [outDir filesep propNames{iProp, 1} 'Histogram-AllCells.pdf']);
                close(figHandle);

                %left/right bias
                [axesHandle, figHandle] = PlotUtil.newAxesHandle();
                lefts = false(size(finSimTfs));
                lefts(1:2:end, 1) = true;
                yData = NaN(nGen, 1);
                for iGen = 0:nGen-1
                    yData(iGen + 1, 1) = nanmean(simData( lefts & finSimTfs & ancestry(:, MultiGenerations.ANCESTRY_COLIDX_GEN) == iGen, iProp));
                    yData(iGen + 1, 2) = nanmean(simData(~lefts & finSimTfs & ancestry(:, MultiGenerations.ANCESTRY_COLIDX_GEN) == iGen, iProp));
                end
                h = plot(axesHandle, 1:nGen-1, yData(2:end, :));
                legend(h, {'Left', 'Right'}, 'Location', 'NorthEastOutside');
                title(axesHandle, [propNames{iProp, 1}])
                xlabel(axesHandle, 'Generation');
                ylabel(axesHandle, propNames{iProp, 2});
                xlim(axesHandle, [0.5 nGen-0.5]);
                saveas(figHandle, [outDir filesep propNames{iProp, 1} 'MeanLeftRight.pdf']);
                close(figHandle);

                %by family
                [axesHandle, figHandle] = PlotUtil.newAxesHandle();
                hold(axesHandle, 'on');
                colors = MultiGenerations.calcRedGreenColors(nCellFirstGen);
                labels = cell(nCellFirstGen, 1);
                h = zeros(nCellFirstGen, 1);
                for i = 1:nCellFirstGen
                    tfs = finSimTfs & ancestry(:, MultiGenerations.ANCESTRY_COLIDX_FAMILY) == i;
                    h(i) = plot(axesHandle, ancestry(tfs, MultiGenerations.ANCESTRY_COLIDX_GEN), ...
                        simData(tfs, iProp), '.', 'Color', colors(i, :));

                    labels{i} = num2str(i);
                end
                legend(h, labels, 'Location', 'NorthEastOutside');
                title(axesHandle, [propNames{iProp, 1} ' - All Cells'])
                xlabel(axesHandle, 'Generation');
                ylabel(axesHandle, propNames{iProp, 2});
                xlim(axesHandle, [-0.5 nGen-0.5]);
                saveas(figHandle, [outDir filesep propNames{iProp, 1} 'Distribution.pdf']);
                close(figHandle);

                [axesHandle, figHandle] = PlotUtil.newAxesHandle();
                hold(axesHandle, 'on');
                colors = MultiGenerations.calcRedGreenColors(nCellFirstGen);
                labels = cell(nCellFirstGen, 1);
                h = zeros(nCellFirstGen, 1);
                for i = 1:nCellFirstGen
                    tfs = finSimTfs & ancestry(:, MultiGenerations.ANCESTRY_COLIDX_FAMILY) == i;
                    gens = ancestry(tfs, MultiGenerations.ANCESTRY_COLIDX_GEN);
                    x = simData(tfs, iProp);
                    avgs = zeros(nGen, 1);
                    for iGen = 0:nGen-1
                        avgs(iGen + 1) = nanmean(x(gens == iGen));
                    end
                    h(i) = plot(axesHandle, 0:nGen-1, avgs, 'Color', colors(i, :));

                    labels{i} = num2str(i);
                end
                legend(h, labels, 'Location', 'NorthEastOutside');
                title(axesHandle, [propNames{iProp, 1} ' - All Cells'])
                xlabel(axesHandle, 'Generation');
                ylabel(axesHandle, propNames{iProp, 2});
                xlim(axesHandle, [-0.25 nGen-0.75]);
                saveas(figHandle, [outDir filesep propNames{iProp, 1} 'MeanByLineage-AllCells.pdf']);
                close(figHandle);

                [axesHandle, figHandle] = PlotUtil.newAxesHandle();
                hold(axesHandle, 'on');
                colors = MultiGenerations.calcRedGreenColors(nCellFirstGen);
                labels = cell(nCellFirstGen, 1);
                h = zeros(nCellFirstGen, 1);
                for i = 1:nCellFirstGen
                    tfs = divSimTfs & ancestry(:, MultiGenerations.ANCESTRY_COLIDX_FAMILY) == i;
                    gens = ancestry(tfs, MultiGenerations.ANCESTRY_COLIDX_GEN);
                    x = simData(tfs, iProp);
                    avgs = zeros(nGen, 1);
                    for iGen = 0:nGen-1
                        avgs(iGen + 1) = nanmean(x(gens == iGen));
                    end
                    h(i) = plot(axesHandle, 0:nGen-1, avgs, 'Color', colors(i, :));

                    labels{i} = num2str(i);
                end
                legend(h, labels, 'Location', 'NorthEastOutside');
                title(axesHandle, [propNames{iProp, 1} ' - Successful Divisions'])
                xlabel(axesHandle, 'Generation');
                ylabel(axesHandle, propNames{iProp, 2});
                xlim(axesHandle, [-0.25 nGen-0.75]);
                saveas(figHandle, [outDir filesep propNames{iProp, 1} 'MeanByLineage-SuccessfulCells.pdf']);
                close(figHandle);
            end

            %% mother-daughter correlations within cell cycle
            deltaGenerations = relations(finSimTfs, finSimTfs, 1);
            removal = relations(finSimTfs, finSimTfs, 2);

            [idxs1, idxs2] = find(deltaGenerations == 1 & removal == 0);

            ensemble = SimulationEnsemble(simBatchDir, cell(0, 1), [1 2], find(finSimTfs));
            simEndTimes = ensemble.stateData.simulationEndTimes;
            states = SimulationEnsemble.load(simBatchDir, {'MetabolicReaction' 'growth'}, [], [], 1, 'extract', find(finSimTfs));
            growth = permute(states.MetabolicReaction.growth, [4 3 1 2]);
            growth(growth == 0) = NaN;
            growth0 = growth(sub2ind(size(growth), (1:numel(simEndTimes))', simEndTimes));

            r = zeros(1, size(growth, 2));
            n = zeros(1, size(growth, 2));
            for i = 1:size(growth, 2)
                tfs = ~isnan(growth(idxs2, i));
                r(i) = corr(growth0(idxs1(tfs)), growth(idxs2(tfs), i), 'type', 'Pearson');
                n(i) = sum(simEndTimes >= i);
            end

            [~, figHandle] = PlotUtil.newAxesHandle();
            clf(figHandle);

            axesHandle = subplot(2, 1, 1);
            plot(axesHandle, (1:size(growth, 2))/3600, r);
            xlabel(axesHandle, 'Time (h)');
            ylabel(axesHandle, 'Corr');

            axesHandle = subplot(2, 1, 2);
            plot(axesHandle, (1:size(growth, 2))/3600, n);
            xlabel(axesHandle, 'Time (h)');
            ylabel(axesHandle, 'No. cells');

            saveas(figHandle, [outDir filesep 'GrowthCorrelationDecay.pdf']);
            close(figHandle);

            %% descendents
            deltaGenerations = relations(:, :, 1);
            removal = relations(:, :, 2);

            colors = MultiGenerations.calcRedGreenColors(nGen - 1);
            for iProp = 1:size(propNames, 1)
                %all relationships
                [~, figHandle] = PlotUtil.newAxesHandle();
                clf(figHandle);

                axesHandle = subplot(2, 1, 1);
                hold(axesHandle, 'on');

                r = zeros(nGen - 1, 1);
                h = zeros(nGen - 1, 1);
                labels = cell(nGen - 1, 1);
                for iGen = 1:nGen-1
                    [idxs1, idxs2] = find(deltaGenerations == iGen & removal == 0);
                    tfs = finSimTfs(idxs1) & finSimTfs(idxs2);
                    idxs1 = idxs1(tfs, :);
                    idxs2 = idxs2(tfs, :);

                    r(iGen) = corr(simData(idxs1, iProp), simData(idxs2, iProp), 'type', 'Pearson');
                    h(iGen) = plot(simData(idxs1, iProp), simData(idxs2, iProp), '.', 'Color', colors(iGen, :));
                    labels{iGen} = sprintf('{\\Delta}g=%d (n=%d)', iGen, numel(idxs1));
                end

                legend(h, labels, 'Location', 'NorthEastOutside');
                title(axesHandle, ['All descendents - ' propNames{iProp, 1}])
                xlabel(axesHandle, 'Parent');
                ylabel(axesHandle, 'Child');

                axesHandle = subplot(2, 1, 2);
                plot(axesHandle, 1:nGen-1, r);
                line([0.5 nGen-0.5], [0 0], 'Color', 0.25 * [1 1 1], 'Parent', axesHandle);
                xlabel(axesHandle, '{\Delta}Generation');
                ylabel(axesHandle, 'Corr');
                xlim(axesHandle, [0.5 nGen-0.5]);
                ylim(axesHandle, [-1 1]);
                set(axesHandle, 'XTick', 1:nGen-1);
                box(axesHandle, 'off');

                saveas(figHandle, [outDir filesep 'Descendents-all-' propNames{iProp, 1} '.pdf']);
                close(figHandle);

                %relationships from founders
                [~, figHandle] = PlotUtil.newAxesHandle();
                clf(figHandle);

                axesHandle = subplot(2, 1, 1);
                hold(axesHandle, 'on');

                r = zeros(nGen - 1, 1);
                h = zeros(nGen - 1, 1);
                labels = cell(nGen - 1, 1);
                for iGen = 1:nGen-1
                    [idxs1, idxs2] = find(deltaGenerations == iGen & removal == 0);
                    tfs = finSimTfs(idxs1) & finSimTfs(idxs2);
                    tfs = tfs & ancestry(idxs1, MultiGenerations.ANCESTRY_COLIDX_GEN) == 0;
                    idxs1 = idxs1(tfs, :);
                    idxs2 = idxs2(tfs, :);

                    r(iGen) = corr(simData(idxs1, iProp), simData(idxs2, iProp), 'type', 'Pearson');
                    h(iGen) = plot(simData(idxs1, iProp), simData(idxs2, iProp), '.', 'Color', colors(iGen, :));
                    labels{iGen} = sprintf('{\\Delta}g=%d (n=%d)', iGen, numel(idxs1));
                end

                legend(h, labels, 'Location', 'NorthEastOutside');
                title(axesHandle, ['Generations from founder - ' propNames{iProp, 1}])
                xlabel(axesHandle, 'Parent');
                ylabel(axesHandle, 'Child');

                axesHandle = subplot(2, 1, 2);
                plot(axesHandle, 1:nGen-1, r);
                line([0.5 nGen-0.5], [0 0], 'Color', 0.25 * [1 1 1], 'Parent', axesHandle);
                xlabel(axesHandle, '{\Delta}Generation');
                ylabel(axesHandle, 'Corr');
                xlim(axesHandle, [0.5 nGen-0.5]);
                ylim(axesHandle, [-1 1]);
                set(axesHandle, 'XTick', 1:nGen-1);
                box(axesHandle, 'off');

                saveas(figHandle, [outDir filesep 'Descendents-FromFounders-' propNames{iProp, 1} '.pdf']);
                close(figHandle);

                %relationships from founders
                [~, figHandle] = PlotUtil.newAxesHandle();
                clf(figHandle);

                axesHandle = subplot(2, 1, 1);
                hold(axesHandle, 'on');

                r = zeros(nGen - 1, 1);
                h = zeros(nGen - 1, 1);
                labels = cell(nGen - 1, 1);
                for iGen = 1:nGen-1
                    [idxs1, idxs2] = find(deltaGenerations == 1 & removal == 0);
                    tfs = finSimTfs(idxs1) & finSimTfs(idxs2);
                    tfs = tfs & ancestry(idxs1, MultiGenerations.ANCESTRY_COLIDX_GEN) == (iGen - 1);
                    idxs1 = idxs1(tfs, :);
                    idxs2 = idxs2(tfs, :);

                    r(iGen) = corr(simData(idxs1, iProp), simData(idxs2, iProp), 'type', 'Pearson');
                    h(iGen) = plot(simData(idxs1, iProp), simData(idxs2, iProp), '.', 'Color', colors(iGen, :));
                    labels{iGen} = sprintf('g_1=%d (n=%d)', iGen, numel(idxs1));
                end

                legend(h, labels, 'Location', 'NorthEastOutside');
                title(axesHandle, propNames{iProp, 1})
                xlabel(axesHandle, 'Parent');
                ylabel(axesHandle, 'Child');

                axesHandle = subplot(2, 1, 2);
                plot(axesHandle, 1:nGen-1, r);
                line([0.5 nGen-0.5], [0 0], 'Color', 0.25 * [1 1 1], 'Parent', axesHandle);
                xlabel(axesHandle, '{\Delta}Generation');
                ylabel(axesHandle, 'Corr');
                xlim(axesHandle, [0.5 nGen-0.5]);
                ylim(axesHandle, [-1 1]);
                set(axesHandle, 'XTick', 1:nGen-1);
                box(axesHandle, 'off');

                saveas(figHandle, [outDir filesep 'Descendents-Children-' propNames{iProp, 1} '.pdf']);
                close(figHandle);
            end

            %% siblings, cousins
            deltaGenerations = relations(:, :, 1);
            removal = relations(:, :, 2);

            colors = MultiGenerations.calcRedGreenColors(nGen - 1);
            for iProp = 1:size(propNames, 1)
                %all cousins
                [~, figHandle] = PlotUtil.newAxesHandle();
                clf(figHandle);

                axesHandle = subplot(2, 1, 1);
                hold(axesHandle, 'on');

                r = zeros(nGen - 1, 1);
                h = zeros(nGen - 1, 1);
                labels = cell(nGen - 1, 1);
                for iGen = 1:nGen-1
                    [idxs1, idxs2] = find(deltaGenerations == 0 & removal == iGen);
                    tfs = finSimTfs(idxs1) & finSimTfs(idxs2);
                    idxs1 = idxs1(tfs, :);
                    idxs2 = idxs2(tfs, :);

                    r(iGen) = corr(simData(idxs1, iProp), simData(idxs2, iProp), 'type', 'Pearson');
                    h(iGen) = plot(simData(idxs1, iProp), simData(idxs2, iProp), '.', 'Color', colors(iGen, :));
                    labels{iGen} = sprintf('{\\Delta}r=%d (n=%d)', iGen, numel(idxs1));
                end

                legend(h, labels, 'Location', 'NorthEastOutside');
                title(axesHandle, ['All cousins - ' propNames{iProp, 1}])
                xlabel(axesHandle, 'Parent');
                ylabel(axesHandle, 'Child');

                axesHandle = subplot(2, 1, 2);
                plot(axesHandle, 1:nGen-1, r);
                line([0.5 nGen-0.5], [0 0], 'Color', 0.25 * [1 1 1], 'Parent', axesHandle);
                xlabel(axesHandle, 'Removal');
                ylabel(axesHandle, 'Corr');
                xlim(axesHandle, [0.5 nGen-0.5]);
                ylim(axesHandle, [-1 1]);
                set(axesHandle, 'XTick', 1:nGen-1);
                box(axesHandle, 'off');

                saveas(figHandle, [outDir filesep 'Cousins-all-' propNames{iProp, 1} '.pdf']);
                close(figHandle);

                %cousins of last generation
                [~, figHandle] = PlotUtil.newAxesHandle();
                clf(figHandle);

                axesHandle = subplot(2, 1, 1);
                hold(axesHandle, 'on');

                r = zeros(nGen - 1, 1);
                h = zeros(nGen - 1, 1);
                labels = cell(nGen - 1, 1);
                for iGen = 1:nGen-1
                    [idxs1, idxs2] = find(deltaGenerations == 0 & removal == iGen);
                    tfs = finSimTfs(idxs1) & finSimTfs(idxs2) & ancestry(idxs1, MultiGenerations.ANCESTRY_COLIDX_GEN) == iGen;
                    idxs1 = idxs1(tfs, :);
                    idxs2 = idxs2(tfs, :);

                    r(iGen) = corr(simData(idxs1, iProp), simData(idxs2, iProp), 'type', 'Pearson');
                    h(iGen) = plot(simData(idxs1, iProp), simData(idxs2, iProp), '.', 'Color', colors(iGen, :));
                    labels{iGen} = sprintf('{\\Delta}r=%d (n=%d)', iGen, numel(idxs1));
                end

                legend(h, labels, 'Location', 'NorthEastOutside');
                title(axesHandle, ['Cousins of last generation - ' propNames{iProp, 1}])
                xlabel(axesHandle, 'Parent');
                ylabel(axesHandle, 'Child');

                axesHandle = subplot(2, 1, 2);
                plot(axesHandle, 1:nGen-1, r);
                line([0.5 nGen-0.5], [0 0], 'Color', 0.25 * [1 1 1], 'Parent', axesHandle);
                xlabel(axesHandle, 'Removal');
                ylabel(axesHandle, 'Corr');
                xlim(axesHandle, [0.5 nGen-0.5]);
                ylim(axesHandle, [-1 1]);
                set(axesHandle, 'XTick', 1:nGen-1);
                box(axesHandle, 'off');

                saveas(figHandle, [outDir filesep 'Cousins-LastGeneration-' propNames{iProp, 1} '.pdf']);
                close(figHandle);

                %first cousins by generation
                [~, figHandle] = PlotUtil.newAxesHandle();
                clf(figHandle);

                axesHandle = subplot(2, 1, 1);
                hold(axesHandle, 'on');

                r = NaN(nGen - 1, 1);
                h = NaN(nGen - 1, 1);
                labels = cell(nGen - 1, 1);
                for iGen = 1:nGen-1
                    [idxs1, idxs2] = find(deltaGenerations == 0 & removal == 1);
                    tfs = finSimTfs(idxs1) & finSimTfs(idxs2);
                    tfs = tfs & ancestry(idxs1, MultiGenerations.ANCESTRY_COLIDX_GEN) == iGen;
                    if ~any(tfs)
                        continue;
                    end
                    idxs1 = idxs1(tfs, :);
                    idxs2 = idxs2(tfs, :);

                    r(iGen) = corr(simData(idxs1, iProp), simData(idxs2, iProp), 'type', 'Pearson');
                    h(iGen) = plot(simData(idxs1, iProp), simData(idxs2, iProp), '.', 'Color', colors(iGen, :));
                    labels{iGen} = sprintf('g=%d (n=%d)', iGen, numel(idxs1));
                end

                legend(h(~isnan(h)), labels(~isnan(h)), 'Location', 'NorthEastOutside');
                title(axesHandle, ['All First Cousins - ' propNames{iProp, 1}])
                xlabel(axesHandle, 'Parent');
                ylabel(axesHandle, 'Child');

                axesHandle = subplot(2, 1, 2);
                plot(axesHandle, 1:nGen-1, r);
                line([0.5 nGen-0.5], [0 0], 'Color', 0.25 * [1 1 1], 'Parent', axesHandle);
                xlabel(axesHandle, '{\Delta}Generation');
                ylabel(axesHandle, 'Corr');
                xlim(axesHandle, [0.5 nGen-0.5]);
                ylim(axesHandle, [-1 1]);
                set(axesHandle, 'XTick', 1:nGen-1);
                box(axesHandle, 'off');

                saveas(figHandle, [outDir filesep 'FirstCousins-all-' propNames{iProp, 1} '.pdf']);
                close(figHandle);
            end

            %% Varying effective DnaA copy number
            MultiGenerations.runEffectiveDnaACopyNumberAnalysis();
        end

        function runEffectiveDnaACopyNumberAnalysis(outDir)
            import edu.stanford.covert.cell.sim.analysis.MultiGenerations;
            import edu.stanford.covert.cell.sim.util.PlotUtil;
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;

            [simBatchDirs, simData, ~, ~, finSimTfs, ~, propNames] = MultiGenerations.cacheSimData();

            nSimBatch = size(simBatchDirs, 1);

            %% create output directory
            if nargin < 1
                outDir = [SimulationDiskUtil.getBaseDir() filesep 'multiGenerations'];
            end
            if ~exist(outDir, 'dir')
                mkdir(outDir)
            end

            %% analyze
            colors = MultiGenerations.calcRedGreenColors(nSimBatch);
            for iProp = 1:size(propNames, 1)
                [axesHandle, figHandle] = PlotUtil.newAxesHandle();
                hold(axesHandle, 'on');

                h = zeros(nSimBatch, 1);
                for iSimBatch = 1:nSimBatch
                    nGen = simBatchDirs{iSimBatch, 2};
                    nCellFirstGen = simBatchDirs{iSimBatch, 3};
                    ancestry = MultiGenerations.calcAncestry(nGen, nCellFirstGen);

                    avgs = zeros(1, nGen);
                    for iGen = 0:nGen - 1
                        avgs(1, iGen + 1) = nanmean(simData{iSimBatch}(finSimTfs{iSimBatch} & ancestry(:, MultiGenerations.ANCESTRY_COLIDX_GEN) == iGen, iProp));
                    end
                    h(iSimBatch) = plot(axesHandle, 0:nGen-1, avgs, 'Color', colors(iSimBatch, :));
                end

                xlim(axesHandle, [-0.5 nGen-0.5]);
                title(axesHandle, ['DnaA Copy Number vs ' propNames{iProp, 1}]);
                xlabel(axesHandle, 'Generation');
                ylabel(axesHandle, propNames{iProp, 2});
                legend(h, cellfun(@num2str, simBatchDirs(:, 4), 'UniformOutput', false), 'Location', 'NorthEastOutside');

                saveas(figHandle, [outDir filesep 'DnaACopyNumberVs' propNames{iProp, 1} '.pdf']);
                close(figHandle);
            end
        end

        function [simBatchDirs, simData, simStartTimes, divSimTfs, finSimTfs, finSimIdxs, propNames, ...
                dnaBndMons, dnaBndCpxs, otherBndCpxs, dnaBndMonTfs, dnaBndCpxTfs, otherBndCpxTfs] = cacheSimData(simBatchDirs)
            import edu.stanford.covert.cell.sim.analysis.MultiGenerations;
            import edu.stanford.covert.cell.sim.util.CachedSimulationObjectUtil;
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;

            if nargin >= 1
                if ~iscell(simBatchDirs)
                    simBatchDirs = {simBatchDirs};
                end
            else
                simBatchDirs = {
                    '2012_11_15_18_48_23'  3 8  3
                    '2012_10_24_00_49_53'  6 8  1
                    };
            end

            nSimBatch = size(simBatchDirs, 1);

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

            %% get data
            simData = cell(nSimBatch, 1);
            simStartTimes = cell(nSimBatch, 1);
            divSimTfs = cell(nSimBatch, 1);
            finSimTfs = cell(nSimBatch, 1);
            finSimIdxs = cell(nSimBatch, 1);
            dnaBndMons = cell(nSimBatch, 1);
            dnaBndCpxs = cell(nSimBatch, 1);
            otherBndCpxs = cell(nSimBatch, 1);
            dnaBndMonTfs = cell(nSimBatch, 1);
            dnaBndCpxTfs = cell(nSimBatch, 1);
            otherBndCpxTfs = cell(nSimBatch, 1);
            for iSimBatch = 1:nSimBatch
                simBatchDir = simBatchDirs{iSimBatch, 1};
                if exist([SimulationDiskUtil.getBaseDir() filesep simBatchDir filesep 'MultiGenerationsData.mat'], 'file')
                    tmp = load([SimulationDiskUtil.getBaseDir() filesep simBatchDir filesep 'MultiGenerationsData.mat']);
                else
                    nGen = simBatchDirs{iSimBatch, 2};
                    nCellFirstGen = simBatchDirs{iSimBatch, 3};

                    [tmp_simData, tmp_simStartTimes, tmp_divSimTfs, tmp_finSimTfs, tmp_finSimIdxs, tmp_propNames, ...
                        tmp_dnaBndMons, tmp_dnaBndCpxs, tmp_otherBndCpxs, tmp_dnaBndMonTfs, tmp_dnaBndCpxTfs, tmp_otherBndCpxTfs] = ...
                        MultiGenerations.getSimData(simBatchDir, nGen, nCellFirstGen, sim);

                    tmp = struct;
                    tmp.simData = tmp_simData;
                    tmp.simStartTimes = tmp_simStartTimes;
                    tmp.divSimTfs = tmp_divSimTfs;
                    tmp.finSimTfs = tmp_finSimTfs;
                    tmp.finSimIdxs = tmp_finSimIdxs;
                    tmp.propNames = tmp_propNames;
                    tmp.dnaBndMons = tmp_dnaBndMons;
                    tmp.dnaBndCpxs = tmp_dnaBndCpxs;
                    tmp.otherBndCpxs = tmp_otherBndCpxs;
                    tmp.dnaBndMonTfs = tmp_dnaBndMonTfs;
                    tmp.dnaBndCpxTfs = tmp_dnaBndCpxTfs;
                    tmp.otherBndCpxTfs = tmp_otherBndCpxTfs;
                    save([SimulationDiskUtil.getBaseDir() filesep simBatchDir filesep 'MultiGenerationsData.mat'], '-struct', 'tmp');
                end
                simData{iSimBatch} = tmp.simData;
                simStartTimes{iSimBatch} = tmp.simStartTimes;
                divSimTfs{iSimBatch} = tmp.divSimTfs;
                finSimTfs{iSimBatch} = tmp.finSimTfs;
                finSimIdxs{iSimBatch} = tmp.finSimIdxs;

                dnaBndMons{iSimBatch} = tmp.dnaBndMons;
                dnaBndCpxs{iSimBatch} = tmp.dnaBndCpxs;
                otherBndCpxs{iSimBatch} = tmp.otherBndCpxs;
                dnaBndMonTfs{iSimBatch} = tmp.dnaBndMonTfs;
                dnaBndCpxTfs{iSimBatch} = tmp.dnaBndCpxTfs;
                otherBndCpxTfs{iSimBatch} = tmp.otherBndCpxTfs;

                propNames = tmp.propNames;
            end
        end

        function [simData, simStartTimes, divSimTfs, finSimTfs, finSimIdxs, propNames, ...
                dnaBndMons, dnaBndCpxs, otherBndCpxs, ...
                dnaBndMonTfs, dnaBndCpxTfs, otherBndCpxTfs] = ...
                getSimData(simBatchDir, nGen, nCellFirstGen, sim)
            import edu.stanford.covert.cell.sim.analysis.MultiGenerations;
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;
            import edu.stanford.covert.cell.sim.util.SimulationEnsemble;
            import edu.stanford.covert.cell.sim.util.SummaryLogger;

            ancestry = MultiGenerations.calcAncestry(nGen, nCellFirstGen);
            nCells = size(ancestry, 1);

            g = sim.gene;
            c = sim.compartment;
            massState = sim.state('Mass');
            met = sim.state('Metabolite');
            pm = sim.state('ProteinMonomer');
            pc = sim.state('ProteinComplex');

            simStats = SummaryLogger.getSimulationStatistics([SimulationDiskUtil.getBaseDir() filesep simBatchDir], 1:nCells);
            simStats(:, SummaryLogger.SIM_STATUS_INDEX_INIT_GROWTH) = ...
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_INIT_GROWTH)' * ...
                massState.cellInitialDryWeight / (1 - massState.fractionWetWeight) * 3600 * 1e15;
            simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_INIT_TIME) = simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_INIT_TIME) / 3600;
            simStats(:, SummaryLogger.SIM_STATUS_INDEX_INIT_REP_TIME) = simStats(:, SummaryLogger.SIM_STATUS_INDEX_INIT_REP_TIME) / 3600;
            simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_TIME) = simStats(:, SummaryLogger.SIM_STATUS_INDEX_REP_TIME) / 3600;
            simStats(:, SummaryLogger.SIM_STATUS_INDEX_CYTOKINESIS_TIME) = simStats(:, SummaryLogger.SIM_STATUS_INDEX_CYTOKINESIS_TIME) / 3600;
            simStats(:, SummaryLogger.SIM_STATUS_INDEX_MASS_DOUBLING_TIME) = simStats(:, SummaryLogger.SIM_STATUS_INDEX_MASS_DOUBLING_TIME) / 3600;
            simStats(:, SummaryLogger.SIM_STATUS_INDEX_FIN_MASS) = simStats(:, SummaryLogger.SIM_STATUS_INDEX_FIN_MASS) * 1e15;
            simStats(:, SummaryLogger.SIM_STATUS_INDEX_INIT_MASS) = simStats(:, SummaryLogger.SIM_STATUS_INDEX_INIT_MASS) * 1e15;

            simStartTimes = zeros(size(ancestry, 1), 1);
            for i = 1:size(ancestry, 1)
                if isnan(ancestry(i, MultiGenerations.ANCESTRY_COLIDX_PARENT))
                    simStartTimes(i) = 0;
                else
                    iParent = ancestry(i, MultiGenerations.ANCESTRY_COLIDX_PARENT);
                    simStartTimes(i) = ...
                        simStartTimes(iParent) +  ...
                        simStats(iParent, SummaryLogger.SIM_STATUS_INDEX_CYTOKINESIS_TIME) + ...
                        1 / 3600;
                end

                if ~isnan(ancestry(i, MultiGenerations.ANCESTRY_COLIDX_CHILD1)) &&  ...
                        simStats(i, SummaryLogger.SIM_STATUS_INDEX_STATUS) ~= SummaryLogger.SIM_STATUS_COMPLETED_WITH_DIVISION
                    simStats(ancestry(i, MultiGenerations.ANCESTRY_COLIDX_CHILD1), SummaryLogger.SIM_STATUS_INDEX_STATUS) = SummaryLogger.SIM_STATUS_DIDNT_START;
                    simStats(ancestry(i, MultiGenerations.ANCESTRY_COLIDX_CHILD2), SummaryLogger.SIM_STATUS_INDEX_STATUS) = SummaryLogger.SIM_STATUS_DIDNT_START;
                end
            end

            divSimTfs = simStats(:, SummaryLogger.SIM_STATUS_INDEX_STATUS) == SummaryLogger.SIM_STATUS_COMPLETED_WITH_DIVISION;
            nonDivSimTfs = simStats(:, SummaryLogger.SIM_STATUS_INDEX_STATUS) == SummaryLogger.SIM_STATUS_COMPLETED_WITHOUT_DIVISION;
            finSimTfs = divSimTfs | nonDivSimTfs;
            finSimIdxs = find(finSimTfs);

            props = {'dnaA_total'; 'dnaA_boxes'; 'rnaPolymerases'; 'ribosomes';
                'amino_acids'; 'ntps'; 'rnas'; 'immatureRnas'; 'matureMonomers'; 'immatureMonomers';
                'matureComplexs'; 'immatureComplexs'; 'atp'; 'adp'; 'amp'};
            ensemble = SimulationEnsemble(simBatchDir, props, [1 2], finSimIdxs);

            rnaPolMonIdxs = pm.matureIndexs(any(pc.proteinComplexComposition(g.mRNAIndexs, pc.rnaPolymeraseIndexs(1), :), 3));
            stateNames = {
                'ProteinMonomer'  'counts'  rnaPolMonIdxs   c.cytosolIndexs
                };
            rnaPolSubunitCnts = SimulationEnsemble.load(simBatchDir, stateNames, 1, 1, 1, 'extract', finSimIdxs);

            stateNames = {
                'ProteinMonomer'  'counts'                pm.boundIndexs  c.cytosolIndexs
                'ProteinComplex'  'counts'                pc.boundIndexs  c.cytosolIndexs
                'Chromosome'      'damagedBases'          ':'  ':'
                'Chromosome'      'superhelicalDensity'   ':'  ':'
                };
            chromState = SimulationEnsemble.load(simBatchDir, stateNames, 1, 1, 1, 'extract', finSimIdxs);

            dnaBndMonTfs = any(any(any(chromState.ProteinMonomer.counts, 2), 3), 4);
            dnaBndCpxTfs = any(any(any(chromState.ProteinComplex.counts, 2), 3), 4) & ~any(any(pc.proteinComplexComposition(g.getIndexs('MG_469'), :, :), 1), 3)';
            otherBndCpxTfs = any(any(any(chromState.ProteinComplex.counts, 2), 3), 4) &  any(any(pc.proteinComplexComposition(g.getIndexs('MG_469'), :, :), 1), 3)';

            dnaBndMons = NaN(numel(finSimIdxs), sum(dnaBndMonTfs));
            dnaBndCpxs = NaN(numel(finSimIdxs), sum(dnaBndCpxTfs));
            otherBndCpxs = NaN(numel(finSimIdxs), sum(otherBndCpxTfs));
            methylations = NaN(numel(finSimIdxs), 1);
            sigmas = NaN(numel(finSimIdxs), 1);

            for i = 1:numel(finSimIdxs)
                dnaBndMons(i, :) = chromState.ProteinMonomer.counts(dnaBndMonTfs, 1, 1, i);
                dnaBndCpxs(i, :) = chromState.ProteinComplex.counts(dnaBndCpxTfs, 1, 1, i);
                otherBndCpxs(i, :) = chromState.ProteinComplex.counts(otherBndCpxTfs, 1, 1, i);
                methylations(i, 1) = full(sum(sum(chromState.Chromosome.damagedBases(:, :, 1, i) == met.m6ADIndexs, 1), 2));
                sigmas(i, 1) = full(chromState.Chromosome.superhelicalDensity(1, 1, 1, i));
            end

            simData = [
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_INIT_GROWTH) ...
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_FIN_MASS) ...
                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) ...
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_CYTOKINESIS_TIME) ...
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_INIT_MASS) ...
                simStats(:, SummaryLogger.SIM_STATUS_INDEX_FIN_MASS) ...
                zeros(size(simStats, 1), size(ensemble.stateData.values, 1) + size(rnaPolSubunitCnts.ProteinMonomer.counts, 1) + 2)
                ];
            simData(finSimTfs, 9:end) = [
                permute(ensemble.stateData.values(:, :, 1, :), [4 1 2 3]) ...
                permute(rnaPolSubunitCnts.ProteinMonomer.counts, [4 1 2 3]) ...
                methylations ...
                sigmas
                ];
            propNames = {
                'InitialGrowth'                   'Growth (fg h^{-1})'
                'FinalGrowth'                     'Growth (fg h^{-1})'
                'ReplicationInitiationDuration'   'Time (h)'
                'ReplicationDuration'             'Time (h)'
                'CytokinesisDuration'             'Time (h)'
                'CellCycleLength'                 'Time (h)'
                'InitialMass'                     'Mass (fg)'
                'FinalMass'                       'Mass (fg)'
                'DnaA'                            'Count'
                'BoundDnaA'                       'Count'
                'RNAPolymerases'                  'Count'
                'Ribosomes'                       'Count'
                'AminoAcids'                      'Count'
                'Ntps'                            'Count'
                'Rnas'                            'Count'
                'ImmatureRnas'                    'Count'
                'MatureMonomers'                  'Count'
                'ImmatureMonomers'                'Count'
                'MatureComplexes'                 'Count'
                'ImmatureComplexes'               'Count'
                'Atp'                             'Count'
                'Adp'                             'Count'
                'Amp'                             'Count'
                'RnaPolDelta'                     'Count'
                'RnaPolAlpha'                     'Count'
                'RnaPolBeta1'                     'Count'
                'RnaPolBeta2'                     'Count'
                'Methylations'                    'Count'
                'SuperhelicalDensity'             'dimensionless'
                };
        end

        %Returns matrix with five colums:
        %- cell index
        %- generation
        %- parent index
        %- child-1 index
        %- child-2 index
        function ancestry = calcAncestry(nGen, nCellFirstGen)
            import edu.stanford.covert.cell.sim.analysis.MultiGenerations;

            nCells = nCellFirstGen * (2^nGen - 1);

            ancestry = NaN(nCells, 5);

            genOffset = 0;
            iCell = 0;
            for iGen = 0:nGen-1
                for iGenCell = 1:nCellFirstGen * 2^iGen
                    iCell = iCell + 1;

                    if iGen == 0
                        iParent = NaN;
                        iFamily = iGenCell;
                    else
                        iParent = genOffset - nCellFirstGen * 2^(iGen-1) + ceil(iGenCell / 2);
                        iFamily = ancestry(iParent, MultiGenerations.ANCESTRY_COLIDX_FAMILY);
                    end

                    if iGen < nGen-1
                        iChild1 = genOffset + nCellFirstGen * 2^iGen + 2 * iGenCell - 1;
                        iChild2 = genOffset + nCellFirstGen * 2^iGen + 2 * iGenCell;
                    else
                        iChild1 = NaN;
                        iChild2 = NaN;
                    end

                    ancestry(iCell, MultiGenerations.ANCESTRY_COLIDX_CELL) = genOffset + iGenCell;
                    ancestry(iCell, MultiGenerations.ANCESTRY_COLIDX_GEN) = iGen;
                    ancestry(iCell, MultiGenerations.ANCESTRY_COLIDX_FAMILY) = iFamily;
                    ancestry(iCell, MultiGenerations.ANCESTRY_COLIDX_PARENT) = iParent;
                    ancestry(iCell, MultiGenerations.ANCESTRY_COLIDX_CHILD1) = iChild1;
                    ancestry(iCell, MultiGenerations.ANCESTRY_COLIDX_CHILD2) = iChild2;
                end
                genOffset = genOffset + nCellFirstGen * 2^iGen;
            end
        end

        function relations = calcAncestryRelations(nGen, nCellFirstGen)
            import edu.stanford.covert.cell.sim.analysis.MultiGenerations;

            nCells = 2^nGen - 1;

            %intergenerational -- difference in generations
            generations = NaN(nCells, 1);
            genOffset = 0;
            for iGen = 0:nGen-1
                iCell = genOffset + (1:2^iGen);
                genOffset = genOffset + 2^iGen;
                generations(iCell) = iGen;
            end
            deltaGeneration = -generations(:, ones(nCells, 1)) + generations(:, ones(nCells, 1))';

            %intragenerational -- number of generations since common ancestor
            ancestry = MultiGenerations.calcAncestry(nGen, 1);
            removal = NaN(nCells, nCells);
            for iCell1 = 1:nCells
                removal(iCell1, iCell1) = 0;
                for iCell2 = iCell1+1:nCells
                    [~, removal(iCell1, iCell2)] = MultiGenerations.calcCommonAncestor(iCell1, iCell2, ancestry);
                    removal(iCell2, iCell1) = removal(iCell1, iCell2);
                end
            end

            tmp = cat(3, deltaGeneration, removal);
            relations = NaN(nCellFirstGen * nCells, nCellFirstGen * nCells, 2); %cell-1 x cell-2 x [generation separation,  removal]
            ancestry = MultiGenerations.calcAncestry(nGen, nCellFirstGen);
            for i = 1:nCellFirstGen
                relations(ancestry(:, MultiGenerations.ANCESTRY_COLIDX_FAMILY) == i, ancestry(:, MultiGenerations.ANCESTRY_COLIDX_FAMILY) == i, :) = tmp;
            end
        end

        function [iCommonAncestor, removal] = calcCommonAncestor(iCell1, iCell2, ancestry)
            import edu.stanford.covert.cell.sim.analysis.MultiGenerations;

            if iCell1 == iCell2
                iCommonAncestor = iCell1;
                removal = 0;
            elseif ancestry(iCell1, MultiGenerations.ANCESTRY_COLIDX_GEN) > ancestry(iCell2, MultiGenerations.ANCESTRY_COLIDX_GEN)
                [iCommonAncestor, removal] = MultiGenerations.calcCommonAncestor(ancestry(iCell1, MultiGenerations.ANCESTRY_COLIDX_PARENT), iCell2, ancestry);
            elseif ancestry(iCell1, MultiGenerations.ANCESTRY_COLIDX_GEN) < ancestry(iCell2, MultiGenerations.ANCESTRY_COLIDX_GEN)
                [iCommonAncestor, removal] = MultiGenerations.calcCommonAncestor(iCell1, ancestry(iCell2, MultiGenerations.ANCESTRY_COLIDX_PARENT), ancestry);
            else
                [iCommonAncestor, removal] = MultiGenerations.calcCommonAncestor(...
                    ancestry(iCell1, MultiGenerations.ANCESTRY_COLIDX_PARENT), ...
                    ancestry(iCell2, MultiGenerations.ANCESTRY_COLIDX_PARENT), ...
                    ancestry);
                removal = removal + 1;
            end
        end

        function colors = calcRedGreenColors(n)
            tmp = (0:n-1)';
            r1 = 1;
            r2 = 0;
            g1 = 0;
            g2 = 1;
            b1 = 0;
            b2 = 0;
            colors = [
                r1*2*max(0, 0.5-tmp/(n-1)) + r2*2*max(0, tmp/(n-1)-0.5) ...
                g1*2*max(0, 0.5-tmp/(n-1)) + g2*2*max(0, tmp/(n-1)-0.5) ...
                b1*2*max(0, 0.5-tmp/(n-1)) + b2*2*max(0, tmp/(n-1)-0.5) ...
                ];
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function max_index = r8vec_max_abs_index ( n, a )

%*****************************************************************************80
%
%% R8VEC_MAX_ABS_INDEX: index of the maximum absolute value in an R8VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the array.
%
%    Input, real A(N), the array.
%
%    Output, integer MAX_INDEX, the index of the entry of
%    maximum absolute value.
%
  if ( n <= 0 )

    max_index = -1;

  else

    max_index = 1;

    for i = 2 : n
      if ( abs ( a(max_index) ) < abs ( a(i) ) )
        max_index = i;
      end
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%ControlCenter

% compile the mex programming function.
mex SearchCycleNode.c

n=30; % The number of nodes in the tree

CostMatric = rand( n );  % randomly generating a cost matrix
for p=1:n
    CostMatric( p,p ) = 0;
end

% 1. search maximum directed spanning tree by unsymmetrical cost matrix "CostMatric", with specified root node "Root".

Root = 3;    % Root of the tree is predefined ahead .
[MaxTree1,MaxCost1] =  DirectedMaximumSpanningTree( CostMatric,Root )
h = view(biograph( MaxTree1 ));

% As for other provided algorithms, you can look at the code , and change
% the part in the similar way.....cheers!!!

author: Sable

repo: mcbench-benchmarks


text:
%kvnthick 'N-1 Compositions of Thickenings'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vnthick.pane file
%
% Parameters:
% InputFile: i 'Input Image ', required: 'input image'
% InputFile: str1 'Struct Element 1', required: 'Structuring Element'
% InputFile: str2 'Struct Element 2', required: 'Structuring Element'
% Integer: n 'Number of Thickenings', default: 5: 'number of thickenings'
% MultiChoice: r 'Struct El. Rotation Step', default: 1: 'Rotation step'
%    Choices are:
%   1: '1 '
%   2: '2 '
%   3: '4 '
% OutputFile: o 'Output Image ', required: 'resulting output image'
%
% Example: o = kvnthick({i, str1, str2}, {'i','';'str1','';'str2','';'n',5;'r',1;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vnthick - N-1 compositions of Thickenings
%
%  DESCRIPTION
% .I vnthick
% will do N successive thickenings of an image.
%
% The output image data type will be the same as that of the highest data type
% of the input image.
%
% Some useful pairs of structuring elements:
%
%                               0 0 0    1 1 1
%   Conditional Segmentation    0 1 0    0 0 0
%                               1 1 1    0 0 0
%
%                               0 0 0    1 0 0
%   Pseudo Convex Hull          0 1 1    1 0 0
%                               0 0 0    1 0 0
%
%                               1 0 0    0 0 0
%   Convex Hull                 1 0 0    0 1 0
%                               1 0 0    0 0 0
%
%                               0 0 0    1 1 1
%   Homotopic Pruning           0 1 0    1 0 1
%                               0 0 0    1 1 1
%
%
%
%  EXAMPLES
% vnthick -i feath.xv -str1 b1.str -str2 b2.str -n 5 -r 1 -o outimage.xv
%
% Will do 5 successive thickenings (vnthick) of the input image "feath.xv" by the structuring element b1.str and b2.str, rotating the struct. els. 1 step after each thinning. The resulting image will be written to  "outimage.xv".
%
%  "SEE ALSO"
% vcast1(1)
%
%  RESTRICTIONS
% .I vnthick
% can be defined for all data types supported by Khoros, but at the moment it has been implemented just for the bit and unsigned char types.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993-1997 Junior Barrera, Roberto Lotufo.  All rights reserved.
%

function varargout = kvnthick(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kvnthick(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'str1', '__input';'str2', '__input';'n', 5;'r', 1;'o', '__output'};
maxval={0,0,0,2,0,0};
minval={0,0,0,2,0,0};
istoggle=[0,0,0,0,0,0];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','InputFile','Integer','MultiChoice','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'vnthick"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function pathS = MC_M(S0,r,d,T,sigma,muj,sigj,lambda,...
    NTime,NSim,NBatches)
% discretization of Merton model
pathS = zeros(NSim,NTime+1,NBatches);

dT = T/NTime;                        % time step

lnS1 = zeros(NSim,NTime+1);          % used in batch
lnS1(:,1) = log(S0*exp(-d*T));       % set S(0) dividend adjusted

omega = (r-d-0.5 * sigma^2) * dT;    % martingale adjustment

for l = 1:NBatches                   % batch loop
    dW = randn(NSim,NTime);          % Gaussians
    Z = randn(NSim,NTime);           % Gaussians
    for i=2:NTime+1                  % time loop
        % Jumps
        P = poissrnd(lambda*dT,NSim,1);        % Poissonians
        lnY = muj*P + sigj*sqrt(P).* Z(:,i-1); % log jumps

        lnS1(:,i) = lnS1(:,i-1) + omega + sigma* sqrt(dT)*dW(:,i-1)...
            + lnY - dT*lambda*(exp(muj+0.5*sigj^2)-1);
    end
    pathS(:,:,l) = exp(lnS1);
end


text:
%  Figure 3.30a      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%

clf;
num=1;
a=10;

zeta =1;
den2=[1/(zeta*a) 1];
den1=[1 2*zeta 1];
den=conv(den1,den2);
t=0:.1:8;
y1=step(num,den,t);

zeta =.7;
den2=[1/(zeta*a) 1];
den1=[1 2*zeta 1];
den=conv(den1,den2);
y2=step(num,den,t);

zeta =.5;
den2=[1/(zeta*a) 1];
den1=[1 2*zeta 1];
den=conv(den1,den2);
y3=step(num,den,t);

axis([0 5 .1 .9])
plot(t,y1,'-',t,y2,'--',t,y3,'-.'),grid
title('Fig. 3.30a Step response with extra pole, \alpha= 10')
xlabel('\omega_n t')
ylabel('y(t)')


text:
function B = imageblocks(img, sz, st)

% imageblocks -- Extract blocks of specified size from image or
%                array of images
%
% Usage:
%       B = imageblocks(img, sz, st)
%
% Input:
%       img         Image from which to extract blocks
%       sz          Block size vector
%       st          Block step vector
%
% Output:
%       B           Array of image blocks
%
%
% Author: Brendt Wohlberg <brendt@lanl.gov>  Modified: 2014-10-20
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'License' file distributed with
% the library.

if nargin < 3,
  st = [1 1];
end

Nbr = sz(1);
Nbc = sz(2);
[Nir Nic Nih] = size(img);
Nib = length(1:st(1):(Nir-Nbr+1))*length(1:st(2):(Nic-Nbc+1))*Nih;
B = zeros(Nbr, Nbc, Nib);

% Loop over all blocks in images, adding each one to block set
n = 1;
for k = 1:Nih,
  for l=1:st(1):(Nir-Nbr+1),
    for m=1:st(2):(Nic-Nbc+1),
      B(:,:,n) = img(l:(l+Nbr-1), m:(m+Nbc-1), k);
      n = n + 1;
    end
  end
end

return

author: thfylsty


text:
function [H] = rotate(D);

% ROTATE returns the homogenous coordinate transformation matrix
% corresponding to a rotation around the x, y and z-axis. The direction of
% the rotation is according to the right-hand rule.
%
% Use as
%   [H] = rotate(R)
% where
%   R		[rx, ry, rz] in degrees
%   H		corresponding homogenous transformation matrix
%
% Note that the order in which the rotations are performs matters. The
% rotation is first done around the z-axis, then the y-axis and finally the
% x-axis.

% Copyright (C) 2000-2005, Robert Oostenveld
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% $Log: rotate.m,v $
% Revision 1.1  2009/01/30 04:02:11  arno
% *** empty log message ***
%
% Revision 1.5  2006/09/12 13:35:28  roboos
% convert input rotation from degrees (according to documentation) into radians (needed for computations)
%
% Revision 1.4  2005/08/15 08:15:33  roboos
% reimplemented the rotate function, which contained an error (the error is in the AIR technical reference)
% changed all functions to be dependent on the rotate, translate and scale function
% all functions now behave consistenly, which also means that they are not compleetly backward compatible w.r.t. the order of the rotations
%
% Revision 1.3  2004/05/19 09:57:07  roberto
% added GPL copyright statement, added CVS log item
%

% convert degrees to radians
R = D*pi/180;

% get the individual angles (in radians)
rx = R(1);
ry = R(2);
rz = R(3);

% precompute the sin/cos values of the angles
cX = cos(rx);
cY = cos(ry);
cZ = cos(rz);
sX = sin(rx);
sY = sin(ry);
sZ = sin(rz);

% according to Roger Woods' http://bishopw.loni.ucla.edu/AIR5/homogenous.html
% it should be this, but I cannot reproduce his rotation matrix
% H = eye(4,4);
% H(1,1) = cZ*cY + sZ*sX*sY;
% H(1,2) = sZ*cY - cZ*sX*sY;
% H(1,3) =            cX*sY;
% H(2,1) = -sZ*cX;
% H(2,2) =  cZ*cX;
% H(2,3) =     sX;
% H(3,1) =  sZ*sX*cY - cZ*sY;
% H(3,2) = -cZ*sX*cY - sZ*sY;
% H(3,3) =             cX*cY;

% instead, the following rotation matrix does work according my
% expectations. It rotates according to the right hand rule and first
% rotates around z, then y and then x axis
H = [
           cZ*cY,          -sZ*cY,              sY,               0
  cZ*sY*sX+sZ*cX, -sZ*sY*sX+cZ*cX,          -cY*sX,               0
 -cZ*sY*cX+sZ*sX,  sZ*sY*cX+cZ*sX,           cY*cX,               0
               0,               0,               0,               1
];

if 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following code can be used to construct the combined rotation matrix
% for either xyz or zyx ordering (using the Matlab symbolic math toolbox)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  syms sX sY sZ cX cY cZ
  % this is for only rotating around x
  Rx = [
    1   0    0   0
    0   cX  -sX  0
    0   sX   cX  0
    0   0    0   1
    ];
  % this is for only rotating around y
  Ry = [
    cY  0   sY  0
    0   1   0   0
    -sY  0   cY  0
    0   0   0   1
    ];
  % this is for only rotating around z
  Rz = [
    cZ -sZ  0   0
    sZ  cZ  0   0
    0   0   1   0
    0   0   0   1
    ];
  % combine them
  Rzyx = Rz * Ry * Rx  % rotate around x, y, then z
  Rxyz = Rx * Ry * Rz  % rotate around z, y, then x
end

author: PatternRecognition

repo: OpenBMI


text:
function [k] = createStimCutMovieMontage(inputMovie,nAlignPts,timeVector,varargin)
	% Creates a montage movie aligned to specific timepoints.
	% Biafra Ahanonu
	% fxn started: 2014.08.13 - broke off from controllerAnalysis script from ~2014.03
	% inputs
		% inputMovie - path to movie file, in cell array, e.g. {'path.h5'}
		% inputAlignPts - vector containing the frames to align to
		% savePathName - path to save output movie, exclude the extension.
	% outputs
		% 2015.11.05

	[k] = ciapkg.video.createStimCutMovieMontage(inputMovie,nAlignPts,timeVector,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function test_suite = test_regression1

%   Run specific demo and save values for comparison.
%
%   See also
%     TEST_ALL, DEMO_REGRESSION1

initTestSuite;

  function testDemo
    % Set random number stream so that test failing isn't because randomness.
    % Run demo & save test values.
    prevstream=setrandstream(0);

    disp('Running: demo_regression1')
    demo_regression1
    Eft_map = Eft_map(1:50);
    Varft_map = Varft_map(1:50);
    Eft_ia = Eft_ia(1:50);
    Varft_ia = Varft_ia(1:50);
    Eft_mc = Eft_mc(1:50);
    Varft_mc = Varft_mc(1:50);
    path = which('test_regression1.m');
    path = strrep(path,'test_regression1.m', 'testValues');
    if ~(exist(path, 'dir') == 7)
      mkdir(path)
    end
    path = strcat(path, '/testRegression1');
    save(path, 'K', 'C', 'w', 'Eft_map', 'Varft_map', ...
      'Eft_ia', 'Varft_ia', 'Eft_mc', 'Varft_mc');

    % Set back initial random stream
    setrandstream(prevstream);
    drawnow;clear;close all

% Test saved values with multiple tests. Covariance matrices and
% optimized parameters are tested with zero tolerance, while mean and
% variances from various approximations are tested with relative tolerance
% (5% for MCMC and 1% for grid and IA)

    function testCovarianceMatrices
        values.real = load('realValuesRegression1.mat','K','C');
        values.test = load(strrep(which('test_regression1.m'), 'test_regression1.m', 'testValues/testRegression1.mat'),'K','C');
        assertElementsAlmostEqual(values.real.K, values.test.K);
        assertElementsAlmostEqual(values.real.C, values.test.C);

    function testOptimizedParameter
        values.real = load('realValuesRegression1.mat','w');
        values.test = load(strrep(which('test_regression1.m'), 'test_regression1.m', 'testValues/testRegression1.mat'),'w');
        assertElementsAlmostEqual(values.real.w, values.test.w, 'relative', 0.1);

    function testPredictedMeanVarianceGrid
        values.real = load('realValuesRegression1.mat','Eft_map','Varft_map');
        values.test = load(strrep(which('test_regression1.m'), 'test_regression1.m', 'testValues/testRegression1.mat'),'Eft_map','Varft_map');
        if length(values.test.Eft_map) > 50
            assertElementsAlmostEqual(values.real.Eft_map(1:50), values.test.Eft_map(1:50), 'relative', 0.1);
            assertElementsAlmostEqual(values.test.Varft_map(1:50), values.real.Varft_map(1:50), 'relative', 0.1);
        else
            assertElementsAlmostEqual(values.test.Eft_map, values.real.Eft_map, 'relative', 0.1);
            assertElementsAlmostEqual(values.test.Varft_map, values.real.Varft_map, 'relative', 0.1);
        end

    function testPredictedMeanVarianceMC
        values.real = load('realValuesRegression1.mat','Eft_mc','Varft_mc');
        values.test = load(strrep(which('test_regression1.m'), 'test_regression1.m', 'testValues/testRegression1.mat'),'Eft_mc','Varft_mc');
        if length(values.test.Eft_mc) > 50
            assertElementsAlmostEqual(mean(mean(values.test.Eft_mc(1:50))), mean(mean(values.real.Eft_mc(1:50))), 'relative', 0.1);
            assertElementsAlmostEqual(mean(values.test.Varft_mc(1:50)), mean(values.real.Varft_mc(1:50)), 'absolute', 0.3);
       else
            assertElementsAlmostEqual(mean(mean(values.test.Eft_mc)), mean(mean(values.real.Eft_mc)), 'relative', 0.1);
            assertElementsAlmostEqual(mean(mean(values.test.Varft_mc)), mean(mean(values.real.Varft_mc)), 'absolute', 0.3);
        end

    function testPredictedMeanVarianceIA
        values.real = load('realValuesRegression1.mat','Eft_ia','Varft_ia');
        values.test = load(strrep(which('test_regression1.m'), 'test_regression1.m', 'testValues/testRegression1.mat'),'Eft_ia','Varft_ia');
        if length(values.test.Eft_ia) > 50
            assertElementsAlmostEqual(mean(values.test.Eft_ia(1:50)), mean(values.real.Eft_ia(1:50)), 'relative', 0.1);
            assertElementsAlmostEqual(mean(values.test.Varft_ia(1:50)), mean(values.real.Varft_ia(1:50)), 'relative', 0.1);
        else
            assertElementsAlmostEqual(mean(values.test.Eft_ia), mean(values.real.Eft_ia), 'relative', 0.1);
            assertElementsAlmostEqual(mean(values.test.Varft_ia), mean(values.real.Varft_ia), 'relative', 0.1);
        end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function obj3 = binaryopp(f, obj1, obj2)

    if isa(obj2, 'mmo'), tmpobj = obj2; obj2 = obj1; obj1 = tmpobj; clear tmpobj; end
    if ~isequal(size(obj1), size(obj2)) && prod(size(obj2)) ~= 1
        error('Matrix dimensions must agree.');
    end
    data1 = memmapfile(obj1.dataFile, 'writable', obj1.writable, 'format', { 'single' obj1.dimensions 'x' });
    if isa(obj2, 'mmo'),
         data2 = memmapfile(obj2.dataFile, 'writable', obj2.writable, 'format', { 'single' obj2.dimensions 'x' });
    else data2 = obj2;
    end

    % make new memory mapped data file (blank)
    % --------------------------------
    newFileName = mmo.getnewfilename;
    fid = fopen(newFileName, 'w');
    s1  = size(obj1);
    ss.type = '()';
    ss.subs(1:length(s1)-1) = { ':' };
    for index = 1:s1(end)
        ss.subs{length(s1)} = index;
        if prod(size(data2)) == 1
             tmpdata = f(subsref(data1.Data.x, ss), data2);
        else tmpdata = f(subsref(data1.Data.x, ss), subsref(data2.Data.x, ss));
        end
        fwrite(fid, tmpdata, 'float');
    end
    fclose(fid);

    % create object
    % -------------
    obj3 = mmo(newFileName, s1, true, obj1.transposed);
    obj3 = updateWorkspace(obj3);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function test_bug2820

% WALLTIME 00:10:00
% MEM 3gb
% DEPENDENCY

% it should be is2Dana && ~is2Dfun, so start with a surface
[ftver, ftpath] = ft_version;
surffile = fullfile(ftpath, 'template', 'anatomy', 'surface_pial_both.mat');
anatomical = ft_read_headshape(surffile);

functional     = [];
functional.dim = [256 256 256];
functional.transform = eye(4);
functional.transform(1,4) = -128.5; % place the center of the cube in the origin
functional.transform(2,4) = -128.5;
functional.transform(3,4) = -128.5;
functional.pow = randn(functional.dim);

cfg = [];
cfg.funparameter = 'pow';
ft_sourceplot(cfg, functional);

% see whether they fit
ft_determine_coordsys(functional, 'interactive', 'no');
ft_plot_mesh(anatomical, 'edgecolor', 'none', 'facecolor', 'skin');
camlight

cfg = [];
cfg.parameter = 'pow';
% cfg.interpmethod  default does not cause  aproblem
interp1 = ft_sourceinterpolate(cfg, functional, anatomical);

cfg = [];
cfg.parameter = 'pow';
cfg.interpmethod = 'project'; % needed to get to the bug
interp2 = ft_sourceinterpolate(cfg, functional, anatomical);

%% the figures should be specled green (on MATLAB 2014b and later)
figure; ft_plot_mesh(interp1, 'vertexcolor', interp1.pow, 'edgecolor', 'none'); camlight; colorbar
figure; ft_plot_mesh(interp2, 'vertexcolor', interp2.pow, 'edgecolor', 'none'); camlight; colorbar

%% use the high-level functionality to get the same figure, this will call ft_sourceinterpolate
cfg = [];
cfg.surffile = surffile;
cfg.funparameter = 'pow';
cfg.method = 'surface';
ft_sourceplot(cfg, functional);

%% another approach to use the high-level functionality to get the same figure, again calling ft_sourceinterpolate internally
cfg = [];
cfg.funparameter = 'pow';
cfg.method = 'surface';
ft_sourceplot(cfg, functional, anatomical);


text:
function [yPred,PInvPred]=discEIFPred(yPrev,PInvPrev,f,FJacob,Q,FHessian)
%%DISCEIFPRED Perform the discrete-time prediction step that comes with
%             the first and second-order extended information filter (EIF).
%             As the extended information filter as presented in [1] has no
%             real propagation step, this just extracts the state, calls
%             discEKFPred, and converts the result back to an information
%             state. This function is useful if one is using an information
%             state because of how the measurement update in an information
%             filter is performed.
%
%INPUTS:    yPrev   The xDimX1 information state at the previous time-step.
%                   The information state is the inverse covariance matrix
%                   times the target state.
%        PInvPrev   The xDimXxDim inverse of the state covariance matrix at
%                   the previous time-step.
%               f   A function handle for the state transition function
%                   that takes the state as its parameter.
%           FJacob  A function handle for calculating the xDim X xDim
%                   Jacobian of f, or the xDim X xDim Jacobian matrix
%                   itself. If an empty matrix is passed, then
%                   FJacob will be found using numerical differentiation
%                   via the numDiff function with default parameters.
%               Q   The xDimX xDim process noise covariance matrix.
%         FHessian  This parameter is only provided if a second-order EKF
%                   is desired. This is either a function handle for the
%                   state transition Hessian hypermatrix, or it is the
%                   state transition Hessian hypermatrix itself. The matrix
%                   is xDim X xDim X xDim. The matrix FH=FHessian(x) is
%                   such that FH(i,j,k) is the second derivative of the
%                   kth element of the vector returned by f with
%                   respect to the ith and jth components of x. The Hessian
%                   matrix is symmetric. If this parameter is omitted, a
%                   first-order filter is used.
%
%OUTPUTS:   yPred    The xDim X 1 predicted information state vector.
%           PInvPred The predicted xDim X xDim inverse state covariance
%                    matrix.
%
%REFERENCES:
%[1] K. P. B. Chandra, D.-W. Gu, and I. Postlethwaite, "Square root
%    cubature information filter," IEEE Sensors Journal, vol. 13, no. 2,
%    pp. 750-758, Feb. 2013.
%
%October 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<6)
    FHessian=[];
end

PPrev=pinv(PInvPrev);
xPrev=PPrev*yPrev;

[xPred, PPred]=discEKFPred(xPrev,PPrev,f,FJacob,Q,FHessian);

PInvPred=pinv(PPred);
yPred=PInvPred*xPred;

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [sys,x0,str,ts]=NL_PID_3fal(t,x,u,flag,r0,h0,r1,h1,h,Bet,A,D)
switch flag,
    case 0,
        [sys,x0,str,ts]=mdlInitializeSizes(h);
    case 2,
        sys=mdlUpdates(x,u,r0,r1,h0,h1,h);
    case 3,
        sys=mdlOutputs(x,Bet,A,D);
    case 4,
        sys=mdlGetTimeOfNextVarHit(t,h);
    case {1,9},
        sys=[];
    otherwise
        error(['Unhandled flag=',num2str(flag)]);
end
function [sys,x0,str,ts]=mdlInitializeSizes(h)
    sizes=simsizes;
    sizes.NumContStates=0;
    sizes.NumDiscStates=5;
    sizes.NumOutputs=2;
    sizes.NumInputs=2;
    sizes.DirFeedthrough=1;
    sizes.NumSampleTimes=1;
    sys=simsizes(sizes);
    x0=[0;0;0;0;0];
    str=[];
    ts=[h 0];
function sys=mdlUpdates(x,u,r0,r1,h0,h1,h)
    fh0=fhan(x(1)-u(1),x(2),r0,h0);
    x(1)=x(1)+h*x(2);
    x(2)=x(2)+h*fh0;
    fh1=fhan(x(3)-u(2),x(4),r1,h1);
    x(3)=x(3)+h*x(4);
    x(4)=x(4)+h*fh1;
    x(5)=x(5)+h*(x(1)-x(3));
    sys=x;
function sys=mdlOutputs(x,Bet,A,D)
    sys(1)=Bet(1)*fal(x(5),A(1),D)+Bet(2)*fal(x(1)-x(3),A(2),D)+Bet(3)*fal(x(2)-x(4),A(3),D);
    sys(2)=x(1);
function sys=mdlGetTimeOfNextVarHit(t,h)
    sys=t+h;
function y=fhan(x1,x2,r,h)
d=r*h;
d0=h*d;
y=x1+h*x2;
a0=sqrt(d^2+8*r*abs(y));
if abs(y)>d0
    a=x2+(a0-d)*sign(y)/2;
else
    a=x2+y/h;
end

if abs(a)>d
    y=-r*sign(a);
else
    y=-r*a/d;
end
function f=fal(e,a,d)
    if abs(e)<d
        f=e*d^(a-1);
    else f=(abs(e))^a*sign(e);
    end

author: TianfaYao

repo: ADRC

save_path: github-repos/MATLAB/TianfaYao-ADRC


text:
function [data_endo,favar]=favar_gensample3(data_endo,favar)

% determine the numbers of variables other than factors
favar.numdata_exfactors=size(favar.data_exfactors,2);

if favar.onestep==1
    % initalise variables
    XZ=favar.XZ;
    nfactorvar=favar.nfactorvar;
    numpc=favar.numpc;
    Lf=favar.l;
    % identify factors in the case of onestep estimation: rotate the
    % factors and the loadings matrix following Bernanke, Boivin, Eliasz (2005)
    Lfy=bear.favar_olssvd(favar.X(:,numpc+1:nfactorvar),data_endo)';% upper KxM block of Ly set to zero
    Lf=[Lf(1:numpc,:);Lfy(:,1:numpc)];
    Ly=[zeros(numpc,favar.numdata_exfactors);Lfy(:,numpc+1:numpc+favar.numdata_exfactors)];

    %transform factors and loadings for LE normalization
    [ql,rl]=qr(Lf');
    Lf=rl;  % do not transpose yet, is upper triangular
    XZ=XZ*ql;
    %need identity in the first K columns of Lf, call them A for now
    A=Lf(:,1:numpc); %index here for factors
    Lf=[eye(numpc),inv(A)*Lf(:,numpc+1:nfactorvar)]';
    favar.XZ_rotated=XZ*A;

    %rotated loadings
    favar.L=[Lf Ly;zeros(favar.numdata_exfactors,numpc),eye(favar.numdata_exfactors)];

    %replace factors with factors rotated in data_endo
    for ii=1:size(favar.variablestrings_factorsonly)
        data_endo(:,favar.variablestrings_factorsonly(ii))=favar.XZ_rotated(:,ii);
    end

    %errors of factor equation (observation equation)
    if favar.numdata_exfactors~=0
        favar.evf=favar.X-favar.data_exfactors*Ly'-favar.XZ_rotated*Lf';
    elseif favar.numdata_exfactors==0 % in this case we have a pure factor model
        favar.evf=favar.X-favar.XZ_rotated*Lf';
    end
    %favar.evf=favar.XY-data_endo*favar.L';
    favar.Sigma=favar.evf'*favar.evf/size(favar.X,1);
    favar.Sigma=diag(diag(favar.Sigma));
    favar.Sigma=diag([diag(favar.Sigma);zeros(favar.numdata_exfactors,1)]);

    % from here proceed with non-standardised data
    favar.data_exfactors=favar.data_exfactors_temp;
    favar.X=favar.X_temp;

    % replace standardised data with not-standardised data in endo
    for ii=1:size(favar.variablestrings_exfactors)
        data_endo(:,favar.variablestrings_exfactors(ii))=favar.data_exfactors(:,ii);
    end

    % state-space representation
    favar.XY=[favar.X,favar.data_exfactors];

elseif favar.onestep==0 % two-step
    if favar.numdata_exfactors==0
        favar.slowfast=0; % this identifiaction is not applicable in a pure factor model
    end
    if favar.slowfast==1 %apply slowfast recursive identification as in BBE (2005)
        % factor roation with slow/fast scheme
        favar.XZ_rotated=bear.favar_facrot(favar.XZ,favar.data_exfactors(:,end),favar.XZ_slow); %end, has eventually to be changed
        %replace factors with factors rotated
        for ii=1:size(favar.variablestrings_factorsonly)
            data_endo(:,favar.variablestrings_factorsonly(ii))=favar.XZ_rotated(:,ii);
        end
    end
    % state-space representation
    favar.XY=[favar.X,favar.data_exfactors];
    % new loadings
    favar.L=(bear.favar_olssvd(favar.XY,data_endo))';
    % manipulate loadings matrix for blocks
    if favar.blocks==1
        for bb=1:favar.nbnames
            Xbindex=favar.blockindex_each{bb,1};
            Xbindex2=favar.blockindex_each{bb,1}==0;
            %
            favar.L(Xbindex,favar.blocks_index{bb,1})=favar.l_block{bb,1};
            favar.L(Xbindex2,favar.blocks_index{bb,1})=0;
            favar.L(1:favar.nfactorvar,favar.variablestrings_exfactors)=0;
            favar.L(favar.nfactorvar+1:end,favar.variablestrings_factorsonly)=0;
        end
        for vv=1:favar.numdata_exfactors
            favar.L(favar.nfactorvar+1:end,favar.variablestrings_exfactors(vv,1))=1;
        end
    end

    %errors of factor equation (observation equation)
    favar.evf=favar.XY-data_endo*favar.L';
    favar.Sigma=favar.evf'*favar.evf/size(favar.XY,1); %
    favar.Sigma=diag(diag(favar.Sigma));
end

% to activate routines in the BVAR framwork and IRFt 4, where we have It-Bu
% x X, Y and L
favar.bvar=1;
%%% this should be equivalent
%regressing the factors from all x on the slow moving factors and the FFR
% % %         favar.Beta = mvregress(data_endo,favar.XZ);
% % %         favar.Beta_exfactors = favar.Beta(favar.variablestrings_exfactors(:,end),:); %we assume here that favar.variablestrings_exfactors, all the variables that are not factors, are
% % %         favar.Frot=favar.XZ-favar.data_endo(:,favar.variablestrings_exfactors(:,end))*favar.Beta_exfactors;

author: european-central-bank

repo: BEAR-toolbox


text:
function [teE, teW] = nfgCompareAllVolumeTrueError(fgG,g_bundleID,g_radius,fibersSel,iBundleSel,phantomDir)
%Compare volume of all fiber matches using trueError program.
%
%   [teE, teW] =
%   nfgCompareAllVolumeTrueError(fgG,g_bundleID,g_radius,fibersSel,iBundleS
%   el,phantomDir)
%
%
%
% NOTES:

% if ieNotDefined('projName'); projName = 'projectome'; end

% XXX THE BELOW TEST MUST BECOME GENERIC
% First lets find the optimum radius for the trueError test using a
% particular bundle
disp('Searching for near optimal test diameter ...');
selID = 4;
sMask = zeros(size(iBundleSel));
gMask = zeros(size(g_bundleID));
for ss=selID
    sMask = sMask | iBundleSel==ss;
    gMask = gMask | g_bundleID==ss;
end
vRadiusTE = 0.01:0.005:0.03;
fgSS = dtiNewFiberGroup();
fgGS = dtiNewFiberGroup();
fgSS.fibers = fibersSel(sMask);
fgGS.fibers = fgG.fibers(gMask);
[imgE] = nfgCompareVolumeTrueError(fgGS, fgSS, g_radius(1), vRadiusTE, phantomDir);
[mIE, iIE] = min(imgE);
te_radius = vRadiusTE(iIE);
% te_radius = 0.014;
disp(['Found optimal diameter = ' num2str(te_radius*2) 'mm.']);

% Calculate arclengths of all bundles found in projectome and in gold
teE = zeros(1,max(g_bundleID));
teW = zeros(1,max(g_bundleID));
fgSS = dtiNewFiberGroup();
fgGS = dtiNewFiberGroup();
for bb=1:max(g_bundleID)
    disp(['Computing volume for bundle ' num2str(bb) ' ...']);
    fgSS.fibers = fibersSel(iBundleSel==bb);
    fgGS.fibers = fgG.fibers(g_bundleID==bb);
    if isempty(fgSS.fibers)
        if ~isempty(fgGS.fibers)
            teE(bb) = 100;
        end
    else
        [teE(bb), teW(bb)] = nfgCompareVolumeTrueError(fgGS, fgSS, g_radius(bb), te_radius, phantomDir);
    end

end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function varargout = drawGrid3d(varargin)
%DRAWGRID3D Draw a 3D grid on the current axis
%
%   drawGrid3d
%   draws a 3D square grid, with origin (0,0,0) and spacing 1 in each
%   direction, with bounds corresponding to the bounds of current axis.
%
%   drawGrid3d(SPACING)
%   where spacing is either a scalar or a [1x3] matrix, specifies the size
%   of the unit cell.
%
%   drawGrid3d(ORIGIN, SPACING)
%   Also specify origin of grid. ORIGIN is a [1x3] array.
%
%   drawGrid3d(..., EDGE)
%   specifies whether function should draw edges touching edges of axis.
%   EDGE is a characheter string, which can be :
%   - 'OPEN' : each line start from one face of window to the opposite
%   face. This results in a 'spiky' grid.
%   - 'CLOSED' (default value) : each line stops at the last visible point
%   of the grid for this line. The result looks like a box (no free spikes
%   around the grid).
%
%   H = drawGrid3d(...);
%   return a vector of handles for each LINE object which was crated.
%

%   ------
%   Author: David Legland
%   e-mail: david.legland@grignon.inra.fr
%   Created: 2005-11-17
%   Copyright 2005 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

%% initialize variables -----

% default values
closed = true;
origin = [0 0 0];
spacing = [1 1 1];

% check if grid is open or not
str = '';
if ~isempty(varargin)
    str = varargin{end};
end
if ischar(str)
    if strncmpi(str, 'open', 4)
        closed = false;
    end
    varargin = varargin(1:end-1);
end

% check origin and grid spacing
if length(varargin)==1
    spacing = varargin{1};
elseif length(varargin)==2
    origin = varargin{1};
    spacing = varargin{2};
end

%% Compute internam data -----

% get axis limits
ax = axis;
x0 = ax(1); x1 = ax(2);
y0 = ax(3); y1 = ax(4);
z0 = ax(5); z1 = ax(6);

% get first and last coordinates of the grid in each direction
dx = spacing(1); dy = spacing(2); dz = spacing(3);
xe = x0 + mod(origin(1) - x0, dx);
xf = x1 - mod(x1 - origin(1), dx);
ye = y0 + mod(origin(2) - y0, dy);
yf = y1 - mod(y1 - origin(2), dy);
ze = z0 + mod(origin(1) - z0, dz);
zf = z1 - mod(z1 - origin(1), dz);

% update first and last coordinate if grid is 'closed'
if closed
    x0 = xe; x1 = xf;
    y0 = ye; y1 = yf;
    z0 = ze; z1 = zf;
end

%% Draw the grid -----

h = [];
%TODO: rewrite code, avoiding loops

% draw lines parallel to x axis
for y = ye:dy:yf
    for z = ze:dz:zf
        h = [h; drawEdge3d([x0 y z x1 y z])]; %#ok<AGROW>
    end
end

% draw lines parallel to y axis
for x = xe:dx:xf
    for z = ze:dz:zf
        h = [h; drawEdge3d([x y0 z x y1 z])]; %#ok<AGROW>
    end
end

% draw lines parallel to z axis
for x = xe:dx:xf
    for y = ye:dy:yf
        h = [h; drawEdge3d([x y z0 x y z1])]; %#ok<AGROW>
    end
end

%% Check output arguments -----

if nargout>0
    varargout{1} = h;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [wz,AA,RA] = Interaction(PopObj,Point)
% Identify the preferred point and all the others are treated as
% non-preferred points

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

    %% Identify the unique and non-dominated solutions
    PopObj = unique(PopObj,'rows');
    PopObj = PopObj(NDSort(PopObj,1)==1,:);

    %% Identify the preferred solution
    ideal = min(PopObj,[],1);
    % Calculate the Tchebycheff function value of each solution on Point
    Fitness = max((PopObj-repmat(ideal,size(PopObj,1),1))./repmat(Point,size(PopObj,1),1),[],2);
    % The one having the minimal function value is treated as the preferred
    % point
    [~,prefer] = min(Fitness);
    AA = PopObj(prefer,:);
    RA = PopObj([1:prefer-1,prefer+1:end],:);

    %% Calculate the weight distribution function value
    RefPoint = max(PopObj,[],1) + 0.1;
    wz = [0,1,1+CalWHV(AA,RefPoint,ones(1,size(AA,1)))/CalWHV(RA,RefPoint,ones(1,size(RA,1)))];
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
num = length(uLINK);

for n=1:num
    uLINK(n).sister = 0;
    uLINK(n).child  = 0;
end

for n=1:num
    mom = uLINK(n).mother;  % My mother
    if mom ~= 0
        if uLINK(mom).child == 0
            uLINK(mom).child = n;  % I am the first child.
        else
            eldest_sister = uLINK(mom).child;  % I have a sister.
            SetYoungestSister(eldest_sister,n);
        end
    end
end

author: s-kajita

repo: IntroductionToHumanoidRobotics


text:
%[2018]-"Tree growth algorithm (TGA): A novel approach for solving
%optimization problems"

function TGA = jTreeGrowthAlgorithm(feat,label,opts)
% Parameters
lb         = 0;
ub         = 1;
thres      = 0.5;
num_tree1  = 3;    % size of first group
num_tree2  = 5;    % size of second group
num_tree4  = 3;    % size of fourth group
theta      = 0.8;  % tree reduction rate of power
lambda     = 0.5;  % control nearest tree

if isfield(opts,'T'), max_Iter = opts.T; end
if isfield(opts,'N'), N = opts.N; end
if isfield(opts,'N1'), num_tree1 = opts.N1; end
if isfield(opts,'N2'), num_tree2 = opts.N2; end
if isfield(opts,'N4'), num_tree4 = opts.N4; end
if isfield(opts,'theta'), theta = opts.theta; end
if isfield(opts,'lambda'), lambda = opts.lambda; end
if isfield(opts,'thres'), thres = opts.thres; end

% Limit number of N4 to N1
if num_tree4 > num_tree1 + num_tree2
  num_tree4 = num_tree1 + num_tree2;
end
% Objective function
fun = @jFitnessFunction;
% Number of dimensions
dim = size(feat,2);
% Initial
X   = zeros(N,dim);
for i = 1:N
	for d = 1:dim
    X(i,d) = lb + (ub - lb) * rand();
	end
end
% Fitness
fit  = zeros(1,N);
fitG = inf;
for i = 1:N
  fit(i) = fun(feat,label,(X(i,:) > thres),opts);
  % Best
  if fit(i) < fitG
    fitG = fit(i);
    Xgb  = X(i,:);
  end
end
% Sort tree from best to worst
[fit, idx] = sort(fit,'ascend');
X          = X(idx,:);
% Initial
dist = zeros(1,num_tree1 + num_tree2);
X1   = zeros(num_tree1,dim);
Xnew = zeros(num_tree4,dim);
Fnew = zeros(1,num_tree4);

curve = zeros(1,max_Iter);
curve(1) = fitG;
t = 2;
% Iterations
while t <= max_Iter
	% {1} Best trees group
  for i = 1:num_tree1
    r1 = rand();
    for d = 1:dim
      % Local search (1)
      X1(i,d) = (X(i,d) / theta) + r1 * X(i,d);
    end
  	% Boundary
    XB = X1(i,:); XB(XB > ub) = ub; XB(XB < lb) = lb;
    X1(i,:) = XB;
    % Fitness
    fitT = fun(feat,label,(X1(i,:) > thres),opts);
    % Greedy selection
    if fitT <= fit(i)
      X(i,:) = X1(i,:);
      fit(i) = fitT;
    end
  end
  % {2} Competitive for light tree group
  X_ori = X;
  for i = num_tree1 + 1 : num_tree1 + num_tree2
    % Neighbor tree
    for j = 1 : num_tree1 + num_tree2
      if j ~= i
        % Compute Euclidean distance (2)
        dist(j) = sqrt(sum((X_ori(j,:) - X_ori(i,:)) .^ 2));
      else
        % Solve same tree problem
        dist(j) = inf;
      end
    end
    % Find 2 trees with shorter distance
    [~, idx] = sort(dist,'ascend');
    T1       = X_ori(idx(1),:);
    T2       = X_ori(idx(2),:);
    % Alpha in [0,1]
    alpha    = rand();
    for d = 1:dim
      % Compute linear combination between 2 shorter tree (3)
      y = lambda * T1(d) + (1 - lambda) * T2(d);
      % Move tree i between 2 adjacent trees (4)
      X(i,d) = X(i,d) + alpha * y;
    end
    % Boundary
    XB = X(i,:); XB(XB > ub) = ub; XB(XB < lb) = lb;
    X(i,:) = XB;
    % Fitness
    fit(i) = fun(feat,label,(X(i,:) > thres),opts);
  end
  % {3} Remove and replace group
  for i = num_tree1 + num_tree2 + 1 : N
    for d = 1:dim
      % Generate new tree by remove worst tree
      X(i,d) = lb + (ub - lb) * rand();
    end
    % Fitness
    fit(i) = fun(feat,label,(X(i,:) > thres),opts);
  end
  % {4} Reproduction group
  for i = 1:num_tree4
    % Random a best tree
    r     = randi([1,num_tree1]);
    Xbest = X(r,:);
    % Mask operator
    mask  = randi([0,1],1,dim);
    % Mask opration between new & best trees
    for d = 1:dim
      % Generate new solution
      Xn = lb + (ub - lb) * rand();
      if mask(d) == 1
        Xnew(i,d) = Xbest(d);
      elseif mask(d) == 0
        % Generate new tree
        Xnew(i,d) = Xn;
      end
    end
    % Fitness
    Fnew(i) = fun(feat,label,(Xnew(i,:) > thres),opts);
  end
  % Sort population get best nPop trees
  XX        = [X; Xnew];
  FF        = [fit, Fnew];
  [FF, idx] = sort(FF,'ascend');
  X         = XX(idx(1:N),:);
  fit       = FF(1:N);
  % Global best
  if fit(1) < fitG
    fitG = fit(1);
    Xgb  = X(1,:);
  end
  curve(t) = fitG;
  fprintf('\nIteration %d Best (TGA)= %f',t,curve(t))
  t = t + 1;
end
% Select features
Pos   = 1:dim;
Sf    = Pos((Xgb > thres) == 1);
sFeat = feat(:,Sf);
% Store results
TGA.sf = Sf;
TGA.ff = sFeat;
TGA.nf = length(Sf);
TGA.c  = curve;
TGA.f  = feat;
TGA.l  = label;
end

author: JingweiToo

repo: Wrapper-Feature-Selection-Toolbox