diff --git a/src/aevol_modify.cpp b/src/aevol_modify.cpp
index ae1fe162e9e62729678dd50e0f148180fa69d7e4..1247561ceb616c5f4723a21980958ce472d0fd2b 100644
--- a/src/aevol_modify.cpp
+++ b/src/aevol_modify.cpp
@@ -25,6 +25,12 @@
// ****************************************************************************
+#include "AeTime.h"
+#include "DnaMutator.h"
+#include "Individual_7.h"
+#include "MutationParams.h"
+#include "ae_enums.h"
+#include "legacy/population/Individual.h"
const char* DEFAULT_PARAM_FILE_NAME = "param.in";
// ============================================================================
@@ -34,6 +40,7 @@ const char* DEFAULT_PARAM_FILE_NAME = "param.in";
#include <cstdio>
#include <cstring>
#include <fstream>
+#include <filesystem>
#include <getopt.h>
@@ -234,11 +241,203 @@ int main(int argc, char* argv[]) {
exp_manager->output_m()->set_big_backup_step(atol(line->words[1]));
}
else if (strcmp(line->words[0], "POPULATION_SIZE") == 0) {
- printf("ERROR in param file \"%s\" on line %" PRId32
- ": the change of population size is not implemented yet\n"
- " for spatially structured populations",
- param_file_name, cur_line);
- exit(EXIT_FAILURE);
+ if (exp_manager->sel()->selection_scope() != aevol::SCOPE_GLOBAL) {
+ Utils::ExitWithDevMsg("changes of population size supports only (for now) global selection scope", param_file_name, cur_line);
+ }
+
+ if (line->nb_words < 2) {
+ Utils::ExitWithDevMsg("population size parameter is missing.", param_file_name, cur_line);
+ } else if (line->nb_words < 3) {
+ Utils::ExitWithDevMsg("sampling type parameter is missing.", param_file_name, cur_line);
+ }
+
+ bool custom_name = false;
+ if (line->nb_words == 4) {
+ if (strncmp(line->words[3], "custom", 6)) {
+ custom_name = true;
+ }
+ }
+
+ if (strncmp(line->words[2], "with", 4) != 0 && strncmp(line->words[2], "no", 2) != 0) {
+ std::string error_msg("unknown sampling type \"");
+ error_msg += line->words[2];
+ error_msg += "\".";
+ Utils::ExitWithUsrMsg(error_msg);
+ }
+
+ int32_t new_pop_size = atoi(line->words[1]);
+ if (floor(sqrt((double)new_pop_size)) != ceil(sqrt((double)new_pop_size))) {
+ Utils::ExitWithUsrMsg("Population size should be a perfect square.");
+ }
+
+ std::string population_changes_filename("population_changes");
+ if (custom_name) {
+ population_changes_filename += "_";
+ population_changes_filename += AeTime::time();
+ }
+ population_changes_filename += ".csv";
+
+ std::filesystem::path population_changes_path(population_changes_filename);
+ if (std::filesystem::exists(population_changes_path)) {
+ std::string error_msg("'");
+ error_msg += population_changes_path;
+ error_msg += "' already exists.";
+ Utils::ExitWithUsrMsg(error_msg);
+ }
+
+ bool with_selection = strncmp(line->words[2], "no", 2) != 0;
+
+ grid_width_ = (int)sqrt((double)new_pop_size);
+ grid_height_ = grid_width_;
+
+ World *old_world = exp_manager->non_const_world();
+ int32_t old_pop_size = old_world->width() * old_world->height();
+
+ World *new_world = new World();
+ new_world->set_prng(old_world->non_const_prng());
+ new_world->set_mut_prng(old_world->non_const_mut_prng());
+ new_world->set_stoch_prng(old_world->non_const_stoch_prng());
+ new_world->InitGrid(grid_width_, grid_height_, old_world->grid(0, 0)->habitat_non_const(), true);
+
+ for (int32_t x = 0; x < grid_width_; x++) {
+ for (int32_t y = 0; y < grid_height_; y++) {
+ int32_t seed = new_world->prng()->random(1000000);
+
+ new_world->grid(x, y)->set_reprod_prng(std::make_unique<JumpingMT>(seed));
+ new_world->grid(x, y)->set_reprod_prng_simd(std::make_unique<JumpingMT>(seed));
+ }
+ }
+
+ DnaFactory *dna_factory = new DnaFactory(DnaFactory_Policy::FIRSTFIT, 32, 5000);
+ FuzzyFactory_7 *fuzzy_factory = new FuzzyFactory_7(exp_manager->exp_s()->get_fuzzy_flavor(), old_pop_size * 4, exp_manager->world()->phenotypic_target_handler()->sampling());
+
+ double *fitness_array = new double[old_pop_size];
+ double *probas = new double[old_pop_size];
+ double fitness_sum = 0.0;
+
+ for (int32_t indiv_id = 0; indiv_id < old_pop_size; indiv_id++) {
+ if (!with_selection) {
+ fitness_array[indiv_id] = 1;
+ fitness_sum += 1;
+ } else {
+ int32_t x = indiv_id % old_world->height();
+ int32_t y = indiv_id / old_world->height();
+
+ Individual *indiv = new Individual(*old_world->grid(x, y)->individual());
+ indiv->Evaluate();
+
+ fitness_array[indiv_id] = indiv->fitness();
+ fitness_sum += indiv->fitness();
+
+ delete indiv;
+ }
+ }
+
+ for (int32_t indiv_id = 0; indiv_id < old_pop_size; indiv_id++) {
+ probas[indiv_id] = fitness_array[indiv_id] / fitness_sum;
+ }
+
+ int32_t *nb_offsprings = new int32_t[old_pop_size];
+ exp_manager->world()->grid(0, 0)->reprod_prng_simd_->multinomial_drawing(nb_offsprings, probas, new_pop_size, old_pop_size);
+
+ int64_t index = 0;
+
+ std::ofstream pop_changes;
+ pop_changes.open(population_changes_path, std::ofstream::trunc);
+ if(!pop_changes.is_open()) {
+ std::ostringstream error_msg("Couldn't open file '");
+ error_msg << population_changes_path;
+ error_msg << "'.";
+ Utils::ExitWithUsrMsg(error_msg.str());
+ }
+ pop_changes << "Generation,old_index,new_index" << std::endl;
+
+ for (int32_t indiv_id = 0; indiv_id < old_pop_size; indiv_id++) {
+ for (int32_t offspring_nb = 0; offspring_nb < nb_offsprings[indiv_id]; offspring_nb++) {
+ pop_changes << AeTime::time() << "," << indiv_id << "," << index << std::endl;
+
+ int32_t x = index % grid_height_;
+ int32_t y = index / grid_height_;
+ int32_t old_x = indiv_id % old_world->height();
+ int32_t old_y = indiv_id / old_world->height();
+
+ Individual* indiv = new Individual(old_world->indiv_at(old_x, old_y),
+ index,
+ new_world->grid(x, y)->mut_prng(),
+ new_world->grid(x, y)->stoch_prng());
+ new_world->PlaceIndiv(indiv, x, y, true);
+
+ index++;
+ }
+ }
+
+ pop_changes.close();
+
+ Individual_7** previous_individuals = new Individual_7*[new_pop_size];
+ DnaMutator** dna_mutator_array = new DnaMutator*[new_pop_size];
+
+ for (int32_t indiv_id = 0; indiv_id < std::min(old_pop_size, new_pop_size); ++indiv_id) {
+ previous_individuals[indiv_id] = exp_manager->exp_m_7_->previous_individuals[indiv_id];
+ dna_mutator_array[indiv_id] = exp_manager->dna_mutator_array_[indiv_id];
+ }
+
+ for (int32_t indiv_id = old_pop_size; indiv_id < new_pop_size; ++indiv_id) {
+ Individual_7* indiv = new Individual_7(previous_individuals[0]->w_max_, dna_factory, fuzzy_factory);
+ indiv->dna_ = dna_factory->get_dna(old_world->grid(0, 0)->individual()->genetic_unit_seq_length(0));
+ indiv->dna_->set_indiv(old_world->grid(0, 0)->individual()->genetic_unit(0).dna(), dna_factory);
+ indiv->dna_->set_indiv(indiv);
+ indiv->parent_id = 0;
+
+ dna_mutator_array[indiv_id] = new DnaMutator(old_world->grid(0, 0)->mut_prng(),
+ indiv->dna_->length(),
+ exp_s->mut_params()->duplication_rate(),
+ exp_s->mut_params()->deletion_rate(),
+ exp_s->mut_params()->translocation_rate(),
+ exp_s->mut_params()->inversion_rate(),
+ exp_s->mut_params()->point_mutation_rate(),
+ exp_s->mut_params()->small_insertion_rate(),
+ exp_s->mut_params()->small_deletion_rate(),
+ exp_s->mut_params()->max_indel_size(),
+ exp_s->min_genome_length(),
+ exp_s->max_genome_length(),
+ indiv_id,
+ 0,
+ 0);
+
+ previous_individuals[indiv_id] = indiv;
+ }
+
+ for (int32_t indiv_id = new_pop_size; indiv_id < old_pop_size; indiv_id++) {
+ delete exp_manager->exp_m_7_->previous_individuals[indiv_id];
+ exp_manager->exp_m_7_->previous_individuals[indiv_id] = nullptr;
+
+ delete exp_manager->dna_mutator_array_[indiv_id];
+ }
+
+ delete[] exp_manager->exp_m_7_->previous_individuals;
+ delete[] exp_manager->dna_mutator_array_;
+ exp_manager->exp_m_7_->previous_individuals = previous_individuals;
+ exp_manager->dna_mutator_array_ = dna_mutator_array;
+
+ delete[] nb_offsprings;
+ delete[] fitness_array;
+ delete[] probas;
+
+ new_world->update_best();
+ exp_manager->set_world(new_world);
+
+ delete old_world;
+
+ std::ofstream log_file;
+ log_file.open("aevol_modify.log", std::ofstream::app);
+ log_file << "Change of population size at generation " << AeTime::time() << " : "
+ << "from " << old_pop_size << " to " << new_pop_size << " "
+ << (with_selection ? "with" : "without")
+ << " selection." << std::endl;
+ log_file.close();
+
+ exp_manager->output_m()->init_tree(exp_manager,exp_manager->output_m()->tree()->tree_step());
+ tree = exp_manager->tree();
}
else if (strcmp(line->words[0], "SELECTION_SCHEME") == 0) {
if (strncmp(line->words[1], "lin", 3) == 0) {