Index: /framspy/README.txt
===================================================================
--- /framspy/README.txt	(revision 1190)
+++ /framspy/README.txt	(revision 1190)
@@ -0,0 +1,41 @@
+______________________________ "framspy" directory _______________________________
+
+This directory contains sources that allow one to interact with the native
+Framsticks library from Python, and to perform evolutionary optimization.
+
+"frams.py" is the most fundamental file. It allows to access Framsticks script
+objects in Python. You will find detailed descriptions in the source of this file.
+
+"frams-test.py" uses "frams.py" and tests the connection between Python
+and the native Framsticks library by performing a number of simple yet diversified
+operations,  so you should run it first to ensure everything works correctly. 
+Again, read the comments in the source for more information.
+
+"FramsticksLib.py" uses "frams.py" and provides a few fundamental building blocks
+for optimization: functions to mutate, crossover, evaluate a solution, etc.
+
+
+There are three independent implementations of evolutionary algorithms in this directory
+and they all rely on "FramsticksLib.py":
+
+1)
+"FramsticksEvolution.py" - the simplest one, everything is in one file. Uses the DEAP
+framework ( https://deap.readthedocs.io/en/master/ ) to perform the optimization.
+Supports multiple criteria, the generational architecture, and numerous parameters
+to customize the evolutionary search. Invalid genotypes are assigned the fitness of "-1"
+and the selection process is responsible for their extinction.
+See "run-deap-examples.cmd".
+
+2)
+The "evolalg_steps" module - the most complicated one, every piece of an algorithm
+is considered a "step" and can be customized. This results in a large number of small
+files (classes), but allows for the highest flexibility in building new algorithms.
+This module is no longer supported, but basic usage scenarios should work correctly.
+See "run-evolalg_steps-examples.cmd".
+
+3)
+The "evolalg" module - the middle ground. Minimalistic yet it allows for easy
+customizations, supports the generational architecture, the island model, quality-diversity
+techniques and numerical optimization (completely independent from Framsticks).
+Invalid genotypes are not accepted, so new ones are generated until they are valid.
+See "run-evolalg-examples.cmd".
Index: /framspy/evolalg/base/experiment_abc.py
===================================================================
--- /framspy/evolalg/base/experiment_abc.py	(revision 1190)
+++ /framspy/evolalg/base/experiment_abc.py	(revision 1190)
@@ -0,0 +1,209 @@
+import argparse
+import json
+import os
+import pickle
+import time
+from abc import ABC, abstractmethod
+
+from ..base.random_sequence_index import RandomIndexSequence
+from ..constants import BAD_FITNESS
+from ..json_encoders import Encoder
+from ..structures.hall_of_fame import HallOfFame
+from ..structures.individual import Individual
+from ..structures.population import PopulationStructures
+
+
+class ExperimentABC(ABC):
+
+    population_structures = None
+    hof = []
+    stats = []
+    current_generation = 0
+    time_elapsed = 0
+    
+
+    def __init__(self, popsize, hof_size, save_only_best=True) -> None:
+        self.hof = HallOfFame(hof_size)
+        self.popsize = popsize
+        self.save_only_best = save_only_best
+
+    def select(self, individuals, tournament_size, random_index_sequence):
+        """Tournament selection, returns the index of the best individual from those taking part in the tournament"""
+        best_index = None
+        for i in range(tournament_size):
+            rnd_index = random_index_sequence.getNext()
+            if best_index is None or individuals[rnd_index].fitness > best_index.fitness:
+                best_index = individuals[rnd_index]
+        return best_index
+
+    def addGenotypeIfValid(self, ind_list, genotype):
+        new_individual = Individual()
+        new_individual.set_and_evaluate(genotype, self.evaluate)
+        if new_individual.fitness is not BAD_FITNESS:  # this is how we defined BAD_FITNESS in evaluate()
+            ind_list.append(new_individual)
+
+    def make_new_population(self, individuals, prob_mut, prob_xov, tournament_size):
+        """'individuals' is the input population (a list of individuals).
+        Assumptions: all genotypes in 'individuals' are valid and evaluated (have fitness set).
+        Returns: a new population of the same size as 'individuals' with prob_mut mutants, prob_xov offspring, and the remainder of clones."""
+
+        newpop = []
+        N = len(individuals)
+        expected_mut = int(N * prob_mut)
+        expected_xov = int(N * prob_xov)
+        assert expected_mut + \
+            expected_xov <= N, "If probabilities of mutation (%g) and crossover (%g) added together exceed 1.0, then the population would grow every generation..." % (
+                prob_mut, prob_xov)
+        ris = RandomIndexSequence(N)
+
+        # adding valid mutants of selected individuals...
+        while len(newpop) < expected_mut:
+            ind = self.select(
+                individuals, tournament_size=tournament_size, random_index_sequence=ris)
+            self.addGenotypeIfValid(newpop, self.mutate(ind.genotype))
+
+        # adding valid crossovers of selected individuals...
+        while len(newpop) < expected_mut + expected_xov:
+            ind1 = self.select(
+                individuals, tournament_size=tournament_size, random_index_sequence=ris)
+            ind2 = self.select(
+                individuals, tournament_size=tournament_size, random_index_sequence=ris)
+            self.addGenotypeIfValid(
+                newpop, self.cross_over(ind1.genotype, ind2.genotype))
+
+        # select clones to fill up the new population until we reach the same size as the input population
+        while len(newpop) < len(individuals):
+            ind = self.select(
+                individuals, tournament_size=tournament_size, random_index_sequence=ris)
+            newpop.append(Individual().copyFrom(ind))
+
+        return newpop
+
+    def save_state(self, state_filename):
+        state = self.get_state()
+        if state_filename is None:
+            return
+        state_filename_tmp = state_filename + ".tmp"
+        try:
+            with open(state_filename_tmp, "wb") as f:
+                pickle.dump(state, f)
+            # ensures the new file was first saved OK (e.g. enough free space on device), then replace
+            os.replace(state_filename_tmp, state_filename)
+        except Exception as ex:
+            raise RuntimeError("Failed to save evolution state '%s' (because: %s). This does not prevent the experiment from continuing, but let's stop here to fix the problem with saving state files." % (
+                state_filename_tmp, ex))
+
+    def load_state(self, state_filename):
+        if state_filename is None:
+            # print("Loading evolution state: file name not provided")
+            return None
+        try:
+            with open(state_filename, 'rb') as f:
+                state = pickle.load(f)
+                self.set_state(state)
+        except FileNotFoundError:
+            return None
+        print("...Loaded evolution state from '%s'" % state_filename)
+        return True
+
+    def get_state_filename(self, save_file_name):
+        return None if save_file_name is None else save_file_name + '_state.pkl'
+
+    def get_state(self):
+        return [self.time_elapsed, self.current_generation, self.population_structures, self.hof, self.stats]
+
+    def set_state(self, state):
+        self.time_elapsed, self.current_generation, self.population_structures, hof_, self.stats = state
+        # sorting: ensure that we add from worst to best so all individuals are added to HOF
+        for h in sorted(hof_, key=lambda x: x.rawfitness):
+            self.hof.add(h)
+
+    def update_stats(self, generation, all_individuals):
+        worst = min(all_individuals, key=lambda item: item.rawfitness)
+        best = max(all_individuals, key=lambda item: item.rawfitness)
+        # instead of single best, could add all individuals in population here, but then the outcome would depend on the order of adding
+        self.hof.add(best)
+        self.stats.append(
+            best.rawfitness if self.save_only_best else best)
+        print("%d\t%d\t%g\t%g" % (generation, len(
+            all_individuals), worst.rawfitness, best.rawfitness))
+
+    def initialize_evolution(self, initialgenotype):
+        self.current_generation = 0
+        self.time_elapsed = 0
+        self.stats = []  # stores the best individuals, one from each generation
+        initial_individual = Individual()
+        initial_individual.set_and_evaluate(initialgenotype, self.evaluate)
+        self.hof.add(initial_individual)
+        self.stats.append(
+            initial_individual.rawfitness if self.save_only_best else initial_individual)
+        self.population_structures = PopulationStructures(
+            initial_individual=initial_individual, archive_size=0, popsize=self.popsize)
+
+    def evolve(self, hof_savefile, generations, initialgenotype, pmut, pxov, tournament_size):
+        file_name = self.get_state_filename(hof_savefile)
+        state = self.load_state(file_name)
+        if state is not None:  # loaded state from file
+            # saved generation has been completed, start with the next one
+            self.current_generation += 1
+            print("...Resuming from saved state: population size = %d, hof size = %d, stats size = %d, archive size = %d, generation = %d/%d" % (len(self.population_structures.population), len(self.hof),
+                                                                                                                                                 len(self.stats),  (len(self.population_structures.archive)), self.current_generation, generations))  # self.current_generation (and g) are 0-based, parsed_args.generations is 1-based
+
+        else:
+            self.initialize_evolution(initialgenotype)
+        time0 = time.process_time()
+        for g in range(self.current_generation, generations):
+            self.population_structures.population = self.make_new_population(
+                self.population_structures.population, pmut, pxov, tournament_size)
+            self.update_stats(g, self.population_structures.population)
+            if hof_savefile is not None:
+                self.current_generation = g
+                self.time_elapsed += time.process_time() - time0
+                self.save_state(file_name)
+        if hof_savefile is not None:
+            self.save_genotypes(hof_savefile)
+        return self.population_structures.population, self.stats
+
+    @abstractmethod
+    def mutate(self, gen1):
+        pass
+
+    @abstractmethod
+    def cross_over(self, gen1, gen2):
+        pass
+
+    @abstractmethod
+    def evaluate(self, genotype):
+        pass
+
+    def save_genotypes(self, filename):
+        """Implement if you want to save finall genotypes,in default implementation this function is run once at the end of evolution"""
+        state_to_save = {
+            "number_of_generations": self.current_generation,
+            "hof": [{"genotype": individual.genotype,
+                     "fitness": individual.rawfitness} for individual in self.hof.hof]}
+        with open(f"{filename}.json", 'w') as f:
+            json.dump(state_to_save, f, cls=Encoder)
+
+    
+    @staticmethod
+    def get_args_for_parser():
+        parser = argparse.ArgumentParser()
+        parser.add_argument('-popsize',type= int, default= 50,
+                            help="Population size, default: 50.")
+        parser.add_argument('-generations',type= int, default= 5,
+                                help="Number of generations, default: 5.")
+        parser.add_argument('-tournament',type= int, default= 5,
+                            help="Tournament size, default: 5.")
+        parser.add_argument('-pmut',type= float, default= 0.7,
+                        help="Probability of mutation, default: 0.7")
+        parser.add_argument('-pxov',type= float, default= 0.2,
+                        help="Probability of crossover, default: 0.2")
+        parser.add_argument('-hof_size',type= int, default= 10,
+                            help="Number of genotypes in Hall of Fame. Default: 10.")
+        parser.add_argument('-hof_savefile',type= str, required= False,
+                                help= 'If set, Hall of Fame will be saved in Framsticks file format (recommended extension *.gen. This also activates saving state (checpoint} file and auto-resuming from the saved state, if this file exists.')
+        parser.add_argument('-save_only_best',type= bool, default= True, required= False,
+                            help="")
+        
+        return parser
Index: /framspy/evolalg/base/experiment_islands_model_abc.py
===================================================================
--- /framspy/evolalg/base/experiment_islands_model_abc.py	(revision 1190)
+++ /framspy/evolalg/base/experiment_islands_model_abc.py	(revision 1190)
@@ -0,0 +1,99 @@
+import time
+from abc import ABC
+from typing import List
+
+from ..structures.individual import Individual
+from ..structures.population import PopulationStructures
+from .experiment_abc import ExperimentABC
+
+
+class ExperimentIslands(ExperimentABC, ABC):
+
+    number_of_populations = 5
+    popsize = 100
+    populations: List[PopulationStructures] = []
+    migration_interval = 10
+
+    def __init__(self, popsize, hof_size, number_of_populations, migration_interval, save_only_best) -> None:
+        super().__init__(popsize=popsize, hof_size=hof_size, save_only_best=save_only_best)
+        self.number_of_populations=number_of_populations
+        self.migration_interval=migration_interval
+
+    def migrate_populations(self):
+        print("Performing base migration")
+        pool_of_all_individuals = []
+        for p in self.populations:
+            pool_of_all_individuals.extend(p.population)
+        print(f"Pool of individuals: {len(pool_of_all_individuals)}")
+        sorted_individuals = sorted(
+            pool_of_all_individuals, key=lambda x: x.rawfitness)
+        print(f"Best indiviudal for new islands:")
+        for i in range(self.number_of_populations):
+            shift = i*self.popsize
+            self.populations[i].population = sorted_individuals[shift:shift+self.popsize]
+            print(i, self.populations[i].population[-1].rawfitness)
+
+    def initialize_evolution(self, initialgenotype):
+        self.current_generation = 0
+        self.time_elapsed = 0
+        self.stats = []  # stores the best individuals, one from each generation
+        initial_individual = Individual()
+        initial_individual.set_and_evaluate(initialgenotype, self.evaluate)
+        self.stats.append(initial_individual.rawfitness)
+        [self.populations.append(PopulationStructures(initial_individual=initial_individual,
+                                                      popsize=self.popsize))
+         for _ in range(self.number_of_populations)]
+
+    def get_state(self):
+        return [self.time_elapsed, self.current_generation, self.populations, self.hof, self.stats]
+
+    def set_state(self, state):
+        self.time_elapsed, self.current_generation, self.populations, hof_, self.stats = state
+        # sorting: ensure that we add from worst to best so all individuals are added to HOF
+        for h in sorted(hof_, key=lambda x: x.rawfitness):
+            self.hof.add(h)
+
+    def evolve(self, hof_savefile, generations, initialgenotype, pmut, pxov, tournament_size):
+        file_name = self.get_state_filename(hof_savefile)
+        state = self.load_state(file_name)
+        if state is not None:  # loaded state from file
+            # saved generation has been completed, start with the next one
+            self.current_generation += 1
+            print("...Resuming from saved state: population size = %d, hof size = %d, stats size = %d, generation = %d/%d" % (len(self.populations[0].population), len(
+                self.hof), len(self.stats), self.current_generation, generations))  # self.current_generation (and g) are 0-based, parsed_args.generations is 1-based
+        else:
+            self.initialize_evolution(initialgenotype)
+        time0 = time.process_time()
+        for g in range(self.current_generation, generations):
+            for p in self.populations:
+                p.population = self.make_new_population(
+                    p.population, pmut, pxov, tournament_size)
+
+            if g % self.migration_interval == 0:
+                print("---------Start of migration-------")
+                self.migrate_populations()
+                print("---------End of migration---------")
+
+            pool_of_all_individuals = []
+            [pool_of_all_individuals.extend(p.population)
+             for p in self.populations]
+            self.update_stats(g, pool_of_all_individuals)
+            if hof_savefile is not None:
+                self.current_generation = g
+                self.time_elapsed += time.process_time() - time0
+                self.save_state(file_name)
+
+        if hof_savefile is not None:
+            self.save_genotypes(hof_savefile)
+
+        return self.hof, self.stats
+
+    @staticmethod
+    def get_args_for_parser():
+        parser = ExperimentABC.get_args_for_parser()
+
+        parser.add_argument("-islands",type=int, default=5,
+                            help="Number of subpopulations (islands)")
+        parser.add_argument("-generations_migration",type=int, default=10,
+                            help="Number of generations separating migration events when genotypes migrate between subpopulations (islands)")
+        return parser
Index: /framspy/evolalg/base/experiment_niching_abc.py
===================================================================
--- /framspy/evolalg/base/experiment_niching_abc.py	(revision 1190)
+++ /framspy/evolalg/base/experiment_niching_abc.py	(revision 1190)
@@ -0,0 +1,201 @@
+import time
+from abc import ABC, abstractmethod
+from tkinter import W
+
+import numpy as np
+from deap import base, tools
+from deap.tools.emo import assignCrowdingDist
+
+from ..constants import BAD_FITNESS
+from ..structures.individual import Individual
+from .experiment_abc import ExperimentABC
+from .remove_diagonal import remove_diagonal
+
+
+class DeapFitness(base.Fitness):
+    weights = (1, 1)
+
+    def __init__(self, *args, **kwargs):
+        super(DeapFitness, self).__init__(*args, **kwargs)
+
+
+class ExperimentNiching(ExperimentABC, ABC):
+    fit: str = "niching"
+    normalize: str = "None"
+    archive_size: int = None
+
+    def __init__(self, fit, normalize, popsize, hof_size, save_only_best=True, knn_niching=5, knn_nslc=10, archive_size=0) -> None:
+        ExperimentABC.__init__(self,popsize=popsize, hof_size=hof_size, save_only_best=save_only_best)
+        self.fit = fit
+        self.knn_niching = knn_niching
+        self.knn_nslc = knn_nslc
+        self.normalize = normalize
+        self.archive_size=archive_size
+        if popsize < self.knn_niching:
+            self.knn_niching = popsize - 2
+        if popsize < self.knn_nslc:
+            self.knn_nslc = popsize - 2
+
+    def transform_indexes(self, i, index_array):
+        return [x+1 if x >= i else x for x in index_array]
+
+    def normalize_dissim(self, dissim_matrix):
+        dissim_matrix = remove_diagonal(np.array(dissim_matrix))
+        if self.normalize == "none":
+            return dissim_matrix
+        elif self.normalize == "max":
+            divide_by = np.max(dissim_matrix)
+        elif self.normalize == "sum":
+            divide_by = np.sum(dissim_matrix)
+        else:
+            raise Exception(f"Wrong normalization method,", self.normalize)
+        if divide_by != 0:
+            return dissim_matrix/divide_by
+        else:
+            return dissim_matrix
+
+    def do_niching(self, population_structures):
+        population_archive = population_structures.population + population_structures.archive
+        dissim_matrix = self.dissimilarity(population_archive)
+        if "knn" not in self.fit:
+            dissim_list = np.mean(self.normalize_dissim(dissim_matrix), axis=1)
+        else:
+            dissim_list = np.mean(np.partition(
+                self.normalize_dissim(dissim_matrix), self.knn_niching)[:, :self.knn_niching], axis=1)
+
+        if "niching" in self.fit:
+            for i, d in zip(population_archive, dissim_list):
+                i.fitness = i.rawfitness * d
+        elif "novelty" in self.fit:
+            for i, d in zip(population_archive, dissim_list):
+                i.fitness = d
+        else:
+            raise Exception("Wrong fit type: ", self.fit,
+                            f" chose correct one or implement new behaviour")
+        population_structures.update_archive(dissim_matrix, population_archive)
+
+    def do_nsga2_dissim(self, population):
+        dissim_matrix = self.dissimilarity(population)
+        dissim_list = np.mean(self.normalize_dissim(dissim_matrix), axis=1)
+        for i, d in zip(population, dissim_list):
+            i.fitness = DeapFitness(tuple((d, i.rawfitness)))
+
+    def do_nslc_dissim(self, population):
+        dissim_matrix = self.dissimilarity(population)
+        normalized_matrix = self.normalize_dissim(dissim_matrix)
+        for i in range(len(normalized_matrix)):
+            temp_dissim = normalized_matrix[i]
+            index_array = np.argpartition(temp_dissim, kth=self.knn_nslc, axis=-1)[:self.knn_nslc]
+            dissim_value = np.mean(np.take_along_axis(
+                temp_dissim, index_array, axis=-1))
+            temp_fitness = population[i].rawfitness
+            population_of_most_different = list(
+                map(population.__getitem__, self.transform_indexes(i, index_array)))
+            temp_ind_fit = sum(
+                [1 for ind in population_of_most_different if ind.rawfitness < temp_fitness])
+            population[i].fitness = DeapFitness(
+                tuple((dissim_value, temp_ind_fit)))
+
+    def make_new_population_nsga2(self, population, prob_mut, prob_xov):
+        N = len(population)
+        expected_mut = int(N * prob_mut)
+        expected_xov = int(N * prob_xov)
+        assignCrowdingDist(population)
+        offspring = tools.selTournamentDCD(population, N)
+
+        def addGenotypeIfValid(ind_list, genotype):
+            new_individual = Individual()
+            new_individual.set_and_evaluate(genotype, self.evaluate)
+            if new_individual.fitness is not BAD_FITNESS:  # this is how we defined BAD_FITNESS in frams_evaluate()
+                ind_list.append(new_individual)
+
+        counter = 0
+
+        def get_indyvidual(pop, c):
+            if c < len(pop):
+                ind = pop[c]
+                c += 1
+                return ind, c
+            else:
+                c = 0
+                ind = pop[c]
+                c += 1
+                return ind, c
+
+        newpop = []
+        while len(newpop) < expected_mut:
+            ind, counter = get_indyvidual(offspring, counter)
+            addGenotypeIfValid(newpop, self.mutate(ind.genotype))
+
+        # adding valid crossovers of selected individuals...
+        while len(newpop) < expected_mut + expected_xov:
+            ind1, counter = get_indyvidual(offspring, counter)
+            ind2, counter = get_indyvidual(offspring, counter)
+            addGenotypeIfValid(newpop, self.cross_over(
+                ind1.genotype, ind2.genotype))
+
+        # select clones to fill up the new population until we reach the same size as the input population
+        while len(newpop) < len(population):
+            ind, counter = get_indyvidual(offspring, counter)
+            newpop.append(Individual().copyFrom(ind))
+
+        pop_offspring = population+newpop
+        print(len(pop_offspring))
+        if self.fit == "nslc":
+            self.do_nslc_dissim(pop_offspring)
+        elif self.fit == "nsga2":
+            self.do_nsga2_dissim(pop_offspring)
+        out_pop = tools.selNSGA2(pop_offspring, len(population))
+        return out_pop
+
+    def evolve(self, hof_savefile, generations, initialgenotype, pmut, pxov, tournament_size):
+        file_name = self.get_state_filename(hof_savefile)
+        state = self.load_state(file_name)
+        if state is not None:  # loaded state from file
+            # saved generation has been completed, start with the next one
+            self.current_generation += 1
+            print("...Resuming from saved state: population size = %d, hof size = %d, stats size = %d, archive size = %d, generation = %d/%d" % (len(self.population_structures.population), len(self.hof),
+                                                                                                                                                 len(self.stats),  (len(self.population_structures.archive)), self.current_generation, generations))  # self.current_generation (and g) are 0-based, parsed_args.generations is 1-based
+        else:
+            self.initialize_evolution(self.genformat, initialgenotype)
+
+        time0 = time.process_time()
+        for g in range(self.current_generation, generations):
+            if self.fit != "raw" and self.fit != "nsga2" and self.fit != "nslc":
+                self.do_niching(self.population_structures)
+
+            if type(self.population_structures.population[0].fitness) == DeapFitness:
+                self.population_structures.population = self.make_new_population_nsga2(
+                    self.population_structures.population, pmut, pxov)
+            else:
+                self.population_structures.population = self.make_new_population(
+                    self.population_structures.population, pmut, pxov, tournament_size)
+
+            self.update_stats(g, self.population_structures.population)
+
+            if hof_savefile is not None:
+                self.current_generation = g
+                self.time_elapsed += time.process_time() - time0
+                self.save_state(file_name)
+        if hof_savefile is not None:
+            self.save_genotypes(hof_savefile)
+        return self.population_structures.population, self.stats
+
+    @staticmethod
+    def get_args_for_parser():
+        parser = ExperimentABC.get_args_for_parser()
+        parser.add_argument("-dissim",type= int, default= "frams",
+                   help="Dissimilarity measure type. Availible -2:emd, -1:lev, 1:frams1 (default}, 2:frams2")
+        parser.add_argument("-fit",type= str, default= "raw",
+                        help="Fitness type, availible  types: niching, novelty, nsga2, nslc and raw (default}")
+        parser.add_argument("-archive",type= int, default= 50,
+                            help="Maximum archive size")
+        parser.add_argument("-normalize",type= str, default= "max",
+                            help="What normalization use for dissimilarity matrix, max (default}, sum and none")
+        parser.add_argument("-knn",type= int, default= 0,
+                        help="Nearest neighbours parameter for local novelty/niching, if knn==0 global is performed.Default:0")
+        return parser 
+        
+    @abstractmethod
+    def dissimilarity(self, population: list):
+        pass
Index: /framspy/evolalg/base/random_sequence_index.py
===================================================================
--- /framspy/evolalg/base/random_sequence_index.py	(revision 1190)
+++ /framspy/evolalg/base/random_sequence_index.py	(revision 1190)
@@ -0,0 +1,16 @@
+import numpy as np
+
+
+class RandomIndexSequence:
+    """A helper function for tournament selection. Allows each individual to take part in a tournament,
+    thus reducing the risk of overlooking individuals by not sampling them due to uncontrolled
+    randomness, as can happen in the naive implementation of tournament selection."""
+
+    def __init__(self, popsize):
+        self.popsize: int = popsize
+        self.permut: list = []
+
+    def getNext(self):
+        if len(self.permut) < 1:
+            self.permut = list(np.random.permutation(self.popsize))
+        return self.permut.pop()
Index: /framspy/evolalg/base/remove_diagonal.py
===================================================================
--- /framspy/evolalg/base/remove_diagonal.py	(revision 1190)
+++ /framspy/evolalg/base/remove_diagonal.py	(revision 1190)
@@ -0,0 +1,11 @@
+import numpy as np
+
+
+def remove_diagonal(arrays):
+    numpy_arrays = np.array(arrays)
+    n_of_arrays = numpy_arrays.shape[0]
+    strided = np.lib.stride_tricks.as_strided
+    stride_0, stride_1 = numpy_arrays.strides
+    return strided(numpy_arrays.ravel()[1:],
+                   shape=(n_of_arrays-1, n_of_arrays),
+                   strides=(stride_0+stride_1, stride_1)).reshape(n_of_arrays, -1)
Index: /framspy/evolalg/constants.py
===================================================================
--- /framspy/evolalg/constants.py	(revision 1190)
+++ /framspy/evolalg/constants.py	(revision 1190)
@@ -0,0 +1,1 @@
+BAD_FITNESS = None
Index: /framspy/evolalg/frams_base/experiment_frams.py
===================================================================
--- /framspy/evolalg/frams_base/experiment_frams.py	(revision 1190)
+++ /framspy/evolalg/frams_base/experiment_frams.py	(revision 1190)
@@ -0,0 +1,126 @@
+from ..base.experiment_abc import ExperimentABC
+from ..constants import BAD_FITNESS
+from ..structures.individual import Individual
+from ..structures.population import PopulationStructures
+from ..utils import ensureDir
+
+
+class ExperimentFrams(ExperimentABC):
+    def __init__(self, hof_size, popsize, frams_lib, optimization_criteria, genformat, save_only_best=True, constraints={}) -> None:
+        ExperimentABC.__init__(self, hof_size=hof_size, popsize=popsize, save_only_best=save_only_best)
+        self.optimization_criteria = optimization_criteria
+        self.frams_lib = frams_lib
+        self.constraints = constraints
+        self.genformat = genformat
+
+    def frams_getsimplest(self, genetic_format, initial_genotype):
+        return initial_genotype if initial_genotype is not None else self.frams_lib.getSimplest(genetic_format)
+
+    def genotype_within_constraint(self, genotype, dict_criteria_values, criterion_name, constraint_value):
+        REPORT_CONSTRAINT_VIOLATIONS = False
+        if constraint_value is not None:
+            actual_value = dict_criteria_values[criterion_name]
+            if actual_value > constraint_value:
+                if REPORT_CONSTRAINT_VIOLATIONS:
+                    print('Genotype "%s" assigned low fitness because it violates constraint "%s": %s exceeds threshold %s' % (
+                        genotype, criterion_name, actual_value, constraint_value))
+                return False
+        return True
+
+    def check_valid_constraints(self, genotype, default_evaluation_data):
+        valid = True
+        valid &= self.genotype_within_constraint(
+            genotype, default_evaluation_data, 'numparts', self.constraints.get('max_numparts'))
+        valid &= self.genotype_within_constraint(
+            genotype, default_evaluation_data, 'numjoints', self.constraints.get('max_numjoints'))
+        valid &= self.genotype_within_constraint(
+            genotype, default_evaluation_data, 'numneurons', self.constraints.get('max_numneurons'))
+        valid &= self.genotype_within_constraint(
+            genotype, default_evaluation_data, 'numconnections', self.constraints.get('max_numconnections'))
+        valid &= self.genotype_within_constraint(
+            genotype, default_evaluation_data, 'numgenocharacters', self.constraints.get('max_numgenochars'))
+        return valid
+
+    def evaluate(self, genotype):
+        data = self.frams_lib.evaluate([genotype])
+        # print("Evaluated '%s'" % genotype, 'evaluation is:', data)
+        valid = True
+        try:
+            first_genotype_data = data[0]
+            evaluation_data = first_genotype_data["evaluations"]
+            default_evaluation_data = evaluation_data[""]
+            fitness = [default_evaluation_data[crit] for crit in self.optimization_criteria][0]
+        # the evaluation may have failed for an invalid genotype (such as X[@][@] with "Don't simulate genotypes with warnings" option) or for some other reason
+        except (KeyError, TypeError) as e:
+            valid = False
+            print('Problem "%s" so could not evaluate genotype "%s", hence assigned it fitness: %s' % (
+                str(e), genotype, BAD_FITNESS))
+        if valid:
+            default_evaluation_data['numgenocharacters'] = len(genotype)  # for consistent constraint checking below
+            valid = self.check_valid_constraints(genotype, default_evaluation_data) 
+        if not valid:
+            fitness = BAD_FITNESS
+        return fitness
+        
+
+    def mutate(self, gen1):
+        return self.frams_lib.mutate([gen1])[0]
+
+    def cross_over(self, gen1, gen2):
+        return self.frams_lib.crossOver(gen1, gen2)
+
+    def initialize_evolution(self, initialgenotype):
+        self.current_generation = 0
+        self.time_elapsed = 0
+        self.stats = []  # stores the best individuals, one from each generation
+        initial_individual = Individual()
+        initial_individual.set_and_evaluate(self.frams_getsimplest(
+            '1' if self.genformat is None else self.genformat, initialgenotype), self.evaluate)
+        self.hof.add(initial_individual)
+        self.stats.append(
+            initial_individual.rawfitness if self.save_only_best else initial_individual)
+        self.population_structures = PopulationStructures(
+            initial_individual=initial_individual, popsize=self.popsize)
+
+    def save_genotypes(self, filename):
+        from framsfiles import writer as framswriter
+        with open(filename, "w") as outfile:
+            for ind in self.hof:
+                keyval = {}
+                # construct a dictionary with criteria names and their values
+                for i, k in enumerate(self.optimization_criteria):
+                    # .values[i]  # TODO it would be better to save in Individual (after evaluation) all fields returned by Framsticks, and get these fields here, not just the criteria that were actually used as fitness in evolution.
+                    keyval[k] = ind.rawfitness
+                # Note: prior to the release of Framsticks 5.0, saving e.g. numparts (i.e. P) without J,N,C breaks re-calcucation of P,J,N,C in Framsticks and they appear to be zero (nothing serious).
+                outfile.write(framswriter.from_collection(
+                    {"_classname": "org", "genotype": ind.genotype, **keyval}))
+                outfile.write("\n")
+        print("Saved '%s' (%d)" % (filename, len(self.hof)))
+
+    @staticmethod
+    def get_args_for_parser():
+        parser = ExperimentABC.get_args_for_parser()
+        parser.add_argument('-path',type= ensureDir, required= True,
+                        help= 'Path to Framsticks CLI without trailing slash.')
+        parser.add_argument('-lib',type= str, required= False,
+                        help= 'Library name. If not given, "frams-objects.dll" or "frams-objects.so" is assumed depending on the platform.')
+        parser.add_argument('-sim',type= str, required= False, default= "eval-allcriteria.sim",
+                        help="The name of the .sim file with settings for evaluation, mutation, crossover, and similarity estimation. If not given, \"eval-allcriteria.sim\" is assumed by default. Must be compatible with the \"standard-eval\" expdef. If you want to provide more files, separate them with a semicolon ';'.")
+
+        parser.add_argument('-genformat',type= str, required= False,
+                            help= 'Genetic format for the simplest initial genotype, for example 4, 9, or B. If not given, f1 is assumed.')
+        parser.add_argument('-initialgenotype',type= str, required= False,
+                                    help= 'The genotype used to seed the initial population. If given, the -genformat argument is ignored.')
+        parser.add_argument('-opt',required=True, help='optimization criteria: vertpos, velocity, distance, vertvel, lifespan, numjoints, numparts, numneurons, numconnections (or other as long as it is provided by the .sim file and its .expdef).')
+
+        parser.add_argument('-max_numparts',type= int, default= None,
+                                help="Maximum number of Parts. Default: no limit")
+        parser.add_argument('-max_numjoints',type= int, default= None,
+                                help="Maximum number of Joints. Default: no limit")
+        parser.add_argument('-max_numneurons',type= int, default= None,
+                                help="Maximum number of Neurons. Default: no limit")
+        parser.add_argument('-max_numconnections',type= int, default= None,
+                                    help="Maximum number of Neural connections. Default: no limit")
+        parser.add_argument('-max_numgenochars',type= int, default= None,
+                                    help="Maximum number of characters in genotype (including the format prefix, if any}. Default: no limit")
+        return parser
Index: /framspy/evolalg/frams_base/experiment_frams_islands.py
===================================================================
--- /framspy/evolalg/frams_base/experiment_frams_islands.py	(revision 1190)
+++ /framspy/evolalg/frams_base/experiment_frams_islands.py	(revision 1190)
@@ -0,0 +1,36 @@
+
+from ..base.experiment_islands_model_abc import ExperimentIslands
+from ..frams_base.experiment_frams import ExperimentFrams
+from ..structures.individual import Individual
+from ..structures.population import PopulationStructures
+from ..utils import merge_two_parsers
+
+
+class ExperimentFramsIslands(ExperimentIslands, ExperimentFrams):
+    def __init__(self, frams_lib, optimization_criteria, hof_size,
+                 popsize, constraints, genformat,
+                 number_of_populations, migration_interval, save_only_best) -> None:
+        ExperimentFrams.__init__(self, frams_lib=frams_lib, optimization_criteria=optimization_criteria,
+                                 hof_size=hof_size, popsize=popsize,
+                                 genformat=genformat, save_only_best=save_only_best, constraints=constraints)
+
+        self.number_of_populations = number_of_populations
+        self.migration_interval = migration_interval
+
+    def initialize_evolution(self, initialgenotype):
+        self.current_generation = 0
+        self.time_elapsed = 0
+        # stores the best individuals, one from each generation across all populations
+        self.stats = []
+        initial_individual = Individual()
+        initial_individual.set_and_evaluate(self.frams_getsimplest(
+            '1' if self.genformat is None else self.genformat, initialgenotype), self.evaluate)
+        self.stats.append(initial_individual.rawfitness)
+        self.populations= [PopulationStructures(initial_individual=initial_individual,
+                                              popsize=self.popsize)
+                         for _ in range(self.number_of_populations)]
+    @staticmethod
+    def get_args_for_parser():
+        parser1 = ExperimentFrams.get_args_for_parser()
+        parser2 = ExperimentIslands.get_args_for_parser()
+        return merge_two_parsers(parser1, parser2)
Index: /framspy/evolalg/frams_base/experiment_frams_niching.py
===================================================================
--- /framspy/evolalg/frams_base/experiment_frams_niching.py	(revision 1190)
+++ /framspy/evolalg/frams_base/experiment_frams_niching.py	(revision 1190)
@@ -0,0 +1,54 @@
+from ..base.experiment_niching_abc import ExperimentNiching
+from ..structures.individual import Individual
+from ..structures.population import PopulationStructures
+from ..utils import merge_two_parsers
+from .experiment_frams import ExperimentFrams
+
+
+class ExperimentFramsNiching(ExperimentFrams, ExperimentNiching):
+    def __init__(self, frams_lib, optimization_criteria, hof_size, popsize, constraints, normalize, dissim, fit, genformat, archive_size, save_only_best) -> None:
+        ExperimentFrams.__init__(self, hof_size=hof_size,
+                                 popsize=popsize,
+                                 frams_lib=frams_lib,
+                                 constraints=constraints,
+                                 optimization_criteria=optimization_criteria,
+                                 genformat=genformat,
+                                 save_only_best=save_only_best
+                                 )
+        ExperimentNiching.__init__(self, hof_size=hof_size,
+                                   popsize=popsize,
+                                   fit=fit,
+                                   normalize=normalize,
+                                   save_only_best=save_only_best,
+                                   archive_size=archive_size
+                                   )
+        self.dissim = dissim
+        
+
+    def initialize_evolution(self, genformat, initialgenotype):
+        self.current_generation = 0
+        self.time_elapsed = 0
+        self.stats = []  # stores the best individuals, one from each generation
+        initial_individual = Individual()
+        initial_individual.set_and_evaluate(self.frams_getsimplest(
+            '1' if genformat is None else genformat, initialgenotype), self.evaluate)
+        self.hof.add(initial_individual)
+        self.stats.append(
+            initial_individual.rawfitness if self.save_only_best else initial_individual)
+        self.population_structures = PopulationStructures(
+            initial_individual=initial_individual, archive_size=self.archive_size, popsize=self.popsize)
+        if self.fit == "nsga2":
+            self.do_nsga2_dissim(self.population_structures.population)
+        if self.fit == "nslc":
+            self.do_nslc_dissim(self.population_structures.population)
+
+    def dissimilarity(self, population):
+        return self.frams_lib.dissimilarity([i.genotype for i in population], self.dissim)
+
+
+    @staticmethod
+    def get_args_for_parser():
+        p1 = ExperimentFrams.get_args_for_parser()
+        p2 = ExperimentNiching.get_args_for_parser()
+        return merge_two_parsers(p1, p2)
+
Index: /framspy/evolalg/json_encoders.py
===================================================================
--- /framspy/evolalg/json_encoders.py	(revision 1190)
+++ /framspy/evolalg/json_encoders.py	(revision 1190)
@@ -0,0 +1,12 @@
+import json
+
+import numpy as np
+
+
+class Encoder(json.JSONEncoder):
+    def default(self, obj):
+        if isinstance(obj, np.ndarray):
+            return obj.tolist()
+        if isinstance(obj, np.integer):
+            return int(obj)
+        return super(Encoder, self).default(obj)
Index: /framspy/evolalg/numerical_example/numerical_example.py
===================================================================
--- /framspy/evolalg/numerical_example/numerical_example.py	(revision 1190)
+++ /framspy/evolalg/numerical_example/numerical_example.py	(revision 1190)
@@ -0,0 +1,20 @@
+import numpy as np
+
+from ..base.experiment_abc import ExperimentABC
+from ..structures.hall_of_fame import HallOfFame
+
+
+class ExperimentNumerical(ExperimentABC):
+    def __init__(self, hof_size, popsize, save_only_best) -> None:
+        ExperimentABC.__init__(self,popsize=popsize,
+                               hof_size=hof_size,
+                               save_only_best=save_only_best)
+
+    def mutate(self, gen1):
+        return list(gen1 + np.random.randint(-10, 10, len(gen1)))
+
+    def cross_over(self, gen1, gen2):
+        return gen1
+
+    def evaluate(self, genotype):
+        return 1/sum([x*x for x in genotype])
Index: /framspy/evolalg/numerical_example/numerical_islands_example.py
===================================================================
--- /framspy/evolalg/numerical_example/numerical_islands_example.py	(revision 1190)
+++ /framspy/evolalg/numerical_example/numerical_islands_example.py	(revision 1190)
@@ -0,0 +1,22 @@
+import numpy as np
+
+from ..base.experiment_islands_model_abc import ExperimentIslands
+from ..structures.hall_of_fame import HallOfFame
+
+
+class ExperimentNumericalIslands(ExperimentIslands):
+    def __init__(self, hof_size, popsize, number_of_populations, migration_interval, save_only_best) -> None:
+        ExperimentIslands.__init__(self,popsize=popsize,
+                                hof_size=hof_size,
+                                number_of_populations=number_of_populations,
+                                migration_interval=migration_interval,
+                                save_only_best=save_only_best)
+
+    def mutate(self, gen1):
+        return gen1 + np.random.randint(-10, 10, len(gen1))
+
+    def cross_over(self, gen1, gen2):
+        return gen1
+
+    def evaluate(self, genotype):
+        return 1/sum([x*x for x in genotype])
Index: /framspy/evolalg/run_frams_islands.py
===================================================================
--- /framspy/evolalg/run_frams_islands.py	(revision 1190)
+++ /framspy/evolalg/run_frams_islands.py	(revision 1190)
@@ -0,0 +1,44 @@
+from FramsticksLib import FramsticksLib
+
+from .frams_base.experiment_frams_islands import ExperimentFramsIslands
+
+
+def main():
+    # random.seed(123)  # see FramsticksLib.DETERMINISTIC below, set to True if you want full determinism
+    # must be set before FramsticksLib() constructor call
+    FramsticksLib.DETERMINISTIC = False
+    parsed_args =ExperimentFramsIslands.get_args_for_parser().parse_args()
+    print("Argument values:", ", ".join(
+        ['%s=%s' % (arg, getattr(parsed_args, arg)) for arg in vars(parsed_args)]))
+    # multiple criteria not supported here. If needed, use FramsticksEvolution.py
+    opt_criteria = parsed_args.opt.split(",")
+    framsLib = FramsticksLib(
+        parsed_args.path, parsed_args.lib, parsed_args.sim.split(";"))
+    constrains = {"max_numparts": parsed_args.max_numparts,
+                  "max_numjoints": parsed_args.max_numjoints,
+                  "max_numneurons": parsed_args.max_numneurons,
+                  "max_numconnections": parsed_args.max_numconnections,
+                  "max_numgenochars": parsed_args.max_numgenochars,
+                  }
+
+    print('Best individuals:')
+    experiment = ExperimentFramsIslands(frams_lib=framsLib,
+                                        optimization_criteria=opt_criteria,
+                                        hof_size=parsed_args.hof_size,
+                                        constraints=constrains,
+                                        genformat=parsed_args.genformat,
+                                        popsize=parsed_args.popsize,
+                                        migration_interval=parsed_args.generations_migration,
+                                        number_of_populations=parsed_args.islands,
+                                        save_only_best=parsed_args.save_only_best)
+
+    experiment.evolve(hof_savefile=parsed_args.hof_savefile,
+                      generations=parsed_args.generations,
+                      initialgenotype=parsed_args.initialgenotype,
+                      pmut=parsed_args.pmut,
+                      pxov=parsed_args.pxov,
+                      tournament_size=parsed_args.tournament)
+
+
+if __name__ == "__main__":
+    main()
Index: /framspy/evolalg/run_frams_niching.py
===================================================================
--- /framspy/evolalg/run_frams_niching.py	(revision 1190)
+++ /framspy/evolalg/run_frams_niching.py	(revision 1190)
@@ -0,0 +1,46 @@
+from FramsticksLib import FramsticksLib
+
+from .frams_base.experiment_frams_niching import ExperimentFramsNiching
+
+
+def main():
+    # random.seed(123)  # see FramsticksLib.DETERMINISTIC below, set to True if you want full determinism
+    # must be set before FramsticksLib() constructor call
+    FramsticksLib.DETERMINISTIC = False
+    parsed_args = ExperimentFramsNiching.get_args_for_parser().parse_args()
+    print("Argument values:", ", ".join(
+        ['%s=%s' % (arg, getattr(parsed_args, arg)) for arg in vars(parsed_args)]))
+    # multiple criteria not supported here. If needed, use FramsticksEvolution.py
+    opt_criteria = parsed_args.opt.split(",")
+    framsLib = FramsticksLib(
+        parsed_args.path, parsed_args.lib, parsed_args.sim.split(";"))
+    constrains = {"max_numparts": parsed_args.max_numparts,
+                  "max_numjoints": parsed_args.max_numjoints,
+                  "max_numneurons": parsed_args.max_numneurons,
+                  "max_numconnections": parsed_args.max_numconnections,
+                  "max_numgenochars": parsed_args.max_numgenochars,
+                  }
+
+    print('Best individuals:')
+    experiment = ExperimentFramsNiching(frams_lib=framsLib,
+                                        optimization_criteria=opt_criteria,
+                                        hof_size=parsed_args.hof_size,
+                                        constraints=constrains,
+                                        normalize=parsed_args.normalize,
+                                        dissim=parsed_args.dissim,
+                                        fit=parsed_args.fit,
+                                        genformat=parsed_args.genformat,
+                                        popsize=parsed_args.popsize,
+                                        archive_size=parsed_args.archive,
+                                        save_only_best=parsed_args.save_only_best)
+
+    experiment.evolve(hof_savefile=parsed_args.hof_savefile,
+                      generations=parsed_args.generations,
+                      initialgenotype=parsed_args.initialgenotype,
+                      pmut=parsed_args.pmut,
+                      pxov=parsed_args.pxov,
+                      tournament_size=parsed_args.tournament)
+
+
+if __name__ == "__main__":
+    main()
Index: /framspy/evolalg/run_numerical_example.py
===================================================================
--- /framspy/evolalg/run_numerical_example.py	(revision 1190)
+++ /framspy/evolalg/run_numerical_example.py	(revision 1190)
@@ -0,0 +1,28 @@
+from .numerical_example.numerical_example import ExperimentNumerical
+
+
+def main():
+    parsed_args = ExperimentNumerical.get_args_for_parser().parse_args()
+    print("Argument values:", ", ".join(
+        ['%s=%s' % (arg, getattr(parsed_args, arg)) for arg in vars(parsed_args)]))
+
+    initialgenotype = [100, 100, 100, 100]
+    print('Best individuals:')
+    experiment = ExperimentNumerical(
+        hof_size=parsed_args.hof_size,
+        popsize=parsed_args.popsize,
+        save_only_best=parsed_args.save_only_best)
+
+    hof, stats = experiment.evolve(hof_savefile=parsed_args.hof_savefile,
+                                   generations=parsed_args.generations,
+                                   initialgenotype=initialgenotype,
+                                   pmut=parsed_args.pmut,
+                                   pxov=parsed_args.pxov,
+                                   tournament_size=parsed_args.tournament)
+    print(len(hof))
+    for ind in hof:
+        print(ind.genotype, ind.rawfitness)
+
+
+if __name__ == "__main__":
+    main()
Index: /framspy/evolalg/run_numerical_islands_example.py
===================================================================
--- /framspy/evolalg/run_numerical_islands_example.py	(revision 1190)
+++ /framspy/evolalg/run_numerical_islands_example.py	(revision 1190)
@@ -0,0 +1,31 @@
+from .numerical_example.numerical_islands_example import \
+    ExperimentNumericalIslands
+
+
+def main():
+    parsed_args = ExperimentNumericalIslands.get_args_for_parser().parse_args()
+    print("Argument values:", ", ".join(
+        ['%s=%s' % (arg, getattr(parsed_args, arg)) for arg in vars(parsed_args)]))
+
+    initialgenotype = [100, 100, 100, 100]
+    print('Best individuals:')
+    experiment = ExperimentNumericalIslands(
+        hof_size=parsed_args.hof_size,
+        popsize=parsed_args.popsize,
+        migration_interval=parsed_args.generations_migration,
+        number_of_populations=parsed_args.islands,
+        save_only_best=parsed_args.save_only_best)
+
+    hof, stats = experiment.evolve(hof_savefile=parsed_args.hof_savefile,
+                                   generations=parsed_args.generations,
+                                   initialgenotype=initialgenotype,
+                                   pmut=parsed_args.pmut,
+                                   pxov=parsed_args.pxov,
+                                   tournament_size=parsed_args.tournament)
+    print(len(hof))
+    for ind in hof:
+        print(ind.genotype, ind.rawfitness)
+
+
+if __name__ == "__main__":
+    main()
Index: /framspy/evolalg/structures/hall_of_fame.py
===================================================================
--- /framspy/evolalg/structures/hall_of_fame.py	(revision 1190)
+++ /framspy/evolalg/structures/hall_of_fame.py	(revision 1190)
@@ -0,0 +1,23 @@
+class HallOfFame:
+    """A simple function that keeps the specified number of individuals, adding only better ones.
+    The list remains sorted from best to worst."""
+
+    def __init__(self, hofsize):
+        self.hofsize = hofsize
+        self.hof = []
+
+    def __iter__(self):
+        return iter(self.hof)
+
+    def __len__(self):
+        return len(self.hof)
+
+    def add(self, individual):
+        if len(self.hof) < 1:  # empty hof?
+            self.hof.append(individual)  # then add the first individual
+        else:  # we have some individuals in hof?
+            # only add if the new one is better than the first stored in hof
+            if individual.rawfitness > self.hof[0].rawfitness:
+                self.hof.insert(0, individual)  # add as first
+        while len(self.hof) > self.hofsize:  # exceeded desired hof capacity?
+            self.hof.pop()  # delete last (=worst)
Index: /framspy/evolalg/structures/individual.py
===================================================================
--- /framspy/evolalg/structures/individual.py	(revision 1190)
+++ /framspy/evolalg/structures/individual.py	(revision 1190)
@@ -0,0 +1,43 @@
+from ..constants import BAD_FITNESS
+
+
+class Individual:
+    only_positive_fitness = True
+
+    def __init__(self):
+        self.genotype = None
+        # used for stats. This is raw fitness value, None = not evaluated or invalid genotype
+        self.rawfitness: float = None
+        # used in selection and can be modified e.g. by diversity or niching techniques
+        self.fitness: float = None
+
+    def copyFrom(self, individual):  # "copying constructor"
+        self.genotype = individual.genotype
+        self.rawfitness = individual.rawfitness
+        self.fitness = individual.fitness
+        return self
+
+    def copy(self):
+        new_copy = Individual()
+        new_copy.genotype = self.genotype
+        new_copy.rawfitness = self.rawfitness
+        new_copy.fitness = self.fitness
+        return new_copy
+
+    def set_and_evaluate(self, genotype, evaluate):
+        self.genotype = genotype
+        fitness = evaluate(genotype)
+        if fitness is not BAD_FITNESS:  # BAD_FITNESS indicates that the genotype was not valid or some other problem occurred during evaluation
+            if self.only_positive_fitness:
+                if fitness < 0:
+                    fitness = 0
+        self.fitness = self.rawfitness = fitness
+
+    def __str__(self):
+        try:
+            return "%g\t%g\t'%s'" % (self.rawfitness, self.fitness, self.genotype)
+        except:
+            return "%g\t'%s'" % (self.rawfitness, self.genotype)
+
+    def __gt__(self, other):
+        return self.rawfitness > other.rawfitness
Index: /framspy/evolalg/structures/population.py
===================================================================
--- /framspy/evolalg/structures/population.py	(revision 1190)
+++ /framspy/evolalg/structures/population.py	(revision 1190)
@@ -0,0 +1,40 @@
+import random
+
+import numpy as np
+
+from ..base.remove_diagonal import remove_diagonal
+
+
+class PopulationStructures:
+    def __init__(self,  initial_individual, archive_size=0, popsize=100) -> None:
+        self.population_size = popsize
+        self.population = [initial_individual.copy()
+                           for _ in range(self.population_size)]
+        self.archive = []
+        self.archive_size = archive_size
+
+    def change_individual(self, dissim):
+        no_diagonal = remove_diagonal(dissim.copy())
+        minval = np.min(no_diagonal)
+        to_remove_pair = list(set(np.where(no_diagonal == minval)[0]))
+        last_index = len(dissim)-1
+        if last_index in to_remove_pair:
+            return last_index
+        else:
+            return random.choice(to_remove_pair)
+
+    def update_archive(self, dissim_matrix, population_archive):
+        if self.archive_size < 1:
+            return
+        current_archive_ind = [i for i in range(
+            self.population_size, len(dissim_matrix))]
+        for i in range(self.population_size):
+            current_archive_ind = current_archive_ind + [i]
+            temp_dissim_matrix = dissim_matrix[current_archive_ind,
+                                               ][:, current_archive_ind]
+            if len(current_archive_ind) > self.archive_size:
+                to_remove = self.change_individual(temp_dissim_matrix)
+                current_archive_ind.pop(to_remove)
+
+        self.archive = np.array(population_archive)[
+            current_archive_ind].tolist()
Index: /framspy/evolalg/tests/test_diferent_settings.py
===================================================================
--- /framspy/evolalg/tests/test_diferent_settings.py	(revision 1190)
+++ /framspy/evolalg/tests/test_diferent_settings.py	(revision 1190)
@@ -0,0 +1,239 @@
+import argparse
+import sys
+
+from FramsticksLib import FramsticksLib
+
+from ..frams_base.experiment_frams_niching import ExperimentFramsNiching
+from ..frams_base.experiment_frams_islands import ExperimentFramsIslands
+from ..numerical_example.numerical_example import ExperimentNumerical
+from ..numerical_example.numerical_islands_example import ExperimentNumericalIslands
+
+
+from ..utils import ensureDir
+
+SETTINGS_TO_TEST_NUMERIC = {
+    'hof_size': [0, 10],
+    'popsize': [8],
+    'archive': [8],
+    'pmut': [0.7],
+    'pxov': [0.2],
+    'tournament': [5],
+    'initialgenotype':[[100, 100, 100, 100], [-100,-100]]
+}
+
+SETTINGS_TO_TEST_NUMERIC_ISLAND = {
+    'hof_size': [0, 10],
+    'popsize': [8],
+    'archive': [8],
+    'pmut': [0.7],
+    'pxov': [0.2],
+    'tournament': [5],
+    'migration_interval': [1,5],
+    'number_of_populations':[1,5],
+    'initialgenotype':[[100, 100, 100, 100], [-100,-100]]
+}
+
+SETTINGS_TO_TEST_FRAMS_NICHING = {
+    'opt': ['velocity', 'vertpos'],
+    'max_numparts': [None],
+    'max_numjoints': [20],
+    'max_numneurons': [20],
+    'max_numconnections': [None],
+    'max_numgenochars': [20],
+    'hof_size': [0, 10],
+    'normalize': ['none', 'max', 'sum'],
+    'dissim': [-2, -1, 1, 2],
+    'fit': ['niching', 'novelty', 'nsga2', 'nslc', 'raw'],
+    'genformat': ['1'],
+    'popsize': [8],
+    'archive': [8],
+    'initialgenotype': [None],
+    'pmut': [0.7],
+    'pxov': [0.2],
+    'tournament': [5]
+}
+
+SETTINGS_TO_TEST_FRAMS_ISLANDS = {
+    'opt': ['velocity', 'vertpos'],
+    'max_numparts': [None],
+    'max_numjoints': [20],
+    'max_numneurons': [20],
+    'max_numconnections': [None],
+    'max_numgenochars': [20],
+    'hof_size': [0, 10],
+    'migration_interval': [1,5],
+    'number_of_populations':[1,5],
+    'genformat': ['1'],
+    'popsize': [8],
+    'initialgenotype': [None],
+    'pmut': [0.7],
+    'pxov': [0.2],
+    'tournament': [5]
+}
+
+def test_run_experiment_numerical(params):
+    # multiple criteria not supported here. If needed, use FramsticksEvolution.py
+
+    experiment = ExperimentNumerical(
+                                    hof_size=params['hof_size'],
+                                    popsize=params['popsize'],
+                                    save_only_best=True,)
+
+    experiment.evolve(hof_savefile=None,
+                      generations=2,
+                      initialgenotype=params['initialgenotype'],
+                      pmut=params['pmut'],
+                      pxov=params['pxov'],
+                      tournament_size=params['tournament'])
+
+
+def test_run_experiment_numerical_islands(params):
+    # multiple criteria not supported here. If needed, use FramsticksEvolution.py
+
+    experiment = ExperimentNumericalIslands(hof_size=params['hof_size'],
+                                    popsize=params['popsize'],
+                                    save_only_best=True,
+                                    migration_interval=params['migration_interval'],
+                                    number_of_populations=params['number_of_populations'])
+
+
+    experiment.evolve(hof_savefile=None,
+                      generations=2,
+                      initialgenotype=params['initialgenotype'],
+                      pmut=params['pmut'],
+                      pxov=params['pxov'],
+                      tournament_size=params['tournament'])
+
+def test_run_experiment_frams_niching(params):
+    # multiple criteria not supported here. If needed, use FramsticksEvolution.py
+    opt_criteria = params['opt'].split(",")
+    framsLib = FramsticksLib(
+        parsed_args.path, parsed_args.lib, parsed_args.sim.split(";"))
+    constrains = {"max_numparts": params['max_numparts'],
+                  "max_numjoints": params['max_numjoints'],
+                  "max_numneurons": params['max_numneurons'],
+                  "max_numconnections": params['max_numconnections'],
+                  "max_numgenochars": params['max_numgenochars'],
+                  }
+
+    print('Best individuals:')
+    experiment = ExperimentFramsNiching(frams_lib=framsLib,
+                                        optimization_criteria=opt_criteria,
+                                        hof_size=params['hof_size'],
+                                        constraints=constrains,
+                                        normalize=params['normalize'],
+                                        dissim=params['dissim'],
+                                        fit=params['fit'],
+                                        genformat=params['genformat'],
+                                        popsize=params['popsize'],
+                                        archive_size=params['archive'],
+                                        save_only_best=True)
+
+    experiment.evolve(hof_savefile=None,
+                      generations=2,
+                      initialgenotype=params['initialgenotype'],
+                      pmut=params['pmut'],
+                      pxov=params['pxov'],
+                      tournament_size=params['tournament'])
+
+def test_run_experiment_frams_island(params):
+    # multiple criteria not supported here. If needed, use FramsticksEvolution.py
+    opt_criteria = params['opt'].split(",")
+    framsLib = FramsticksLib(
+        parsed_args.path, parsed_args.lib, parsed_args.sim.split(";"))
+    constrains = {"max_numparts": params['max_numparts'],
+                  "max_numjoints": params['max_numjoints'],
+                  "max_numneurons": params['max_numneurons'],
+                  "max_numconnections": params['max_numconnections'],
+                  "max_numgenochars": params['max_numgenochars'],
+                  }
+
+    print('Best individuals:')
+    experiment = ExperimentFramsIslands(frams_lib=framsLib,
+                                        optimization_criteria=opt_criteria,
+                                        hof_size=params['hof_size'],
+                                        constraints=constrains,
+                                        genformat=params['genformat'],
+                                        popsize=params['popsize'],
+                                        migration_interval=params['migration_interval'],
+                                        number_of_populations=params['number_of_populations'],
+                                        save_only_best=True)
+
+    experiment.evolve(hof_savefile=None,
+                      generations=2,
+                      initialgenotype=params['initialgenotype'],
+                      pmut=params['pmut'],
+                      pxov=params['pxov'],
+                      tournament_size=params['tournament'])
+
+
+def parseArguments():
+    parser = argparse.ArgumentParser(
+        description='Run this program with "python -u %s" if you want to disable buffering of its output.' % sys.argv[0])
+    parser.add_argument('-path', type=ensureDir, required=True,
+                        help='Path to Framsticks CLI without trailing slash.')
+    parser.add_argument('-lib', required=False,
+                        help='Library name. If not given, "frams-objects.dll" or "frams-objects.so" is assumed depending on the platform.')
+    parser.add_argument('-sim', required=False, default="eval-allcriteria.sim",
+                        help="The name of the .sim file with settings for evaluation, mutation, crossover, and similarity estimation. If not given, \"eval-allcriteria.sim\" is assumed by default. Must be compatible with the \"standard-eval\" expdef. If you want to provide more files, separate them with a semicolon ';'.")
+
+    return parser.parse_args()
+
+
+def get_params_sets(settings):
+    params_sets = []
+    for k in settings.keys():
+        temp_param_set = []
+        for value in settings[k]:
+            if params_sets:
+                for exsiting_set in params_sets:
+                    copy_of_set = exsiting_set.copy()
+                    copy_of_set[k] = value
+                    temp_param_set.append(copy_of_set)
+            else:
+                temp_param_set.append({k: value})
+        params_sets = temp_param_set
+    return params_sets
+
+
+def cover_to_test(params, run_exp):
+    try:
+        run_exp(params)
+        return [1, None]
+    except Exception as e:
+        return [0, f"Experiment {run_exp.__name__} with params:{params} failied with the stack:{e}"]
+
+
+def run_tests():
+    result = []
+    print("TESTING NUMERICAL")
+    params_sets = get_params_sets(SETTINGS_TO_TEST_NUMERIC)
+    print(f"Starting executing {len(params_sets)} experiments")
+    result.extend([cover_to_test(params,test_run_experiment_numerical) for params in params_sets])
+
+    print("TESTING NUMERICAL ISLANDS")
+    params_sets = get_params_sets(SETTINGS_TO_TEST_NUMERIC_ISLAND)
+    print(f"Starting executing {len(params_sets)} experiments")
+    result.extend([cover_to_test(params,test_run_experiment_numerical_islands) for params in params_sets])
+
+    print("TESTING FRAMS NICHING")
+    params_sets = get_params_sets(SETTINGS_TO_TEST_FRAMS_NICHING)
+    print(f"Starting executing {len(params_sets)} experiments")
+    result.extend([cover_to_test(params, test_run_experiment_frams_niching) for params in params_sets])
+
+    print("TESTING FRAMS ISLANDS")
+    params_sets = get_params_sets(SETTINGS_TO_TEST_FRAMS_ISLANDS)
+    print(f"Starting executing {len(params_sets)} experiments")
+    result.extend([cover_to_test(params,test_run_experiment_frams_island) for params in params_sets])
+
+    Successful = sum([1 for r in result  if r[0] == 1])
+    Failed = sum([1 for r in result if r[0] == 0])
+    print(f"{Successful} out of {len(result)} passed")
+    if Successful < Failed:
+        print("Experiments that failed :{Failed}")
+        [print(r[1]) for r in result if r[1]]
+
+
+if __name__ == "__main__":
+    parsed_args = parseArguments()
+    run_tests()
Index: /framspy/evolalg/utils.py
===================================================================
--- /framspy/evolalg/utils.py	(revision 1190)
+++ /framspy/evolalg/utils.py	(revision 1190)
@@ -0,0 +1,56 @@
+import argparse
+import os
+
+
+def ensureDir(string):
+    if os.path.isdir(string):
+        return string
+    else:
+        raise NotADirectoryError(string)
+
+def merge_two_parsers(p1,p2):
+    """
+    This function is a modification of  argparse _ActionsContainer._add_container_actions
+    that allows for merging duplicates 
+    """
+    # collect groups by titles
+    title_group_map = {}
+    for group in p1._action_groups:
+        if group.title in title_group_map:
+            msg = ('cannot merge actions - two groups are named %r')
+            raise ValueError(msg % (group.title))
+        title_group_map[group.title] = group
+
+    # map each action to its group
+    group_map = {}
+    for group in p2._action_groups:
+
+        # if a group with the title exists, use that, otherwise
+        # create a new group matching the p2's group
+        if group.title not in title_group_map:
+            title_group_map[group.title] = p1.add_argument_group(
+                title=group.title,
+                description=group.description,
+                conflict_handler=group.conflict_handler)
+
+        # map the actions to their new group
+        for action in group._group_actions:
+            group_map[action] = title_group_map[group.title]
+
+    # add p2's mutually exclusive groups
+    # NOTE: if add_mutually_exclusive_group ever gains title= and
+    # description= then this code will need to be expanded as above
+    for group in p2._mutually_exclusive_groups:
+        mutex_group = p1.add_mutually_exclusive_group(
+            required=group.required)
+
+        # map the actions to their new mutex group
+        for action in group._group_actions:
+            group_map[action] = mutex_group
+    # add all actions to this p2 or their group
+    for action in p2._actions:
+        try:
+            group_map.get(action, p1)._add_action(action)
+        except:
+            print("Warning:",action.option_strings, "is a duplicate" )
+    return p1
Index: /framspy/run-evolalg-examples.cmd
===================================================================
--- /framspy/run-evolalg-examples.cmd	(revision 1190)
+++ /framspy/run-evolalg-examples.cmd	(revision 1190)
@@ -0,0 +1,23 @@
+rem To learn about all available options of each .py algorithm below, add "-h" to its parameters.
+rem Use the source code of the examples as a starting point for your customizations.
+rem Example usage:
+
+set DIR_WITH_FRAMS_LIBRARY=.................
+
+
+
+
+rem evolution with niching
+python -m evolalg.run_frams_niching  -path %DIR_WITH_FRAMS_LIBRARY%  -sim eval-allcriteria.sim;deterministic.sim;sample-period-longest.sim -opt velocity -genformat 0 -dissim 1 -fit knn_niching -archive 50 -max_numparts 15 -max_numneurons 15 -max_numjoints 30 -max_numconnections 30  -max_numgenochars 10000 -popsize 50 -generations 10 -normalize none -hof_savefile HoF-niching.gen 
+
+rem a generic island model example
+python -m evolalg.run_frams_islands  -path %DIR_WITH_FRAMS_LIBRARY%  -islands=10 -generations_migration=5 -sim eval-allcriteria.sim;only-body.sim;deterministic.sim;sample-period-2.sim -opt vertpos -genformat 1 -max_numparts 15 -max_numneurons 15 -max_numjoints 5 -max_numconnections 5  -max_numgenochars 100 -popsize 20 -generations 15  -hof_savefile HoF-islands.gen
+
+rem a combination of various parameters for thorough testing, see source of 'test_diferent_settings.py'
+rem python -m evolalg.tests.test_diferent_settings  -path %DIR_WITH_FRAMS_LIBRARY%  -sim eval-allcriteria.sim
+
+
+
+rem numerical examples (unrelated to Framsticks)
+python -m evolalg.run_numerical_example
+python -m evolalg.run_numerical_islands_example
