Skip to content
Snippets Groups Projects
Select Git revision
  • 9f88329ebdf1302ab5c3332ab688e5888ad337c7
  • master default
  • dev
  • ci/windows
  • feature/installer
  • release/1.0.2
  • feature/tox
  • pydoc
  • jenkins
  • issue10
  • json-config
  • v2.0.0
  • v2.0.0-beta.5
  • v1.0.2
  • v2.0.0-beta.4
  • v2.0.0-beta.3
  • v1.0.1
  • v1.0.0
18 results

tutorial.rst

Blame
  • Program.cs 62.35 KiB
    using MathNet.Numerics.LinearAlgebra;
    using MathNet.Numerics.Random;
    using OxyPlot;
    using OxyPlot.Axes;
    using OxyPlot.Series;
    using System;
    using System.Collections.Generic;
    using System.IO;
    using System.Linq;
    using System.Threading.Tasks;
    
    namespace IVMCTrim
    {
        class Program
        {
            // This is the entry point.
            // Consult readme.md for additional information about how to run this program.
            // Consult fjn1g13@soton.ac.uk if you don't understand something.
            // This code is not exemplar C#; rather, the idea is that it should be fairly easy for anyone to reimplement the model, but it should also reproduce the experiments from the paper.
            static void Main()
            {
                // none of our matrices are big enough to warrant multi-threading
                MathNet.Numerics.Control.UseSingleThread();
    
                // runs a bunch of experiments with the 'standard' configuration, only varying Z
                // seeds should match those used in the paper
                Run("Z000", Configurations.Standard(Z: 0.0), seed: 42); // (no dual-peaked environments)
                Run("Z050", Configurations.Standard(Z: 0.5), seed: 542);
                Run("Z075", Configurations.Standard(Z: 0.75), seed: 1842);
                Run("Z095", Configurations.Standard(Z: 0.95), seed: 2042); // Corresponds to Figures 3&4
                Run("Z000PerfectG", Configurations.PerfectG(Z: 0.0), seed: 42); // (we align G with the environment automatically and disable mutations in G, eliminating any oportunity to experiance an evolvability benefit)
    
                // prepare some experiments with unusual configurations
                Run("FixedHigh", Configurations.FixedHigh, seed: 42); // Corresponds to Figure 6
                Run("Z095L2", Configurations.StandardL2(Z: 0.95), seed: 2042); // (example with super-additive cost function (quadratic))
                Run("Z095MMSO", Configurations.StandardMMSO(Z: 0.95), seed: 2042); // (example with sub-additive cost function)
                Run("Z050E204", Configurations.Standard204(Z: 0.95), seed: 42, reportingPeriodMultiplier: 10); // (example with more modules; this one takes a while to run)
    
                // Vary λ spread (will take many days...)
                RunSpread("VaryL", new[] { 0.00, 0.02, 0.04, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16, 0.18, 2.00 }, λ => Configurations.Standard(Z: 0.95, λ: λ), λ => (100 * λ).ToString("000"), 42, 40);
    
                // Vary Z spread (will take many many days...)
                RunSpread("VaryZ", new[] { 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95 }, z => Configurations.Standard(Z: z), z => (100 * z).ToString("000"), 42, 40);
            }
    
            /// <summary>
            /// Runs a single experiment with the given parameters and seed producing output into the given output directory
            /// </summary>
            public static void Run(string outputDirectory, Parameters parameters, int seed, int reportingPeriodMultiplier = 1)
            {
                Console.WriteLine($"Running into {outputDirectory} with seed {seed}");
                System.IO.Directory.CreateDirectory(outputDirectory);
    
                // initialise reporter
                var reporter = Recording.CreateEndOfEpisodeReporter(outputDirectory, parameters, reportingPeriodMultiplier, out var recorders);
    
                // initialise random source
                var ctx = new ExecutionContext(new MersenneTwister(seed, false));
    
                // run the actual experiment
                var runonSolveFrequency = Model.RunExperiment(ctx, parameters, reporter);
    
                // produce output from recorders
                Recording.ProduceOutput(outputDirectory, parameters, recorders, runonSolveFrequency);
            }
    
            /// <summary>
            /// Runs many experiments in parallel with the same parameters and sequential seeds.
            /// </summary>
            public static void RunRepeats(string outputDirectory, Parameters parameters, int baseSeed, int repeats, int startFrom = 0, int reportingPeriodMultiplier = 1)
            {
                var tasks = new List<Action>();
    
                baseSeed += startFrom;
    
                for (int r = startFrom; r < repeats + startFrom; r++)
                {
                    // increment seeds
                    int seed = baseSeed++;
    
                    // stuff it into a sub-dir
                    var dir = System.IO.Path.Combine(outputDirectory, $"r{r}");
    
                    // run
                    tasks.Add(new Action(() => Run(dir, parameters, seed, reportingPeriodMultiplier)));
                }
    
                // run them all in parallel
                Parallel.ForEach(tasks, new ParallelOptions { MaxDegreeOfParallelism = 8 }, x => x());
            }
    
            /// <summary>
            /// Run a spread of experiments with different parameters.
            /// </summary>
            public static void RunSpread<T>(string outputDirectory, T[] spread, Func<T, Parameters> parameterFunc, Func<T, string> formatter, int baseSeed, int repeats, int startFrom = 0, int reportingPeriodMultiplier = 1) where T : IComparable<T>
            {
                // perform runs
                foreach (var t in spread)
                {
                    var subDir = Path.Combine(outputDirectory, formatter(t));
                    var parameters = parameterFunc(t);
                    RunRepeats(subDir, parameters, baseSeed, repeats, startFrom, reportingPeriodMultiplier);
    
                    baseSeed += 100;
                }
    
                // generate runon solve frequency plot
                PlotRunonSolveFrequencies(outputDirectory);
            }
    
            public static void PlotRunonSolveFrequencies(string directory)
            {
                var runonSolveFrequenciesTable = new Dictionary<string, double[]>();
    
                foreach (var subDir in Directory.GetDirectories(directory))
                {
                    var subName = Path.GetFileName(subDir);
                    var list = new List<double>();
    
                    foreach (var runDir in Directory.GetDirectories(subDir))
                    {
                        var runName = Path.GetFileName(runDir);
                        var fileContents = File.ReadAllText(Path.Combine(runDir, "runonSolveFrequency.txt"));
                        list.Add(double.Parse(fileContents));
                    }
    
                    runonSolveFrequenciesTable.Add(subName, list.ToArray());
                }
    
                var runonSolveFrequencyPlot = new PlotModel();
                var xaxis = new CategoryAxis() { Position = AxisPosition.Bottom };
                var yaxis = new LinearAxis() { Position = AxisPosition.Left, Title = "Mean Frequency of Max Fitness Phenotype" };
                runonSolveFrequencyPlot.Axes.Add(xaxis);
                runonSolveFrequencyPlot.Axes.Add(yaxis);
                var columns = new ColumnSeries();
                foreach (var t in runonSolveFrequenciesTable.Keys.OrderBy(x => x))
                {
                    xaxis.Labels.Add(t);
                    columns.Items.Add(new ColumnItem(runonSolveFrequenciesTable[t].Average()));
                }
                runonSolveFrequencyPlot.Series.Add(columns);
                runonSolveFrequencyPlot.ExportPdf(Path.Combine(directory, "meanRunonSolveFrequency.pdf"), 1.0, false);
            }
        }
    
        /*
         * Herebelow is the model and all that.
         */
    
        /// <summary>
        /// High-level model operations.
        /// </summary>
        public static class Model
        {
            /// <summary>
            /// Represents an operation to perform at the end of an episode, for reporting/recording purposes.
            /// </summary>
            public delegate void EndOfEpisodeReporter(IvmcEnvironment environment, IndividualEvaluation individualEvaluation, int episode);
    
            /// <summary>
            /// Represents an operation to perform at the end of an episode, for reporting/recording purposes.
            /// </summary>
            public delegate void EndOfEvolutionaryStepReporter(IvmcEnvironment environment, IndividualEvaluation individualEvaluation, int episode, int evolutionaryStep);
    
            /// <summary>
            /// Runs a single experiment with the given parameters in the given context, reporting the completion of each episode.
            /// </summary>
            /// <returns>The runon frequency.</returns>
            public static double RunExperiment(ExecutionContext ctx, Parameters parameters, EndOfEpisodeReporter endOfEpisode)
            {
                // prepare the initial individual
                int N = parameters.ModuleCount * parameters.ModuleSize;
                var individual = Individual.CreateZeroGenome(N);
                individual.DevelopInplace(ctx, parameters);
    
                // report initial state
                endOfEpisode?.Invoke(null, new IndividualEvaluation(individual, default), 0);
    
                // prepare the environment
                var environment = new IvmcEnvironment(parameters.ModuleSize, parameters.ModuleCount);
    
                // run many episode
                for (int episode = 1; episode <= parameters.TotalEpisodes; episode++)
                {
                    // run one episode
                    var result = Model.RunEpisode(episode, ctx, parameters, environment, ref individual, null);
    
                    // report state at end of episode
                    endOfEpisode?.Invoke(environment, result, episode);
                }
    
                // compute runon stuff
                var runonFrequency = MeasureRunonSolveFrequency(ctx, parameters, individual);
                return runonFrequency;
            }
    
            /// <summary>
            /// Runs an episode with the given parameters, environment, and individual.
            /// The individual may be modified or replaced with a different instance.
            /// </summary>
            public static IndividualEvaluation RunEpisode(int episode, ExecutionContext ctx, Parameters parameters, IvmcEnvironment environment, ref Individual individual, EndOfEvolutionaryStepReporter endOfEvolutionaryStepReporter)
            {
                // perform start of episode action (e.g. randomise environment)
                parameters.EpisodeStartOperation(ctx, parameters, individual, environment);
    
                // reset check (for compatibilty with M4M, not important for model)
                ctx.Rand.NextDouble();
    
                // evaluate the current individual
                var evaluation = Evaluate(ctx, parameters, environment, individual);
                var result = new IndividualEvaluation(individual, evaluation);
    
                endOfEvolutionaryStepReporter?.Invoke(environment, result, episode, 0);
    
                // initialise an candidate individual
                var candidate = Individual.CreateZeroGenome(individual.Size);
    
                for (int evolutionaryStep = 0; evolutionaryStep < parameters.EvolutionaryStepsPerEpisode; evolutionaryStep++)
                {
                    // mutate the current individual into the candidate individual
                    Individual.Mutate(ctx, parameters, individual, candidate);
    
                    // evaluate the candidate
                    var candidateEvaluation = Evaluate(ctx, parameters, environment, candidate);
    
                    result = new IndividualEvaluation(individual, evaluation);
    
                    // keep the candidate if he is best
                    if (candidateEvaluation.Fitness > evaluation.Fitness)
                    {
                        // swap
                        (individual, candidate) = (candidate, individual);
    
                        evaluation = candidateEvaluation;
                    }
    
                    // report
                    endOfEvolutionaryStepReporter?.Invoke(environment, new IndividualEvaluation(individual, evaluation), episode, evolutionaryStep + 1);
                }
    
                // NOTE NOTE: we return the penultimate individual for compatibility with M4M
                return result;
            }
    
            /// <summary>
            /// Evaluates the fitness of the given individual in the given environment with the given parameters.
            /// </summary>
            /// <remarks>
            /// Eq2 in the Model section under Evolution
            /// </remarks>
            public static Evaluation Evaluate(ExecutionContext ctx, Parameters parameters, IvmcEnvironment environment, Individual individual)
            {
                var benefit = environment.EvaluateBenefit(individual.P);
                var cost = parameters.CostFunction(individual.B);
                var fitness = benefit - cost * parameters.CostCoefficient;
    
                return new Evaluation(fitness, benefit, cost);
            }
    
            /// <summary>
            /// Performs a run-on, and counts the proportion of solves.
            /// </summary>
            /// <remarks>
            /// So long as the variation in the environment has a unbiased distribution, then this is an unbiased runon.
            /// </remarks>
            public static double MeasureRunonSolveFrequency(ExecutionContext ctx, Parameters parameters, Individual individual, int runonEpisodeCount = 100)
            {
                individual = individual.Clone();
    
                int[] wins = new int[parameters.ModuleCount + 1];
    
                var environment = new IvmcEnvironment(parameters.ModuleSize, parameters.ModuleCount);
    
                for (int episode = 0; episode < runonEpisodeCount; episode++)
                {
                    var result = Model.RunEpisode(episode, ctx, parameters, environment, ref individual, null);
                    Recorders.CountSolves(environment, result.Inidividual, result.Evaluation, wins);
                }
    
                return (double)wins.Last() / runonEpisodeCount;
            }
    
            /// <summary>
            /// Runs a few episodes, recording the fitness per evolutionary step with the given sample period.
            /// </summary>
            /// <returns>
            /// A list of <see cref="DataPoint"/> comprising the evolutionary step and fitness.
            /// </returns>
            public static List<DataPoint> RunTracee(Parameters parameters, Individual individual, int traceeEpisodeCount = 4, int samplePeriod = 100)
            {
                individual = individual.Clone();
    
                // run tracee
                var ctx = new ExecutionContext(new MersenneTwister(3)); // common seed (chosen so that it has a mix of cL and c0)
                var env = new IvmcEnvironment(parameters.ModuleSize, parameters.ModuleCount);
    
                var fitnessSamples = new List<DataPoint>();
    
                void sample(IvmcEnvironment _env, IndividualEvaluation stepIndividualEvaluation, int _episode, int evolutionaryStep)
                {
                    if (evolutionaryStep % samplePeriod == 0)
                        fitnessSamples.Add(new DataPoint(_episode * parameters.EvolutionaryStepsPerEpisode + evolutionaryStep, stepIndividualEvaluation.Evaluation.Fitness));
                }
    
                for (int traceeEpisodes = 0; traceeEpisodes < traceeEpisodeCount; traceeEpisodes++)
                {
                    Model.RunEpisode(traceeEpisodes, ctx, parameters, env, ref individual, sample);
                }
    
                return fitnessSamples;
            }
        }
    
        public static class Configurations
        {
            /// <summary>
            /// Prepares a 'standard' 4x4 set of parameters with the given Z value.
            /// </summary>
            public static Parameters Standard(double Z, double λ = 0.1)
            {
                return new Parameters()
                {
                    DualPeakedEnvironmentalInstanceFrequency = Z,
                    CostCoefficient = λ,
                };
            }
    
            /// <summary>
            /// Prepares a 'standard' 4x4 set of parameters with the given Z value and an L2 cost function.
            /// </summary>
            public static Parameters StandardL2(double Z, double λ = 0.02)
            {
                return new Parameters()
                {
                    DualPeakedEnvironmentalInstanceFrequency = Z,
                    CostCoefficient = λ,
                    CostFunction = CostFunctions.JudgeCostL2,
                };
            }
    
            // TODO: add L1And2 if I need to
    
            /// <summary>
            /// Prepares a 'standard' 4x4 set of parameters with the given Z value and a sub-additive cost function.
            /// </summary>
            public static Parameters StandardMMSO(double Z, double λ = 0.1)
            {
                return new Parameters()
                {
                    DualPeakedEnvironmentalInstanceFrequency = Z,
                    CostCoefficient = λ,
                    CostFunction = CostFunctions.JudgeCostMMSO,
                };
            }
    
            /// <summary>
            /// Prepares a 'standard' 20x4 set of parameters with the given Z value.
            /// </summary>
            public static Parameters Standard204(double Z)
            {
                return new Parameters()
                {
                    DualPeakedEnvironmentalInstanceFrequency = Z,
                    ModuleCount = 20,
                    ModuleSize = 4,
                    TotalEpisodes = 500000,
                    BMutationRate = 0.0002,
                    CostCoefficient = 8,
                    EvolutionaryStepsPerEpisode = 2000,
                };
            }
    
            /// <summary>
            /// Prepares a 4x4 set of parameters with the given Z value,
            /// where G is aligned with the environment at the start of each episode.
            /// </summary>
            public static Parameters PerfectG(double Z)
            {
                return new Parameters()
                {
                    DualPeakedEnvironmentalInstanceFrequency = Z,
                    EpisodeStartOperation = EpisodeStartOperations.VaryEnvironmentWithPerfectG,
                    MutateBProbability = 1, // disable mutations in G
                };
            }
    
            /// <summary>
            /// Prepares a 4x4 set of parameters with the given Z value,
            /// where G and the environment are set high at the start of the environment.
            /// </summary>
            public static Parameters FixedHigh
            {
                get
                {
                    return new Parameters()
                    {
                        EpisodeStartOperation = EpisodeStartOperations.FixEnvironmentAndG,
                        TotalEpisodes = 30,
                    };
                }
            }
        }
    
        /// <summary>
        /// The parameters for an experiment.
        /// </summary>
        public class Parameters
        {
            // default should be equivalent to the 'main' figure in the poaper
            // m4m gen=composedense bprob=0.5 BEx=true lambda=0.1 Z=0.95 regfunc=l1 ch=1.0 cl=0.7 c0=-1.0 epochs=400000 mb=1E-3 judgemode=Split K=1000 environmentpackage=ivmc:4 transmatmutator=singlecell
    
            public int ModuleCount { get ; set; } = 4;
    
            public int ModuleSize { get ; set; } = 4;
    
            public double BMutationRate { get; set; } = 0.001;
            public double MB => BMutationRate;
    
            public bool AllowGMutation { get; set; } = true;
    
            public int NumberOfDevelopmentalTimeSteps { get; set; } = 10;
            public int T => NumberOfDevelopmentalTimeSteps;
    
            public double DualPeakedEnvironmentalInstanceFrequency { get; set; } = 0.95;
            public double Z => DualPeakedEnvironmentalInstanceFrequency;
    
            public double CostCoefficient { get; set; } = 0.1;
            public double λ => CostCoefficient;
    
            public CostFunction CostFunction { get; set; } = CostFunctions.JudgeCostL1;
            public CostFunction ϕ => CostFunction; // ϕ(B) should equal Sum_ij(φ|B_ij|)
    
            public double DecayRate { get; set; } = 0.2;
            public double τ => DecayRate;
    
            public Func<double, double> SquashFunction { get; set; } = x => Math.Tanh(x * 0.5);
            public Func<double, double> σ => SquashFunction;
    
            public double HighPeakBenefit { get; set; } = 1.0;
            public double CH => HighPeakBenefit;
    
            public double LowPeakBenefit { get; set; } = 0.7;
            public double CL => LowPeakBenefit;
    
            public double NoPeakBenefit { get; set; } = -1;
            public double C0 => NoPeakBenefit;
    
            public int EvolutionaryStepsPerEpisode { get; set; } = 1000;
            public int K => EvolutionaryStepsPerEpisode;
    
            public double MutateBProbability { get; set; } = 0.5;
            public double RB => MutateBProbability;
    
            public int TotalEpisodes { get; set; } = 400_000;
    
            public EpisodeStartOperation EpisodeStartOperation { get; set; } = EpisodeStartOperations.VaryEnvironment;
        }
    
        /// <summary>
        /// A function that determines the connection cost of a grn.
        /// See <see cref="CostFunctions"/> for standard implementations.
        /// </summary>
        public delegate double CostFunction(Matrix<double> grn);
    
        /// <summary>
        /// An operation to be performed at the start of each episode.
        /// See <see cref="EpisodeStartOperations"/> for standard implementations.
        /// </summary>
        public delegate void EpisodeStartOperation(ExecutionContext ctx, Parameters parameters, Individual individual, IvmcEnvironment environment);
    
        /// <summary>
        /// Represents an individual, with genotype and phenotype.
        /// </summary>
        public class Individual
        {
            private Individual(Vector<double> G, Matrix<double> B)
            {
                if (G.Count != B.RowCount || B.RowCount != B.ColumnCount)
                    throw new ArgumentException("Mismatched dimensions");
    
                this.G = G;
                this.B = B;
                P = CreateVector.Dense<double>(Size);
            }
    
            /// <summary>
            /// Creates an individual with all zeros in the genome.
            /// Does not (can not) develop the individual.
            /// </summary>
            public static Individual CreateZeroGenome(int size)
            {
                return new Individual(CreateVector.Dense<double>(size), CreateMatrix.Dense<double>(size, size));
            }
    
            /// <summary>
            /// The embryonic gene expression of the individual.
            /// </summary>
            public Vector<double> G { get; }
            
            /// <summary>
            /// The regulatory matrix of the individual.
            /// </summary>
            public Matrix<double> B { get; }
            
            /// <summary>
            /// The phenotype vector of the individual.
            /// Computed by <see cref="DevelopInplace(ExecutionContext, Parameters)"/> from <see cref="G"/> and <see cref="B"/>.
            /// </summary>
            public Vector<double> P { get; }
    
            /// <summary>
            /// The number of 'expressed' genes.
            /// </summary>
            public int Size => G.Count;
    
            /// <summary>
            /// Produces a mutant type from the given <paramref name="parent"/>.
            /// This is done by cloning the parent's genome, performing a mutation, and then developing the  <paramref name="offspring"/> individual.
            /// These operations are perform in-place in the offspring individual.
            /// </summary>
            /// <remarks>
            /// This implements B-Exclusive (BEx) and SingleCell mutation for B, and 'binary' mutation for G.
            /// These are the main departures from the models of Kouvaris et al. and Kounois et al.
            /// </remarks>
            public static void Mutate(ExecutionContext ctx, Parameters parameters, Individual parent, Individual offspring)
            {
                var rand = ctx.Rand;
    
                if (parent.Size != offspring.Size)
                    throw new ArgumentException("Parent and offspring dimensions do not match.");
    
                // clone
                parent.B.CopyTo(offspring.B);
                parent.G.CopyTo(offspring.G);
    
                // mutate
                bool mutateB = rand.NextBool(parameters.MutateBProbability);
                if (mutateB)
                {
                    // pick a random cell
                    int r = rand.Next(parent.Size);
                    int c = rand.Next(parent.Size);
    
                    // draw delta
                    double delta = rand.NextSignedUniform(parameters.BMutationRate);
    
                    // additive update (unclamped)
                    offspring.B[r, c] = delta + offspring.B[r, c];
                }
                else if (parameters.AllowGMutation) // mutate G
                {
                    // pick a random gene
                    int r = rand.Next(parent.Size);
    
                    // assign randomly (NOTE: half of these mutations do nothing at all)
                    offspring.G[r] = rand.NextBool(0.5) ? 1 : -1;
                }
    
                // develop
                offspring.DevelopInplace(ctx, parameters);
            }
    
            /// <summary>
            /// Develops the individuals' phenotype from it's genotype.
            /// </summary>
            /// <remarks>
            /// Eq1 in the Model section under Development
            /// </remarks>
            public void DevelopInplace(ExecutionContext ctx, Parameters parameters)
            {
                // t=0
                var state = G.Clone();
    
                // t=1..T
                for (int t = 1; t <= parameters.NumberOfDevelopmentalTimeSteps; t++)
                {
                    // y_i(t+1) = y_i(t)*(1-τ) + σ(sum_j(B_{ij}*y_j(t)))
                    state = state * (1 - parameters.DecayRate) + B.Multiply(state).Map(parameters.SquashFunction);
                }
    
                // normalise to [-1, +1], and place result in our phenotype
                state.Multiply(parameters.DecayRate, P);
            }
    
            /// <summary>
            /// Clones this individual.
            /// </summary>
            public Individual Clone()
            {
                var clone = new Individual(G.Clone(), B.Clone());
                P.CopyTo(clone.P);
                return clone;
            }
        }
    
        /// <summary>
        /// Represents an individual and its evaluation.
        /// </summary>
        public struct IndividualEvaluation
        {
            public IndividualEvaluation(Individual inidividual, Evaluation evaluation)
            {
                Inidividual = inidividual ?? throw new ArgumentNullException(nameof(inidividual));
                Evaluation = evaluation;
            }
    
            public Individual Inidividual { get; }
            public Evaluation Evaluation { get; }
        }
    
        /// <summary>
        /// Represents the benefit (b), cost (c), and combined fitness (f) of an individual.
        /// f = b - λc
        /// </summary>
        public struct Evaluation
        {
            public Evaluation(double fitness, double benefit, double cost)
            {
                Fitness = fitness;
                Benefit = benefit;
                Cost = cost;
            }
    
            public double Fitness { get; }
            public double Benefit { get; }
            public double Cost { get; }
        }
    
        /// <summary>
        /// Represents an IVMC changing environment.
        /// </summary>
        public class IvmcEnvironment
        {
            public IvmcEnvironment(int moduleSize, int moduleCount)
            {
                ModuleSize = moduleSize;
                ModuleCount = moduleCount;
    
                PlusCoefs = new double[ModuleCount];
                MinusCoefs = new double[ModuleCount];
            }
    
            /// <summary>
            /// Gets the size of each module in the environment.
            /// </summary>
            public int ModuleSize { get; }
    
            /// <summary>
            /// Gets the number of modules in the environment.
            /// </summary>
            public int ModuleCount { get; }
    
            /// <summary>
            /// Gets the number of elements in the module.
            /// </summary>
            public int Size => ModuleCount * ModuleSize;
    
            public double[] PlusCoefs { get; }
            public double[] MinusCoefs { get; }
    
            /// <summary>
            /// Randomises the benefit coefficients.
            /// </summary>
            public void Randomise(ExecutionContext ctx, Parameters parameters)
            {
                var rand = ctx.Rand;
    
                for (int i = 0; i < ModuleCount; i++)
                {
                    bool dualPeaked = rand.NextBool(parameters.Z);
                    bool minusHigh = rand.NextBoolean();
    
                    PlusCoefs[i] = !minusHigh
                        ? parameters.CH
                        : (dualPeaked ? parameters.CL : parameters.C0);
                    MinusCoefs[i] = minusHigh
                        ? parameters.CH
                        : (dualPeaked ? parameters.CL : parameters.C0);
                }
            }
    
            /// <summary>
            /// Determines the benefit award to given phenotype by the environment in its current state.
            /// </summary>
            /// <remarks>
            /// Eq4 in the Model section under Environment
            /// </remarks>
            public double EvaluateBenefit(Vector<double> phenotype)
            {
                double benefitAcc = 0;
    
                for (int i = 0; i < ModuleCount * ModuleSize; i += ModuleSize)
                {
                    double acc = 0;
                    for (int j = i; j < i + ModuleSize; j++)
                        acc += phenotype[j];
    
                    acc /= ModuleSize; // mean expression within module
    
                    double cplus = PlusCoefs[i / ModuleSize];
                    double cminus = MinusCoefs[i / ModuleSize];
                    double moduleFitness = Math.Max(acc * cplus, -acc * cminus);
    
                    benefitAcc += moduleFitness;
                }
    
                return benefitAcc / ModuleCount; // mean benefit of all modules
            }
        }
    
        /// <summary>
        /// Common reset operations functions.
        /// </summary>
        public static class EpisodeStartOperations
        {
            /// <summary>
            /// Varies the environment every episode.
            /// </summary>
            public static void VaryEnvironment(ExecutionContext ctx, Parameters parameters, Individual individual, IvmcEnvironment environment)
            {
                environment.Randomise(ctx, parameters);
            }
    
            /// <summary>
            /// Varies the environment every episode,
            /// and sets G as per the ideal phenotype for the new environmental conditions.
            /// </summary>
            public static void VaryEnvironmentWithPerfectG(ExecutionContext ctx, Parameters parameters, Individual individual, IvmcEnvironment environment)
            {
                VaryEnvironment(ctx, parameters, individual, environment);
    
                var perfectP = environment.PerfectP();
                perfectP.CopyTo(individual.G);
    
                individual.DevelopInplace(ctx, parameters);
            }
    
            /// <summary>
            /// Sets the environment to all D+ environment conditions,
            /// and sets G to be all 1s.
            /// </summary>
            public static void FixEnvironmentAndG(ExecutionContext ctx, Parameters parameters, Individual individual, IvmcEnvironment environment)
            {
                for (int i = 0; i < environment.ModuleCount; i++)
                {
                    environment.PlusCoefs[i] = parameters.HighPeakBenefit;
                    environment.MinusCoefs[i] = parameters.LowPeakBenefit;
                }
    
                individual.G.Clear();
                individual.G.Add(1, individual.G);
    
                individual.DevelopInplace(ctx, parameters);
            }
        }
    
        /// <summary>
        /// Common cost functions.
        /// </summary>
        public static class CostFunctions
        {
            /// <summary>
            /// Linear (L1) cost function.
            /// Corresponds to φ(x) = x
            /// </summary>
            /// <remarks>
            /// Eq3 in the Model section under Evolution
            /// </remarks>
            public static double JudgeCostL1(Matrix<double> grn)
            {
                double acc = 0.0;
    
                for (int i = 0; i < grn.RowCount; i++)
                {
                    for (int j = 0; j < grn.ColumnCount; j++)
                    {
                        acc += Math.Abs(grn[i, j]);
                    }
                }
    
                return acc / ((double)grn.ColumnCount * grn.RowCount);
            }
    
            /// <summary>
            /// Quadratic (L2) cost function.
            /// Corresponds to φ(x) = x^2
            /// </summary>
            public static double JudgeCostL2(Matrix<double> grn)
            {
                double acc = 0.0;
    
                for (int i = 0; i < grn.RowCount; i++)
                {
                    for (int j = 0; j < grn.ColumnCount; j++)
                    {
                        acc += grn[i, j] * grn[i, j];
                    }
                }
    
                return acc / ((double)grn.ColumnCount * grn.RowCount);
            }
    
            /// <summary>
            /// An arbitrary sub-additive cost function.
            /// </summary>
            public static double JudgeCostMMSO(Matrix<double> grn)
            {
                // sum of a linear component and a michaelis-mention style function
                // gradient is 1 at 0
    
                double a = 1.0;
                double b = 0.2; // non-linear (sub-additive) contribution
                double c = 0.8; // linear contribution
    
                double acc = 0.0;
    
                for (int i = 0; i < grn.RowCount; i++)
                {
                    for (int j = 0; j < grn.ColumnCount; j++)
                    {
                        var d = Math.Abs(grn[i, j]);
                        acc += c * d + b * d / (d + a);
                    }
                }
    
                return acc / ((double)grn.ColumnCount * grn.RowCount);
            }
        }
    
        /// <summary>
        /// Analysis helpers.
        /// </summary>
        public static class Analysis
        {
            /// <summary>
            /// Determines the 'perfect' phenotype for an environment in its present conditions.
            /// </summary>
            public static Vector<double> PerfectP(this IvmcEnvironment environment)
            {
                var p = CreateVector.Dense<double>(environment.Size);
    
                for (int i = 0; i < environment.ModuleCount; i++)
                {
                    var val = environment.PlusCoefs[i] > environment.MinusCoefs[i] ? 1 : -1;
    
                    for (int j = 0; j < environment.ModuleSize; j++)
                    {
                        p[i * environment.ModuleSize + j] = val;
                    }
                }
    
                return p;
            }
    
            /// <summary>
            /// Computes the degree of hierarchy of sequencial block modules of the given size
            /// </summary>
            public static IEnumerable<double> ComputeBlockModuleHierarchy(Matrix<double> dtm, int moduleSize)
            {
                for (int i = 0; i < dtm.ColumnCount; i += moduleSize)
                {
                    yield return ComputeModuleHierarchy(dtm, Enumerable.Range(i, moduleSize).ToList());
                }
            }
    
            /// <summary>
            /// Computes the degree of hierarchy of sequencial block modules of the given size
            /// </summary>
            public static IEnumerable<double> ComputeBlockModuleIndependence(Matrix<double> dtm, int moduleSize)
            {
                for (int i = 0; i < dtm.ColumnCount; i += moduleSize)
                {
                    yield return ComputeModuleIndependence(dtm, Enumerable.Range(i, moduleSize).ToList());
                }
            }
    
            /// <summary>
            /// Computes the degree of hierarchy of a module containing 2 or more traits
            /// degree of hierarchy = (max_col_weight / total_weight), scaled between 0 (dense/nothing) and 1 (ideal hierarchy)
            /// Returns 0 if there are no weights
            /// </summary>
            public static double ComputeModuleHierarchy(Matrix<double> dtm, IReadOnlyList<int> moduleTraits)
            {
                double[] columnSums = moduleTraits.Select(j => moduleTraits.Sum(i => Math.Abs(dtm[i, j]))).ToArray();
    
                int count = columnSums.Length;
                if (count < 2)
                    throw new ArgumentException(nameof(moduleTraits), "A module must comprise 2 or more traits");
    
                double max = columnSums.Max();
                if (max == 0) // avoid div0: this is not a hierarchy
                    return 0.0;
    
                double total = columnSums.Sum();
                double hratio = max / total;
    
                if (double.IsNaN(hratio) || double.IsInfinity(hratio))
                    return double.NaN; // probably ought throw...
    
                // re-scale hratio from [1/n, 1] to [0, 1]
                int n = count;
                return (hratio - (1.0 / n)) / (1.0 - (1.0 / n));
            }
    
            /// <summary>
            /// Computes the independence of a module
            /// module_independence = (weight_within_module / weight_in_module_row), on scale between 0 (no control) and 1 (total independence)
            /// Returns NaN if there are no weights
            /// </summary>
            public static double ComputeModuleIndependence(Matrix<double> dtm, IReadOnlyList<int> moduleTraits)
            {
                double withinSum = 0.0;
                double rowSum = 0.0;
    
                foreach (var i in moduleTraits)
                {
                    foreach (var j in moduleTraits)
                    {
                        withinSum += Math.Abs(dtm[i, j]);
                    }
    
                    for (int j = 0; j < dtm.ColumnCount; j++)
                    {
                        rowSum += Math.Abs(dtm[i, j]);
                    }
                }
    
                var iratio = withinSum / rowSum;
    
                if (double.IsNaN(iratio) || double.IsInfinity(iratio))
                    return double.NaN;
    
                return iratio;
            }
        }
    
        /*
         * Herebelow are things that have nothing to do with the model.
         */
    
        public static class Recording
        {
            public static Model.EndOfEpisodeReporter CreateEndOfEpisodeReporter(string outputDirectory, Parameters parameters, int reportingPeriodMultiplier, out Recorders recorders)
            {
                string path(string f) => System.IO.Path.Combine(outputDirectory, f);
    
                // prepare the recorders
                var _recorders = recorders = new Recorders(parameters.ModuleSize, parameters.ModuleCount, reportingPeriodMultiplier);
    
                // endOfEpisode reporting actions
                void endOfEpisode(IvmcEnvironment environment, IndividualEvaluation individualEvaluation, int episode)
                {
                    var individual = individualEvaluation.Inidividual;
                    var evaluation = individualEvaluation.Evaluation;
    
                    // inform recorders
                    _recorders.Sample(environment, individual, evaluation, episode);
    
                    if (episode % (reportingPeriodMultiplier * 1000) == 0)
                    {
                        // report progress
                        Console.WriteLine($"{outputDirectory}: {episode}/{parameters.TotalEpisodes} ({(double)episode / parameters.TotalEpisodes:P})\tf={evaluation.Fitness}");
                    }
    
                    if (episode % (reportingPeriodMultiplier * 10000) == 0)
                    {
                        // produce grn plot
                        var grnPlot = PlotHelpers.IntegralCartesianAxes($"Episode {episode}", "row", "col");
                        grnPlot.Axes.Add(PlotHelpers.Gray());
                        grnPlot.Series.Add(PlotHelpers.PlotMatrix(individual.B));
                        grnPlot.ExportPdf(path($"episode{episode}.pdf"), 0.5, true);
    
                        // run tracee
                        var fitnessSamples = Model.RunTracee(parameters, individual);
    
                        // plot tracee
                        var traceePlot = PlotHelpers.LinearAxes($"Episode {episode}", "Evolutionary Step", "Fitness");
                        traceePlot.Axes[0].MajorGridlineStyle = LineStyle.Solid;
                        traceePlot.Axes[1].Minimum = -0.02;
                        traceePlot.Axes[1].Maximum = 1.02;
                        var fitnessLine = new LineSeries() { Color = OxyColors.Red, MarkerType = MarkerType.Circle };
                        fitnessLine.Points.AddRange(fitnessSamples);
                        traceePlot.Series.Add(fitnessLine);
                        traceePlot.ExportPdf(path($"tracee{episode}.pdf"), 0.5, false, 0.5);
                    }
                }
    
                return endOfEpisode;
            }
    
            public static void ProduceOutput(string outDir, Parameters parameters, Recorders recorders, double runonSolveFrequency)
            {
                string path(string f) => System.IO.Path.Combine(outDir, f);
    
                const double plotSize = 0.5;
    
                Console.WriteLine($"{outDir} runon frequency: {runonSolveFrequency}");
                System.IO.File.WriteAllText(path("runonSolveFrequency.txt"), runonSolveFrequency.ToString("R"));
    
                // data export and mindless plotting
                var rcsData = recorders.GetRcsDataFlat();
                TrajectoryHelpers.SaveTrajectories(path("rcs_t.dat"), rcsData, recorders.Grns.SamplePeriod);
    
                var rcsPlot = PlotHelpers.LinearAxes("Regulation Coefficients", "Episodes", "Connection Strength");
                rcsPlot.Series.AddRange(PlotHelpers.PlotTrajectoryLines2D(recorders.GetRcsData(), recorders.Grns.SamplePeriod, OxyColors.DarkCyan, OxyColors.LightCyan, OxyColors.DarkGreen, OxyColors.LightGreen));
                rcsPlot.ExportPdf(path("rcs_t.pdf"), plotSize, false);
    
                var hierarchyPlot = PlotHelpers.LinearAxes("Evolution of Hierarchy", "Episodes", "Degree of Hierarchy");
                hierarchyPlot.Series.AddRange(
                    PlotHelpers.PlotTrajectoryLines1D(
                        recorders.Grns.Samples.Select(grn => Analysis.ComputeBlockModuleHierarchy(grn, parameters.ModuleSize).ToList()).ToArray().Transpose(),
                        recorders.Grns.SamplePeriod,
                        OxyColors.DarkSlateGray, OxyColors.SlateGray));
                hierarchyPlot.ExportPdf(path("hierarchy_t.pdf"), plotSize, false);
    
                var moduleIndependencePlot = PlotHelpers.LinearAxes("Evolution of ModuleIndependence", "Episodes", "Degree of Hierarchy");
                moduleIndependencePlot.Series.AddRange(
                    PlotHelpers.PlotTrajectoryLines1D(
                        recorders.Grns.Samples.Select(grn => Analysis.ComputeBlockModuleIndependence(grn, parameters.ModuleSize).ToList()).ToArray().Transpose(),
                        recorders.Grns.SamplePeriod,
                        OxyColors.Teal, OxyColors.DarkTurquoise));
                moduleIndependencePlot.ExportPdf(path("moduleIndependence_t.pdf"), plotSize, false);
    
                var fitnessData = recorders.GetFitnessData();
                TrajectoryHelpers.SaveTrajectories(path("fitness_t.dat"), fitnessData, recorders.Evaluations.SamplePeriod);
    
                var fitnessPlot = PlotHelpers.LinearAxes("Fitness", "Episodes", "Fitness");
                fitnessPlot.Series.AddRange(PlotHelpers.PlotTrajectoryScatters1D(fitnessData.Take(1).ToArray().DownsampleTrajectories(100), recorders.Evaluations.SamplePeriod * 100, OxyColors.Red, OxyColors.DarkOrange));
                fitnessPlot.ExportPdf(path("fitness_t.pdf"), plotSize, false);
    
                var meanFitnessPlot = PlotHelpers.LinearAxes("Mean Fitness", "Episodes", "Mean Fitness");
                meanFitnessPlot.Series.AddRange(PlotHelpers.PlotTrajectoryScatters1D(fitnessData.Take(1).ToArray().MeanDownsampleTrajectories(100), recorders.Evaluations.SamplePeriod * 100, OxyColors.Red, OxyColors.DarkOrange));
                meanFitnessPlot.ExportPdf(path("meanfitness_t.pdf"), plotSize, false);
    
                var solveData = recorders.GetSolveData().Map2D(x => (double)x);
                TrajectoryHelpers.SaveTrajectories(path("ivmcsolves_t.dat"), solveData, recorders.Solves.SamplePeriod);
    
                var solvePlot = PlotHelpers.LinearAxes("Solves", "Episodes", "Solve Frequency");
                solvePlot.Series.AddRange(PlotHelpers.PlotTrajectoryScatters1D(solveData.Map2D(x => x / recorders.Solves.SamplePeriod).TakeLast(1).ToArray(), recorders.Solves.SamplePeriod, OxyColors.IndianRed, OxyColors.DarkRed));
                solvePlot.ExportPdf(path("solves_t.pdf"), plotSize, false);
    
                var switchData = recorders.GetSwitchData().Map2D(x => (double)x);
                TrajectoryHelpers.SaveTrajectories(path("ivmcswitches_t.dat"), switchData, recorders.Switches.SamplePeriod);
    
                var switchPlot = PlotHelpers.LinearAxes("Switches", "Episodes", "Switch Frequency");
                switchPlot.Series.AddRange(PlotHelpers.PlotTrajectoryScatters1D(switchData.Map2D(x => x / recorders.Switches.SamplePeriod).TakeLast(1).ToArray(), recorders.Switches.SamplePeriod, OxyColors.DeepPink, OxyColors.HotPink));
                switchPlot.ExportPdf(path("switches_t.pdf"), plotSize, false);
            }
        }
    
        /// <summary>
        /// Records something over the course of an experiment.
        /// </summary>
        public class Recorder<T>
        {
            public Recorder(Func<IvmcEnvironment, Individual, Evaluation, T> sampler, int samplePeriod)
            {
                Sampler = sampler ?? throw new ArgumentNullException(nameof(sampler));
                SamplePeriod = samplePeriod;
            }
    
            public Func<IvmcEnvironment, Individual, Evaluation, T> Sampler { get; }
            public int SamplePeriod { get; }
    
            public List<T> Samples { get; } = new List<T>();
    
            public void Sample(IvmcEnvironment environment, Individual individual, Evaluation Evaluation, int episode)
            {
                if (episode % SamplePeriod == 0)
                    Samples.Add(Sampler(environment, individual, Evaluation));
            }
        }
    
        /// <summary>
        /// Counts someting over the course of an experiment.
        /// </summary>
        public class CountRecorder
        {
            public CountRecorder(int sampleSize, Action<IvmcEnvironment, Individual, Evaluation, int[]> counter, int samplePeriod)
            {
                Counter = counter;
                SamplePeriod = samplePeriod;
                SampleSize = sampleSize;
            }
    
            public Action<IvmcEnvironment, Individual, Evaluation, int[]> Counter { get; }
            public int SampleSize { get; }
            public int SamplePeriod { get; }
    
            private int[] CurrentCount = null;
            public List<int[]> Samples { get; } = new List<int[]>();
    
            public void Sample(IvmcEnvironment environment, Individual individual, Evaluation Evaluation, int episode)
            {
                if (CurrentCount == null)
                {
                    CurrentCount = new int[SampleSize];
                }
    
                Counter(environment, individual, Evaluation, CurrentCount);
    
                if (episode % SamplePeriod == 0)
                {
                    Samples.Add(CurrentCount);
                    CurrentCount = new int[SampleSize];
                }
            }
        }
    
        /// <summary>
        /// Records information about an experiment.
        /// </summary>
        /// <remarks>
        /// These correspond to the feedback stuff in M4M.
        /// </remarks>
        public class Recorders
        {
            public Recorders(int moduleSize, int moduleCount, int samplePeriodMultiplier = 1)
            {
                Evaluations = new Recorder<Evaluation>((t, i, j) => j, EvaluationSamplePeriod * samplePeriodMultiplier);
                Grns = new Recorder<Matrix<double>>((t, i, j) => i.B.Clone(), GrnSamplePeriod * samplePeriodMultiplier);
                Ps = new Recorder<Vector<double>>((t, i, j) => i.P.Clone(), OtherSamplePeriod * samplePeriodMultiplier);
                Solves = new CountRecorder(moduleCount + 1, CountSolves, OtherSamplePeriod * samplePeriodMultiplier);
                Switches = new CountRecorder(moduleCount + 1, CountSwitches(), OtherSamplePeriod * samplePeriodMultiplier);
            }
    
            public static Action<IvmcEnvironment, Individual, Evaluation, int[]> CountSwitches()
            {
                Vector<double> last = null;
    
                void count(IvmcEnvironment environment, Individual individual, Evaluation Evaluation, int[] switches)
                {
                    if (last == null)
                    {
                        last = individual.P.Clone();
                        return;
                    }
    
                    for (int i = 0; i < environment.ModuleCount; i++)
                    {
                        var val = environment.PlusCoefs[i] > environment.MinusCoefs[i] ? 1 : -1;
                        bool isDualPeaked = environment.PlusCoefs[i] >= 0 && environment.MinusCoefs[i] >= 0;
    
                        bool allOk = true;
                        for (int j = 0; j < environment.ModuleSize; j++)
                        {
                            var idx = i * environment.ModuleSize + j;
                            bool isSwitch = Math.Sign(individual.P[idx]) == -Math.Sign(last[idx]);
    
                            // NOTE: this doesn't check if we were correct (per M4M)
                            if (!isDualPeaked || !isSwitch)
                            {
                                allOk = false;
                                break;
                            }
                        }
    
                        if (allOk)
                        {
                            switches[i]++;
                            switches[environment.ModuleCount]++;
                        }
                    }
    
                    individual.P.CopyTo(last);
                }
    
                return count;
            }
    
            public static void CountSolves(IvmcEnvironment environment, Individual individual, Evaluation Evaluation, int[] wins)
            {
                if (environment == null)
                    return; // pre-start
    
                bool allAllOk = true;
                for (int i = 0; i < environment.ModuleCount; i++)
                {
                    // is a + or - type environmental instance?
                    var val = environment.PlusCoefs[i] > environment.MinusCoefs[i] ? 1 : -1;
    
                    // does the phenotype agree with this module?
                    bool allOk = true;
                    for (int j = 0; j < environment.ModuleSize; j++)
                    {
                        var idx = i * environment.ModuleSize + j;
                        if (Math.Sign(individual.P[idx]) != val)
                        {
                            allOk = false;
                            break;
                        }
                    }
    
                    // did we get the module correct?
                    if (allOk)
                        wins[i]++;
                    else
                        allAllOk = false;
                }
    
                // were all modules correct?
                if (allAllOk)
                    wins[environment.ModuleCount]++;
            }
    
            private const int GrnSamplePeriod = 500;
            private const int EvaluationSamplePeriod = 1;
            private const int OtherSamplePeriod = 500;
    
            public Recorder<Evaluation> Evaluations { get; }
            public Recorder<Matrix<double>> Grns { get; }
            public Recorder<Vector<double>> Ps { get; }
            public CountRecorder Solves { get; }
            public CountRecorder Switches { get; }
    
            public void Sample(IvmcEnvironment environment, Individual individual, Evaluation Evaluation, int episode)
            {
                Evaluations.Sample(environment, individual, Evaluation, episode);
                Grns.Sample(environment, individual, Evaluation, episode);
                Ps.Sample(environment, individual, Evaluation, episode);
                Solves.Sample(environment, individual, Evaluation, episode);
                Switches.Sample(environment, individual, Evaluation, episode);
            }
    
            public double[][] GetFitnessData()
            {
                return new[] {
                    Evaluations.Samples.Select(j => j.Fitness).ToArray(),
                    Evaluations.Samples.Select(j => j.Benefit).ToArray(),
                    Evaluations.Samples.Select(j => j.Fitness - j.Benefit).ToArray() // cost * lambda
                    };
            }
    
            public int[][] GetSolveData()
            {
                return Solves.Samples.Transpose();
            }
    
            public int[][] GetSwitchData()
            {
                return Switches.Samples.Transpose();
            }
    
            public double[][] GetRcsDataFlat()
            {
                var N = Grns.Samples[0].RowCount;
                return Enumerable.Range(0, N * N).Select(i => Grns.Samples.Select(grn => grn[i / N, i % N]).ToArray()).ToArray();
            }
    
            public double[][][] GetRcsData()
            {
                var N = Grns.Samples[0].RowCount;
                return Enumerable.Range(0, N).Select(i => Enumerable.Range(0, N).Select(j => Grns.Samples.Select(grn => grn[i, j]).ToArray()).ToArray()).ToArray();
            }
        }
    
        // C&P'd from M4M
        public static class TrajectoryHelpers
        {
            public static TOut[][] Map2D<TIn, TOut>(this TIn[][] arr, Func<TIn, TOut> map)
            {
                return arr.Select(a => a.Select(map).ToArray()).ToArray();
            }
    
            public static T[][] Transpose<T>(this IReadOnlyList<IReadOnlyList<T>> samples)
            {
                var l = samples.Count;
                var n = samples[0].Count;
    
                var res = new T[n][];
    
                for (int i = 0; i < n; i++)
                    res[i] = new T[l];
    
                for (int i = 0; i < l; i++)
                    for (int j = 0; j < n; j++)
                        res[j][i] = samples[i][j];
    
                return res;
            }
    
            public static void SaveTrajectories(String filename, double[][] trajectories, int samplePeriod)
            {
                using (System.IO.FileStream fs = System.IO.File.Open(filename, System.IO.FileMode.Create))
                {
                    fs.WriteInt32(trajectories.Length); // length
                    fs.WriteInt32(trajectories[0].Length); // width
    
                    for (int i = 0; i < trajectories.Length; i++)
                    {
                        double[] curr = trajectories[i];
    
                        for (int j = 0; j < curr.Length; j++)
                        {
                            fs.WriteFloat64(curr[j]);
                        }
                    }
    
                    if (samplePeriod >= 0)
                    {
                        fs.WriteInt32(-1);
                        fs.WriteInt32(samplePeriod);
                    }
                }
            }
    
            public static double[][] LoadTrajectories(String filename, out int samplePeriod)
            {
                using (System.IO.FileStream fs = System.IO.File.Open(filename, System.IO.FileMode.Open))
                {
                    int length = fs.ReadInt32();
                    int width = fs.ReadInt32();
    
                    double[][] trajectories = new double[length][];
    
                    for (int i = 0; i < trajectories.Length; i++)
                    {
                        double[] curr = trajectories[i] = new double[width];
    
                        for (int j = 0; j < curr.Length; j++)
                        {
                            curr[j] = fs.ReadFloat64();
                        }
                    }
    
                    if (fs.Position < fs.Length)
                    { // there is more
                        int signaller = fs.ReadInt32();
                        samplePeriod = fs.ReadInt32();
                    }
                    else
                    {
                        // (legacy support)
                        samplePeriod = -1;
                    }
    
                    return trajectories;
                }
            }
        }
    
        public static class PlotHelpers
        {
            public static void AddRange<T>(this ElementCollection<T> col, IEnumerable<T> things)
                where T : PlotElement
            {
                foreach (var thing in things)
                    col.Add(thing);
            }
    
            public static void ExportPdf(this PlotModel plot, string filename, double size, bool forceAspect, double relativeYSize = 1)
            {
                PdfExporter exporter = new PdfExporter();
    
                plot.ResetAllAxes();
                plot.InvalidatePlot(true);
                ((IPlotModel)plot).Update(true);
    
                double width = 297 / 25.4 * 72 * size; // A4
                double height = 210 / 25.4 * 72 * size * relativeYSize;
    
                exporter.Width = width;
                exporter.Height = height;
    
                void doPlot()
                {
                    plot.ResetAllAxes();
                    plot.InvalidatePlot(true);
    
                    using (var fs = new System.IO.FileStream(filename, FileMode.Create))
                    {
                        exporter.Export(plot, fs);
                    }
                }
    
                void doForceAspect()
                {
                    if (plot.PlotArea.Width > plot.PlotArea.Height)
                    {
                        width = width - plot.PlotArea.Width + plot.PlotArea.Height;
                    }
                    else
                    {
                        height = height + plot.PlotArea.Width - plot.PlotArea.Height;
                    }
    
                    exporter.Width = width;
                    exporter.Height = height;
                }
    
                doPlot();
    
                if (forceAspect)
                {
                    doForceAspect();
                    doPlot();
                }
            }
    
            public static LinearColorAxis Gray()
            {
                return new LinearColorAxis() { Position = AxisPosition.Right, Palette = OxyPalettes.Gray(100) };
            }
    
            public static PlotModel IntegralCartesianAxes(string title, string x, string y)
            {
                var plot = new PlotModel() { Title = title, PlotType = PlotType.Cartesian };
                var xaxis = new LinearAxis() { Key = "x", Title = x, Position = AxisPosition.Bottom, MinimumMinorStep = 1, MinimumMajorStep = 1 };
                var yaxis = new LinearAxis() { Key = "y", Title = y, Position = AxisPosition.Left, MinimumMinorStep = 1, MinimumMajorStep = 1 };
                plot.Axes.Add(xaxis);
                plot.Axes.Add(yaxis);
                return plot;
            }
    
            public static PlotModel LinearAxes(string title, string x, string y)
            {
                var plot = new PlotModel() { Title = title };
                var xaxis = new LinearAxis() { Key = "x", Title = x, Position = AxisPosition.Bottom, MaximumPadding = 0.02, MinimumPadding = 0.02 };
                var yaxis = new LinearAxis() { Key = "y", Title = y, Position = AxisPosition.Left, MaximumPadding = 0.02, MinimumPadding = 0.02 };
                plot.Axes.Add(xaxis);
                plot.Axes.Add(yaxis);
                return plot;
            }
    
            public static HeatMapSeries PlotMatrix(Matrix<double> matrix)
            {
                double[,] darr = matrix.Transpose().ToArray();
    
                var hms = new HeatMapSeries()
                {
                    Interpolate = false,
                    Data = darr,
                    X0 = 0,
                    X1 = matrix.ColumnCount - 1,
                    Y0 = 0,
                    Y1 = matrix.ColumnCount - 1,
                };
    
                return hms;
            }
    
            public static LineSeries[] PlotTrajectoryLines1D(double[][] trajectories, int samplePeriod, OxyColor c0, OxyColor c1)
            {
                var lines = new LineSeries[trajectories.Length];
    
                for (int i = 0; i < trajectories.Length; i++)
                {
                    var ls = lines[i] = new LineSeries()
                    {
                        Color = OxyColor.Interpolate(c0, c1, i / (double)Math.Max(1, trajectories.Length - 1)),
                    };
    
                    ls.Points.AddRange(trajectories[i].Select((s, x) => new DataPoint(x * samplePeriod, s)));
                }
    
                return lines;
            }
    
            public static ScatterSeries[] PlotTrajectoryScatters1D(double[][] trajectories, int samplePeriod, OxyColor c0, OxyColor c1)
            {
                var scatters = new ScatterSeries[trajectories.Length];
    
                for (int i = 0; i < trajectories.Length; i++)
                {
                    var ss = scatters[i] = new ScatterSeries()
                    {
                        MarkerFill = OxyColor.Interpolate(c0, c1, i / (double)Math.Max(1, trajectories.Length - 1)),
                        MarkerType = MarkerType.Diamond,
                        MarkerSize = 2,
                    };
    
                    ss.Points.AddRange(trajectories[i].Select((s, x) => new ScatterPoint(x * samplePeriod, s)));
                }
    
                return scatters;
            }
    
            public static LineSeries[] PlotTrajectoryLines2D(double[][][] trajectories, int samplePeriod, OxyColor c00, OxyColor c10, OxyColor c01, OxyColor c11)
            {
                var lines = new List<LineSeries>();
    
                for (int i = 0; i < trajectories.Length; i++)
                {
                    for (int j = 0; j < trajectories[i].Length; j++)
                    {
                        var iz = (double)i / Math.Max(1, trajectories.Length - 1);
                        var jz = (double)j / Math.Max(1, trajectories[i].Length - 1);
    
                        var ls = new LineSeries()
                        {
                            Color = OxyColor.Interpolate(OxyColor.Interpolate(c00, c10, iz), OxyColor.Interpolate(c01, c11, iz), jz),
                        };
    
                        lines.Add(ls);
    
                        ls.Points.AddRange(trajectories[i][j].Select((s, x) => new DataPoint(x * samplePeriod, s)));
                    }
                }
    
                return lines.ToArray();
            }
        }
    
        public static class SamplingExtentions
        {
            public static double[][] DownsampleTrajectories(this IReadOnlyList<IReadOnlyList<double>> trajectories, int samplePeriod)
            {
                return trajectories.Select(t => t == null ? null : Downsample(t, samplePeriod).ToArray()).ToArray();
            }
    
            public static IEnumerable<T> Downsample<T>(this IEnumerable<T> sequence, int samplePeriod)
            {
                int c = 0;
    
                foreach (var thing in sequence)
                {
                    if (c == 0)
                    {
                        yield return thing;
                        c = samplePeriod;
                    }
    
                    c--;
                }
            }
    
            public static double[][] MeanDownsampleTrajectories(this IReadOnlyList<IReadOnlyList<double>> trajectories, int resolution)
            {
                return trajectories.Select(t => t == null ? null : MeanDownsample(t, resolution)).ToArray();
            }
    
            public static double[] MeanDownsample(IReadOnlyList<double> trajectory, int resolution)
            {
                double[] res = new double[trajectory.Count / resolution];
    
                for (int i = 0; i < res.Length; i++)
                {
                    double total = 0.0;
                    int c = 0;
                    for (int j = i * resolution; j < i * resolution + resolution; j++)
                    {
                        if (j >= 0 && j < trajectory.Count)
                        {
                            total += trajectory[j];
                            c++;
                        }
                    }
    
                    res[i] = total / c;
                }
    
                return res;
            }
        }
    
        /// <summary>
        /// Provides a random source for a single thread of execution.
        /// </summary>
        public class ExecutionContext
        {
            public ExecutionContext(Random rand)
            {
                Rand = rand ?? throw new ArgumentNullException(nameof(rand));
            }
    
            public Random Rand { get; }
        }
    
        /// <summary>
        /// Implements big-endian (ish) primitive operations for compatibility purposes.
        /// </summary>
        /// <remarks>
        /// Compatible with netstate, but simplified implementation.</remarks>
        public static class StreamExtentions
        {
            /// <summary>
            /// Modifies input.
            /// </summary>
            private static void WriteEndianedBytes(this Stream s, byte[] barr)
            {
                if (BitConverter.IsLittleEndian)
                    Array.Reverse(barr);
                s.Write(barr);
            }
    
            private static byte[] ReadEndianedBytes(this Stream s, int count)
            {
                var barr = new byte[count];
                s.Read(barr);
                if (BitConverter.IsLittleEndian)
                    Array.Reverse(barr);
                return barr;
            }
    
            public static void WriteInt32(this Stream s, int i)
            {
                s.WriteEndianedBytes(BitConverter.GetBytes(i));
            }
    
            public static void WriteFloat64(this Stream s, double d)
            {
                // not sure if doubles have an endianness... but they do now
                s.WriteEndianedBytes(BitConverter.GetBytes(d));
            }
    
            public static int ReadInt32(this Stream s)
            {
                return BitConverter.ToInt32(s.ReadEndianedBytes(sizeof(int)));
            }
    
            public static double ReadFloat64(this Stream s)
            {
                // not sure if doubles have an endianness... but they do now
                return BitConverter.ToDouble(s.ReadEndianedBytes(sizeof(double)));
            }
        }
    
        /// <summary>
        /// Helper methods for random sampling.
        /// </summary>
        public static class RandomExtensions
        {
            /// <summary>
            /// Returns a random number in the range [-maxAbs, maxAbs).
            /// </summary>
            public static double NextSignedUniform(this Random rand, double maxAbs)
            {
                return (rand.NextDouble() * 2.0 - 1.0) * maxAbs;
            }
    
            /// <summary>
            /// Returns true or false with equal probability.
            /// </summary>
            /// <remarks>
            /// Provided for consistency with M4M.
            /// RandomBoolean in MathNet prouduces the opposite result.
            /// </remarks>
            public static bool NextBool(this Random rand, double probability)
            {
                return rand.NextDouble() < probability;
            }
        }
    }