Select Git revision
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);
}
}
}