Skip to content
Snippets Groups Projects
Select Git revision
  • 45191b47d838d14cfdbaa8d4123fb1e9e22a1d34
  • master default protected
2 results

Program.cs

Blame
  • Program.cs 54.63 KiB
    using MathNet.Numerics.LinearAlgebra;
    using MathNet.Numerics.Random;
    using OxyPlot;
    using OxyPlot.Axes;
    using OxyPlot.Series;
    using System;
    using System.Collections.Concurrent;
    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 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 from looking at it without taking too long to run.
            static void Main(string[] args)
            {
                // none of our matrices are big enough to warrent multi-threading
                MathNet.Numerics.Control.UseSingleThread();
    
                // non-essential stuff so we can run in parallel if we don't have time to wait
                var tasks = new List<Task>();
                void CreateRun(string outDir, Parameters parameters, int seed, int reportingPeriodMultiplier = 1)
                {
                    tasks.Add(new Task(() => Run(outDir, parameters, seed, reportingPeriodMultiplier)));
                }
    
                // prepare a bunch of experiments with the 'standard' configuration, only varying Z
                // seeds should match those used in the paper
                CreateRun("Z000", Configurations.Standard(Z: 0.0), seed: 42);
                CreateRun("Z050", Configurations.Standard(Z: 0.5), seed: 542);
                CreateRun("Z075", Configurations.Standard(Z: 0.75), seed: 1842);
                CreateRun("Z095", Configurations.Standard(Z: 0.95), seed: 2042);
                CreateRun("Z000PerfectG", Configurations.PerfectG(Z: 0.0), seed: 42);
    
                // prepare some experiments with unusual configurations
                CreateRun("FixedHigh", Configurations.FixedHigh, seed: 42);
                CreateRun("Z050E204", Configurations.Standard204(Z: 0.95), seed: 42, reportingPeriodMultiplier: 10); // this one takes a while
                CreateRun("Z095L2", Configurations.StandardL2(Z: 0.95), seed: 2042);
                CreateRun("Z095MMSO", Configurations.StandardMMSO(Z: 0.95), seed: 2042);
    
                // run the experiments
                bool runInParallel = true;
                if (runInParallel)
                {
                    // parallel: queue all the experiments are once (will be run subject to the thread pool)
                    foreach (var t in tasks)
                        t.Start();
                    Task.WaitAll(tasks.ToArray());
                }
                else
                {
                    // sequential: run the experiments one at a time
                    foreach (var t in tasks)
                        t.RunSynchronously();
                }
            }
    
            /// <summary>
            /// Runs a single experiment with the given parameters and seed into the given output directory
            /// </summary>
            public static void Run(string outDir, Parameters parameters, int seed, int reportingPeriodMultiplier = 1)
            {
                Console.WriteLine($"Running into {outDir}");
                System.IO.Directory.CreateDirectory(outDir);
                string path(string f) => System.IO.Path.Combine(outDir, f);
    
                var ctx = new ExecutionContext(new MersenneTwister(seed, false));
    
                // prepare the recorders
                var recorders = new Recorders(parameters.ModuleSize, parameters.ModuleCount, reportingPeriodMultiplier);
    
                // endOfEpoch reporting actions
                void endOfEpoch(IvmcEnvironment environment, Individual individual, Judgement judgement, int epoch)
                {
                    // inform recorders
                    recorders.Sample(environment, individual, judgement, epoch);
    
                    if (epoch % (reportingPeriodMultiplier * 1000) == 0)
                    {
                        // report progress
                        Console.WriteLine($"{outDir}: {epoch}/{parameters.TotalEpochs} ({(double)epoch / parameters.TotalEpochs:P})\tf={judgement.Fitness}");
    
                        // produce grn plot
                        var grnPlot = PlotHelpers.IntegralCartesianAxes($"Epoch {epoch}", "row", "col");
                        grnPlot.Axes.Add(PlotHelpers.Gray());
                        grnPlot.Series.Add(PlotHelpers.PlotMatrix(individual.B));
                        grnPlot.ExportPdf(path($"epoch{epoch}.pdf"), 0.5, true);
                    }
                }
    
                // run the actual experiment
                Model.RunExperiment(ctx, parameters, endOfEpoch);
    
                // produce output from recorders
                Recording.ProduceOutput(outDir, parameters, recorders);
            }
        }
    
        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,
                };
            }
    
            /// <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,
                    TotalEpochs = 500000,
                    BMutationRate = 0.0002,
                    CostCoefficient = 8,
                    GenerationsPerEpoch = 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 epoch.
            /// </summary>
            public static Parameters PerfectG(double Z)
            {
                return new Parameters()
                {
                    DualPeakedEnvironmentalInstanceFrequency = Z,
                    EpochStartOperation = EpochStartOperations.VaryEnvironmentWithPerfectG,
                    MutateBProbability = 1,
                };
            }
    
            /// <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()
                    {
                        EpochStartOperation = EpochStartOperations.FixEnvironmentAndG,
                        TotalEpochs = 30,
                    };
                }
            }
        }
    
        /*
         * 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 epoch, for reporting purposes.
            /// </summary>
            public delegate void EndOfEpochReporter(IvmcEnvironment environment, Individual individual, Judgement judgement, int epoch);
    
            /// <summary>
            /// Runs a single experiment with the given parameters in the given context, reporting the completion of each epoch.
            /// </summary>
            public static void RunExperiment(ExecutionContext ctx, Parameters parameters, EndOfEpochReporter endOfEpoch)
            {
                // prepare the initial individual
                int N = parameters.ModuleCount * parameters.ModuleSize;
                var individual = Individual.CreateZeroGenome(N);
                individual.DevelopInplace(ctx, parameters);
    
                // report initial state
                endOfEpoch(null, individual, default, 0);
    
                // prepare the environment
                var environment = new IvmcEnvironment(parameters.ModuleSize, parameters.ModuleCount);
    
                // run many epochs
                for (int epoch = 1; epoch <= parameters.TotalEpochs; epoch++)
                {
                    // run one epoch
                    var judgement = Model.RunEpoch(epoch, ctx, parameters, environment, ref individual);
    
                    // report state at end of epoch
                    endOfEpoch(environment, individual, judgement, epoch);
                }
            }
    
            /// <summary>
            /// Judges the given individual in the given environment with the given parameters.
            /// </summary>
            /// <remarks>
            /// Eq2 in the Model section under Evolution
            /// </remarks>
            public static Judgement Judge(ExecutionContext ctx, Parameters parameters, IvmcEnvironment environment, Individual individual)
            {
                var benefit = environment.JudgeBenefit(individual.P);
                var cost = parameters.CostFunction(individual.B);
                var fitness = benefit - cost * parameters.CostCoefficient;
    
                return new Judgement(fitness, benefit, cost);
            }
    
            /// <summary>
            /// Runs an epoch with the given parameters, environment, and individual.
            /// The individual may be modified or replaced with a different instance.
            /// </summary>
            public static Judgement RunEpoch(int epoch, ExecutionContext ctx, Parameters parameters, IvmcEnvironment environment, ref Individual individual)
            {
                parameters.EpochStartOperation(ctx, parameters, individual, environment);
    
                // reset check (compat)
                ctx.Rand.NextDouble();
    
                var judgement = Judge(ctx, parameters, environment, individual);
    
                var candidate = Individual.CreateZeroGenome(individual.Size);
    
                for (int generation = 0; generation < parameters.GenerationsPerEpoch; generation++)
                {
                    Individual.Mutate(ctx, parameters, individual, candidate);
                    var candidateJudgement = Judge(ctx, parameters, environment, candidate);
    
                    if (candidateJudgement.Fitness > judgement.Fitness)
                    {
                        // swap
                        (individual, candidate) = (candidate, individual);
    
                        judgement = candidateJudgement;
                    }
                }
    
                return judgement;
            }
        }
    
        /// <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=120000 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 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 GenerationsPerEpoch { get; set; } = 1000;
            public int K => GenerationsPerEpoch;
    
            public double MutateBProbability { get; set; } = 0.5;
            public double RB => MutateBProbability;
    
            public int TotalEpochs { get; set; } = 120_000;
    
            public EpochStartOperation EpochStartOperation { get; set; } = EpochStartOperations.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 epoch.
        /// See <see cref="EpochStartOperations"/> for standard implementations.
        /// </summary>
        public delegate void EpochStartOperation(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 model 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)
            {
                // borrow temporaries (probably not worth pooling these... but whatever)
                using (var _update = ctx.VectorPool.Borrow(Size, out var update))
                using (var _state = ctx.VectorPool.Borrow(Size, out var state))
                {
                    // t=0
                    G.CopyTo(state);
    
                    // 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))
                        B.Multiply(state, update);
                        state.Multiply(1 - parameters.DecayRate, state);
                        update.MapInplace(parameters.SquashFunction);
                        state.Add(update, state);
                    }
    
                    // normalise to [-1, +1], and place result in our phenotype
                    state.Multiply(parameters.DecayRate, P);
                }
            }
        }
    
        /// <summary>
        /// Represents the benefit (b), cost (c), and combined fitness (f) of an individual.
        /// f = b - λc
        /// </summary>
        public struct Judgement
        {
            public Judgement(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 JudgeBenefit(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
    
                    // Equivalent to Eq5 in the Model section under Environment for the coefficients we use
                    double moduleFitness = acc > 0
                        ? acc * PlusCoefs[i / ModuleSize]
                        : -acc * MinusCoefs[i / ModuleSize];
    
                    benefitAcc += moduleFitness;
                }
    
                return benefitAcc / ModuleCount; // mean fitness of all modules
            }
        }
    
        /// <summary>
        /// Common reset operations functions.
        /// </summary>
        public static class EpochStartOperations
        {
            /// <summary>
            /// Varies the environment every epoch.
            /// </summary>
            public static void VaryEnvironment(ExecutionContext ctx, Parameters parameters, Individual individual, IvmcEnvironment environment)
            {
                environment.Randomise(ctx, parameters);
            }
    
            /// <summary>
            /// Varies the environment every epoch,
            /// 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 (LASSO) 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 (Ridge) 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 void ProduceOutput(string outDir, Parameters parameters, Recorders recorders)
            {
                string path(string f) => System.IO.Path.Combine(outDir, f);
    
                // 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", "Epoch", "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"), 0.5, false);
    
                var hierarchyPlot = PlotHelpers.LinearAxes("Evolution of Hierarchy", "Epoch", "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"), 0.5, false);
    
                var moduleIndependencePlot = PlotHelpers.LinearAxes("Evolution of ModuleIndependence", "Epoch", "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"), 0.5, false);
    
                var fitnessData = recorders.GetFitnessData();
                TrajectoryHelpers.SaveTrajectories(path("fitness_t.dat"), fitnessData, recorders.Judgements.SamplePeriod);
    
                var fitnessPlot = PlotHelpers.LinearAxes("Fitness", "Epoch", "Fitness");
                fitnessPlot.Series.AddRange(PlotHelpers.PlotTrajectoryScatters1D(fitnessData.Take(1).ToArray().DownsampleTrajectories(100), recorders.Judgements.SamplePeriod * 100, OxyColors.Red, OxyColors.DarkOrange));
                fitnessPlot.ExportPdf(path("fitness_t.pdf"), 0.5, false);
    
                var meanFitnessPlot = PlotHelpers.LinearAxes("Mean Fitness", "Epoch", "Mean Fitness");
                fitnessPlot.Series.AddRange(PlotHelpers.PlotTrajectoryScatters1D(fitnessData.Take(1).ToArray().MeanDownsampleTrajectories(100), recorders.Judgements.SamplePeriod * 100, OxyColors.Red, OxyColors.DarkOrange));
                fitnessPlot.ExportPdf(path("meanfitness_t.pdf"), 0.5, false);
    
                var solveData = recorders.GetSolveData().Map2D(x => (double)x);
                TrajectoryHelpers.SaveTrajectories(path("ivmcsolves_t.dat"), solveData, recorders.Solves.SamplePeriod);
    
                var solvePlot = PlotHelpers.LinearAxes("Solves", "Epoch", "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"), 0.5, false);
    
                var switchData = recorders.GetSwitchData().Map2D(x => (double)x);
                TrajectoryHelpers.SaveTrajectories(path("ivmcswitches_t.dat"), switchData, recorders.Switches.SamplePeriod);
    
                var switchPlot = PlotHelpers.LinearAxes("Switches", "Epoch", "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"), 0.5, false);
            }
        }
    
        /// <summary>
        /// Records something over the course of an experiment.
        /// </summary>
        public class Recorder<T>
        {
            public Recorder(Func<IvmcEnvironment, Individual, Judgement, T> sampler, int samplePeriod)
            {
                Sampler = sampler ?? throw new ArgumentNullException(nameof(sampler));
                SamplePeriod = samplePeriod;
            }
    
            public Func<IvmcEnvironment, Individual, Judgement, T> Sampler { get; }
            public int SamplePeriod { get; }
    
            public List<T> Samples { get; } = new List<T>();
    
            public void Sample(IvmcEnvironment environment, Individual individual, Judgement judgement, int epoch)
            {
                if (epoch % SamplePeriod == 0)
                    Samples.Add(Sampler(environment, individual, judgement));
            }
        }
    
        /// <summary>
        /// Counts someting over the course of an experiment.
        /// </summary>
        public class CountRecorder
        {
            public CountRecorder(int sampleSize, Action<IvmcEnvironment, Individual, Judgement, int[]> counter, int samplePeriod)
            {
                Counter = counter;
                SamplePeriod = samplePeriod;
                SampleSize = sampleSize;
            }
    
            public Action<IvmcEnvironment, Individual, Judgement, 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, Judgement judgement, int epoch)
            {
                if (CurrentCount == null)
                {
                    CurrentCount = new int[SampleSize];
                }
    
                Counter(environment, individual, judgement, CurrentCount);
    
                if (epoch % 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)
            {
                Judgements = new Recorder<Judgement>((t, i, j) => j, JudgementSamplePeriod * 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);
            }
    
            private static Action<IvmcEnvironment, Individual, Judgement, int[]> CountSwitches()
            {
                Vector<double> last = null;
    
                void count(IvmcEnvironment environment, Individual individual, Judgement judgement, 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;
            }
    
            private static void CountSolves(IvmcEnvironment environment, Individual individual, Judgement judgement, 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 = 100;
            private const int JudgementSamplePeriod = 1;
            private const int OtherSamplePeriod = 100;
    
            public Recorder<Judgement> Judgements { 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, Judgement judgement, int epoch)
            {
                Judgements.Sample(environment, individual, judgement, epoch);
                Grns.Sample(environment, individual, judgement, epoch);
                Ps.Sample(environment, individual, judgement, epoch);
                Solves.Sample(environment, individual, judgement, epoch);
                Switches.Sample(environment, individual, judgement, epoch);
            }
    
            public double[][] GetFitnessData()
            {
                return new[] {
                    Judgements.Samples.Select(j => j.Fitness).ToArray(),
                    Judgements.Samples.Select(j => j.Benefit).ToArray(),
                    Judgements.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)
            {
                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;
    
                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 };
                var yaxis = new LinearAxis() { Key = "y", Title = y, Position = AxisPosition.Left };
                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 and vector pooling for a single thread of execution.
        /// </summary>
        public class ExecutionContext
        {
            public ExecutionContext(Random rand)
            {
                Rand = rand ?? throw new ArgumentNullException(nameof(rand));
                VectorPool = new Pool<int, Vector<double>>(CreateVector.Dense<double>);
            }
    
            public Random Rand { get; }
            public Pool<int, Vector<double>> VectorPool { 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;
            }
        }
    
        /// <summary>
        /// Provides a pool of resuable objects.
        /// </summary>
        /// <typeparam name="TParams">The parameters that describe an object.</typeparam>
        /// <typeparam name="TResult">The type of objects pooled.</typeparam>
        public class Pool<TParams, TResult> where TParams : IEquatable<TParams>
        {
            public Pool(Func<TParams, TResult> factory)
            {
                Factory = factory ?? throw new ArgumentNullException(nameof(factory));
            }
    
            private Func<TParams, TResult> Factory { get; }
            private ConcurrentDictionary<TParams, ConcurrentBag<TResult>> Table { get; }
                = new ConcurrentDictionary<TParams, ConcurrentBag<TResult>>();
    
            private ConcurrentBag<TResult> FetchBag(TParams p)
            {
                return Table.GetOrAdd(p, _ => new ConcurrentBag<TResult>());
            }
    
            public Borrowed<TParams, TResult> Borrow(TParams p, out TResult res)
            {
                var bag = FetchBag(p);
                if (!bag.TryTake(out res))
                    res = Factory(p);
                return new Borrowed<TParams, TResult>(this, p, res = Factory(p));
            }
    
            public void Release(Borrowed<TParams, TResult> borrowed)
            {
                FetchBag(borrowed.Param).Add(borrowed.Value);
            }
        }
    
        public class Borrowed<TParams, TResult> : IDisposable where TParams : IEquatable<TParams>
        {
            public Borrowed(Pool<TParams, TResult> pool, TParams param, TResult value)
            {
                Pool = pool ?? throw new ArgumentNullException(nameof(pool));
                Param = param;
                Value = value;
            }
    
            private Pool<TParams, TResult> Pool;
            public TParams Param { get; private set; }
            public TResult Value { get; private set; }
    
            public void Dispose()
            {
                Pool?.Release(this);
                Pool = null;
                Param = default(TParams);
                Value = default(TResult);
            }
        }
    }