diff --git a/src/libaevol/BitSet_SIMD.h b/src/libaevol/BitSet_SIMD.h
index 43fed401539ba1a5b027574e777830d80225ce3b..8d37fb459d70d18d651b0b7604ccbe66326ffdf2 100644
--- a/src/libaevol/BitSet_SIMD.h
+++ b/src/libaevol/BitSet_SIMD.h
@@ -430,17 +430,21 @@ class BitSet_SIMD {
     char* to_char() {
 #ifdef _DYNAMIC_CUSTOM_BITSET
       char* ret = new char[length_+1];
+
       for (int32_t i = 0; i < length_; i++) {
         if (get(i))
           ret[i] = '1';
         else
           ret[i] = '0';
       }
+
       ret[length_] = '\0';
       return ret;
 #elif _STATIC_BITSET
       char *ret = (char*)data_.to_string().c_str();
       return ret;
+#else
+      return nullptr;
 #endif
     }
 
diff --git a/src/libaevol/CMakeLists.txt b/src/libaevol/CMakeLists.txt
index 3e34681e85486cbb3b732f002d7ba0e44eb2335f..e0440fa22f27b931d6578e86d149300a23db8aa1 100644
--- a/src/libaevol/CMakeLists.txt
+++ b/src/libaevol/CMakeLists.txt
@@ -177,7 +177,9 @@ if ( ${with-raevol} )
           raevol/Protein_R.h
           raevol/Protein_R.cpp
           raevol/Rna_R.h
-          raevol/Rna_R.cpp)
+          raevol/Rna_R.cpp 
+	  #"SIMD_Abstract Metadata.cpp" "SIMD_Abstract Metadata.h" SIMD_Map_Metadata.cpp SIMD_Map_Metadata.h
+	  )
 endif ()
 
 add_library(aevol ${aevol_sources})
diff --git a/src/libaevol/DnaMutator.h b/src/libaevol/DnaMutator.h
index 968a569f5a7f971d6ebe028710e1f1387779450c..b862c0f8b4f3dcf398b4894a95287900bbfc5f5c 100644
--- a/src/libaevol/DnaMutator.h
+++ b/src/libaevol/DnaMutator.h
@@ -44,7 +44,7 @@ class DnaMutator {
 
     bool hasMutate() {return hasMutate_;}
 
-    bool setMutate(bool mutate) {hasMutate_ = mutate;}
+    void setMutate(bool mutate) {hasMutate_ = mutate;}
 
     int x_,y_;
 
diff --git a/src/libaevol/ExpManager.cpp b/src/libaevol/ExpManager.cpp
index b7497b533c94c719bbc03b7510a8adf882fd7985..765d7fd47e93afb410d018d56ece9fecd1a0d359 100644
--- a/src/libaevol/ExpManager.cpp
+++ b/src/libaevol/ExpManager.cpp
@@ -324,7 +324,7 @@ void ExpManager::step_to_next_generation() {
 #ifdef __PERF_LOG__
 auto     t1 = high_resolution_clock::now();
 
-  auto     s_t1 = high_resolution_clock::now();
+  //auto     s_t1 = high_resolution_clock::now();
 #endif
 
 #ifdef __CUDACC__
diff --git a/src/libaevol/ExpManager_X11.cpp b/src/libaevol/ExpManager_X11.cpp
index 73325f8e9bf3febde943f6e4a09039af78b67f78..8d1110f188b249d7199ea43759dc0db83755997f 100644
--- a/src/libaevol/ExpManager_X11.cpp
+++ b/src/libaevol/ExpManager_X11.cpp
@@ -1043,16 +1043,16 @@ void ExpManager_X11::refresh_window(int8_t win_number) {
 #endif
 
       // Display all the phenotypes (blue)
-      for (const auto& indiv: indivs())
+#ifndef __REGUL
+
+        for (const auto& indiv: indivs())
       {
-        #ifndef __REGUL
         display(cur_win, *(indiv->phenotype()), BLUE);
         if (indiv->allow_plasmids())
         {
           display(cur_win, *(indiv->genetic_unit(0).phenotypic_contribution()), YELLOW);
           display(cur_win, *(indiv->genetic_unit(1).phenotypic_contribution()), GREEN);
         }
-        #else
         /*Individual_R_X11* indiv_r = dynamic_cast<Individual_R_X11*>(indiv);
 
         display(cur_win, *(indiv_r->get_phenotype()), BLUE);
@@ -1061,8 +1061,8 @@ void ExpManager_X11::refresh_window(int8_t win_number) {
           display(cur_win, *(indiv_r->get_genetic_unit(0).get_phenotypic_contribution()), YELLOW);
           display(cur_win, *(indiv_r->get_genetic_unit(1).get_phenotypic_contribution()), GREEN);
         }*/
-        #endif
       }
+#endif
 
       // Display best indiv's phenotype (white)
       #ifndef __REGUL
diff --git a/src/libaevol/GeneticUnit.cpp b/src/libaevol/GeneticUnit.cpp
index 435c51a7532049dcf9d5f1fe5f0f4cea9e3e8caf..abe7bb50dd73627438d1e0e085b56b669be4e896 100644
--- a/src/libaevol/GeneticUnit.cpp
+++ b/src/libaevol/GeneticUnit.cpp
@@ -1494,24 +1494,24 @@ bool GeneticUnit::is_promoter(Strand strand, int32_t pos, int8_t& dist) const {
   const char* genome = dna_->data();
   int32_t len = dna_->length();
 
-  float dist_a[22] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+  int8_t dist_a[22] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
   int32_t pos_a[22];
   if (strand == LEADING) {
     //~ printf("LEADING\n");
     if (pos + PROM_SIZE < len) {
-#pragma vector always
+//#pragma vector always
       for (int32_t i  = 0; i < PROM_SIZE; i++)
         pos_a[i] = pos + i;
     } else {
-#pragma vector always
+//#pragma vector always
       for (int32_t i  = 0; i < PROM_SIZE; i++)
         pos_a[i] = (pos + i) % len;
 
     }
 
-    #pragma vector always
+//    #pragma vector always
     for (int16_t i = 0; i < PROM_SIZE; i++) {
-      dist_a[i] = genome[pos_a[i]] != PROM_SEQ[i] ? 1 : 0;
+      dist_a[i] = genome[pos_a[i]] != PROM_SEQ[i] ? (int8_t ) 1 : (int8_t ) 0;
     }
   }
   else // (strand == LAGGING)
@@ -1529,9 +1529,9 @@ bool GeneticUnit::is_promoter(Strand strand, int32_t pos, int8_t& dist) const {
       pos_a[i] = (pos - i) >= 0 ? (pos - i) % len : (len - abs_i) % len;
     }
 
-    #pragma vector always
+  //  #pragma vector always
     for (int16_t i = 0; i < PROM_SIZE; i++) {
-      dist_a[i] = genome[pos_a[i]] == PROM_SEQ[i] ? 1 : 0;
+      dist_a[i] = genome[pos_a[i]] == PROM_SEQ[i] ? (int8_t ) 1 : (int8_t )  0;
     }
   }
 
diff --git a/src/libaevol/HabitatFactory.cpp b/src/libaevol/HabitatFactory.cpp
index d7ddbbfd505ffdb4c5626acec76cd563a67e8e41..05f4fd7f2b279840dbc40df6869d9762d9366fbf 100644
--- a/src/libaevol/HabitatFactory.cpp
+++ b/src/libaevol/HabitatFactory.cpp
@@ -53,11 +53,7 @@ HabitatFactory::create_unique_habitat(gzFile backup_file,
                                                                    bool share_phenotypic_target) {
 
 #ifndef __REGUL
-        #if __cplusplus == 201103L
-  return make_unique<Habitat> (habitat, share_phenotypic_target);
-#else
-  return std::make_unique<Habitat> (habitat, share_phenotypic_target);
-#endif
+
 #else
 #if __cplusplus == 201103L
         return make_unique<Habitat_R> (dynamic_cast<Habitat_R&>(habitat), share_phenotypic_target);
diff --git a/src/libaevol/ParamLoader.cpp b/src/libaevol/ParamLoader.cpp
index efb605aab5b1df3cd0a489519307159322d19972..fe7dec55f01c54a0d009fcf648eb6c6fc55e41a0 100644
--- a/src/libaevol/ParamLoader.cpp
+++ b/src/libaevol/ParamLoader.cpp
@@ -951,7 +951,7 @@ void ParamLoader::interpret_line(ParameterLine * line, int32_t cur_line)
     #ifdef __REGUL
       // le premier chiffre est l'indice d'environment en convention humaine ( le premier a 1)
       // On vérifie que cet indice n'est pas trop élevé ni négatif pour éviter les crash
-      if ( atoi(line->words[1]) - 1 < _env_gaussians_list.size() && atoi(line->words[1]) > 0)
+      if ( atoi(line->words[1]) - 1 < (int)_env_gaussians_list.size() && atoi(line->words[1]) > 0)
       {
         (_env_gaussians_list.at( atoi(line->words[1]) - 1)).push_back
         ( Gaussian(  atof( line->words[2] ), atof( line->words[3] ), atof( line->words[4] ) ) );
@@ -1416,7 +1416,7 @@ void ParamLoader::interpret_line(ParameterLine * line, int32_t cur_line)
     {
       // le premier chiffre est l'indice d'environment en convention humaine ( le premier a 1)
       // On vérifie que cet indice n'est pas trop élevé ni négatif pour éviter les crash
-      if ( atoi(line->words[1]) - 1 < _env_signals_list.size() && atoi(line->words[1]) > 0)
+      if ( atoi(line->words[1]) - 1 < (int)_env_signals_list.size() && atoi(line->words[1]) > 0)
       {
         (_env_signals_list.at( atoi(line->words[1]) - 1)).push_back(atoi(line->words[2]) - 1);
       }
diff --git a/src/libaevol/Protein.cpp b/src/libaevol/Protein.cpp
index f5ed23d2ddd28bb37524c44474f176385eb4b8cd..b454c140695ca74e912363d1b00f56c3b1615934 100644
--- a/src/libaevol/Protein.cpp
+++ b/src/libaevol/Protein.cpp
@@ -83,7 +83,10 @@ Protein::Protein(GeneticUnit* gen_unit, const Protein &model)
   // Copy the list of amino-acids
 
   // TODO vld: check if deep copy needed (at least, in case of heritage yes)
-  AA_list_ = model.AA_list_;
+  //AA_list_ = model.AA_list_;
+          for (auto aa_c : model.AA_list_) {
+            AA_list_.push_back(new Codon(aa_c->value()));
+          }
 
   // Copy triangle parameters
   mean_   = model.mean_;
diff --git a/src/libaevol/SIMD_Individual.cpp b/src/libaevol/SIMD_Individual.cpp
index e5d38ea5242979e61922d01bfbeb9289a13c8d53..ff899aa0ce4d0a81facabb282ac734eb0cc5ffbb 100644
--- a/src/libaevol/SIMD_Individual.cpp
+++ b/src/libaevol/SIMD_Individual.cpp
@@ -54,10 +54,10 @@ SIMD_Individual::SIMD_Individual(ExpManager* exp_m) {
   }
 
   dna_size = new int[exp_m_->nb_indivs()];
-  int x, y;
+  //int x, y;
   for (int indiv_id = 0; indiv_id < exp_m_->nb_indivs(); indiv_id++) {
-    x = indiv_id / exp_m_->world()->height();
-    y = indiv_id % exp_m_->world()->height();
+    //x = indiv_id / exp_m_->world()->height();
+    //y = indiv_id % exp_m_->world()->height();
 
 
     dna_size[indiv_id] = internal_simd_struct[indiv_id]->dna_->length();
@@ -970,7 +970,7 @@ void SIMD_Individual::start_stop_RNA() {
   //int x, y;
 
   //
-  ExpManager* exp_m = exp_m_;
+  //ExpManager* exp_m = exp_m_;
 
   //#pragma omp parallel for collapse(2) default(shared)
 //#pragma omp parallel for
@@ -1000,8 +1000,8 @@ void SIMD_Individual::start_stop_RNA() {
         bool is_terminator_lag = internal_simd_struct[indiv_id]->dna_->bitset_->is_terminator(
             false, dna_pos);
 #else
-      int x = indiv_id / exp_m->world()->height();
-      int y = indiv_id % exp_m->world()->height();
+      //int x = indiv_id / exp_m->world()->height();
+      //int y = indiv_id % exp_m->world()->height();
 
       int len = dna_size[indiv_id];
       if (len >= PROM_SIZE) {
@@ -1207,11 +1207,11 @@ void SIMD_Individual::start_stop_RNA() {
 
 
     void SIMD_Individual::start_stop_RNA(int indiv_id) {
-        int nb_indiv = exp_m_->nb_indivs();
+        //int nb_indiv = exp_m_->nb_indivs();
         //int x, y;
 
         //
-        ExpManager* exp_m = exp_m_;
+        //ExpManager* exp_m = exp_m_;
 
         //#pragma omp parallel for collapse(2) default(shared)
 //#pragma omp parallel for
@@ -1240,8 +1240,8 @@ void SIMD_Individual::start_stop_RNA() {
         bool is_terminator_lag = internal_simd_struct[indiv_id]->dna_->bitset_->is_terminator(
             false, dna_pos);
 #else
-                int x = indiv_id / exp_m->world()->height();
-                int y = indiv_id % exp_m->world()->height();
+                //int x = indiv_id / exp_m->world()->height();
+                //int y = indiv_id % exp_m->world()->height();
 
                 int len = dna_size[indiv_id];
                 if (len >= PROM_SIZE) {
@@ -1738,7 +1738,7 @@ void SIMD_Individual::opt_prom_compute_RNA() {
 
     void SIMD_Individual::opt_prom_compute_RNA(int indiv_id) {
 
-        int nb_indiv = exp_m_->nb_indivs();
+        //int nb_indiv = exp_m_->nb_indivs();
 
 
             if (exp_m_->dna_mutator_array_[indiv_id]->hasMutate()) {
@@ -1769,7 +1769,7 @@ void SIMD_Individual::opt_prom_compute_RNA() {
 #ifdef WITH_FINETASKLOOP
 #pragma omp taskloop grainsize(rna_grain_size)
 #endif
-                for (int prom_idx = 0; prom_idx< internal_simd_struct[indiv_id]->promoters.size(); prom_idx++) {
+                for (int prom_idx = 0; prom_idx<(int) internal_simd_struct[indiv_id]->promoters.size(); prom_idx++) {
 
                     if (internal_simd_struct[indiv_id]->promoters[prom_idx] != nullptr) {
                         int rna_idx = prom_idx;
@@ -2427,8 +2427,8 @@ void SIMD_Individual::compute_RNA() {
 //#pragma omp task firstprivate(indiv_id, rna_idx) depend(out: internal_simd_struct[indiv_id])
             {
                 if (internal_simd_struct[indiv_id]->rnas[rna_idx]->is_init_) {
-                    int x = indiv_id / exp_m_->world()->height();
-                    int y = indiv_id % exp_m_->world()->height();
+                    //int x = indiv_id / exp_m_->world()->height();
+                    //int y = indiv_id % exp_m_->world()->height();
 
                     int c_pos = internal_simd_struct[indiv_id]->rnas[rna_idx]->begin;
 
@@ -2562,8 +2562,8 @@ void SIMD_Individual::start_protein() {
 //#pragma omp task firstprivate(indiv_id, rna_idx) depend(inout: internal_simd_struct[indiv_id])
           {
             if (internal_simd_struct[indiv_id]->rnas[rna_idx]->is_init_) {
-              int x = indiv_id / exp_m_->world()->height();
-              int y = indiv_id % exp_m_->world()->height();
+              //int x = indiv_id / exp_m_->world()->height();
+              //int y = indiv_id % exp_m_->world()->height();
 
               int c_pos = internal_simd_struct[indiv_id]->rnas[rna_idx]->begin;
 
@@ -2713,8 +2713,8 @@ void SIMD_Individual::compute_protein(int indiv_id) {
                                      rnas[rna_idx]->start_prot.size(); protein_idx++) {
 //#pragma omp task firstprivate(indiv_id, rna_idx, protein_idx) depend(inout: internal_simd_struct[indiv_id])
                             {
-                                int x = indiv_id / exp_m_->world()->height();
-                                int y = indiv_id % exp_m_->world()->height();
+                                //int x = indiv_id / exp_m_->world()->height();
+                                //int y = indiv_id % exp_m_->world()->height();
 
                                 int start_protein_pos = internal_simd_struct[indiv_id]->
                                         rnas[rna_idx]->leading_lagging == 0 ?
@@ -3055,8 +3055,8 @@ void SIMD_Individual::compute_protein() {
                      rnas[rna_idx]->start_prot.size(); protein_idx++) {
 //#pragma omp task firstprivate(indiv_id, rna_idx, protein_idx) depend(in: internal_simd_struct[indiv_id])
               {
-                int x = indiv_id / exp_m_->world()->height();
-                int y = indiv_id % exp_m_->world()->height();
+                //int x = indiv_id / exp_m_->world()->height();
+                //int y = indiv_id % exp_m_->world()->height();
 
                 int start_protein_pos = internal_simd_struct[indiv_id]->
                     rnas[rna_idx]->leading_lagging == 0 ?
@@ -3384,8 +3384,8 @@ void SIMD_Individual::translate_protein(double w_max) {
 //#pragma omp task firstprivate(indiv_id, protein_idx) depend(inout: internal_simd_struct[indiv_id])
           {
             if (internal_simd_struct[indiv_id]->proteins[protein_idx]->is_init_) {
-              int x = indiv_id / exp_m_->world()->height();
-              int y = indiv_id % exp_m_->world()->height();
+              //int x = indiv_id / exp_m_->world()->height();
+              //int y = indiv_id % exp_m_->world()->height();
 
               int c_pos = internal_simd_struct[indiv_id]->proteins[protein_idx]->protein_start, t_pos;
               int end_pos = internal_simd_struct[indiv_id]->proteins[protein_idx]->protein_end;
@@ -3650,8 +3650,8 @@ void SIMD_Individual::translate_protein(double w_max) {
 //#pragma omp task firstprivate(indiv_id, protein_idx) depend(inout: internal_simd_struct[indiv_id])
             {
                 if (internal_simd_struct[indiv_id]->proteins[protein_idx]->is_init_) {
-                    int x = indiv_id / exp_m_->world()->height();
-                    int y = indiv_id % exp_m_->world()->height();
+                    //int x = indiv_id / exp_m_->world()->height();
+                    //int y = indiv_id % exp_m_->world()->height();
 
                     int c_pos = internal_simd_struct[indiv_id]->proteins[protein_idx]->protein_start, t_pos;
                     int end_pos = internal_simd_struct[indiv_id]->proteins[protein_idx]->protein_end;
@@ -4544,7 +4544,7 @@ void SIMD_Individual::run_a_step(double w_max, double selection_pressure,bool op
         for (int indiv_id = 0; indiv_id < (int) exp_m_->nb_indivs(); indiv_id++) {
             //if (indiv_id == 0)
             //    printf("%d -- usage %d -- \n", indiv_id, internal_simd_struct[indiv_id]->usage_count_);
-            int usage_cpt = 0;
+            //int usage_cpt = 0;
 
 #pragma omp critical
             {
@@ -4574,7 +4574,7 @@ void SIMD_Individual::run_a_step(double w_max, double selection_pressure,bool op
 //#pragma omp single
 //#pragma omp taskloop
         for (int indiv_id = 0; indiv_id < (int) exp_m_->nb_indivs(); indiv_id++) {
-            int usage_cpt = 0;
+            //int usage_cpt = 0;
 
 #pragma omp critical
             {
@@ -4866,7 +4866,7 @@ void SIMD_Individual::check_individual(int i, int x, int y) {
         }
     }
 
-    int prot_cpt_a=0,prot_cpt_b=0;
+    int prot_cpt_b=0;
     idx = 0;
     for (auto prot : exp_m_->world()->grid(x, y)->old_one->protein_list()) {
         bool found = false;
@@ -4943,14 +4943,8 @@ void SIMD_Individual::check_result() {
     float i_fit_2 = roundf(fit_2 * 100);
 
 
-    int prot_size = 0;
-
+    int prot_size = (int) exp_m_->world()->grid(x, y)->individual()->protein_list().size();
 
-    for (auto prot : exp_m_->world()->grid(x, y)->individual()->protein_list()) {
-
-        prot_size++;
-      //}
-    }
 
     if (i_fit_1 != i_fit_2 && dna_size[i] > 300)
     //if (i == 268) {
@@ -4998,7 +4992,7 @@ void SIMD_Individual::check_result() {
 */
 
       idx = 0;
-        int prot_cpt_a=0,prot_cpt_b=0;
+        int prot_cpt_b=0;
 
       for (auto prot : exp_m_->world()->grid(x, y)->individual()->protein_list()) {
           bool found = false;
@@ -5183,7 +5177,7 @@ Internal_SIMD_Struct::~Internal_SIMD_Struct() {
  * We need some index for the promoter optimization
  */
 void Internal_SIMD_Struct::rebuild_index() {
-  if (count_prom > promoters.size()/2) {
+  if (count_prom > (int)promoters.size()/2) {
     /**
      * Do the reindexation process
      */
@@ -5276,7 +5270,7 @@ void Internal_SIMD_Struct::insert_promoters_at(std::vector<std::list<promoterStr
 
     // Insert the promoters in the individual's RNA list
     for (auto& to_insert: promoters_to_insert[strand]) {
-        int prev_pos = to_insert->pos;
+//        int prev_pos = to_insert->pos;
       // Update promoter position
       to_insert->pos = Utils::mod(to_insert->pos + pos, dna_->length());
       if (strand == LEADING) {
diff --git a/src/libaevol/Selection.cpp b/src/libaevol/Selection.cpp
index c0eb3a58f8f363cce5cb24b44563c7c857bee354..4548c4d81e28a55dfe44d3bcd6ebc933feb99517 100644
--- a/src/libaevol/Selection.cpp
+++ b/src/libaevol/Selection.cpp
@@ -419,7 +419,7 @@ void Selection::step_to_next_generation() {
 #pragma omp target teams distribute parallel for schedule(static,1)
 #endif
 #endif
-  for (int i = 0; i < to_evaluate.size(); i++) {
+  for (int i = 0; i < (int) to_evaluate.size(); i++) {
 #ifdef __REGUL
     if ((dynamic_cast<PhenotypicTargetHandler_R*>(&to_evaluate[i]->grid_cell()->habitat().
         phenotypic_target_handler_nonconst())->hasChanged()) ||
diff --git a/src/libaevol/StatRecord.cpp b/src/libaevol/StatRecord.cpp
index 3457fcb443ad0f0dfb1f913b09285b1ffa6d0619..6478172756a60f5be6507fbb7498ae16da816a36 100644
--- a/src/libaevol/StatRecord.cpp
+++ b/src/libaevol/StatRecord.cpp
@@ -164,7 +164,7 @@ StatRecord::StatRecord(const StatRecord &model)
   Individual_R* indiv_r = dynamic_cast<Individual_R*>(indiv);
 
   for (auto& rna: indiv_r->get_rna_list_coding()) {
-    for (unsigned int i = 0; i < ((Rna_R*)rna)->nb_influences(); i++) {
+    for (int i = 0; i < ((Rna_R*)rna)->nb_influences(); i++) {
       //compute the activity
       if (((Rna_R*)rna)->_enhancing_coef_list[i] > 0)
       {
diff --git a/src/libaevol/World.cpp b/src/libaevol/World.cpp
index 12a0c0fbb0a61a8f62616960605c3b6a1c3ab250..a06c8546d5c7c69087324e9ccbe50e3125aceaa9 100644
--- a/src/libaevol/World.cpp
+++ b/src/libaevol/World.cpp
@@ -123,7 +123,11 @@ void World::InitGrid(int16_t width, int16_t height,
       if (share_phenotypic_target)
         grid_[x][y] =
                 new GridCell(x, y,
+#ifdef __REGUL
                              HabitatFactory::create_unique_habitat(dynamic_cast<Habitat_R&>(habitat),share_phenotypic_target),
+#else
+                             HabitatFactory::create_unique_habitat(habitat,share_phenotypic_target),
+#endif
                              NULL,std::make_shared<JumpingMT>(mut_seed),
                              std::make_shared<JumpingMT>(stoch_seed));
     }
@@ -606,12 +610,12 @@ Individual* World::indiv_by_id(int32_t id) const {
   Individual* indiv = grid_1d_[id]->individual();
   // When the population isn't mixed at all, the individual with id n is in
   // grid_1d_[n]. Try this first...
-  if (indiv->id() == id)
+  if ((int32_t)indiv->id() == id)
     return indiv;
   // ... If it isn't, do a basic search
   int32_t nb_indivs = width_ * height_;
   for (int32_t i = 0 ; i < nb_indivs ; i++) {
-    if (grid_1d_[i]->individual()->id() == id)
+    if ((int32_t ) grid_1d_[i]->individual()->id() == id)
       return grid_1d_[i]->individual();
   }
   return nullptr;
diff --git a/src/libaevol/raevol/Individual_R.cpp b/src/libaevol/raevol/Individual_R.cpp
index bbb3484068d7139067c4bf2c86ecbc07405fa997..50ca577210646f2ae4a355dc94a775a573a1d49b 100644
--- a/src/libaevol/raevol/Individual_R.cpp
+++ b/src/libaevol/raevol/Individual_R.cpp
@@ -146,7 +146,7 @@ Individual_R::Individual_R(ExpManager* exp_m, gzFile backup_file) : Individual(
     int32_t nb_inherited_proteins = 0;
     gzread( backup_file, &nb_inherited_proteins,  sizeof(nb_inherited_proteins) );
 
-    printf("Nb prot %d\n",nb_inherited_proteins);
+    //printf("Nb prot %d\n",nb_inherited_proteins);
 
     for ( int16_t i = 0 ; i < nb_inherited_proteins ; i++ )
     {
@@ -665,11 +665,11 @@ void Individual_R::make_rna_list( void )
 
   // Parse the newly created RNA list and copy the coding RNAs in _rna_list_coding.
   for (const auto& gen_unit: genetic_unit_list_) {
-    GeneticUnit* genu = const_cast<GeneticUnit*>(&gen_unit);
+    //GeneticUnit* genu = const_cast<GeneticUnit*>(&gen_unit);
     // Create proxies
     const auto& rna_list = gen_unit.rna_list();
-    const auto& lead = rna_list[LEADING];
-    const auto& lagg = rna_list[LAGGING];
+    //const auto& lead = rna_list[LEADING];
+    //const auto& lagg = rna_list[LAGGING];
 
     // append pointers to rna material to local _rna_list
     for (auto& strand: {LEADING, LAGGING})
@@ -771,7 +771,7 @@ void Individual_R::create_csv(char *directory_name) {
 
   fprintf(drawingfile, "Generation,Protein,Concentration\n");
 
-  int16_t life_time =  exp_m_->exp_s()->get_nb_indiv_age();
+  //int16_t life_time =  exp_m_->exp_s()->get_nb_indiv_age();
   //set the concentrations of proteins to their initial value
   double* concentrations = new double[protein_list_.size()]; // initialise le tableau de concentrations.
   //  int16_t prot_index = 0;
@@ -782,7 +782,7 @@ void Individual_R::create_csv(char *directory_name) {
   }
 
 
-  std::set<int>* eval = exp_m_->exp_s()->get_list_eval_step();
+  //std::set<int>* eval = exp_m_->exp_s()->get_list_eval_step();
   // i is thus the age of the individual
   for (int16_t i = 1; i <= exp_m_->exp_s()->get_nb_indiv_age(); i++) {
     //Set the concentration of signals for this age
@@ -808,7 +808,8 @@ void Individual_R::create_csv(char *directory_name) {
     int proti = 0;
     //printf("Age[%d] : ",i);
     for (const auto& prot : protein_list_) {
-      if (((Protein_R*)prot)->is_signal()) printf("Protein List %ld -> %f (signal: %ld) (%p)\n",((Protein_R*)prot)->get_id(),prot->concentration(),((Protein_R*)prot)->is_signal(),prot);
+      if (((Protein_R*)prot)->is_signal()) printf("Protein List %ld -> %f (signal: %d) (%p)\n",
+                                                  ((Protein_R*)prot)->get_id(),prot->concentration(),((Protein_R*)prot)->is_signal(),prot);
     
 
       // morceau ajouté pour colorer les protéines en fonctions de leur paramètres
@@ -820,7 +821,7 @@ void Individual_R::create_csv(char *directory_name) {
     update_phenotype();
 
     for (auto prot : signal_list) {
-      printf("after %ld prot %ld : %f\n",i,prot.second->get_id(),prot.second->concentration());
+      printf("after %d prot %ld : %f\n",i,prot.second->get_id(),prot.second->concentration());
     }
 
     printf("Creating the EPS file with the phenotype of the chosen individual at step %d... ",i);
@@ -831,7 +832,7 @@ void Individual_R::create_csv(char *directory_name) {
   }
 
     for (const auto& prot : protein_list_) {
-      if (((Protein_R*)prot)->is_signal()) printf("Protein List %d -> %f (signal: %d)\n",((Protein_R*)prot)->get_id(),prot->concentration(),((Protein_R*)prot)->is_signal() );
+      if (((Protein_R*)prot)->is_signal()) printf("Protein List %ld -> %f (signal: %d)\n",((Protein_R*)prot)->get_id(),prot->concentration(),((Protein_R*)prot)->is_signal() );
     }
 
   delete[] concentrations;
diff --git a/src/libaevol/raevol/Individual_R_X11.cpp b/src/libaevol/raevol/Individual_R_X11.cpp
index 0d467c4e3ac84e2dd5b1fce3764d26dd6c95eca1..bcc3841b199dcb8e535727ce908ddb4c1f094470 100644
--- a/src/libaevol/raevol/Individual_R_X11.cpp
+++ b/src/libaevol/raevol/Individual_R_X11.cpp
@@ -359,8 +359,8 @@ void Individual_R_X11::display_regulation( X11Window* win )
   // NB : As we want OriC to be at the "top" of the circle and the orientation
   //      to be clockwise, the drawing angle (theta) will be given as
   //      (90 - alpha), alpha being the "classical" trigonometric angle
-  int16_t alpha_rna_first,alpha_prot_first; // Angles of first and last transcribed bases from OriC (degrees)
-  int16_t theta_rna_first, theta_prot_first; // Transposed angles on the trigonometric circle (degrees)
+  //int16_t alpha_rna_first,alpha_prot_first; // Angles of first and last transcribed bases from OriC (degrees)
+  //int16_t theta_rna_first, theta_prot_first; // Transposed angles on the trigonometric circle (degrees)
   // Same as above with precision = 1/64 degree
   int16_t alpha_rna_first_64,alpha_prot_first_64;
   int16_t theta_rna_first_64,theta_prot_first_64;
@@ -374,7 +374,7 @@ void Individual_R_X11::display_regulation( X11Window* win )
     //  Draw each regulation link
     // ---------------
 
-    for (unsigned int i = 0; i < rna->nb_influences(); i++) {
+    for (int i = 0; i < rna->nb_influences(); i++) {
       //if (rna->_protein_list[i] != nullptr) {
         //compute the activity
         if (rna->_enhancing_coef_list[i] > 0) {
@@ -416,7 +416,7 @@ void Individual_R_X11::display_regulation( X11Window* win )
               max_merged_activator_activity)
             max_merged_activator_activity = merged_activity;
           if (merged_activity <
-              min_activator_activity)
+                  min_merged_activator_activity)
             min_merged_activator_activity = merged_activity;
         }
         if (merged_activity < 0) {
@@ -446,52 +446,52 @@ void Individual_R_X11::display_regulation( X11Window* win )
     // Alpha : angles from OriC (in degrees)
     // Theta : angles on the trigonometric circle (in degrees)
     // nb_sect : "length" in degrees of the arc to be drawn
-    alpha_rna_first = (int16_t) round(
-        360 * ((double) rna->first_transcribed_pos() / (double) genome_length));
-    theta_rna_first = std::fmod(90 - alpha_rna_first, 360);
+    //alpha_rna_first = (int16_t) round(
+    //    360 * ((double) rna->first_transcribed_pos() / (double) genome_length));
+    //theta_rna_first = std::fmod(90 - alpha_rna_first, 360);
 
     // These are the same as above but with a higher precision (1/64 degrees)
     alpha_rna_first_64 = (int16_t) round(64 * 360 *
                                          ((double) rna->first_transcribed_pos() /
                                           (double) genome_length));
-    theta_rna_first_64 = std::fmod(64 * 90 - alpha_rna_first_64, 64 * 360);
+    theta_rna_first_64 = (int16_t ) std::fmod(64 * 90 - alpha_rna_first_64, 64 * 360);
 
-    pos_rna_x = (win->width() / 2.0) +
-                (cos((theta_rna_first_64 / (64 * 180.0) * M_PI)) * diam / 2.0);
-    pos_rna_y = (win->height() / 2.0) -
-                (sin((theta_rna_first_64 / (64 * 180.0) * M_PI)) * diam / 2.0);
+    pos_rna_x = (int16_t ) ((win->width() / 2.0) +
+                (cos((theta_rna_first_64 / (64 * 180.0) * M_PI)) * diam / 2.0));
+    pos_rna_y = (int16_t ) ((win->height() / 2.0) -
+                (sin((theta_rna_first_64 / (64 * 180.0) * M_PI)) * diam / 2.0));
 
     // ---------------
     //  Draw each regulation link
     // ---------------
-    for (unsigned int i = 0; i < rna->_nb_influences; i++) {
+    for (int i = 0; i < rna->_nb_influences; i++) {
       Protein_R* prot = rna->_protein_list[i];
 
       if (!(prot->is_signal())) {
-        alpha_prot_first = (int16_t) round(360 *
-                                           ((double) prot->first_translated_pos() /
-                                            (double) genome_length));
-        theta_prot_first = std::fmod(90 - alpha_prot_first, 360);
+        //alpha_prot_first = (int16_t) round(360 *
+        //                                   ((double) prot->first_translated_pos() /
+        //                                    (double) genome_length));
+        //theta_prot_first = std::fmod(90 - alpha_prot_first, 360);
 
         alpha_prot_first_64 = (int16_t) round(64 * 360 *
                                               ((double) prot->first_translated_pos() /
                                                (double) genome_length));
-        theta_prot_first_64 = std::fmod(64 * 90 - alpha_prot_first_64,
+        theta_prot_first_64 = (int16_t ) std::fmod(64 * 90 - alpha_prot_first_64,
                                         64 * 360);
 
-        pos_prot_x = (win->width() / 2.0) +
+        pos_prot_x = (int16_t ) ((win->width() / 2.0) +
                      (cos((theta_prot_first_64 / (64 * 180.0) * M_PI)) *
                       diam /
-                      2.0);
-        pos_prot_y = (win->height() / 2.0) -
+                      2.0));
+        pos_prot_y = (int16_t ) ((win->height() / 2.0) -
                      (sin((theta_prot_first_64 / (64 * 180.0) * M_PI)) *
                       diam /
-                      2.0);
+                      2.0));
       }
       else {
         nb_signals += 1;
-        pos_prot_x = (win->width() / 10.0) * (((Protein_R*)prot)->get_id()+1);
-        pos_prot_y = (win->height() * 0.9);
+        pos_prot_x = (int16_t ) ((win->width() / 10.0) * (((Protein_R*)prot)->get_id()+1));
+        pos_prot_y = (int16_t ) (win->height() * 0.9);
       }
 
       // compute the color of the link
@@ -501,21 +501,19 @@ void Individual_R_X11::display_regulation( X11Window* win )
       //printf("Merged influence of RNA %ld with prot %ld is %f (%f %f)\n",rna->get_id(),rna->_protein_list[i]->get_id(),
       //       merged_influence,rna->_enhancing_coef_list[i],rna->_operating_coef_list[i]);
       if (merged_influence > 0) {
-            //printf("ONE %lf %lf %d\n", merged_influence, max_merged_activator_activity,(int)((255 * merged_influence) / max_merged_activator_activity));
-            //printf("COLOR : #%02x%02x%02x\n", 0,(int)((255 * merged_influence) / max_merged_activator_activity),0);
-            sprintf(color, "#%02x%02x%02x", 0,
-                    (int) ((255 * merged_influence) /
-                           max_merged_activator_activity), 0);
-
-        }
-        else {
-            //printf("TWO %lf %lf %d\n", merged_influence, max_merged_activator_activity,(int)((255 * merged_influence) / max_merged_activator_activity));
-            //printf("COLOR : #%02x%02x%02x\n", (int)((255 * merged_influence) / max_merged_operator_activity),0,0);
-            sprintf(color, "#%02x%02x%02x", (int) ((255 * merged_influence) /
-                                                   max_merged_operator_activity),
-                    0, 0);
+        //printf("ONE %lf %lf %d\n", merged_influence, max_merged_activator_activity,(int)((255 * merged_influence) / max_merged_activator_activity));
+        //printf("COLOR : #%02x%02x%02x\n", 0,(int)((255 * merged_influence) / max_merged_activator_activity),0);
+        sprintf(color, "#%02x%02x%02x", 0,
+                (int) ((255 * merged_influence) /
+                       max_merged_activator_activity), 0);
+      } else {
+        //printf("TWO %lf %lf %d\n", merged_influence, max_merged_activator_activity,(int)((255 * merged_influence) / max_merged_activator_activity));
+        //printf("COLOR : #%02x%02x%02x\n", (int)((255 * merged_influence) / max_merged_operator_activity),0,0);
+        sprintf(color, "#%02x%02x%02x", (int) ((255 * merged_influence) /
+                                               max_merged_operator_activity),
+                0, 0);
 
-        }
+      }
 
       /*if (rna->_protein_list[i]->is_signal())
         printf("is influenced by signal %ld : %f %f %f %s\n",rna->_protein_list[i]->get_id(),
@@ -574,14 +572,14 @@ void Individual_R_X11::display_phenotype( X11Window* win, const Habitat_R& habit
 
   double dist_temp = 0;
   char* color = new char[8];
-  char* color2 = NULL;
+  //char* color2 = NULL;
   int nb_eval = 0;
   strcpy( color, "#FFFFFF" );
 
   //variable qui sert à stocker du texte à afficher
   char display_string[40];
-  int16_t nb_prot = 0;
-  int16_t life_time = exp_m()->exp_s()->get_nb_indiv_age();
+  //int16_t nb_prot = 0;
+  int16_t life_time = (int16_t ) exp_m()->exp_s()->get_nb_indiv_age();
 
   //set the concentrations of proteins to their initial value
   /*double* concentrations = new double[protein_list_.size()]; // initialise le tableau de concentrations.
@@ -591,8 +589,8 @@ void Individual_R_X11::display_phenotype( X11Window* win, const Habitat_R& habit
   }*/
 
   // compute steps
-  double x_step = 0.8 * win->width() / (double)(life_time * exp_m()->exp_s()->get_nb_degradation_step());
-  double y_step = 0.7 * win->height();
+  //double x_step = 0.8 * win->width() / (double)(life_time * exp_m()->exp_s()->get_nb_degradation_step());
+  //double y_step = 0.7 * win->height();
 
   // Go from an evaluation date to the next
   //int8_t compteur_env   = 0;
@@ -696,7 +694,7 @@ void Individual_R_X11::display_phenotype( X11Window* win, const Habitat_R& habit
   //Draw the evaluation result
   sprintf( display_string, "Mean dist_to_target =  %lf",dist_temp/(double)nb_eval);
 
-  printf("AT %d ID %d Mean dist %lf\n",AeTime::time(),id(),dist_temp/(double)nb_eval);
+  //printf("AT %d ID %d Mean dist %lf\n",AeTime::time(),id(),dist_temp/(double)nb_eval);
   for(Protein_R* prot1 : habitat.signals()) {
     prot1->set_concentration(0.0);
   }
diff --git a/src/libaevol/raevol/PhenotypicTargetHandler_R.cpp b/src/libaevol/raevol/PhenotypicTargetHandler_R.cpp
index 4f3a1347d101fb31058a4fea3b8c75e139f26a42..657eac8038387c4a8797e0f10fe5bf1cc7388bec 100644
--- a/src/libaevol/raevol/PhenotypicTargetHandler_R.cpp
+++ b/src/libaevol/raevol/PhenotypicTargetHandler_R.cpp
@@ -154,7 +154,7 @@ void PhenotypicTargetHandler_R::ApplyVariation() {
 
 //      printf("Nb env targers %d\n",nb_old_env);
 
-      for (int i = 0; i < phenotypic_targets_.size(); i++) {
+      for (int i = 0; i < (int) phenotypic_targets_.size(); i++) {
 
 //        printf("Target at %d is %d\n",i,phenotypic_targets_[i]->get_id());
         list_of_old_target_id[i] = phenotypic_targets_[i]->get_id();
@@ -221,7 +221,7 @@ void PhenotypicTargetHandler_R::ApplyVariation() {
 
       //printf("Computing env changed or not\n");
 
-      if (nb_old_env == phenotypic_targets_.size()) {
+      if (nb_old_env == (int) phenotypic_targets_.size()) {
         i=0;
         for (PhenotypicTarget_R* target : phenotypic_targets_) {
           if (list_of_old_target_id[i] != target->get_id()) {
@@ -268,7 +268,7 @@ void PhenotypicTargetHandler_R::print_geometric_areas() {
   }
   printf("\n");
 
-  for (int16_t i = 0; i < phenotypic_target_models_.size() ; i++) {
+  for (int16_t i = 0; i < (int) phenotypic_target_models_.size() ; i++) {
     area = phenotypic_target_models_.at(i)->fuzzy()->get_geometric_area();
     printf("Entire geometric area of the phenotypic target %d: %f\n", i, area);
     printf("Signal of env %d is ",i);
@@ -414,7 +414,7 @@ void PhenotypicTargetHandler_R::load(gzFile backup_file) {
   int16_t nb_env = 0;
   gzread(backup_file, &nb_env, sizeof(nb_env));
   id = 0;
-  PhenotypicTarget_R* env_to_add = NULL;
+  //PhenotypicTarget_R* env_to_add = NULL;
 
   phenotypic_targets_.resize(nb_env);
 
@@ -514,7 +514,7 @@ void PhenotypicTargetHandler_R::BuildPhenotypicTargetModel( int16_t id) {
 
   // Compute areas (total and by feature)
   phenotypic_target->ComputeArea();
-  double area = phenotypic_target->fuzzy()->get_geometric_area();
+  //double area = phenotypic_target->fuzzy()->get_geometric_area();
     //printf("Entire geometric area of the phenotypic target %d: %f\n", id, area);
 
   //Add its signals to the env
@@ -557,7 +557,7 @@ void PhenotypicTargetHandler_R::InitPhenotypicTargets(int16_t nb_indiv_age) {
   phenotypic_targets_.clear();
   //printf("Taille de l'habitat après le clear dans initialize... : %d\n", phenotypic_targets_.size());
   //phenotypic_targets_.resize(nb_indiv_age);
-  PhenotypicTarget_R* env_to_add;
+  //PhenotypicTarget_R* env_to_add;
   phenotypic_targets_.resize(nb_indiv_age);
 
   for (int i = 0; i < nb_indiv_age; ++i) {
@@ -572,12 +572,12 @@ void PhenotypicTargetHandler_R::InitPhenotypicTargets(int16_t nb_indiv_age) {
 }
 
 void PhenotypicTargetHandler_R::addEnv( int time, int16_t env_id ) {
-  assert(env_id >= 0 && env_id <= phenotypic_target_models_.size());
+  assert(env_id >= 0 && env_id <= (int) phenotypic_target_models_.size());
   phenotypic_targets_[time] = phenotypic_target_models_.at(env_id);
 }
 
 void PhenotypicTargetHandler_R::changeEnv( int16_t ind, int16_t env_id ) {
-  assert(env_id >= 0 && env_id <= phenotypic_target_models_.size());
+  assert(env_id >= 0 && env_id <= (int) phenotypic_target_models_.size());
   phenotypic_targets_.at(ind) = phenotypic_target_models_.at(env_id);
 }
 
diff --git a/src/libaevol/raevol/Rna_R.cpp b/src/libaevol/raevol/Rna_R.cpp
index 7d63ed2b99c82588554a61b030b370b9e33d11af..bb20761f0cba3884268c7b46dced3f30e25159ca 100644
--- a/src/libaevol/raevol/Rna_R.cpp
+++ b/src/libaevol/raevol/Rna_R.cpp
@@ -407,7 +407,7 @@ ProteinConcentration Rna_R::affinity_with_protein( int32_t index, Protein *prote
     }
 
 
-    ProteinConcentration (* binding_matrix)[MAX_QUADON][MAX_CODON] = &(gen_unit_->exp_m()->exp_s()->_binding_matrix);
+    //ProteinConcentration (* binding_matrix)[MAX_QUADON][MAX_CODON] = &(gen_unit_->exp_m()->exp_s()->_binding_matrix);
 #ifdef __SIMD
 #pragma omp simd
 #endif
diff --git a/src/post_treatments/ancestor_stats.cpp b/src/post_treatments/ancestor_stats.cpp
index 12cad904ff40e91eda2c263650cef149aa5754a9..3d57a13c1423307490887dbc0bf2fccfdf335127 100644
--- a/src/post_treatments/ancestor_stats.cpp
+++ b/src/post_treatments/ancestor_stats.cpp
@@ -126,9 +126,9 @@ int main(int argc, char* argv[]) {
                               "for per grid-cell phenotypic target");
   auto phenotypicTargetHandler =
       exp_manager->world()->phenotypic_target_handler();
-  if (not (phenotypicTargetHandler->var_method() == NO_VAR))
-    Utils::ExitWithUsrMsg("sorry, ancestor stats has not yet been implemented "
-                              "for variable phenotypic targets");
+    /*if (not (phenotypicTargetHandler->var_method() == NO_VAR))
+      Utils::ExitWithUsrMsg("sorry, ancestor stats has not yet been implemented "
+                                "for variable phenotypic targets");*/
 
   int64_t backup_step = exp_manager->backup_step();
 
@@ -211,22 +211,96 @@ int main(int argc, char* argv[]) {
   
   
   
-    std::ofstream fitmeta;
-    fitmeta.open("ancestor_composed_fitness.csv",std::ofstream::trunc);
-    fitmeta<<"Generation,EnvID,Composed,Fitness"<<std::endl;
+
 
   // ==================================================
   //  Prepare the initial ancestor and write its stats
   // ==================================================
   GridCell* grid_cell = new GridCell(lineage_file, exp_manager, nullptr);
   auto* indiv = grid_cell->individual();
+#ifdef __REGUL
+    std::ofstream fitmeta;
+    fitmeta.open("ancestor_composed_fitness.csv",std::ofstream::trunc);
+    fitmeta<<"Generation,EnvID,Composed,Fitness"<<std::endl;
+
+    Individual_R *best = dynamic_cast<Individual_R *>(indiv);
+    best->clear_everything_except_dna_and_promoters();
+    //best->do_transcription_translation_folding();
+    //printf("---------> stoch env at %d\n",t);
+    exp_manager->world()->ApplyHabitatVariation();
+
+    //printf("---------> Evaluate indiv at %d\n",t);
+    //best->evaluated_ = false;
+    best->Evaluate();
+    best->compute_statistical_data();
+    printf("Initial fitness     = %e\n", best->fitness());
+    printf("Initial genome size = %" PRId32 "\n", best->total_genome_size());
+
+    double base_metaerror = 0;
+    double base_fitness = 0;
+    int nb_iteration = 100;
+
+    for (int i = 0; i < nb_iteration; i++) {
+        exp_manager->world()->ApplyHabitatVariation();
+        Individual_R* cloned = new Individual_R(*best);
+        cloned->set_grid_cell(exp_manager->world()->grid(0,0));
+
+        cloned->clear_everything_except_dna_and_promoters();
+        cloned->Evaluate();
+
+#pragma omp atomic
+        base_metaerror += cloned->dist_to_target_by_feature(METABOLISM);
+
+#pragma omp atomic
+        base_fitness += cloned->fitness();
+
+#pragma omp critical
+        fitmeta << t0 << ",RANDOM," << i << "," << cloned->dist_to_target_by_feature(METABOLISM) << ","
+                << cloned->fitness() << std::endl;
+
+        //printf("Iteration RANDOM at %d : %d/%d : %lf %e\n", t, i, nb_iteration,
+        //       cloned->dist_to_target_by_feature(METABOLISM), cloned->fitness());
+        delete cloned;
+    }
+
+    best->set_fitness(base_fitness / (double) nb_iteration);
+    best->set_metaerror(base_metaerror / (double) nb_iteration);
+
+    best->compute_statistical_data();
+    best->compute_non_coding();
+    mystats->write_statistics_of_this_indiv(t0,best,nullptr);
+
+    int nb_phenotypic_target_models = (int) dynamic_cast<PhenotypicTargetHandler_R *>(exp_manager->world()->
+            phenotypic_target_handler())->phenotypic_target_models_.size();
+
+    for (int target_id = 0; target_id < nb_phenotypic_target_models; target_id++) {
+        dynamic_cast<PhenotypicTargetHandler_R *>(exp_manager->world()->phenotypic_target_handler())->set_single_env(
+                target_id);
+
+        Individual_R* cloned = new Individual_R(*best);
+        cloned->set_grid_cell(exp_manager->world()->grid(0,0));
+
+        cloned->clear_everything_except_dna_and_promoters();
+        cloned->Evaluate();
+
+        fitmeta << t0 << ",TARGET," << target_id << "," << cloned->dist_to_target_by_feature(METABOLISM) << ","
+                << cloned->fitness() << std::endl;
+
+        //printf("PhenotypicTarget at %d : %d/%d : %lf %e\n", t, target_id, nb_phenotypic_target_models,
+        //       cloned->dist_to_target_by_feature(METABOLISM), cloned->fitness());
+
+        delete cloned;
+    }
+#else
   indiv->Evaluate();
   indiv->compute_statistical_data();
   indiv->compute_non_coding();
 
   mystats->write_statistics_of_this_indiv(t0,indiv, nullptr);
+#endif
 
   // Additional outputs
+#ifndef __REGUL
   write_environment_stats(t0, phenotypicTargetHandler, env_output_file);
   write_terminators_stats(t0, indiv, term_output_file);
   if(phenotypicTargetHandler->phenotypic_target().nb_segments() > 1)
@@ -234,6 +308,7 @@ int main(int argc, char* argv[]) {
     write_zones_stats(t0, indiv, phenotypicTargetHandler, zones_output_file);
   }
   write_operons_stats(t0, indiv, operons_output_file);
+#endif
 
   if (verbose) {
     printf("Initial fitness     = %f\n", indiv->fitness());
@@ -393,7 +468,7 @@ int main(int argc, char* argv[]) {
     }
 
     // 3) All the mutations have been replayed, we can now evaluate the new individual
-    
+
         indiv->evaluated_ = false;
         indiv->Evaluate();
 
@@ -429,8 +504,6 @@ int main(int argc, char* argv[]) {
         int16_t grid_width  = exp_manager->world()->width();
         int16_t grid_height = exp_manager->world()->height();
 
-        int32_t selection_scope_x_ = exp_manager->sel()->selection_scope_x();
-        int32_t selection_scope_y_ = exp_manager->sel()->selection_scope_y();
 
         fitness_sum_local_tab_ = new double*[fitness_function_scope_x_*fitness_function_scope_y_];
         for (int tab_id = 0; tab_id < fitness_function_scope_x_*fitness_function_scope_y_; tab_id++)
@@ -469,28 +542,100 @@ int main(int argc, char* argv[]) {
         double composed_fitness = 0;
         for (int env_id = 0; env_id < number_of_phenotypic_target_models; env_id++) {
             composed_fitness +=  dynamic_cast<Individual_R*>(exp_manager->world()->indiv_at(cur_x, cur_y))->fitness(env_id) / fitness_sum_tab_[env_id];
-            fitmeta<<t0<<","<<env_id<<","<<"0"<<","<<dynamic_cast<Individual_R*>(exp_manager->world()->indiv_at(cur_x, cur_y))->fitness(env_id) / fitness_sum_tab_[env_id]<<std::endl;
+            fitmeta<<time()<<","<<env_id<<","<<"0"<<","<<dynamic_cast<Individual_R*>(exp_manager->world()->indiv_at(cur_x, cur_y))->fitness(env_id) / fitness_sum_tab_[env_id]<<std::endl;
         }
         composed_fitness/=number_of_phenotypic_target_models;
-        fitmeta<<t0<<","<<"-1"<<","<<"1"<<","<<composed_fitness<<std::endl;
+        fitmeta<<time()<<","<<"-1"<<","<<"1"<<","<<composed_fitness<<std::endl;
     } else if (fitness_function_ == FITNESS_LOCAL_SUM) {
         double composed_fitness = 0;
         for (int env_id = 0; env_id < number_of_phenotypic_target_models; env_id++) {
             composed_fitness +=  dynamic_cast<Individual_R*>(exp_manager->world()->indiv_at(cur_x, cur_y))->fitness(env_id) / fitness_sum_local_tab_[0][env_id];
-            fitmeta<<t0<<","<<env_id<<","<<"0"<<","<<dynamic_cast<Individual_R*>(exp_manager->world()->indiv_at(cur_x, cur_y))->fitness(env_id) / fitness_sum_local_tab_[0][env_id]<<std::endl;
+            fitmeta<<time()<<","<<env_id<<","<<"0"<<","<<dynamic_cast<Individual_R*>(exp_manager->world()->indiv_at(cur_x, cur_y))->fitness(env_id) / fitness_sum_local_tab_[0][env_id]<<std::endl;
         }
         composed_fitness/=number_of_phenotypic_target_models;
-        fitmeta<<t0<<","<<"-1"<<","<<"1"<<","<<composed_fitness<<std::endl;
+        fitmeta<<time()<<","<<"-1"<<","<<"1"<<","<<composed_fitness<<std::endl;
     }
 #endif
 
+#ifdef __REGUL
+    if (fitness_function_ == FITNESS_EXP) {
+        Individual_R *best = dynamic_cast<Individual_R *>(indiv);
+        best->clear_everything_except_dna_and_promoters();
+        //best->do_transcription_translation_folding();
+        //printf("---------> stoch env at %d\n",t);
+        exp_manager->world()->ApplyHabitatVariation();
+
+        //printf("---------> Evaluate indiv at %d\n",t);
+        //best->evaluated_ = false;
+        best->Evaluate();
+        best->compute_statistical_data();
+
+        double base_metaerror = 0;
+        double base_fitness = 0;
+        int nb_iteration = 100;
+
+        for (int i = 0; i < nb_iteration; i++) {
+            exp_manager->world()->ApplyHabitatVariation();
+            Individual_R* cloned = new Individual_R(*best);
+            cloned->set_grid_cell(exp_manager->world()->grid(0,0));
+
+            cloned->clear_everything_except_dna_and_promoters();
+            cloned->Evaluate();
+
+#pragma omp atomic
+            base_metaerror += cloned->dist_to_target_by_feature(METABOLISM);
+
+#pragma omp atomic
+            base_fitness += cloned->fitness();
+
+#pragma omp critical
+            fitmeta << time() << ",RANDOM," << i << "," << cloned->dist_to_target_by_feature(METABOLISM) << ","
+                    << cloned->fitness() << std::endl;
+
+            //printf("Iteration RANDOM at %d : %d/%d : %lf %e\n", t, i, nb_iteration,
+            //       cloned->dist_to_target_by_feature(METABOLISM), cloned->fitness());
+            delete cloned;
+        }
+
+        best->set_fitness(base_fitness / (double) nb_iteration);
+        best->set_metaerror(base_metaerror / (double) nb_iteration);
+
+        best->compute_statistical_data();
+        best->compute_non_coding();
+        mystats->write_statistics_of_this_indiv(time(),best,rep);
+
+        int nb_phenotypic_target_models = (int) dynamic_cast<PhenotypicTargetHandler_R *>(exp_manager->world()->
+                phenotypic_target_handler())->phenotypic_target_models_.size();
+
+        for (int target_id = 0; target_id < nb_phenotypic_target_models; target_id++) {
+            dynamic_cast<PhenotypicTargetHandler_R *>(exp_manager->world()->phenotypic_target_handler())->set_single_env(
+                    target_id);
+
+            Individual_R* cloned = new Individual_R(*best);
+            cloned->set_grid_cell(exp_manager->world()->grid(0,0));
+
+            cloned->clear_everything_except_dna_and_promoters();
+            cloned->Evaluate();
+
+            fitmeta << time() << ",TARGET," << target_id << "," << cloned->dist_to_target_by_feature(METABOLISM) << ","
+                    << cloned->fitness() << std::endl;
+
+            //printf("PhenotypicTarget at %d : %d/%d : %lf %e\n", t, target_id, nb_phenotypic_target_models,
+            //       cloned->dist_to_target_by_feature(METABOLISM), cloned->fitness());
+
+            delete cloned;
+        }
+    }
+
+#else
     indiv->Reevaluate();
     indiv->compute_statistical_data();
     indiv->compute_non_coding();
 
     mystats->write_statistics_of_this_indiv(time(),indiv, rep);
-
+#endif
     // Additional outputs
+#ifndef __REGUL
     write_environment_stats(time(), phenotypicTargetHandler, env_output_file);
     write_terminators_stats(time(), indiv, term_output_file);
     if(phenotypicTargetHandler->phenotypic_target().nb_segments() > 1) {
@@ -498,6 +643,7 @@ int main(int argc, char* argv[]) {
                         zones_output_file);
     }
     write_operons_stats(time(), indiv, operons_output_file);
+#endif
 
     if (verbose) printf(" OK\n");
 
diff --git a/src/post_treatments/lineage.cpp b/src/post_treatments/lineage.cpp
index c84c817caf78a601ca8a84f12fc8ddf9f57f4738..04ed81286617bea7935e0eb1542a774c5990800d 100644
--- a/src/post_treatments/lineage.cpp
+++ b/src/post_treatments/lineage.cpp
@@ -141,11 +141,7 @@ int main(int argc, char** argv) {
   // etc.
   //
 
-  #ifdef __REGUL
-    sprintf(tree_file_name,"tree/tree_" TIMESTEP_FORMAT ".rae", t_end);
-  #else
     sprintf(tree_file_name,"tree/tree_" TIMESTEP_FORMAT ".ae", t_end);
-  #endif
 
   tree = new Tree(exp_manager, tree_file_name);
 
@@ -172,12 +168,15 @@ int main(int argc, char** argv) {
       // No index nor rank was given in the command line.
       // By default, we construct the lineage of the best individual, the rank of which
       // is simply the number of individuals in the population.
-      final_indiv_rank = exp_manager->nb_indivs();
+        reports[t_end - t0 - 1] =
+                new ReplicationReport(*(tree->report_by_index(t_end, exp_manager->world()->indiv_at(0,0)->id())));
+    } else {
+        reports[t_end - t0 - 1] =
+                new ReplicationReport(*(tree->report_by_rank(t_end, final_indiv_rank)));
     }
 
     // Retrieve the replication report of the individual of interest (at t_end)
-    reports[t_end - t0 - 1] =
-        new ReplicationReport(*(tree->report_by_rank(t_end, final_indiv_rank)));
+    final_indiv_rank = reports[t_end - t0 - 1]->rank();
     final_indiv_index = reports[t_end - t0 - 1]->id();
 
     indices[t_end - t0]  = final_indiv_index;
@@ -195,15 +194,9 @@ int main(int argc, char** argv) {
   // =======================
   char output_file_name[101];
 
-  #ifdef __REGUL
-    snprintf(output_file_name, 100,
-        "lineage-b" TIMESTEP_FORMAT "-e" TIMESTEP_FORMAT "-i%" PRId32 "-r%" PRId32 ".rae",
-        t0, t_end, final_indiv_index, final_indiv_rank);
-  #else
     snprintf(output_file_name, 100,
         "lineage-b" TIMESTEP_FORMAT "-e" TIMESTEP_FORMAT "-i%" PRId32 "-r%" PRId32 ".ae",
         t0, t_end, final_indiv_index, final_indiv_rank);
-  #endif
 
   gzFile lineage_file = gzopen(output_file_name, "w");
   if (lineage_file == nullptr) {
@@ -245,11 +238,7 @@ int main(int argc, char** argv) {
       // Change the tree file
       delete tree;
 
-      #ifdef __REGUL
-        sprintf(tree_file_name,"tree/tree_" TIMESTEP_FORMAT ".rae", t);
-      #else
         sprintf(tree_file_name,"tree/tree_" TIMESTEP_FORMAT ".ae", t);
-      #endif
       AeTime::set_time(AeTime::time()-tree_step);
 
       tree = new Tree(exp_manager, tree_file_name);
@@ -341,8 +330,8 @@ int main(int argc, char** argv) {
     // Write the replication report of the ancestor for current generation
     if (verbose) {
       printf("Writing the replication report for t= %" PRId64
-             " (built from indiv %" PRId32 " at t= %" PRId64 ")\n",
-             t, indices[i], t-1);
+             " (built from indiv %" PRId32 " %" PRId32" at t= %" PRId64 ")\n",
+             t, indices[i], indices[i+1], t-1);
     }
     of << t << " : " << reports[i]->parent_id() << "  " << reports[i]->id() << std::endl;
     reports[i]->write_to_tree_file(lineage_file);