Stats (#8)

* refactoring: towards dependency injection and IOC

* completed extraction of simulation's business logic to Simulation object

* create Benchmark object

* working fitness plotting service

* use newborns as timeline; add individuals distance plotter

* legend for messy plot

* Abstract stat service

* Add benchmark option for parallel execution

* working aggregate data gathered at the end of parallel simulations execution

.gitignore

@@ -1,4 +1,5 @@
 venv/
 __pycache__/
 .DS_Store
-*.svg
+*.svg
+*.mat

src/config.yaml

@@ -11,7 +11,7 @@ genetic_algo_tuning:
   population_size: 300
   turnover_rate: 0.5  # rate of new individuals generated from crossover
   iterations_number: 400
-  mutation_rate: 0.5  # rate of individuals that will be mutated at each iteration
+  mutation_rate: 0.1  # rate of individuals that will be mutated at each iteration
   cut_points: 1
   use_niches: true
   world_width: 20

src/core/niche.py

@@ -8,7 +8,7 @@ class NichePopulation(Population):
     def __init__(self):
         super().__init__(Config.GeneticAlgoTuning.worldHeight * Config.GeneticAlgoTuning.worldWidth)
 
-        self.world: np.ndarray = np.array(self.population).reshape(
+        self.world: np.ndarray = np.array(self.individuals).reshape(
             Config.GeneticAlgoTuning.worldHeight,
             Config.GeneticAlgoTuning.worldWidth
         )

src/core/population.py

@@ -9,8 +9,9 @@
 
 class Population:
   def __init__(self, pop_size: int = Config.GeneticAlgoTuning.populationSize):
-    self.population = [Gene() for _ in range(pop_size)]
+    self.individuals = [Gene() for _ in range(pop_size)]
     self.generationNumber = 0
+    self.newbornsCounter = 0
     self.king = Gene()
 
   def extractParent(self) -> Gene:
@@ -24,17 +25,17 @@ def extractParentIndex(self) -> Gene:
     Extracts a random parent from the population, with no regard to gene fitness.
     """
 
-    return choice(self.population)
+    return choice(self.individuals)
 
   def extractParentFitness(self) -> Gene:
     """
     Extracts a parent from the population. The likelihood of extraction is proportional to gene fitness.
     """
 
-    fitness = [g.fitness() for g in self.population]  # Can be cached, thus optimized
+    fitness = [g.fitness() for g in self.individuals]  # Can be cached, thus optimized
     fitness = [f + 500 if f > float("-inf") else 0 for f in fitness] # Increment trick
 
-    return choices(self.population, weights=fitness)[0]
+    return choices(self.individuals, weights=fitness)[0]
 
 
   def selectParents(self) -> List[Tuple[Gene]]:
@@ -43,55 +44,57 @@ def selectParents(self) -> List[Tuple[Gene]]:
     half the dimension of the current population
     """
 
-    pop_size = len(self.population)
+    pop_size = len(self.individuals)
     parents_num = pop_size // 2
 
     parents = [
       (self.extractParent(pop_size), self.extractParent(pop_size)) for _ in range(parents_num)
     ]
 
 
-  def generations(self) -> List[Gene]:
+  def generations(self) -> Tuple[List[Gene], int]:
     for _ in range(Config.GeneticAlgoTuning.iterationsNumber):
       self.generateOffspring()
       self.mutate()
       self.fight()
 
-      self.fitnessMean = np.mean([g.fitnessCached for g in self.population])
-      self.fitnessStdDev = np.std([g.fitnessCached for g in self.population])
+      self.fitnessMean = np.mean([g.fitnessCached for g in self.individuals])
+      self.fitnessStdDev = np.std([g.fitnessCached for g in self.individuals])
       logging.info(
         f"\nFitness:\n"
         f"\tMean: {self.fitnessMean:.4f}\n"
         f"\tSd: {self.fitnessStdDev:.4f}\n"
-        f"Population size: {len(self.population)}"
+        f"Population size: {len(self.individuals)}"
       )
 
       self.king = \
-        self.population[0] if self.population[0].fitnessCached > self.king.fitnessCached else self.king
+        self.individuals[0] if self.individuals[0].fitnessCached > self.king.fitnessCached else self.king
 
       if self.fitnessStdDev <= np.finfo(np.float32).eps:
         return
 
       self.generationNumber += 1
-      yield self.population, self.generationNumber
+      yield self.individuals, self.generationNumber
 
 
   def fight(self):
     """This step filters out non-valid individuals and
     keeps only some of them according to the turnover rate
     """
-    oldGenerationSize = len(self.population)
+    oldGenerationSize = len(self.individuals)
 
-    self.population = list(
+    self.individuals = list(
       filter(
       lambda x: x.isValid(),
-      self.population
+      self.individuals
     ))
 
-    logging.warning(f"Killed {oldGenerationSize - len(self.population)} ({(oldGenerationSize - len(self.population)) / oldGenerationSize * 100:.1f}%) genes")
+    self.killedGenes = oldGenerationSize - len(self.individuals)
+    self.killedGenesRatio = self.killedGenes / oldGenerationSize * 100
+    logging.warning(f"Killed {self.killedGenes} ({self.killedGenesRatio:.1f}%) genes")
 
     survivedGenesNumber = ceil(Config.GeneticAlgoTuning.turnoverRate * Config.GeneticAlgoTuning.populationSize)
-    self.population = sorted(self.population, reverse=True)[ : survivedGenesNumber]
+    self.individuals = sorted(self.individuals, reverse=True)[ : survivedGenesNumber]
 
   def crossover(self, mother: Gene, father: Gene):
     cutpointIdx = randrange(Config.GeneEncoding.segmentsNumber)
@@ -108,20 +111,22 @@ def crossover(self, mother: Gene, father: Gene):
 
   def generateOffspring(self):
     newGenerationSize = floor((1.0 - Config.GeneticAlgoTuning.turnoverRate) * Config.GeneticAlgoTuning.populationSize)
-    oldGenerationSize = len(self.population)
+    oldGenerationSize = len(self.individuals)
 
     for _ in range(newGenerationSize // 2):
       momGene = self.extractParent()
       dadGene = self.extractParent()
 
       newGene1, newGene2 = self.crossover(momGene, dadGene)
 
-      self.population.append(newGene1)
-      self.population.append(newGene2)
+      self.individuals.append(newGene1)
+      self.individuals.append(newGene2)
+
+      self.newbornsCounter += 2
 
   def mutate(self):
-    toMutateSize = ceil(Config.GeneticAlgoTuning.mutationRate * len(self.population))
-    genesToMutate = sample(self.population, k = toMutateSize)
+    toMutateSize = ceil(Config.GeneticAlgoTuning.mutationRate * len(self.individuals))
+    genesToMutate = sample(self.individuals, k = toMutateSize)
 
     for gene in genesToMutate:
       mutationAngles = (np.random.rand(Config.GeneEncoding.segmentsNumber) - 0.5) \

src/core/simulation.py

@@ -0,0 +1,43 @@
+import logging
+from typing import List
+from services.plotters import *
+from services.persistence import *
+from services.statistics import *
+
+class Simulation:
+  def __init__(self, population: Population):
+    self.population = population
+    self.plotterServices: List[IPlotterService] = list()
+    self.persistenceServices: List[IPersistenceService] = list()
+    self.statServices: List[IStatService] = list()
+
+  def withService(self, service: Service) -> Any:
+    if isinstance(service, IPlotterService):
+      target = self.plotterServices
+    elif isinstance(service, IPersistenceService):
+      target = self.persistenceServices
+    elif isinstance(service, IStatService):
+      target = self.statServices
+
+    target.append(service)
+
+    return self
+
+  def runServices(self) -> None:
+    for plotter in self.plotterServices:
+      plotter.plot(self.population)
+
+    for saver in self.persistenceServices:
+      saver.save(self.population)
+
+    for stater in self.statServices:
+      stater.stat(self.population)
+
+  def run(self, *_) -> None:
+    generation, epoch = next(self.population.generations())
+
+    logging.info(f"Epoch: {epoch}")
+    logging.debug(generation)
+    logging.info(f"Best gene (fitness={generation[0].fitness():.2f}):\n{generation[0]}")
+
+    self.runServices()

src/poc.py

@@ -3,56 +3,23 @@
 from typing import Callable
 import matplotlib.pyplot as plt
 import matplotlib.animation as animation
+from rf.radiation import RadiationPattern
+from services.plotters import *
+from services.persistence import *
+from services.statistics import *
 from utils.amenities import plotPathAndRad, saveSvg
 from core.config import Config
 from core.population import Population
+from core.simulation import Simulation
+from multiprocessing import Pool
 from functools import partial
+from scipy.io import savemat
 
 CONFIG_FILENAME = "config.yaml"
 
-def simulationStep(
-  pop: Population,
-  doPlot: bool,
-  *_,
-  **kwargs) -> None:
-  try:
-    generation, epoch = next(pop.generations())
-  except StopIteration:
-    quit()
-
-  logging.info(f"Epoch: {epoch}")
-  logging.debug(generation)
-  logging.info(f"Best gene (fitness={generation[0].fitness():.2f}):\n{generation[0]}")
-
-  if doPlot:
-    plotPathAndRad(
-      title = f"Epoch: {epoch} -\n"
-        f"Best fitness: {generation[0].fitness():.2f} -\n"
-        f"Mean fitness: {pop.fitnessMean:.2f} -\n"
-        f"Sd fitness: {pop.fitnessStdDev:.2f}",
-      polychain = generation[0].getCartesianCoords(),  # Plot only best performing individual
-      radiationSagittal = generation[0].getRadiationPatternSagittal(),
-      radiationFrontal = generation[0].getRadiationPatternFrontal(),
-      groundPlaneDistance = generation[0].groundPlaneDistance,
-      axes = (kwargs.pop("shapeAxes"), kwargs.pop("radiationAxesSag"), kwargs.pop("radiationAxesFront"))
-    )
-
-  outputDirectory = kwargs.pop("outputDirectory")
-  if outputDirectory is not None:
-    with open(os.path.join(outputDirectory, f"gen{epoch}.svg"), "w") as outFile:
-      saveSvg(outFile, generation, kwargs.pop("withBoundaries"))
-
-
-def buildSimulation(doPlot: bool, *_, **kwargs) -> Callable[[Population, bool], None]:
-  pop = Population()
-
+def main(doPlot: bool, outdir: str, withBoundaries: bool, statService: IStatService, instanceNumber: int = 0):
   signal.signal(signal.SIGINT, lambda *_: quit())
 
-  return partial(simulationStep, pop, doPlot, **kwargs)
-
-
-def main(doPlot: bool, outdir: str, withBoundaries: bool):
-
   logging.basicConfig(
     level=logging.INFO,
     format="%(asctime)s - %(levelname)s: %(message)s",
@@ -61,27 +28,53 @@ def main(doPlot: bool, outdir: str, withBoundaries: bool):
 
   with open(CONFIG_FILENAME, "r") as f:
     Config.loadYaml(f)
-
-  if doPlot:
-    fig = plt.figure()
-    shape = fig.add_subplot(1, 3, 1)
-    radPatternSag = fig.add_subplot(1, 3, 2, projection='polar')
-    radPatternFront = fig.add_subplot(1, 3, 3, projection='polar')
-    simulation = buildSimulation(
-      doPlot,
-      shapeAxes = shape,
-      radiationAxesSag = radPatternSag,
-      radiationAxesFront = radPatternFront,
-      outputDirectory = outdir,
-      withBoundaries = withBoundaries
-    )
-    anim = animation.FuncAnimation(fig, simulation, interval=10, cache_frame_data=False)
-    plt.show()
-
+
+  if outdir is None:
+    persistenceServiceClass = StubPlotterService
   else:
-    simulation = buildSimulation(doPlot, outputDirectory = outdir, withBoundaries = withBoundaries)
-    while True:
-      simulation()
+    if withBoundaries:
+      persistenceServiceClass = MiniatureWithBoundariesPersistenceService
+    else:
+      persistenceServiceClass = MiniaturePersistenceService
+
+  persistenceService = persistenceServiceClass(outdir)
+
+  PLOT_ROWS = 2
+  PLOT_COLS = 3
+  fig = plt.figure(f"Simulation {instanceNumber}")
+  shape = fig.add_subplot(PLOT_ROWS, PLOT_COLS, 1)
+  radPatternSag = fig.add_subplot(PLOT_ROWS, PLOT_COLS, 2, projection='polar')
+  radPatternFront = fig.add_subplot(PLOT_ROWS, PLOT_COLS, 3, projection='polar')
+  fitnessGraph = fig.add_subplot(PLOT_ROWS, PLOT_COLS, 4)
+  killedGraph = fig.add_subplot(PLOT_ROWS, PLOT_COLS, 5)
+  distanceGraph = fig.add_subplot(PLOT_ROWS, PLOT_COLS, 6)
+  fig.tight_layout()
+
+  statService.withGraphers(
+    FitnessPlotter(fitnessGraph),
+    KilledGenesPlotter(killedGraph),
+    EuclideanDistancePlotter(distanceGraph)
+  )
+
+  pop = Population()
+  sim = Simulation(pop) \
+    .withService(PlanarShapePlotter(shape)) \
+    .withService(RadiationPatternPlotter(radPatternFront, Gene.getRadiationPatternFrontal)) \
+    .withService(RadiationPatternPlotter(radPatternSag, Gene.getRadiationPatternSagittal)) \
+    .withService(persistenceService) \
+    .withService(statService)
+
+  try:
+    if doPlot:
+      anim = animation.FuncAnimation(fig, sim.run, interval = 10, cache_frame_data = False)
+      plt.show()
+    else:
+      while True:
+        sim.run()
+  except StopIteration:
+    return statService.valuesDict
+
+  return statService.valuesDict
 
 
 if __name__ == "__main__":
@@ -104,6 +97,30 @@ def main(doPlot: bool, outdir: str, withBoundaries: bool):
     default=False, action="store_true"
   )
 
+  parser.add_argument(
+    "-bm", "--benchmark-instances", help="Number of benchmark's paralell simulations",
+    type=int, default=1
+  )
+
   args = parser.parse_args()
 
-  main(args.plot, args.outdir, args.with_boundaries)
+  statServices = [StatService(join("results", f"stats{i}.mat")) for i in range(args.benchmark_instances)]
+
+  parallelMain = partial(main, args.plot, args.outdir, args.with_boundaries)
+  with Pool(args.benchmark_instances) as p:
+    statsDicts = p.starmap(parallelMain, [(statService, i) for i, statService in enumerate(statServices)])
+
+  outStats = {}
+  minLength = min([len(d['timeline']) for d in statsDicts])
+  for statName in statsDicts[0].keys():
+    rawValues = None
+
+    for d in statsDicts:
+      if rawValues is not None:
+        rawValues = np.vstack((rawValues, d[statName][:minLength]))
+      else:
+        rawValues = np.array(d[statName][:minLength])
+
+    outStats[statName] = np.mean(rawValues, axis=0)
+
+  savemat(join("results", "aggregate_stats.mat"), outStats)

src/services/__init__.py

@@ -0,0 +1 @@
+__all__ = ["persistence", "plotters"]