Mentions légales du service

Skip to content
Snippets Groups Projects
create_eps.cpp 43.7 KiB
Newer Older
// ****************************************************************************
//          Aevol - An in silico experimental evolution platform
// ****************************************************************************
// Copyright: See the AUTHORS file provided with the package or <www.aevol.fr>
// Web: http://www.aevol.fr/
// E-mail: See <http://www.aevol.fr/contact/>
// Original Authors : Guillaume Beslon, Carole Knibbe, David Parsons
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 2 of the License, or
// (at your option) any later version.
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.
Vincent Liard's avatar
Vincent Liard committed
// ****************************************************************************
// ============================================================================
//                                   Includes
// ============================================================================
#include <cinttypes>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cerrno>
#include <getopt.h>
#include <sys/stat.h>  // for the permission symbols used with mkdir

#include "aevol.h"
using namespace aevol;
//  * option --color ?
//  * Raevol-specific output (EPS file with the network) ?

// Command-line option variables
static int64_t timestep = -1;
static int32_t indiv_index = -1;
static int32_t indiv_rank = -1;
// Helper functions
void interpret_cmd_line_options(int argc, char* argv[]);
// The height of each triangle is proportional to the product c*m,
// where c is the concentration of the protein and m its intrinsic efficacy
// (depending on its aminoacid sequence).
// In the case of Raevol, the concentration used here is the final one,
// i.e. the one reached after all the time steps of the lifetime.
// If a coding sequence has several promoters, only one triangle is drawn.
void draw_triangles(Individual* indiv, const PhenotypicTarget& target,
                    char* directoryName);


// In the case of Raevol, the profile is drawn using the final concentrations
// of the proteins, i.e. the ones reached after all the time steps of the
// lifetime.
void draw_pos_neg_profiles(Individual* indiv, const PhenotypicTarget& target,
                           char* directoryName);


// In the case of Raevol, the phenotype is drawn using the final concentrations
// of the proteins, i.e. the ones reached after all the time steps of the
// lifetime.
void draw_phenotype(Individual* indiv, const PhenotypicTarget& target,
                    char* directoryName);
// The chromosome is drawn as a circle. Coding sequences are represented by
// arcs starting at the start codon and ending at the stop codon.
// Coding sequences that are on the leading (resp. lagging) strand are
// drawn outside (resp. inside) the circle. If a coding sequence has several
// promoters, only one arc is drawn.
void draw_genetic_unit_with_CDS(GeneticUnit* gen_unit, char* directoryName);
Bérénice Batut's avatar
Bérénice Batut committed

// The chromosome is drawn as a circle. Transcribed sequences are represented
// by arcs starting at the first transcribed position and ending at the last
// transcribed position. mRNAs that are on the leading (resp. lagging) strand
// are drawn outside (resp. inside) the circle.
// mRNAs that include at least one coding sequence are black,
// the others are gray.
void draw_genetic_unit_with_mRNAs(GeneticUnit* gen_unit, char* directoryName);
int main(int argc, char* argv[]) {
  interpret_cmd_line_options(argc, argv);
  printf("Creating eps files for generation %" PRId64 "...\n", timestep);
  // =================================================================
  //                       Read the backup file
  // =================================================================
  Individual*  indiv;
Bérénice Batut's avatar
Bérénice Batut committed
  // Load the simulation
  ExpManager* exp_manager = new ExpManager();
  exp_manager->load(timestep, true, false);
  if (indiv_index == -1 && indiv_rank == -1) {
    indiv = exp_manager->best_indiv();
    // TODO <david.parsons@inria.fr> tmp disabled
//     if (indiv_rank != -1) {
//       indiv = new Individual(*exp_manager->indiv_by_rank(indiv_rank), false);
//       indiv = new Individual(*exp_manager->indiv_by_id(indiv_index), false);
Bérénice Batut's avatar
Bérénice Batut committed
  }
  // The constructor of the exp_manager has read the genomes of the individuals
  // and located their promoters, but has not performed the translation nor the
  // phenotype computation. We must do it now.
  // However, as the individuals in the backups are sorted, we don't need to
  // evaluate all the individuals, only the one we are interested in
  indiv->Evaluate();
  // =================================================================
  // =================================================================
  snprintf(directory_name, 63, "analysis-generation_" TIMESTEP_FORMAT,
           timestep);

  // Check whether the directory already exists and is writable
  if (access(directory_name, F_OK) == 0) {
//       struct stat status;
//       stat(directory_name, &status);
//       if (status.st_mode & S_IFDIR) cout << "The directory exists." << endl;
//       else cout << "This path is a file." << endl;

    if (access(directory_name, X_OK | W_OK) != 0) {
      fprintf(stderr, "Error: cannot enter or write in directory %s.\n", directory_name);
    // Create the directory with permissions : rwx r-x r-x
    if (mkdir(directory_name,
              S_IRWXU | S_IRGRP | S_IXGRP | S_IROTH | S_IXOTH) != 0) {
      fprintf(stderr, "Error: cannot create directory %s.\n", directory_name);
  // =================================================================
  // =================================================================

  GeneticUnit*  indiv_main_genome = &indiv->genetic_unit_nonconst(0);
  printf("Creating the EPS file with the triangles of the chosen individual... ");
  draw_triangles(indiv, indiv->phenotypic_target(), directory_name);
  printf("OK\n");
  printf("Creating the EPS file with the positive and negatives profiles of the chosen individual... ");
  fflush(stdout);
  draw_pos_neg_profiles(indiv, indiv->phenotypic_target(), directory_name);
  printf("OK\n");
  printf("Creating the EPS file with the phenotype of the chosen individual... ");
  draw_phenotype(indiv, indiv->phenotypic_target(), directory_name);
  printf("OK\n");
  printf("Creating the EPS file with the CDS of the chosen individual... ");
  draw_genetic_unit_with_CDS(indiv_main_genome, directory_name);
  printf("OK\n");
  printf("Creating the EPS file with the mRNAs of the chosen individual... ");
  draw_genetic_unit_with_mRNAs(indiv_main_genome, directory_name);
  printf("OK\n");
Bérénice Batut's avatar
Bérénice Batut committed
  delete exp_manager;
void draw_triangles(Individual* indiv, const PhenotypicTarget& target,
                    char* directoryName) {
  const uint8_t bbsize = 200;  // a4 paper: 595*842
  double margin = 0.1;
  double scalex = 0.8*(1 - 2*margin);
  double scaley = 0.4*(1 - 2*margin);

  snprintf(filename, 127, "%s/best_triangles.eps", directoryName);
  FILE * drawingfile = fopen(filename, "w");
  fprintf(drawingfile, "%%!PS-Adobe-3.0 EPSF-3.0\n");
  fprintf(drawingfile, "%%%%BoundingBox: 0 0 %d %d\n", bbsize, bbsize);
  fprintf(drawingfile, "%d %d scale\n", bbsize, bbsize);
  fprintf(drawingfile, "%d %d 8 [100 0 0 -100 0 100]\n",bbsize, bbsize);
  fprintf(drawingfile, "{currentfile 3 100 mul string readhexstring pop} bind\n");


  // -----------------------------
  //  paint neutral zones in grey
  // -----------------------------
  if (target.nb_segments() > 1)
    int16_t nb_segments = target.nb_segments();
    PhenotypicSegment** segments = target.segments();
    for (int16_t i = 0 ; i < nb_segments ; i++)
      if (segments[i]->feature == NEUTRAL)
        fprintf(drawingfile, "%lf 0 moveto\n", margin + scalex * segments[i]->start);
        fprintf(drawingfile, "%lf 1 lineto\n", margin + scalex * segments[i]->start);
        fprintf(drawingfile, "%lf 1 lineto\n", margin + scalex * segments[i]->stop);
        fprintf(drawingfile, "%lf 0 lineto\n", margin + scalex * segments[i]->stop);
        fprintf(drawingfile, "closepath\n");
        fprintf(drawingfile, "0.8 setgray\n");
        fprintf(drawingfile, "fill\n");
        fprintf(drawingfile, "0 setgray\n");
  // -----------
  //  draw axis
  // -----------

  double arrowsize = 0.03;
  double arrowangle = 3.14/6;

  fprintf(drawingfile, "0.001 setlinewidth\n");
  fprintf(drawingfile, "0 0 0 setrgbcolor\n");
  fprintf(drawingfile, "%lf %lf moveto\n", margin/2, 0.5);
  fprintf(drawingfile, "%lf %lf lineto\n", 1-margin, 0.5);
  fprintf(drawingfile, "stroke\n");
  fprintf(drawingfile, "%lf %lf moveto\n", 1-margin, 0.5);
  fprintf(drawingfile, "%lf %lf lineto\n", 1-margin-arrowsize*cos(arrowangle), 0.5 + arrowsize*sin(arrowangle));
  fprintf(drawingfile, "stroke\n");
  fprintf(drawingfile, "%lf %lf moveto\n", 1-margin, 0.5);
  fprintf(drawingfile, "%lf %lf lineto\n", 1-margin-arrowsize*cos(arrowangle), 0.5 - arrowsize*sin(arrowangle));
  fprintf(drawingfile, "stroke\n");
  fprintf(drawingfile, "%lf %lf moveto\n", margin, margin/2);
  fprintf(drawingfile, "%lf %lf lineto\n", margin, 1-margin);
  fprintf(drawingfile, "%lf %lf moveto\n", margin, 1-margin);
  fprintf(drawingfile, "%lf %lf lineto\n", margin-arrowsize*sin(arrowangle), 1 - margin - arrowsize*cos(arrowangle));
  fprintf(drawingfile, "stroke\n");
  fprintf(drawingfile, "%lf %lf moveto\n", margin, 1-margin);
  fprintf(drawingfile, "%lf %lf lineto\n", margin+arrowsize*sin(arrowangle), 1 - margin - arrowsize*cos(arrowangle));
  fprintf(drawingfile, "stroke\n");
  fprintf(drawingfile, "[0.02 0.02] 0 setdash\n");
  fprintf(drawingfile, "%lf %lf moveto\n", margin, 0.5 + 1.0*scaley);
  fprintf(drawingfile, "%lf %lf lineto\n", 1-margin, 0.5 + 1.0*scaley);
  fprintf(drawingfile, "stroke\n");
  // ----------------
  //  draw triangles
  // ----------------
  fprintf(drawingfile,"[ ] 0 setdash\n");
  for (auto& gu: indiv->genetic_unit_list_nonconst()) { // should use const version
    for (const auto& prot: gu.protein_list(LEADING)) {
      h = prot.height() * prot.concentration();
      fprintf(drawingfile, "%lf %lf moveto\n", margin, 0.5);
      fprintf(drawingfile, "%lf %lf lineto\n", margin + scalex*(prot.mean() - prot.width()), 0.5);
      fprintf(drawingfile, "%lf %lf lineto\n", margin + scalex*(prot.mean()), 0.5 + scaley*(h));
      fprintf(drawingfile, "%lf %lf lineto\n", margin + scalex*(prot.mean() + prot.width()), 0.5);
      fprintf(drawingfile, "%lf %lf moveto\n", margin + scalex*(1), 0.5);
      fprintf(drawingfile, "stroke\n");
    for (const auto& prot: gu.protein_list(LAGGING)) {
      h = prot.height() * prot.concentration();
      fprintf(drawingfile, "%lf %lf moveto\n", margin, 0.5);
      fprintf(drawingfile, "%lf %lf lineto\n", margin + scalex*(prot.mean() - prot.width()), 0.5);
      fprintf(drawingfile, "%lf %lf lineto\n", margin + scalex*(prot.mean()), 0.5 + scaley*(h));
      fprintf(drawingfile, "%lf %lf lineto\n", margin + scalex*(prot.mean() + prot.width()), 0.5);
      fprintf(drawingfile, "%lf %lf moveto\n", margin + scalex*(1), 0.5);
      fprintf(drawingfile, "stroke\n");
  fprintf(drawingfile,"%%%%EOF\n");
void draw_pos_neg_profiles(Individual* indiv, const PhenotypicTarget& target,
                           char* directoryName) {
  const uint8_t bbsize = 200;  // a4 paper: 595*842
  double margin = 0.1;
  double scale = 0.8*(1 - 2*margin);
  snprintf(filename, 127, "%s/best_pos_neg_profiles.eps", directoryName);
  FILE * drawingfile = fopen(filename, "w");
  fprintf(drawingfile, "%%!PS-Adobe-3.0 EPSF-3.0\n");
  fprintf(drawingfile, "%%%%BoundingBox: 0 0 %d %d\n", bbsize, bbsize);
  fprintf(drawingfile, "%d %d scale\n", bbsize, bbsize);
  fprintf(drawingfile, "%d %d 8 [100 0 0 -100 0 100]\n",bbsize, bbsize);
  fprintf(drawingfile, "{currentfile 3 100 mul string readhexstring pop} bind\n");
  // -----------------------------
  //  paint neutral zones in grey
  // -----------------------------
  if (target.nb_segments() > 1)
    int16_t nb_segments = target.nb_segments();
    PhenotypicSegment** segments = target.segments();
    for (int16_t i = 0 ; i < nb_segments ; i++)
      if (segments[i]->feature == NEUTRAL)
        fprintf(drawingfile, "%lf 0 moveto\n", margin + scale * segments[i]->start);
        fprintf(drawingfile, "%lf 1 lineto\n", margin + scale * segments[i]->start);
        fprintf(drawingfile, "%lf 1 lineto\n", margin + scale * segments[i]->stop);
        fprintf(drawingfile, "%lf 0 lineto\n", margin + scale * segments[i]->stop);
        fprintf(drawingfile, "closepath\n");
        fprintf(drawingfile, "0.8 setgray\n");
        fprintf(drawingfile, "fill\n");
        fprintf(drawingfile, "0 setgray\n");
  // -----------
  //  draw axis
  // -----------

  double arrowsize = 0.03;
  double arrowangle = 3.14/6;

  fprintf(drawingfile, "0.001 setlinewidth\n");
  fprintf(drawingfile, "0 0 0 setrgbcolor\n");
  fprintf(drawingfile, "%lf %lf moveto\n", margin/2, 0.5);
  fprintf(drawingfile, "%lf %lf lineto\n", 1-margin, 0.5);
  fprintf(drawingfile, "stroke\n");
  fprintf(drawingfile, "%lf %lf moveto\n", 1-margin, 0.5);
  fprintf(drawingfile, "%lf %lf lineto\n", 1-margin-arrowsize*cos(arrowangle), 0.5 + arrowsize*sin(arrowangle));
  fprintf(drawingfile, "stroke\n");
  fprintf(drawingfile, "%lf %lf moveto\n", 1-margin, 0.5);
  fprintf(drawingfile, "%lf %lf lineto\n", 1-margin-arrowsize*cos(arrowangle), 0.5 - arrowsize*sin(arrowangle));
  fprintf(drawingfile, "stroke\n");
  fprintf(drawingfile, "%lf %lf moveto\n", margin, margin/2);
  fprintf(drawingfile, "%lf %lf lineto\n", margin, 1-margin);
  fprintf(drawingfile, "%lf %lf moveto\n", margin, 1-margin);
  fprintf(drawingfile, "%lf %lf lineto\n", margin-arrowsize*sin(arrowangle), 1 - margin - arrowsize*cos(arrowangle));
  fprintf(drawingfile, "stroke\n");
  fprintf(drawingfile, "%lf %lf moveto\n", margin, 1-margin);
  fprintf(drawingfile, "%lf %lf lineto\n", margin+arrowsize*sin(arrowangle), 1 - margin - arrowsize*cos(arrowangle));
  fprintf(drawingfile, "stroke\n");


  // -----------------------
  //  draw positive profile
  // -----------------------

  fprintf(drawingfile,"[ ] 0 setdash\n");
  fprintf(drawingfile, "0.002 setlinewidth\n");
  fprintf(drawingfile, "%lf %lf moveto\n", margin, 0.5);
  if (indiv->exp_m()->exp_s()->get_fuzzy_flavor() == 0)
    for (const auto& p: ((Fuzzy*)indiv->phenotype_activ())->points())
      fprintf(drawingfile, "%lf %lf lineto\n", margin + scale * p.x, 0.5 + scale * p.y);
  else
    for (int i=0; i < ((HybridFuzzy*)indiv->phenotype_activ())->get_pheno_size(); i++) {
      int xi = (int) ( i / ((HybridFuzzy*)indiv->phenotype_activ())->get_pheno_size());
      fprintf(drawingfile, "%lf %lf lineto\n", margin +
                                               scale * xi, 0.5 + scale *
                                               ((HybridFuzzy*) indiv->phenotype_activ())->points()[i]);
  fprintf(drawingfile, "stroke\n" );
  // -----------------------
  //  draw negative profile
  // -----------------------

  fprintf(drawingfile,"[ ] 0 setdash\n");
  fprintf(drawingfile, "0.002 setlinewidth\n");
  fprintf(drawingfile, "%lf %lf moveto\n", margin, 0.5);
  if (indiv->exp_m()->exp_s()->get_fuzzy_flavor() == 0)
    for (const auto& p: ((Fuzzy*)indiv->phenotype_inhib())->points())
      fprintf( drawingfile, "%lf %lf lineto\n", margin + scale * p.x, 0.5 + scale * p.y);
  else
    for (int i=0; i < ((HybridFuzzy*)indiv->phenotype_inhib())->get_pheno_size(); i++) {
      int xi = (int) ( i / ((HybridFuzzy*)indiv->phenotype_inhib())->get_pheno_size());
      fprintf(drawingfile, "%lf %lf lineto\n", margin +
                                               scale * xi, 0.5 + scale *
                                                                 ((HybridFuzzy*) indiv->phenotype_inhib())->points()[i]);
  fprintf( drawingfile, "stroke\n" );
  fprintf(drawingfile,"%%%%EOF\n");
void draw_phenotype(Individual* indiv, const PhenotypicTarget& target,
                    char* directoryName) {
  const uint8_t bbsize = 200;  // a4 paper: 595*842
  double margin = 0.1;
  double scale = 0.8*(1 - 2*margin);
  snprintf(filename, 127, "%s/best_phenotype.eps", directoryName);
  FILE * drawingfile = fopen(filename, "w");
  if (drawingfile == NULL)
    {
      fprintf(stderr, "Error: could not create the file %s\n", filename);
      return;
    }
  fprintf(drawingfile, "%%!PS-Adobe-3.0 EPSF-3.0\n");
  fprintf(drawingfile, "%%%%BoundingBox: 0 0 %d %d\n", bbsize, bbsize);
  fprintf(drawingfile, "%d %d scale\n", bbsize, bbsize);
  fprintf(drawingfile, "%d %d 8 [100 0 0 -100 0 100]\n",bbsize, bbsize);
  fprintf(drawingfile, "{currentfile 3 100 mul string readhexstring pop} bind\n");


  // -----------------------------
  //  paint neutral zones in grey
  // -----------------------------
  if (target.nb_segments() > 1)
    int16_t nb_segments = target.nb_segments();
    PhenotypicSegment** segments = target.segments();
    for (int16_t i = 0 ; i < nb_segments ; i++)
      if (segments[i]->feature == NEUTRAL)
        fprintf(drawingfile, "%lf 0 moveto\n", margin + scale * segments[i]->start);
        fprintf(drawingfile, "%lf 1 lineto\n", margin + scale * segments[i]->start);
        fprintf(drawingfile, "%lf 1 lineto\n", margin + scale * segments[i]->stop);
        fprintf(drawingfile, "%lf 0 lineto\n", margin + scale * segments[i]->stop);
        fprintf(drawingfile, "closepath\n");
        fprintf(drawingfile, "0.8 setgray\n");
        fprintf(drawingfile, "fill\n");
        fprintf(drawingfile, "0 setgray\n");
  // -----------
  //  draw axis
  // -----------

  double arrowsize = 0.03;
  double arrowangle = 3.14/6;

  fprintf(drawingfile, "0.001 setlinewidth\n");
  fprintf(drawingfile, "0 0 0 setrgbcolor\n");
  fprintf(drawingfile, "%lf %lf moveto\n", margin/2, margin);
  fprintf(drawingfile, "%lf %lf lineto\n", 1-margin, margin);
  fprintf(drawingfile, "stroke\n");
  fprintf(drawingfile, "%lf %lf moveto\n", 1-margin, margin);
  fprintf(drawingfile, "%lf %lf lineto\n", 1-margin-arrowsize*cos(arrowangle), margin + arrowsize*sin(arrowangle));
  fprintf(drawingfile, "stroke\n");
  fprintf(drawingfile, "%lf %lf moveto\n", 1-margin, margin);
  fprintf(drawingfile, "%lf %lf lineto\n", 1-margin-arrowsize*cos(arrowangle), margin - arrowsize*sin(arrowangle));
  fprintf(drawingfile, "stroke\n");
  fprintf(drawingfile, "%lf %lf moveto\n", margin, margin/2);
  fprintf(drawingfile, "%lf %lf lineto\n", margin, 1-margin);
  fprintf(drawingfile, "%lf %lf moveto\n", margin, 1-margin);
  fprintf(drawingfile, "%lf %lf lineto\n", margin-arrowsize*sin(arrowangle), 1 - margin - arrowsize*cos(arrowangle));
  fprintf(drawingfile, "stroke\n");
  fprintf(drawingfile, "%lf %lf moveto\n", margin, 1-margin);
  fprintf(drawingfile, "%lf %lf lineto\n", margin+arrowsize*sin(arrowangle), 1 - margin - arrowsize*cos(arrowangle));
  fprintf(drawingfile, "stroke\n");
  fprintf(drawingfile, "[0.02 0.02] 0 setdash\n");
  fprintf(drawingfile, "%lf %lf moveto\n", margin, margin + 1*scale);
  fprintf(drawingfile, "%lf %lf lineto\n", 1-margin, margin + 1*scale);
  fprintf(drawingfile, "stroke\n");
  // ----------------
  //  draw phenotype
  // ----------------
  fprintf( drawingfile,"[ ] 0 setdash\n" );
  fprintf( drawingfile, "0.002 setlinewidth\n" );
  fprintf( drawingfile, "%lf %lf moveto\n", margin, margin);
  if (indiv->exp_m()->exp_s()->get_fuzzy_flavor() == 0)
    for (const auto& p: ((Fuzzy*)indiv->phenotype_activ())->points())
      fprintf(drawingfile, "%lf %lf lineto\n", margin + scale * p.x, margin + scale * p.y);
  else
    for (int i=0; i < ((HybridFuzzy*)indiv->phenotype_activ())->get_pheno_size(); i++) {
      int xi = (int) ( i / ((HybridFuzzy*)indiv->phenotype_activ())->get_pheno_size());
      fprintf(drawingfile, "%lf %lf lineto\n", margin +
                                             scale * xi, margin + scale *
                                                               ((HybridFuzzy*) indiv->phenotype_activ())->points()[i]);
  fprintf( drawingfile, "stroke\n" );
  // ------------------
  //  draw environment
  // ------------------
  fprintf( drawingfile,"[ ] 0 setdash\n" );
  fprintf( drawingfile, "0.001 setlinewidth\n" );
  fprintf( drawingfile, "%lf %lf moveto\n", margin, margin);
  if (indiv->exp_m()->exp_s()->get_fuzzy_flavor() == 0)
    for (const auto& p: ((Fuzzy*)target.fuzzy())->points())
      fprintf(drawingfile, "%lf %lf lineto\n", margin + scale * p.x, margin + scale * p.y);
  else
  for (int i=0; i < ((HybridFuzzy*)target.fuzzy())->get_pheno_size(); i++) {
    int xi = (int) ( i / ((HybridFuzzy*)target.fuzzy())->get_pheno_size());
    fprintf(drawingfile, "%lf %lf lineto\n", margin +
                                             scale * xi, margin + scale *
                                                                  ((HybridFuzzy*) target.fuzzy())->points()[i]);
  fprintf( drawingfile, "stroke\n" );
  fprintf(drawingfile,"%%%%EOF\n");
  fclose(drawingfile);
void draw_genetic_unit_with_CDS(GeneticUnit* gen_unit, char* directoryName) {
  const uint8_t bbsize = 200;  // a4 paper: 595*842
  int32_t gen_length = (gen_unit->dna())->length();
  double r = 0.35;
  double scale = 2*M_PI*r/gen_length;

  snprintf(filename, 127, "%s/best_genome_with_CDS.eps", directoryName);
  FILE * drawingfile = fopen(filename, "w");
  fprintf(drawingfile, "%%!PS-Adobe-3.0 EPSF-3.0\n");
  fprintf(drawingfile, "%%%%BoundingBox: 0 0 %d %d\n", bbsize, bbsize);
  fprintf(drawingfile, "%d %d scale\n", bbsize, bbsize);
  fprintf(drawingfile, "%d %d 8 [100 0 0 -100 0 100]\n",bbsize, bbsize);
  fprintf(drawingfile, "{currentfile 3 100 mul string readhexstring pop} bind\n");
  // -----------
  //  chromosome
  // -----------
  fprintf(drawingfile, "0.001 setlinewidth\n");
  fprintf(drawingfile, "0 0 0 setrgbcolor\n");
  fprintf(drawingfile, "%lf %lf %lf 0 360 arc\n", 0.5, 0.5, r); // arcn = clockwise arc
  fprintf(drawingfile, "stroke\n");

  // -----------
  //  scale
  // -----------

  double scalesize = 0.15;
  fprintf(drawingfile, "%lf %lf moveto\n", 0.5-scalesize/2, 0.5);
  fprintf(drawingfile, "%lf %lf lineto\n", 0.5+scalesize/2, 0.5);
  fprintf(drawingfile, "stroke\n");
  fprintf(drawingfile, "/Helvetica findfont\n");
  fprintf(drawingfile, "0.035 scalefont\n");
  fprintf(drawingfile, "setfont\n");
  fprintf(drawingfile, "newpath\n");
  fprintf(drawingfile, "%lf %lf moveto\n", 0.5-scalesize/3, 0.52);
  fprintf(drawingfile, "(scale : %.0lf bp) show\n", scalesize/scale);

  // -----------
  //  genes
  // -----------

  int32_t first;
  int32_t last;
  int8_t  layer = 0;
  int8_t  outmost_layer = 1;
  int16_t nb_sect;
  int16_t rho;
  int16_t angle;
  bool    sectors_free;
  int16_t alpha_first;
  int16_t alpha_last;
  int16_t theta_first;
  int16_t theta_last;

  // To handle gene overlaps, we remember where we have aldready drawn
  // something, with a precision of 1 degree. We handle up to 100 layers:
  //  - 50 layers inside the circle (lagging strand),
  //  - 50 layers outside the cricle (leading strand).
  // At first, only one layer is created, we create new layers when we
  // need them.
  bool* occupied_sectors[2][50];
  occupied_sectors[LEADING][0] = new bool[360];
  occupied_sectors[LAGGING][0] = new bool[360];
  for (int16_t angle = 0 ; angle < 360 ; angle++)
  {
    occupied_sectors[LEADING][0][angle] = false;
    occupied_sectors[LAGGING][0][angle] = false;
  }


  // printf("LEADING\n");
  for (const auto& prot: gen_unit->protein_list(LEADING)) {
    first = prot.first_translated_pos();
    last = prot.last_translated_pos();
    // h = prot.height() * prot.concentration();
    alpha_first   = (int16_t) round((double)(360 * first) / (double)gen_length);  //  == sect1 == alphaB
    alpha_last    = (int16_t) round((double)(360 * last)  / (double)gen_length);  //  == sect2 == alphaA
    theta_first   = Utils::mod(90 - alpha_first, 360);  //  == tetaB
    theta_last    = Utils::mod(90 - alpha_last, 360);  //   == tetaA
    if (theta_first == theta_last) theta_first = Utils::mod(theta_first + 1, 360);
    nb_sect = Utils::mod(theta_first - theta_last + 1, 360);


    // Outside the circle, look for the inmost layer that has all the sectors between
    // theta_first and theta_last free
    layer = 0;
    sectors_free = false;
    while (! sectors_free)
    {
      sectors_free = true;
      for (rho = 0 ; rho < nb_sect ; rho++)
        if (occupied_sectors[LEADING][layer][Utils::mod(theta_first - rho, 360)])
        {
          sectors_free = false;
          break;
        }
      }
      {
        break; // All the needed sectors are free on the current layer
      }
      else
      {
        layer++;

        if ((layer >= outmost_layer) && (layer < 49))
        {
          // Add a new layer (actually, 2 layers, to maintain the symmetry)
          occupied_sectors[LEADING][outmost_layer] = new bool[360];
          occupied_sectors[LAGGING][outmost_layer] = new bool[360];
          for (int16_t angle = 0 ; angle < 360 ; angle++)
            occupied_sectors[LEADING][outmost_layer][angle] = false;
            occupied_sectors[LAGGING][outmost_layer][angle] = false;
          outmost_layer++;
          break; // A new layer is necessarily free, no need to look further
        }
        {
          // We shall not create a 51th layer, the CDS will be drawn on the
          // layer, probably over another CDS
          break;
        }
    // printf("f %d, l %d, af %d, al %d, tf %d, tl %d, nbsect %d, layer %d\n", first, last, alpha_first, alpha_last, theta_first, theta_last, nb_sect, layer);

    // Mark sectors to be drawn as occupied
    for (rho = 0 ; rho < nb_sect ; rho++)
      occupied_sectors[LEADING][layer][Utils::mod(theta_first - rho, 360)] = true;
    // Mark flanking sectors as occupied
    occupied_sectors[LEADING][layer][Utils::mod(theta_first + 1, 360)] = true;
    occupied_sectors[LEADING][layer][Utils::mod(theta_first - nb_sect, 360)] = true;
    fprintf(drawingfile, "0.018 setlinewidth\n");
    // fprintf(drawingfile, "%lf %lf %lf setrgbcolor\n",  1-(0.8*h/max_height + 0.2), 1-(0.8*h/max_height + 0.2),1-(0.8*h/max_height + 0.2));
    layer++; // index starting at 0 but needed to start at 1
      fprintf(drawingfile, "newpath\n");
      fprintf(drawingfile, "%lf %lf %lf %d %d arc\n", 0.5, 0.5, r + layer*0.02, theta_last, 360);
      fprintf(drawingfile, "stroke\n");
      fprintf(drawingfile, "%lf %lf %lf %d %d arc\n", 0.5, 0.5, r + layer*0.02, 0, theta_first);
      fprintf(drawingfile, "stroke\n");
      fprintf(drawingfile, "newpath\n");
      fprintf(drawingfile, "%lf %lf %lf %d %d arc\n", 0.5, 0.5, r + layer*0.02, theta_last, theta_first);
      fprintf(drawingfile, "stroke\n");
  // printf("LAGGING\n");
  for (const auto& prot: gen_unit->protein_list(LAGGING)) {
    first = prot.first_translated_pos();
    last = prot.last_translated_pos();
    // h = prot.height() * prot.concentration();
    alpha_first   = (int16_t) round((double)(360 * first) / (double)gen_length);
    alpha_last    = (int16_t) round((double)(360 * last)  / (double)gen_length);
    theta_first   = Utils::mod(90 - alpha_first, 360);
    theta_last    = Utils::mod(90 - alpha_last, 360);
    if (theta_first == theta_last) theta_last = Utils::mod(theta_last + 1, 360);
    nb_sect = Utils::mod(theta_last - theta_first + 1, 360);


    // Inside the circle, look for the inmost layer that has all the sectors between
    // theta_first and theta_last free
    layer = 0;
    sectors_free = false;
    while (! sectors_free)
      for (rho = 0 ; rho < nb_sect ; rho++)
        if (occupied_sectors[LAGGING][layer][Utils::mod(theta_first + rho, 360)])
      {
        break; // All the needed sectors are free on the current layer
      }
      else
      {
        layer++;

        if ((layer >= outmost_layer) && (layer < 49))
        {
          // Add a new layer (actually, 2 layers, to maintain the symmetry)
          occupied_sectors[LEADING][outmost_layer] = new bool[360];
          occupied_sectors[LAGGING][outmost_layer] = new bool[360];
          for (angle = 0 ; angle < 360 ; angle++)
          {
            occupied_sectors[LEADING][outmost_layer][angle] = false;
            occupied_sectors[LAGGING][outmost_layer][angle] = false;
          }
          outmost_layer++;
          break; // A new layer is necessarily free, no need to look further
        }
        {
          // We shall not create a 51th layer, the CDS will be drawn on the
          // layer, probably over another CDS
          break;
        }
      }
    }

    // printf("f %d, l %d, af %d, al %d, tf %d, tl %d, nbsect %d, layer %d\n", first, last, alpha_first, alpha_last, theta_first, theta_last, nb_sect, layer);

    // Mark sectors to be drawn as occupied
    for (rho = 0 ; rho < nb_sect ; rho++)
      occupied_sectors[LAGGING][layer][Utils::mod(theta_first + rho, 360)] = true;
    }
    // Mark flanking sectors as occupied
    occupied_sectors[LAGGING][layer][Utils::mod(theta_first - 1, 360)] = true;
    occupied_sectors[LAGGING][layer][Utils::mod(theta_first + nb_sect, 360)] = true;
    fprintf(drawingfile, "0.018 setlinewidth\n");
    // fprintf(drawingfile, "%lf %lf %lf setrgbcolor\n",  1-(0.8*h/max_height + 0.2), 1-(0.8*h/max_height + 0.2),1-(0.8*h/max_height + 0.2));
    layer++; // index starting at 0 but needed to start at 1

      fprintf(drawingfile, "newpath\n");
      fprintf(drawingfile, "%lf %lf %lf %d %d arc\n", 0.5, 0.5, r - layer*0.02, theta_first, 360);
      fprintf(drawingfile, "stroke\n");
      fprintf(drawingfile, "%lf %lf %lf %d %d arc\n", 0.5, 0.5, r - layer*0.02, 0, theta_last);
      fprintf(drawingfile, "stroke\n");
      fprintf(drawingfile, "newpath\n");
      fprintf(drawingfile, "%lf %lf %lf %d %d arc\n", 0.5, 0.5, r - layer*0.02, theta_first, theta_last);
      fprintf(drawingfile, "stroke\n");
  fprintf(drawingfile,"showpage\n");
  fprintf(drawingfile,"%%%%EOF\n");
  for (layer = 0 ; layer < outmost_layer ; layer++)
  {
    delete occupied_sectors[LEADING][layer];
    delete occupied_sectors[LAGGING][layer];
  }

}


void draw_genetic_unit_with_mRNAs(GeneticUnit* gen_unit, char* directoryName) {
  const uint8_t bbsize = 200;  // a4 paper: 595*842
  int32_t gen_length = (gen_unit->dna())->length();
  double r = 0.35;
  double scale = 2*M_PI*r/gen_length;

  snprintf(filename, 127, "%s/best_genome_with_mRNAs.eps", directoryName);
  FILE * drawingfile = fopen(filename, "w");
  fprintf(drawingfile, "%%!PS-Adobe-3.0 EPSF-3.0\n");
  fprintf(drawingfile, "%%%%BoundingBox: 0 0 %d %d\n", bbsize, bbsize);
  fprintf(drawingfile, "%d %d scale\n", bbsize, bbsize);
  fprintf(drawingfile, "%d %d 8 [100 0 0 -100 0 100]\n",bbsize, bbsize);
  fprintf(drawingfile, "{currentfile 3 100 mul string readhexstring pop} bind\n");
  // -----------
  //  chromosome
  // -----------
  fprintf(drawingfile, "0.001 setlinewidth\n");
  fprintf(drawingfile, "0 0 0 setrgbcolor\n");
  fprintf(drawingfile, "%lf %lf %lf 0 360 arc\n", 0.5, 0.5, r); // arcn = clockwise arc
  fprintf(drawingfile, "stroke\n");

  // -----------
  //  scale
  // -----------

  double scalesize = 0.15;
  fprintf(drawingfile, "%lf %lf moveto\n", 0.5-scalesize/2, 0.5);
  fprintf(drawingfile, "%lf %lf lineto\n", 0.5+scalesize/2, 0.5);
  fprintf(drawingfile, "stroke\n");
  fprintf(drawingfile, "/Helvetica findfont\n");
  fprintf(drawingfile, "0.035 scalefont\n");
  fprintf(drawingfile, "setfont\n");
  fprintf(drawingfile, "newpath\n");
  fprintf(drawingfile, "%lf %lf moveto\n", 0.5-scalesize/3, 0.52);
  fprintf(drawingfile, "(scale : %.0lf bp) show\n", scalesize/scale);

  // -----------
  //  mRNAs
  // -----------

  int32_t first;
  int32_t last;
  int8_t  layer = 0;
  int8_t  outmost_layer = 1;
  int16_t nb_sect;
  int16_t rho;
  int16_t angle;
  bool    sectors_free;
  int16_t alpha_first;
  int16_t alpha_last;
  int16_t theta_first;
  int16_t theta_last;

  // To handle gene overlaps, we remember where we have aldready drawn
  // something, with a precision of 1 degree. We handle up to 100 layers:
  //  - 50 layers inside the circle (lagging strand),
  //  - 50 layers outside the cricle (leading strand).
  // At first, only one layer is created, we create new layers when we
  // need them.
  bool* occupied_sectors[2][50];
  occupied_sectors[LEADING][0] = new bool[360];
  occupied_sectors[LAGGING][0] = new bool[360];
  for (int16_t angle = 0 ; angle < 360 ; angle++)
  {
    occupied_sectors[LEADING][0][angle] = false;
    occupied_sectors[LAGGING][0][angle] = false;
  }



  for (const auto& rna: gen_unit->rna_list()[LEADING]) {
    first = rna.first_transcribed_pos();
    last = rna.last_transcribed_pos();
    alpha_first   = (int16_t) round((double)(360 * first) / (double)gen_length);  //  == sect1 == alphaB
    alpha_last    = (int16_t) round((double)(360 * last)  / (double)gen_length);  //  == sect2 == alphaA
    theta_first   = Utils::mod(90 - alpha_first, 360);  //  == tetaB
    theta_last    = Utils::mod(90 - alpha_last, 360);  //   == tetaA
    if (theta_first == theta_last) theta_first = Utils::mod(theta_first + 1, 360);
    nb_sect = Utils::mod(theta_first - theta_last + 1, 360);
    // Outside the circle, look for the inmost layer that has all the sectors between
    // theta_first and theta_last free
    layer = 0;
    sectors_free = false;
    while (! sectors_free)
    {
      sectors_free = true;
      for (rho = 0 ; rho < nb_sect ; rho++)
        if (occupied_sectors[LEADING][layer][Utils::mod(theta_first - rho, 360)])
          sectors_free = false;
          break;
        }
      }
      {
        break; // All the needed sectors are free on the current layer
      }
      else
      {
        layer++;

        if ((layer >= outmost_layer) && (layer < 49))
        {
          // Add a new layer (actually, 2 layers, to maintain the symmetry)
          occupied_sectors[LEADING][outmost_layer] = new bool[360];
          occupied_sectors[LAGGING][outmost_layer] = new bool[360];
          for (int16_t angle = 0 ; angle < 360 ; angle++)
          {
            occupied_sectors[LEADING][outmost_layer][angle] = false;
            occupied_sectors[LAGGING][outmost_layer][angle] = false;
          }
          outmost_layer++;
          break; // A new layer is necessarily free, no need to look further
        }
        {
          // We shall not create a 51th layer, the CDS will be drawn on the
          // layer, probably over another CDS
          break;

    // Mark sectors to be drawn as occupied
    for (rho = 0 ; rho < nb_sect ; rho++)
      occupied_sectors[LEADING][layer][Utils::mod(theta_first - rho, 360)] = true;
    // Mark flanking sectors as occupied
    occupied_sectors[LEADING][layer][Utils::mod(theta_first + 1, 360)] = true;
    occupied_sectors[LEADING][layer][Utils::mod(theta_first - nb_sect, 360)] = true;