add support for population size changes in aevol_modify

See merge request aevol/aevol!52

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) {