#!usr/local/bin/perl
use 5.010;
use strict;
use warnings;
use List::Util qw[min max];
use List::Util qw(sum);
use Data::Dumper qw (Dumper);
use Math::Random::Secure qw(irand);
use List::MoreUtils qw(uniq);
use List::MoreUtils qw(any);
use Statistics::Basic qw(:all);
use experimental qw(switch);
use POSIX;
use File::chdir;
use Tie::IxHash;
use Math::Complex;


# PARAMETRI ALGORITMO GENETICO
my $sizepop=50;
my $PercentualeElitimo = 0.1; # ELITISMO AL 10%
my $NumeroGenerazioni = 100;
my $pc = 0.7;           # probabilità di crossing over
my $pm = 0.01;         # probabilità di mutazione


#****************************************************************************************
#                                    FILE TO WRITE                                      #
#****************************************************************************************

#  ---------------------------------------------------->  1) SumUpFile.csv
# Creo un file in cui voglio tener traccia di:
#	> $e: Numero del cromosoma
#	> $G: Numero della generazione a cui il cromosoma appartiene
#	> la sua fitness
#	> la sua objective function
#	> fenotipo del cromosoma (covariate inserite nel modello -> cromosoma)
#	> termine di minimizzazione
my @argsNew = ("touch", "SumUpFile.csv");
system(@argsNew) == 0 or die "Error occurred: is impossible to create a new file for model informations\n";
open my $SUMUPFILE, '>', "SumUpFile.csv", or die $!;

#  ----------------------------------------------------> 2) WSumUpFile.csv
# Creo un file in cui voglio tener traccia dei soli cromosomi peggiori di ogni generazione che vengono sostituiti con elitismo
#	> $e: Numero del cromosoma
#	> $G: Numero della generazione a cui il cromosoma appartiene
#	> la sua objective function
#	> fenotipo del cromosoma (covariate inserite nel modello -> cromosoma)
#	> il run è andato a buon fine?
my @argsNews = ("touch", "Worst_SumUpFile.csv");
system(@argsNews) == 0 or die "Error occurred: is impossible to create a new file for model informations\n";
open my $WSUMUPFILE, '>', "Worst_SumUpFile.csv", or die $!;

#  ----------------------------------------------------> 3) InformazioniFitness.csv
# Creo un file in cui voglio tener traccia dei valori delle fitness:
#   > Fitness massima
#   > Fitness minima
#   > Fitness media
#   > Varianza dei valori della fitness e deviazione standard
my @args = ("touch", "InformazioniFitness.csv");
system(@args) == 0 or die "Error occurred: is impossible to create a new file for fitness information\n";
open my $FILEFITNESS ,'>' ,"InformazioniFitness.csv" , or die $!;
print $FILEFITNESS("Generazione, FitnessMassima,FitnessMimina,FitnessMedia,VarianzaFitness,DeviazioneStandard\n");

#****************************************************************************************
#                                  LETTURA FILE .MOD                                    #
#****************************************************************************************

#print("Please insert NONMEM file name you want to study (name.mod):\t");
#my $FileNONMEMName= <STDIN>;
#chomp $FileNONMEMName;

(my $lineRef, my $maxTheta, my $ParamRef, my $indexPK, my $indexOmega, my $indexPROBLEM, my $FileDataName) = readNONMEM('MODELLO.mod');
my @line = @$lineRef;
my @param = @$ParamRef;

my $MAXTheta = $maxTheta;


#****************************************************************************************
#                               LETTURA FILE CONFIGURAZIONE                             #
#****************************************************************************************

open my $FILE ,'<' ,'FileConfigurazione.scm' , or die 'Il file di configurazione non è stato trovato';
my @configFile = <$FILE>;
close FILE;

# Elimino i commenti dal file di configurazione
dropComment(\@configFile,'FileConfigurazione.scm');

(my $indVS, my $indTR, my $refContVar, my $refCategoricalVar, my $refV)=subConfigFile(\@configFile);

my @categoricalVariables = @{$refCategoricalVar};
my @continuousVariables = @{$refContVar};

my @Covariate = (@continuousVariables, @categoricalVariables);
my @v = @$refV;

#****************************************************************************************
#                               Hash covariata -> parametro                             #
#****************************************************************************************

my $refHash = hashCovariataParametro(\@continuousVariables, \@categoricalVariables, \@configFile, $indTR, $indVS, $refV);
my %hash = %$refHash;

my @covariate_to_print = sort { $hash{$a} <=> $hash{$b} } keys(%hash);

my @to_print=();
for(my $i = 0; $i<= $#covariate_to_print; $i++)
{
    my $param = $hash{$covariate_to_print[$i]};
    my @paramVect = @$param;
    
    for(my $j = 0; $j<=$#paramVect; $j++)
    {
        @to_print = (@to_print,$covariate_to_print[$i].$paramVect[$j]);
    }
}
# $toprint contiene le relazioni covaraiata-parametro che devo scrivere nel file SUMUPFILE
my $toprint = join(',', @to_print);
my $statusCromosom = 'created';


print $SUMUPFILE ("#Cromosoma, #Generazione, FitnessValue, ObjectiveFunction, TermOptimization, $toprint, StatusCromosom\n");
print $WSUMUPFILE ("#Cromosoma, #Generazione, FitnessValue, ObjectiveFunction, TermOptimization, $toprint\n");

#****************************************************************************************
#                                    Find valid states                                  #
#****************************************************************************************

my $refValidStates = ValidStates(\@configFile, $indVS, \@v);
my @ValidStatesContinuous = @$refValidStates;

#****************************************************************************************
#                                    LETTURA data.csv                                   #
#****************************************************************************************

# LETTURA DATA FILE
open FILE ,'<' ,$FileDataName , or die 'Il file.cvs non è stato trovato';
my @wholeCsv=<FILE>;
close FILE;

# Nomi delle colonne del file dei dati
my $names = $wholeCsv[0];
my @splitNames = split(/,|;|\s/, $names); # i dati possono essere separata da ; o , o uno carattere di tabulazione

# Valori di ogni colonna
my @valuess = @wholeCsv[1..$#wholeCsv];
my $i = 0;
foreach my $riga(@valuess)
{
    $riga =~ s/;|,|\s/ /g;
    $valuess[$i] = $riga;
    $i++;
}

# Estrapolo gli indici delle colonne delle variabili categoriche, continue e ID dalla matrice dei dati
(my $indexContinuousRef , my $indexCategoricalRef, my $indexID) = subInfoDataCsv(\@splitNames,\@continuousVariables,\@categoricalVariables);
my %indexContinuous = %{$indexContinuousRef};
my %indexCategorical = %{$indexCategoricalRef};

#  ------------------------------- MATRICE DEI DATI -------------------------------------
my @splitValues;
my @ID;
my %matrix;

# Riempio la matrice di dati - righe i - colonne j
for(my $i=1; $i<=$#valuess+1; $i++)
{
    @splitValues = split(/;/, $valuess[$i-1]);@splitValues = split(/ /, $valuess[$i-1]);
    for(my $j=0; $j<=$#splitNames; $j++)
    {
        $matrix{$i}{$j}= $splitValues[$j];
        $matrix{0}{$j} = $splitNames[$j];   # prima riga, tutte le colonne : nomi delle variabili
    }
    
    $ID[$i-1] =$matrix{$i}{$indexID};       # colonna delle ID
}

# I dati per singolo paziente sono ripetuti: non posso prendere le età su tutte le righe, ma l'età per paziente


#  ---------------------------- MEDIANE DEI DATI CONTINUI  -------------------------------

(my $medianCVariablesRef, my $minRef, my $maxRef) = medianaDatiContinui(\@continuousVariables, \@ID,\%matrix, \%indexContinuous);
my %medianCVariables = %$medianCVariablesRef;   # ---> hash: covariatacontinua -> valore mediano
my %minimoCVariables = %$minRef;                # ---> hash: covariatacontinua -> valore minimo
my %massimoCVariables = %$maxRef;               # ---> hash: covariatacontinua -> valore massimo


#  --------------------------- FREQUENZE DEI DATI CATEGORICI -----------------------------

# Devo trovare valore più frequente nelle variabili categoriche
# 1) trovare unique dei valori
# 2) vedere quali tra questi è il più frequente

(my $FreqTotRef, my $FreqCovRef, my $refUnique)=FrequenzaDatiCategorici(\@categoricalVariables, \@ID, \%matrix, \%indexCategorical);
my %FreqTot = %{$FreqTotRef}; # CovariataCategorica -> Tutte le categorie
my %FreqCov = %{$FreqCovRef}; # CovariataCategorica -> Categoria più frequente
my %VarUnique = %{$refUnique}; # CovariataCategorica -> Numero di livelli -1


#****************************************************************************************
#                                DEFINIZIONE DEL GENOMA                                 #
#****************************************************************************************

# Il genoma dipende da quali sono le relazioni che l'utente vuole studiare
my $numberofStates = $#ValidStatesContinuous+1;
my $numberofBits = CEIL(log2($numberofStates));
my @AmGeni = genome($numberofStates);
my @AmGeniCat = ('0', '1');
my @WholeGenes = (@AmGeni, @AmGeniCat);
my $nAmGeni = $#AmGeni+1;
#print("Geni ammessi per la codifica: @AmGeni\n\n");


#****************************************************************************************
#                          POPOLAZIONE INIZIALE - GENERAZIONE 0                         #
#****************************************************************************************

print("\n\n*********** [GENERAZIONE 0] *********** \n");

my %FitnessHash;
my %popolazione;

my $t = tie %FitnessHash, 'Tie::IxHash';
my $tt = tie %popolazione, 'Tie::IxHash';

my $e; # Conta gli indivdui
my $nrelazioniContinue;
my $NumeroCromosomiElitismo = int($PercentualeElitimo*$sizepop);
my $TermMinimization;
my $infoElitismo;
my $FitnessValue;

my @cromosoma;
my @vecttoprintComplete;
my @ObjFunVect = ();
my @FitnessVect = ();
my @TermMinimizationVect =();
my $toprintVect;

my $ObjFun;

my @arg1 = ("mkdir", "Generazione0");system(@arg1); # crea cartella Generazione 0
my @arg2 = ("cp", "MODELLO.mod", "Generazione0");  system(@arg2); # copia il file del modello
my @arg3 = ("cp", "DataFile.csv", "Generazione0"); system(@arg3); # copia il file dei dati
$CWD ="Generazione0";

my $G=0;

for($e=0; $e<$sizepop; $e++)
{
    my %relationstart;
    
	my @gene = (); 				
    my @Vect_to_print = ();
    @cromosoma = ();
    
	$maxTheta = $MAXTheta;
    
    # Creo un file NONMEM per ogni cromosoma
    my @args = ("touch", "Cromosoma$e.mod");
    system(@args) == 0 or die "system @args failed: $?";
    
    # Apro il file in scrittura '>'
    open my $FILE ,'>' ,"Cromosoma$e.mod" , or die $!;
    
    # Ricopio nel nuovo file le prime righe del modello
    $line[$indexPROBLEM] = substr($line[$indexPROBLEM], index($line[$indexPROBLEM], '$PROBLEM'), -1);
    print $FILE("@line[$indexPROBLEM..$indexPK]\n");
    
    #****************************************************************************************
    #                             1. WRITE NEW NONMEM FILE                                  #
    #****************************************************************************************
    #print("\nI'm writing NONMEM file for the new chromosome ... \n\n");
    (my $RelationStartRef, my $CromosomaRef, $maxTheta, $nrelazioniContinue, my $refhash) =  write_NONMEMfile(\@param, \@AmGeni, \@AmGeniCat, \%relationstart, $maxTheta, \%medianCVariables, $FILE, \@ValidStatesContinuous, \@cromosoma, \%FreqCov, \%FreqTot, \@gene, $nAmGeni, $indexPK, $indexOmega, \@line, $G, \@Covariate,\%hash, \@continuousVariables, \@categoricalVariables, \%VarUnique, \%minimoCVariables, \%massimoCVariables);
    %relationstart = %$RelationStartRef;
    @cromosoma = @$CromosomaRef;
	
    my %FenotipoHash = %{$refhash};
    
    for (my $b=0; $b<= $#to_print; $b++)
    {
        @Vect_to_print = (@Vect_to_print, $FenotipoHash{$to_print[$b]});   # vettore fenotipo scelto per un cromosoma -> es 5 0 1 1 0
    }
    my $vecttoprint = join(',', @Vect_to_print);						   # 5, 0, 1, 1, 0
    $vecttoprintComplete[$e] = join(',', @Vect_to_print); 				   # Vettore fenotipo di tutti i cromosomi della generazione 0
    
    close $FILE;
    
    #****************************************************************************************
    #                                  2. RUN NONMEM                                        #
    #****************************************************************************************
    #print("\nI'm running NONMEM for the estimation...\n\n");
    runNONMEM($e);
    
    #****************************************************************************************
    #                     3.  FIND FITNESS in ResultDirectory                               #
    #****************************************************************************************
    #print("\nI'm finding the fitness value for the $e chomosome...\n\n");
    (my $result, my $ObjFunction) = findFitness($e, $maxTheta, $MAXTheta, \%FitnessHash);
    $FitnessHash{$e} = $result;
    $FitnessValue = $result;
   
    # Find if Hessian
   ($ObjFunction, $FitnessValue, $TermMinimization) = FindHessian($e, $ObjFunction, $FitnessValue);

    print $SUMUPFILE ("$e, $G, $FitnessValue, $ObjFunction, $TermMinimization, $vecttoprint, $statusCromosom \n");
    
    @TermMinimizationVect = (@TermMinimizationVect, $TermMinimization);
    @FitnessVect = (@FitnessVect, $FitnessValue);	# Vettore che contiene tutte le fitness della generazione 0
    @ObjFunVect = (@ObjFunVect, $ObjFunction); 		# Vettore che contiene tutte funzioni obj della generazione 0
    
    #****************************************************************************************
    #                             4.  SAVE POPoLATION                                       #
    #****************************************************************************************
    for(my $K=0; $K<=$#cromosoma;$K++)
    {
        $popolazione{$e}[$K] = $cromosoma[$K];
    }
    
} #close for sizePop
$CWD ="..";

my @nGenesVect= @{$popolazione{0}};
my $nGene = $#nGenesVect;                       # Numero di geni nel cromosoma

#****************************************************************************************
#               Selezione migliori cromosomi (con fitness maggiore)                       #
#****************************************************************************************

# I cromosomi vengono ordinati in modo decrescente secondo il loro valore di fitness
# I cromosomi così ordinati sono immagazinati nell'Hash $FitnessOrder ($indiceCromosoma -> Valore della fitness)
# Seleziono allora gli indici dei primi x Cromosomi e li salvo nella lista @BestIndex

my %FitnessOrder;

# Sort in acending order
foreach my $cromosom ( sort {$FitnessHash{$a} <=> $FitnessHash{$b}} keys %FitnessHash)
{
    $FitnessOrder{"$cromosom"} = $FitnessHash{$cromosom};
};

# Prendo gli indici dei primo 10 cromosomi con fitness mmaggiore
my @Index = keys %FitnessOrder;
my @BestIndex = @Index[0..2]; 

my @FitnessBestIndex = values %FitnessOrder;
my @FitnessBest = @FitnessBestIndex[0..2];

# Informazioni relative alla fitness della prima generazioni
my @Fitnesssss = values %FitnessHash;
my @FitnessVera = ();
for (my $i=0; $i<=$#FitnessVect; $i++)
{
	
	if($FitnessVect[$i] != 0)
	{
		@FitnessVera = (@FitnessVera, $FitnessVect[$i]);
	}
}
my $variance1 = varianza(@FitnessVera);
my $minFit = min(@FitnessVera);
my $maxFit = max(@FitnessVera);
my $meanFit = media(@FitnessVera);
my $devFit = sqrt($variance1);

print $FILEFITNESS("$G, $maxFit,$minFit,$meanFit,$variance1, $devFit\n");

#****************************************************************************************
#                              COUNT GENES X POSIZIONE                                  #
#****************************************************************************************

# Si vuole calcolare la frequenza dei geni per posizione
# Es. Qual è la frequenza del gene 000 in posizione 1? e in posizione 2? ...
my @crom;
my %countResults;
my %count;

for(my $i=0; $i<=$#BestIndex; $i++)
{
    @crom = ();
    
    # Estraggo gli individui migliori dalla popolazione
    for(my $j=0; $j<=$nGene; $j++)
    {
        my @cromx = $popolazione{$BestIndex[$i]}[$j];
        @crom = (@crom, @cromx);
    }
    
    my $countResultsRef= count_x_Position($i, $numberofStates, \@crom, \@WholeGenes, \%count);
    %countResults = %{$countResultsRef};
    %count = %countResults;
}

#****************************************************************************************
#                        GENI PIÙ FREQUENTI IN OGNI POSIZIONE                           #
#****************************************************************************************
# Si vuole conoscere qual è il gene più frequente per ogni posizione nel cromosoma
my @geneMax = GenesMostFrequent_x_Position(\@WholeGenes, \%count);
#print("\nGeni più frequenti per posizione @geneMax\n\n");

#****************************************************************************************
#                               Selezione + Crossing Over                               #
#****************************************************************************************

my @offspring1;
my @offspring2;
my @ObjectiveBest;
my @FitBest;
my @TermBest;
my @Vect_to_printBest;

# Voglio creare vettori che contengano le informazioni contenute all'interno del file SumUpFile.csv
# In particolare devo creare:
#	> @Tabella: contiene tutti i fenotipi all'interno del file
#	> @ObjTabella: contiene tutti i valori della funzione oggettivo del file
# 	> @FitnessTab: contiene tutti i valori della fitness del file
#	> @TermMinimizationTab: contiene tutti i termini di minimizzazione del file
# Quindi all'inizio devo contenere i valori della generazione 0. Successivamente accodo i valori di tutte le altre generazioni
my @Tabella= @vecttoprintComplete; 
my @ObjTabella = @ObjFunVect;
my @FitnessTab = @FitnessVect;
my @TermMinimizationTab = @TermMinimizationVect;

my %NewPopolazione;
my %NewPop;
my $np = tie %NewPop, 'Tie::IxHash';

my $countsOffspring;

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# 					              		GENERAZIONI SUCCESSIVE											 %
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# Segue la formazioni delle generazioni successive. I passaggi da fare sono i seguenti:
#	1) Elitismo1: scelta dei migliori cromosomi della generazione precedente
# 	2) Parent Selection
# 	3) Crossing over
#	4) Mutazione
#	5) Completati i cromosomi, posso scrivere un file NONMEM per ogni cromosoma e farne il run se necessario
# 	6) Elitismo2: scelta dei peggiorni cromosomi della popolazione corrente
#	7) Elitismo3: sostituzione dei peggiorni cromosomi della popolazione corrente con i migliori della popolazione precedente


for (my $G=1; $G<=$NumeroGenerazioni; $G++)
{

    #****************************************************************************************
    #                                 CREO - COPIO FILE                                     #
    #****************************************************************************************
    # Creo una cartella per ogni generazione all'interno della quale copio il file del modello
    # (.mod) e il file dati
    print("\n\n*********** [GENERAZIONE $G] *********** \n");
    
    my @arg1 = ("mkdir", "Generazione$G");              system(@arg1);
    my @arg2 = ("cp", "MODELLO.mod", "Generazione$G");  system(@arg2);
    my @arg3 = ("cp", "DataFile.csv", "Generazione$G"); system(@arg3);
    $CWD ="Generazione$G";
    
    $countsOffspring = 0;
    @Vect_to_printBest = ();
    
    #****************************************************************************************
    #                                 ELITISMO - STEP 1                                     #
    #****************************************************************************************
    # (1) Ordino in modo decrescente le fitness per selezionare cromosomi con fitness migliore
    #****************************************************************************************
    (my $IndiciOrdinatiBestRef,my $FitnessOrdinatiBestRef) = descendHash(\%FitnessHash);
    my @IndiciOrdinatiBest = @$IndiciOrdinatiBestRef;
    my @FitnessOrdiantiBest = @$FitnessOrdinatiBestRef;
    
    # Migliori 10% degli indici
    my @BestIndiciElitismo = @IndiciOrdinatiBest[0..$NumeroCromosomiElitismo];
    
    # Miglior 10% delle fitness
    my @FitnessBestElitismo = @FitnessOrdiantiBest[0..$NumeroCromosomiElitismo];
    
    for (my $obix =0; $obix<=$#BestIndiciElitismo; $obix++)
    {
        $ObjectiveBest[$obix] = $ObjFunVect[$BestIndiciElitismo[$obix]];
        $FitBest[$obix] = $FitnessVect[$BestIndiciElitismo[$obix]];
		$TermBest[$obix] = $TermMinimizationVect[$BestIndiciElitismo[$obix]];
        @Vect_to_printBest = (@Vect_to_printBest, $vecttoprintComplete[$BestIndiciElitismo[$obix]]);
    }
    
    
    for(my $i=0; $i< $sizepop/2; $i++)
    {
        #******************************************************************
        #          PARENT SELECTION - Tournament selection                #
        #******************************************************************
        # Seleziono due genitori per la riproduzione
        #print("[PARENTS]\n");
        #my $ParentCromosomRef1 = RouletteWheel(\%FitnessHash, \%popolazione);
        my $ParentCromosomRef1 = TournamentSelection($sizepop, \%FitnessHash, \%popolazione);
        my @ParentCromosom1 = @{$ParentCromosomRef1};
        
        #my $ParentCromosomRef2 = RouletteWheel(\%FitnessHash, \%popolazione);
        my $ParentCromosomRef2 = TournamentSelection($sizepop, \%FitnessHash, \%popolazione);
        my @ParentCromosom2 = @{$ParentCromosomRef2};
        
        #******************************************************************
        #         CROSSING OVER - OnePoint Crossover Recombination        #
        #******************************************************************
        if(rand(1)<= $pc)
        {
            (my $offspring1Ref, my $offspring2Ref) = onePoinCrossOver(\@ParentCromosom1, \@ParentCromosom2);
            $NewPopolazione{"$countsOffspring"} = join('',@$offspring1Ref);
            $countsOffspring++;
            @offspring2 = @$offspring2Ref;
            $NewPopolazione{"$countsOffspring"} = join('', @$offspring2Ref);
            $countsOffspring++;
        }
        else
        {
            $NewPopolazione{"$countsOffspring"} = join('', @ParentCromosom1);
            $countsOffspring++;
            $NewPopolazione{"$countsOffspring"} = join('', @ParentCromosom2);
            $countsOffspring++;
        }
        
    }
    my @individuoMutato;
    #******************************************************************
    #                            MUTAZIONE                            #
    #******************************************************************
    for(my $i=0; $i<$sizepop; $i++)
    {
        my @individuo = split('', $NewPopolazione{"$i"});
        
        for(my $j=0; $j<=$#individuo; $j++)
        {
            my $randomicNumber = rand(0.1);
            
            if($randomicNumber <= $pm)
            {
                $individuoMutato[$j] = NegazioneBit($individuo[$j]);
            }
            else
            {
                $individuoMutato[$j] = $individuo[$j];
            }
        }
        
        $NewPopolazione{"$i"} = join('', @individuoMutato);
    }
    
    # Splittto il cromosoma e sostituisco con il gene più frequente
    my $NewPopRef = splitNewpop($nrelazioniContinue,$countsOffspring,$nGene,\%NewPopolazione, \@geneMax, \@AmGeni, $numberofBits);
    %NewPop = %{$NewPopRef};
    
	# I vettori che contiengono le infotmazioni sulla fitness, obj fun, ..., devono cambiare ad ogni generazione, quindi li pulisco
    @ObjFunVect = ();
    @FitnessVect = ();
	@TermMinimizationVect = ();
	@vecttoprintComplete = ();
	
    my @statusCromosoms = ();
	
    for(my $e=0; $e<$sizepop; $e++)
    {
        
        my @gene=();
        @cromosoma = ();
        my @Vect_to_print= ();
        my $nAmGeni = $#AmGeni+1;
        my %relationstart;
        my $b;
        my @ret;
        
        # Creo il file e provo a fare questi print nel file
        my @args = ("touch", "Cromosoma$e.mod");
        system(@args) == 0 or die "system @args failed: $?";
        
        # Apro il file in scrittura '>'
        open my $FILE ,'>' ,"Cromosoma$e.mod" , or die $!;
        print $FILE("@line[$indexPROBLEM..$indexPK]\n");
        
        
        #****************************************************************************************
        #                             1. WRITE NEW NONMEM FILE                                  #
        #****************************************************************************************
        # questo non va bene
        my @cromosomaIN = @{$NewPop{"$e"}};
        
        # Scrivo il file di NONMEM relativo ad al cromosoma @cromosomaIN
       (my $RelationStartRef, my $CromosomaRef, $maxTheta, my $nnn, my $refhash ) =  write_NONMEMfile(\@param, \@AmGeni, \@AmGeniCat, \%relationstart, $MAXTheta, \%medianCVariables, $FILE, \@ValidStatesContinuous, \@cromosomaIN, \%FreqCov, \%FreqTot, \@gene, $nAmGeni, $indexPK, $indexOmega, \@line, $G,\@Covariate,\%hash, \@continuousVariables, \@categoricalVariables, \%VarUnique, \%minimoCVariables, \%massimoCVariables);
        %relationstart = %$RelationStartRef;
        
        my %FenotipoHash = %{$refhash};
 
        for (my $b=0; $b<= $#to_print; $b++)
        {
            @Vect_to_print = (@Vect_to_print, $FenotipoHash{$to_print[$b]});
        }
        $vecttoprintComplete[$e] = join(',', @Vect_to_print); #generazioni successive
        @cromosoma = @$CromosomaRef;
        
        close $FILE;
        
        #****************************************************************************************
        #                                  2. RUN NONMEM                                        #
        #****************************************************************************************
        if( any {$_ ~~ $vecttoprintComplete[$e]}@Tabella)
        {
			$statusCromosoms[$e] = 'duplicated';
            
			my $IndiceTrue;
        
            for(my $indiceTabella=0; $indiceTabella<=$#Tabella; $indiceTabella++)
            {
                if($vecttoprintComplete[$e] eq $Tabella[$indiceTabella])
                {
                    $IndiceTrue = $indiceTabella;
                    last;
                }
            }
            
            $FitnessValue = $FitnessTab[$IndiceTrue];
            $ObjFun = $ObjTabella[$IndiceTrue];
            $vecttoprintComplete[$e] = $Tabella[$IndiceTrue];
			#$toprintVect = $Tabella[$IndiceTrue]; # fenotipo
			$TermMinimization = $TermMinimizationTab[$IndiceTrue];
			
			$FitnessHash{$e} = $FitnessTab[$IndiceTrue];
			
        }
        else
        {
			$statusCromosoms[$e] = 'created';
            runNONMEM($e);
			
			#****************************************************************************************
            #                     3.  FIND FITNESS in ResultDirectory                               #
            #****************************************************************************************
            (my $result, $ObjFun) = findFitness($e, $maxTheta, $MAXTheta, \%FitnessHash);
            $FitnessHash{$e} = $result;
            $FitnessValue = $result;
			#$toprintVect = $vecttoprintComplete[$e]; #fenotipo
			($ObjFun, $FitnessValue, $TermMinimization) = FindHessian($e, $ObjFun, $FitnessValue);
		}
        @ObjFunVect = (@ObjFunVect, $ObjFun);
        @FitnessVect = (@FitnessVect, $FitnessValue);
		#@vecttoprintComplete = (@vecttoprintComplete, $toprintVect);
        #@Tabella = (@Tabella,  $vecttoprintComplete[$e]);
		@TermMinimizationVect = (@TermMinimizationVect, $TermMinimization);
    } #close for sizePop (over $e)
    
    #****************************************************************************************
    #                                  ELITISMO - STEP 2                                    #
    #****************************************************************************************
    # (2) Ordino in modo crescente le fitness per selezionare cromosomi con fitness peggiore
    #****************************************************************************************
    # Sort in acending order for ELITISMO
    (my $refIndex, my $refValues)  = ascendHash(\%FitnessHash);
    my @IndexOrdinatiWorst = @$refIndex;
    my @FitnessOrdinatiWorst = @$refValues;
    
    # Prendo gli indici del primo 10% dei cromosomi con fitness maggiore
    
    # Migliori 10% degli indici
    my @WorstIndiciElitismo = @IndexOrdinatiWorst[0..$NumeroCromosomiElitismo];
    
    # Miglior 10% delle fitness
    my @FitnessWorstElitismo = @FitnessOrdinatiWorst[0..$NumeroCromosomiElitismo];
    
    #****************************************************************************************
    #                                  ELITISMO - STEP 3                                    #
    #****************************************************************************************
    # (3) Sostituisco i peggiori cromosomi nella popolazione corrente con i migliori della
    #     popolazione precedente
    #****************************************************************************************
    
    for(my $l=0; $l <= $#BestIndiciElitismo; $l++)
    {
		
		$statusCromosoms[$WorstIndiciElitismo[$l]] = 'elited';
		
        # Nella popolazione corrente (%NewPop) devo sostituire i peggiori cromosomi
        # ovvero quelli di indice (@WorstIndiciElitismo).
        # A questi devo sostituire i migliori cromosomi della popolazione precedente
        # (%popolazione) ovvero quei cromosomi di indice (@BestIndiciElitismo);
		print $WSUMUPFILE("$WorstIndiciElitismo[$l], $G, $FitnessVect[$WorstIndiciElitismo[$l]],$ObjFunVect[$WorstIndiciElitismo[$l]],$TermMinimizationVect[$WorstIndiciElitismo[$l]],$vecttoprintComplete[$WorstIndiciElitismo[$l]]\n");
        $NewPop{$WorstIndiciElitismo[$l]} = $popolazione{$BestIndiciElitismo[$l]};
        $FitnessHash{$WorstIndiciElitismo[$l]} = $FitnessBestElitismo[$l];
        
        $vecttoprintComplete[$WorstIndiciElitismo[$l]] = $Vect_to_printBest[$l];
		$ObjFunVect[$WorstIndiciElitismo[$l]] = $ObjectiveBest[$l]; 
        $FitnessVect[$WorstIndiciElitismo[$l]] = $FitBest[$l];
		$TermMinimizationVect[$WorstIndiciElitismo[$l]] = $TermBest[$l];
    }
    
    @ObjTabella = (@ObjTabella, @ObjFunVect);
    @FitnessTab = (@FitnessTab, @FitnessVect);
	@TermMinimizationTab = (@TermMinimizationTab, @TermMinimizationVect);
	@Tabella = (@Tabella, @vecttoprintComplete);
	
    # POPOLAZIONE finale
    %popolazione = %NewPop;
    
    $CWD ="..";
    
    # Informazioni relative alla fitness delle generazioni successive
    my @Fitness = values(%FitnessHash);
	my @FitnessVera = ();
	for (my $iq=0; $iq<=$#FitnessVect; $iq++)
	{
		if($FitnessVect[$iq] != 0)
		{
			@FitnessVera = (@FitnessVera, $FitnessVect[$iq]);
		}
	}	
    my $MaxValueFit = max(@FitnessVera);
    my $MinValueFit = min(@FitnessVera);
    my $MeanValueFit = media(@FitnessVera);
    my $VarianceValueFit = varianza(@FitnessVera);
	my $DeviazioneFit = sqrt($VarianceValueFit);
 
    print $FILEFITNESS("$G, $MaxValueFit,$MinValueFit,$MeanValueFit,$VarianceValueFit, $DeviazioneFit\n");
    
	for(my $individuoE=0; $individuoE<$sizepop; $individuoE++)
    {
        print $SUMUPFILE ("$individuoE, $G, $FitnessVect[$individuoE], $ObjFunVect[$individuoE] ,$TermMinimizationVect[$individuoE] ,$vecttoprintComplete[$individuoE], $statusCromosoms[$individuoE] \n");
    }

} #close for G
close $FILEFITNESS;
close $SUMUPFILE;
close $WSUMUPFILE;


#**************************************************************************************
#**************************************************************************************
#**************************************************************************************
#                                                                                     *
#                                                                                     *
#                                   SUBRUTINES                                        *
#                                                                                     *
#                                                                                     *
#**************************************************************************************
#**************************************************************************************
#**************************************************************************************
#**************************************************************************************


# **********************************  readNONMEM  ***********************************
sub readNONMEM
{
    my $NomeFile = $_[0];
    
    my @param=();
    my $indexPK = -1;
    my $indexOmega = -1;
    my $indexPROBLEM;
    my $indexDATA;
    my $indexTheta;
    
    my $i=0;
    my $k=0;
    my $j = 0;
    
    my @line;
    my @findteta;
    my $flag = 0;
    my @tetas;
    
    open FILE ,'<' , $NomeFile , or die 'Il file .mod non è stato trovato';
    my @FileNonmem = <FILE>;
    close FILE;
    
    foreach my $riga(@FileNonmem)
    {
        $line[$k] = $riga;
        $line[$k] =~ s/\\//;
        
        # Indice della riga $PROBLEM
        if($riga  =~ m/\$PROBLEM\s*/) {
            $indexPROBLEM = $k;
        }
        
        # Indice della riga $DATA
        if($riga  =~ m/\$DATA\s*/) {
            $indexDATA = $k;
        }
        
        # Indice della riga $PK
        if($riga  =~ m/\$PK\s*/) {
            $indexPK = $k;
        }
        
        # Indice della riga $OMEGA
        if($riga  =~ m/\$OMEGA\s*/)
        {
            $indexOmega = $k; # in questo modo mi prende l'ultimo $THETA
        }
        
        # Indice della riga $THETA
        if($riga  =~ m/\$THETA\s*/)
        {
            $indexTheta = $k;
        }
        
        $k++;
        
        # I parametri su cui effettuare la stima delle covariate vengono riportate nel file come
        # TVparam
        if($riga =~ m/TV.+THETA/)
        {
            $flag = 1;
            
            # Parametri a cui applicare la ricerca delle covariate
            $param[$i] = trim(substr $riga, 2, (index($riga, '=')-2));
            chomp($param[$i]);
            $param[$i]=~s/TV//;
            
            $i++;
        }
        
        my $position = 0;
        my $info;
        
        # Trova i teta
        if($riga =~ m/.+THETA/)
        {
            for (my $indice = 0; $indice< length($riga); $indice++)
            {
                $info = index($riga, 'THETA(', $position+1);
                $position = $info;
                if ($position != -1)
                {
                    $tetas[$j] = substr $riga, $position+6, 1;
                    $j++;
                }
            }
        }
    }
    
    # Teta massimo
    my $maxTheta = max(@tetas);
    
    # Nome file dati
    my @NonNullElement;
    my @lineData = split / /, $line[$indexDATA];
    my $qx=0;
    
    for(my $l = 0; $l<=$#lineData; $l++)
    {
        if (!$lineData[$l] eq '')
        {
            $NonNullElement[$qx] = $lineData[$l];
            $qx++;
        }
    }
    my $FileDataName = $NonNullElement[1];
    
    # Check error
    if($indexPK == -1){die "Error occurred: NO \$PK block was found in file .mod \n";}
    if($indexOmega == -1){die "Error occurred: NO \$OMEGA block was found in file .mod \n";}
    if($flag == 0){die "Error occurred: NO TV - parameters were found in the file.mod\n"}
    
    return(\@line, $maxTheta, \@param, $indexPK, $indexOmega, $indexPROBLEM, $FileDataName);
}

# **********************************  dropComment  ***********************************
# Elimina i commenti dai file ovvero tutte le stringhe che sono precedute dal carattere ";"
sub dropComment
{
    my $refFile = $_[0]; my @FileVect = @$refFile;
    my $nomeFile = $_[1];
    my $q=0;
    
    open my $FILE ,'>' ,$nomeFile , or die 'Error occured in dropComment: Il file di configurazione non è stato trovato';
    
    foreach my $riga(@FileVect)
    {
        my $toCopy = substr($riga, 0, index($riga, ';'));
        my $toDrop = substr($riga, index($riga, ';'), -1);
        $toDrop = " \n";
        $FileVect[$q] = $toCopy.$toDrop;
        print $FILE $FileVect[$q];
        $q++;
    }
    
    close $FILE;
}


# ********************************   subConfigFile ********************************
# Lettura del file di configurazione
sub subConfigFile
{
    my $refFile = $_[0]; my @configFile = @{$refFile};
    
    my $i=0;
    my %ix;
    my $ct; # for categorical covariates
    my $cc; # for continuous covariates
    my $flagContinuous  = -1;
    my $flagCategorical = -1;
    
    my @indiciSezioniVect;
    
    my @search = ('test_relations', 'valid_states');
    
    for(my $j=0; $j<=$#search; $j++)
    {
        $i=0;
        foreach my $riga(@configFile)
        {
            # indici delle sezioni [test relation e valid states]
            if($riga =~ m/\[$search[$j]\]/)
            {
                $ix{$search[$j]} = $i;
            }
            
            # -------------> Covariate categoriche
            if($riga=~ m/categorical_covariates/)
            {
                $ct = substr $riga, (index($riga, '=')+1), -1;
                $ct =~ s/\\//;
                $flagCategorical = 0;
                
            }
            
            # -------------> Covariate continue
            if($riga=~ m/continuous_covariates/)
            {
                $cc = substr $riga, (index($riga, '=')+1), -1;
                $cc =~ s/\\//;
                $flagContinuous = 0;
            }
            $i++;
        }
    }
    
    my @continuousVariables =();
    my @categoricalVariables =();
    if($flagContinuous == 0){  @continuousVariables = split(/,/, $cc);  };
    
    # Trim continuous covariates
    for(my $w= 0 ; $w<= $#continuousVariables; $w++) {   $continuousVariables[$w] = trim($continuousVariables[$w]);  }

    if($flagCategorical == 0)  {@categoricalVariables = split(/,/, $ct);  };
    
    # Trim categorical covariates
    for(my $w= 0 ; $w<= $#categoricalVariables; $w++) { $categoricalVariables[$w] = trim($categoricalVariables[$w]); }
    
    # Check if there are no covariates
    if($flagCategorical == -1 && $flagContinuous == -1){ die "Error Occurred: no covariates were found in configuration file\n" };
    
    @indiciSezioniVect = values(%ix);
    @indiciSezioniVect = sort {$a <=> $b} @indiciSezioniVect; #ascending order
    my $vi = join(' ', @indiciSezioniVect);
    my $indTR;
    my $indVS;
    
    for($i=0; $i<=$#indiciSezioniVect; $i++)
    {
        if(exists $ix{valid_states})
        {
            if($indiciSezioniVect[$i] =~ $ix{valid_states})
            {
                $indVS=$i;
                
            }
        }else{die "Error Occurred: no [valid_states] were found in configuration file!\n"}
        
        if(exists $ix{test_relations})
        {
            if($indiciSezioniVect[$i] =~ $ix{test_relations})
            {
                $indTR=$i;
            }
        }else{die "Error Occurred: no [test_relations] were found in configuration file!\n"}
    }
    
    return($indVS, $indTR, \@continuousVariables, \@categoricalVariables, \@indiciSezioniVect);
}



# ********************************* hashCovariataParametro  ***********************************
sub hashCovariataParametro
{
    # Input
    my $ref1 = $_[0]; my @continuousVariables = @{$ref1};
    my $ref2 = $_[1]; my @categoricalVariables = @{$ref2};
    my $ref3 = $_[2]; my @configFile = @$ref3;
    my $indTR = $_[3];
    my $indVS = $_[4];
    my $ref4 = $_[5]; my @v= @$ref4;
    
    my @Covariate = (@continuousVariables, @categoricalVariables);
    my %hash;
    my @infoIndex;
    my $l = 0;
    
    # Controllo sulle covariate continue
    for(my $i=0; $i<= $#continuousVariables; $i++)
    {
        my $k = 0;
        my $j = 0;
        
        my $covariata = trim($continuousVariables[$i]);
        
        foreach my $riga(@configFile)
        {
            if($j> $v[$indTR] && $j < $v[$indVS] && $riga =~ m/\w\s*\=\w\W*/)
            {
                
                my $parametro = substr($riga, 0, index($riga, '='));
                my $covariate = trim(substr($riga, index($riga, '=')+1, -1));
                $covariate =~ s/\\//;
                
                my @covariateVect = split(/,/ , $covariate);
                #@covariateVect = @{trim(@covariateVect)};
                
                if(any {$_ eq $covariata} @covariateVect)
                {
                    $hash{"$covariata"}[$k] = $parametro;
                    $k++;
                }
            }
            $j++;
        }
    }
    
    # Controllo sulle covariate categoriche
    for(my $i=0; $i<= $#categoricalVariables; $i++)
    {
        my $j = 0;
        my $k = 0;
        my $covariata = $categoricalVariables[$i];
        
        foreach my $riga(@configFile)
        {
            if($j> $v[$indTR] && $j < $v[$indVS] && $riga =~ m/\w\s*\=\w\W*/)
            {
                my $parametro = substr($riga, 0, index($riga, '='));
                my $covariate = trim(substr($riga, index($riga, '=')+1, -1));
                $covariate =~ s/\\//;
                
                my @covariateVect = split(/,/ , $covariate);
                #@covariateVect = @{trim(@covariateVect)};
                if(any {$_ eq $covariata} @covariateVect)
                {
                    $hash{"$covariata"}[$k] = $parametro;
                    $k++;
                }
            }
            $j++;
        }
    }
    return \%hash;
}

# ********************************* valid states ***********************************
sub ValidStates
{
    my $ref1 = $_[0]; my @configFile = @$ref1;
    my $indVS = $_[1];
    my $ref2 = $_[2]; my @v = @$ref2;
    
    my $j = 0;
    my @ValidStatesContinuous;
    
    foreach my $riga(@configFile)
    {
        if($j > $v[$indVS] && $riga =~ m/continuous\=(\d+\W*)/)
        {
            my $relazioni = substr($riga, index($riga, '=')+1, -1);
            $relazioni =~ s/\\//;
            @ValidStatesContinuous = split(/,/,$relazioni);
        }
        $j++;
    }
    for (my $element=0; $element<=$#ValidStatesContinuous; $element++)
    {
        $ValidStatesContinuous[$element]= trim($ValidStatesContinuous[$element]);
    }
    return \@ValidStatesContinuous;
}


# **********************************  subIndoDataCsv **********************************
sub subInfoDataCsv
{
    my $splitNamesRef = $_[0]; my @splitNames = @{$splitNamesRef};
    my $continuousVariablesRef = $_[1]; my @continuousVariables = @{$continuousVariablesRef};
    my $categoricalVariablesRef = $_[2]; my @categoricalVariables = @{$categoricalVariablesRef};
    
    my %indexContinuous;
    my %indexCategorical;
    my $indexID;
    
    for(my $i=0; $i<= $#splitNames; $i++) # sulle colonne
    {
        $splitNames[$i]=~ s/\W//;
        $splitNames[$i]=trim($splitNames[$i]);
        
        # Indici delle colonne che contengono le variabili continue
        for(my $j=0; $j<= $#continuousVariables; $j++)
        {
            my $a = ($continuousVariables[$j]);
            $a=trim($a);
            
            if($splitNames[$i] eq $a)
            {
                $indexContinuous{$a}= $i;
            }
        }
        
        # Indici delle colonne che contengono le variabili categoriche
        for(my $j=0; $j<= $#categoricalVariables; $j++)
        {
            my $b = ($categoricalVariables[$j]);$b = trim($b);
            
            if($b eq $splitNames[$i])
            {
                $indexCategorical{$b}= $i;
            }
        }
        
        # Indice della colonna ID
        if($splitNames[$i] eq 'ID')
        {
            $indexID= $i;
        }
    }
    
    return(\%indexContinuous , \%indexCategorical, $indexID);
}


# ******************************** medianaDatiContinui ********************************
sub medianaDatiContinui
{
    my $ref1 = $_[0]; my @continuousVariables = @{$ref1};
    my $ref2 = $_[1]; my @ID = @{$ref2};
    my $ref3 = $_[2]; my %matrix = %{$ref3};
    my $ref4 = $_[3]; my %indexContinuous = %{$ref4};
    
    my %medianCVariables;  # ---> hash: covariatacontinua -> valore mediano
    my %minimoCVariables;  # ---> hash: covariatacontinua -> valore minimo
    my %massimoCVariables; # ---> hash: covariatacontinua -> valore massimo
    
    my @ValoreCovariata;

    for(my $j=0; $j<=$#continuousVariables; $j++)
    {
        @ValoreCovariata= ();
        my $o=0;
        my $a = $continuousVariables[$j];$a= trim($a);
        
        for(my $i=1; $i< $#ID; $i++)
        {
            if($ID[$i-1] != $ID[$i])
            {
                $ValoreCovariata[$o] = $matrix{$i}{$indexContinuous{$a}};
                $o++;
            }
        }
        $medianCVariables{$a} = median(@ValoreCovariata);
        $minimoCVariables{$a} = min(@ValoreCovariata);
        $massimoCVariables{$a} = max(@ValoreCovariata);
    }

    return (\%medianCVariables, \%minimoCVariables, \%massimoCVariables);
}


# ******************************** FrequenzaDatiCategorici ********************************
sub FrequenzaDatiCategorici
{
    my $ref1 = $_[0]; my @categoricalVariables = @{$ref1};
    my $ref2 = $_[1]; my @ID = @{$ref2};
    my $ref3 = $_[2]; my %matrix = %{$ref3};
    my $ref4 = $_[3]; my %indexCategorical = %{$ref4};
    
    my @ColonnaCategorical;
    my @counts;
    my %FreqCov;
    my %FreqTot;
    my $indexMax;
    my @varfjoin= ();
    my %ValUnique;
    
    for(my $j=0; $j<=$#categoricalVariables; $j++){ # ciclo sulle covariate continue
        
        @ColonnaCategorical= ();
        my $a = $categoricalVariables[$j];$a=~ s/\n//;$a=~ s/\s+//;
        my $o=0;
        
        for(my $i=1; $i< $#ID; $i++)       # ciclo sulle ID dei pazienti
        {
            if($ID[$i-1] != $ID[$i])
            {
                # valori della covariata per ogni paziente
                $ColonnaCategorical[$o] = $matrix{$i}{$indexCategorical{$a}};
                $o++;
                
            }
        }
        @varfjoin = grep {$_ =~ /[\d*]/} @ColonnaCategorical;
        my $varfjoin = join("", @varfjoin);trim($varfjoin);
        
        # valori unici per ogni covariata
        my @unico =(uniq(@ColonnaCategorical));
        # Mi serve per capire quanti valori di teta iniziali aggiungere
        $ValUnique{$categoricalVariables[$j]} = $#unico;
        
        
        for(my $q=0; $q<= $#unico; $q++)
        {
            $FreqTot{$a}[$q] = $unico[$q];# HASH: covariata -> valori unici della covariata
            
            my $count0=0;
            
            for(my $s=0; $s<= $#ColonnaCategorical; $s++)# ciclo sui valori per paziente
            {
                
                if($ColonnaCategorical[$s] == $unico[$q])
                {
                    $count0++;         # conta la frequenza di ogni valore unico
                }
            }
            $counts[$q] = $count0;    # tiene traccia delle frequenze di tutte le covariate
        }
        
        my $m = max(@counts);          # calcola massimo della frequenza in ogni cov
        
        for(my $q=0; $q<= $#unico; $q++)
        {
            if($counts[$q] == $m)
            {
                $indexMax = $q;                    # indice della frequenza massima
            }
        }
        $FreqCov{$a} = $unico[$indexMax]; # HASH: covariata: valore più frequente della stessa covariata
    }
    #print Dumper \%FreqTot;
    #print Dumper \%FreqCov;
    #print("Il valore più frequnete per la covariata FORM è: $FreqCov{HLTH}\n");
    return(\%FreqTot, \%FreqCov, \%ValUnique);
}


# ******************************** DS_continuous ********************************
# Stampa DEFINITION START per le variabili continue
sub DS_continuous
{
    my $gene = $_[0];
    my $covariata = $_[1];
    my $maxTheta = $_[2];
    my $mediana = $_[3]; $mediana =~ s/,/./;
    my $param = $_[4];
    my $FILE = $_[5];
    my $relationsRef = $_[6]; my @relations = @$relationsRef;
    my $AmGeneRef = $_[7]; my @AmGene =  @$AmGeneRef;
    
    my %relazioni;
    my $VectDS_continuous; # elemento parametro-covariata a cui associo un DEFINITION START
    my $r = tie %relazioni, 'Tie::IxHash';
    
    for(my $i=0; $i<=$#relations; $i++)
    {
        $relations[$i] = trim($relations[$i]);
        $relazioni{$AmGene[$i]} =  $relations[$i];
    }
    #print Dumper \%relazioni;
    
    given($gene)
    {
        # **** 1: NON C'è RELAZIONE  PARAMETRO - COVARIATA ***
        #non è necessario scrivere nulla
        
        #    **** 2: RELAZIONE LINEARE ****
        when($relazioni{$gene} eq 2)
        {
            $maxTheta = $maxTheta +1;
            print $FILE(";;; $param$covariata-DEFINITION START\n$param$covariata = (1+THETA($maxTheta)*($covariata-$mediana))\n;;;$param$covariata-DEFINITION END\n\n");
            $VectDS_continuous = "$param$covariata";
            return ($maxTheta, \%relazioni, $VectDS_continuous);
        }
        
        #    **** 3: RELAZIONE PIECE-WISE LINEAR ****
        when($relazioni{$gene} eq 3)
        {
            $maxTheta = $maxTheta +1;
            
            print $FILE(";;; $param$covariata-DEFINITION START\nIF($covariata.LE.$mediana) $param$covariata=(1+THETA($maxTheta)*($covariata-$mediana))\n");
            $maxTheta = $maxTheta +1;
            print $FILE("IF($covariata.GT.$mediana) $param$covariata=(1+THETA($maxTheta)*($covariata-$mediana))\n;;;$param$covariata-DEFINITION END\n\n");
            $VectDS_continuous = "$param$covariata";
            return ($maxTheta, \%relazioni, $VectDS_continuous);
        }
        
        
        #     **** 4: RELAZIONE ESPONENZIALE ****
        when($relazioni{$gene} eq 4)
        {
            $maxTheta = $maxTheta +1;
            print $FILE(";;; $param$covariata-DEFINITION START\n$param$covariata =EXP(THETA($maxTheta)*($covariata - $mediana))\n;;;$param$covariata-DEFINITION END\n\n");
            $VectDS_continuous = "$param$covariata";
            return ($maxTheta, \%relazioni, $VectDS_continuous);
        }
        
        
        #     **** 5: RELAZIONE POTENZA ****
        when($relazioni{$gene} eq 5)
        {
            $maxTheta = $maxTheta +1;
            print $FILE(";;; $param$covariata-DEFINITION START\n$param$covariata = (($covariata/$mediana)**THETA($maxTheta))\n;;;$param$covariata-DEFINITION END\n\n");
            $VectDS_continuous = "$param$covariata";
            return ($maxTheta, \%relazioni, $VectDS_continuous);
        }
    }
}

# ********************************* DS_categorical **********************************
# Stampa DEFINITION START per le variabili categoriche

sub DS_categorical
{
    my $param = $_[0];
    my $cov_categorical = $_[1];
    my $frequenza = $_[2];
    my $maxTheta = $_[3];
    my $gene = $_[4];
    my $FILE = $_[5];
    my $FreqTotRef = $_[6]; my %FreqTot = %{$FreqTotRef};
    my $FreqCovRef = $_[7]; my %FreqCov = %{$FreqCovRef};
    
    my $s;
    my @categoria;
    my @categorie=@{$FreqTot{$cov_categorical}}; # hash covariata -> categorie. Es. { HLTH -> [0 , 1] }
    my $VectDS_categorical;
    
    #    **** 1: RELAZIONE LINEARE ****
    
    print $FILE(";;; $param$cov_categorical-DEFINITION START\n");
    print $FILE("IF($cov_categorical.EQ.$frequenza)$param$cov_categorical= 1 ; Most common \n");
    
    for($s=0; $s<= $#categorie; $s++)
    {
        $categoria[$s] = $FreqTot{$cov_categorical}[$s];
        $categoria[$s] = trim($categoria[$s]);
        
        if($categoria[$s] != $FreqCov{$cov_categorical}) #FreqCov contiene la categoria più frequente
        {
            $maxTheta++;
            print $FILE("IF($cov_categorical.EQ.$categoria[$s])$param$cov_categorical = (1 + THETA($maxTheta))  \n");
        }
    }
    print $FILE(";;; $param$cov_categorical-DEFINITION END\n\n");
    $VectDS_categorical = "$param$cov_categorical";
    
    #    **** 0: NESSUNA RELAZIONE ****
    # non c'è bisogno di inserire alcuna relazione nel modello
    
    return($maxTheta, $VectDS_categorical);
}

# ****************************************** genome **************************************
sub genome
{
    my $numberofStates = $_[0];
    
    my $numberofBits = log2($numberofStates);
    
    my @AmGeni;
    
    given($numberofBits)
    {
        when($numberofBits==1)
        {
            @AmGeni = ('0', '1');
        }
        when($numberofBits>1 & $numberofBits<2)
        {
            @AmGeni = ('00', '01', '10');
        }
        when($numberofBits==2)
        {
            @AmGeni = ('00', '01', '10', '11');
        }
        when(log2($numberofStates)==log2(5))
        {
            @AmGeni = ('000', '001', '010', '011', '100');
        }
    }
    return @AmGeni;
}

# ******************************** write_NONMEMfile ********************************

sub write_NONMEMfile
{
    my $ref1 = $_[0]; my @param = @$ref1;
    my $ref3 = $_[1]; my @AmGeni = @$ref3;
    my $ref4 = $_[2]; my @AmGeniCat = @$ref4;
    my $ref5 = $_[3]; my %relationstart = %$ref5;
    my $maxTheta = $_[4];
    my $ref6 = $_[5]; my %medianCVariables = %$ref6;
    my $FILE = $_[6];
    my $ref7 =$_[7]; my @ValidStatesContinuous = @$ref7;
    my $ref8 = $_[8]; my @cromosomaIN = @$ref8;
    my $ref9 = $_[9]; my %FreqCov = %$ref9;
    my $ref10 = $_[10]; my %FreqTot = %$ref10;
    my $ref11 = $_[11]; my @gene = @$ref11;
    my $nAmGeni = $_[12];
    my $indexPK = $_[13];
    my $indexOmega = $_[14];
    my $ref12 = $_[15]; my @line = @$ref12;
    my $G = $_[16];
    my $ref13 = $_[17]; my @Covariate = @$ref13;
    my $ref14 = $_[18]; my %hash = %$ref14;
    my $ref15 = $_[19]; my @continuousVariables = @$ref15;
    my $ref16 = $_[20]; my @categoricalVariables = @$ref16;
    my $ref17 = $_[21]; my %VarUnique = %$ref17;
    my $ref18 = $_[22]; my %minimoCVariables = %$ref18;
    my $ref19 = $_[23]; my %massimoCVariables = %$ref19;

    my $nrelazioniContinue;
    my $k=0;
    my $nelementi;
    
    my %ParamCov_Gene;
    my %relazioni;
    my @VectDSContinuous;
    my $RelDSContinuous;
    my %paramCov_Relazione;
    
    # Inizializza cromosoma
    my @cromosoma =();
    my $nAmGeniCat = $#AmGeniCat+1;
    
    #****************************************************************************************
    #                         DEFINITION START COVARIATE CONTINUE                           #
    #****************************************************************************************
    for(my $j=0; $j<=$#Covariate; $j++)                           # Start iterations over covariates
    {
        my $covariata = $Covariate[$j];
        $covariata = trim($covariata);
        my $indiceCovariata = 0;
        my $refContinuous = trim(@continuousVariables);
        @continuousVariables = @$refContinuous;
        
        if (exists $hash{$covariata} && any{$_ eq $covariata} @continuousVariables)
        {
            my $refCov = $hash{$covariata};
            my @valore1 = @{$refCov};
            my $relazioniRef;
            
            #  Calcola quanti parametri sono legati alla covariata
            $nelementi = $#valore1+1;
            
            #@VectDSContinuous =();
            for(my $i=0; $i<$nelementi; $i++) # iterazioni sui parametri
            {
                # Seleziona un parametro
                my $parametro = $hash{$covariata}[$i];
                
                #****************************************************************************************
                #                                    find gene                                          #
                #****************************************************************************************
                if ($G==0)  #Genarazione 0: popolazione iniziale casuale
                {
                    my $gene_casuale = irand($nAmGeni);                         # Pesco un indice casuale da 1 a nAmGeni
                    $gene[$i] = $AmGeni[$gene_casuale];                         # Pesco un gene casuale dal pool costrutio
                    
                }
                else       # Altre generazioni non inizializzate in modo casuale
                {
                    $gene[$i] = $cromosomaIN[$k];
                    $k++;
                }
                my $name = "$parametro$covariata";
                $ParamCov_Gene{$name} = $gene[$i];
                
                #****************************************************************************************
                #                     find relations for continuous covariates                          #
                #****************************************************************************************
                unless ($gene[$i] eq $AmGeni[0]) # se il gene è uguale a 0 - 00 - 000 non c'è nessuna relazione
                {
                    #$relationstart : $parametro -> $covariate
                    if(exists $relationstart{$parametro})
                    {
                        my @covariatePerParametro= @{$relationstart{$parametro}};
                        $indiceCovariata = $#covariatePerParametro;
                        $indiceCovariata++;
                    }
                    else
                    {
                        $indiceCovariata = 0;
                    }
                    
                    $relationstart{$parametro}[$indiceCovariata] = $covariata;
                    
                    # DS_continuous
                    ($maxTheta, $relazioniRef, $RelDSContinuous) = DS_continuous($gene[$i],$covariata,$maxTheta,$medianCVariables{$covariata} ,$parametro, $FILE, \@ValidStatesContinuous, \@AmGeni);
                     %relazioni = %$relazioniRef;
                    @VectDSContinuous = (@VectDSContinuous, $RelDSContinuous);
                    $paramCov_Relazione{$covariata.$parametro} = $relazioni{$gene[$i]};
                    
                    
                }
                if ($gene[$i] eq $AmGeni[0])
                {
                    $paramCov_Relazione{$covariata.$parametro} = 0;
                }
                
                
                # Accodo il gene al cromosoma
                @cromosoma = (@cromosoma, $gene[$i]);
            }
            
        }
        $nrelazioniContinue = $#cromosoma;
        
        #****************************************************************************************
        #                          DEFINITION START COVARIATE DISCRETE                          #
        #****************************************************************************************
        #my $ix;
        #my $ixx = 0;
        #print Dumper \%hash;
        if (exists $hash{$covariata} && any{$_ eq $covariata} @categoricalVariables)
        {
            my $refCov = $hash{$covariata};
            my @valore1 = @{$refCov};
            my $relazioniRef;
            
            #  Calcola quanti parametri sono legati alla covariata
            my $nelementi = $#valore1+1;
            
            for(my $i=0; $i<$nelementi; $i++)                 # ciclo sul numero dei parametri
            {
                my $parametro = $hash{$covariata}[$i];
                my $ixx = 0;
                #****************************************************************************************
                #                                    find gene                                          #
                #****************************************************************************************
                if($G==0)
                {
                    my $gene_casuale = irand($nAmGeniCat);
                    $gene[$i] = $AmGeniCat[$gene_casuale];# gene
                }
                else
                {
                    $gene[$i] = $cromosomaIN[$k];
                    $k++;
                }
                
                unless($gene[$i] eq '0')
                {
                    
                    #****************************************************************************************
                    #                     find relations for categorical covariates                         #
                    #****************************************************************************************
                
                    ($maxTheta, my $reldscate)= DS_categorical($parametro,$covariata, $FreqCov{$covariata}, $maxTheta, $gene[$i],$FILE, \%FreqTot, \%FreqCov);
                   @VectDSContinuous = (@VectDSContinuous, $reldscate);
                    
                    if (exists $relationstart{$parametro})
                    {
                        my @relazionecovparam = @{$relationstart{$parametro}};
                        my $ix = $#relazionecovparam;
                        $ix = $ix + 1;
                        $relationstart{$parametro}[$ix]= $covariata;
						
                    }
                    else
                    {
                        $relationstart{$parametro}[$ixx]= $covariata;
                        $ixx++;
                    }
                    
                }
                
                $paramCov_Relazione{$covariata.$parametro} = $gene[$i];
    
                # Continua cromosoma
                @cromosoma = (@cromosoma, $gene[$i]);
            }
        }
    }
   
    my @c= @cromosoma;
    #print Dumper \%paramCov_Relazione;
    #print("Cromosoma: @cromosoma\n");
    #****************************************************************************************
    #                              2. PRINT RELATIONS START                                 #
    #****************************************************************************************
    printRelationStart(\@param, \%relationstart, $FILE);
    
    #****************************************************************************************
    #                            3. PRINT NEW MODEL RELATIONS                               #
    #****************************************************************************************
    printNewModelRelations(\%relationstart, \@line, $FILE, $indexPK, $indexOmega);
    
    #****************************************************************************************
    #                              4. PRINT INITIAL THETA                                   #
    #****************************************************************************************
    printInitialTheta(\@param, $FILE, \@categoricalVariables, \@continuousVariables, \%VarUnique, \%ParamCov_Gene, \%relazioni, \%minimoCVariables, \%massimoCVariables, \%medianCVariables, \@VectDSContinuous);
    
    print $FILE("@line[$indexOmega..$#line]\n");
    return (\%relationstart, \@c, $maxTheta, $nrelazioniContinue, \%paramCov_Relazione);
}

# ************************************* printRelationStart **************************************

sub printRelationStart
{
    my $ref1= $_[0]; my @param = @{$ref1};
    my $ref2 = $_[1]; my %relationstart = %{$ref2};
    my $FILE = $_[2];
    #print Dumper \%relationstart;
    
    for(my $j=0; $j<=$#param; $j++)
    {
        if (exists $relationstart{$param[$j]})
        {
            
            my @refe = @{$relationstart{$param[$j]}};
            
            
            print $FILE(";;;$param[$j]-RELATION START\n$param[$j]COV=");
            
            for(my $i=0; $i<= $#refe; $i++)
            {
                if($i != $#refe)
                {
                    print $FILE("$param[$j]$relationstart{$param[$j]}[$i]*");
                }
                else
                {
                    print $FILE("$param[$j]$relationstart{$param[$j]}[$i]\n");
                }
            }
            print $FILE(";;;$param[$j]-RELATION END\n\n");
        }
    }
}

# ************************************* printNewModelRelations **************************************

sub printNewModelRelations
{
    my $ref1 = $_[0]; my %relationstart = %$ref1;
    my $ref2= $_[1]; my @line = @$ref2;
    my $FILE = $_[2];
    my $indexPK = $_[3];
    my $indexOmega = $_[4];
    
    my @valhash = keys(%relationstart); # Parametri
    
    for(my $i=0; $i<= $#valhash; $i++) # Ciclo sui parametri
    {
        my $parametro = $valhash[$i];
        
        for (my $riga = 0; $riga<=$#line; $riga++)
        {
            if($line[$riga] =~ m/TV$parametro=THETA\(\d+\)/)
            {
                # Rewrite relation
                my @rewrite = ("$line[$riga]\n","TV$parametro=$parametro","COV*TV$parametro\n");
                $line[$riga] = join("", @rewrite);
            }
        }
    }
    
    # stampa il modello
    print $FILE("@line[$indexPK+1..$indexOmega-1]");
}

# *********************************** printInitialTheta **************************************

sub printInitialTheta
{
    # Per ogni relazione covariata - parametro devo assegnare i valori iniziali del teta corrispettivo
    # In particolare in PsN si presentano 3 valori: (limite inferiore, valore iniziale, limite superiore)
    # Il valore minimo può assumere valori negativi perché indica l'effetto che la covariata ha sul parametro
    
    # N.B. La cosa importante è che i theta vengano riportati nello stesso ordine con cui i theta vengano riportati nel modello
    
    my $refParam = $_[0]; my @parametri = @{$refParam};
    my $FILE = $_[1];
    my $refCategorical = $_[2]; my @categoricalVariables = @$refCategorical;
    my $refContinuos = $_[3];my @continuousVariables = @$refContinuos;
    my $cardinalitaCategorical = $_[4]; my %VarUnique = %{$cardinalitaCategorical};
    my $refParamCovGene = $_[5]; my %ParamCov_Gene = %$refParamCovGene;
    my $refRelazioni = $_[6]; my %Relazioni = %$refRelazioni;
    my $minref = $_[7]; my %minimoCVariables = %$minref;
    my $maxref = $_[8]; my %massimoCVariables = %$maxref;
    my $medianref = $_[9]; my %medianCVariables = %$medianref;
    my $dsref= $_[10]; my @VectDSCont= @$dsref;
    
    my $k=0;
    # Le iterazioni sono condotte sui parametri perché i parametri vengono collezionati leggendo il file. In questo modo mantengo l'ordine dei tet
 
    my @covariate = (@categoricalVariables, @continuousVariables);
    
    for(my $ix=0; $ix<=$#VectDSCont; $ix++)
    {
        for(my $i=0; $i<= $#covariate; $i++)
        {
            #print("Covariata: @covariata");
            my $covariata = $covariate[$i];
            for(my $j=0; $j<=$#parametri; $j++)
            {
                my $name = "$parametri[$j]$covariate[$i]";
                
                if($name eq $VectDSCont[$ix])
                {
                    ######################## Controllo se la covariata è CATEGORICA ########################
                    if(any {$_ eq $covariate[$i]} @categoricalVariables)
                    {
                        #print("   > La covariata è categorica!\n");
                        # Le covariate categoriche possono avere due o più livelli. È necessario aggiungere un teta iniziale per i diversi livelli che lo richiedono
                        for(my $q = 0; $q < $VarUnique{$covariate[$i]}; $q++)
                        {
                            my $NumeroVariabile = $q+1;
                            print $FILE ("\$THETA (-1,-0.001,5); $parametri[$j]$covariate[$i]$NumeroVariabile\n");
                        }
                    }
 
                    ######################## Controllo se la covariata è CONTINUA ########################
                    if(any {$_ eq $covariate[$i]} @continuousVariables)
                    {
                        #print("   > La covariata è continua!\n");
                        # Trovo il gene associato alla coppia covariata - parametro
                        my $gene = $ParamCov_Gene{$name};
 
                        my $relationParamCov = $Relazioni{$gene};
                        #print("   > Relazione parametro covariata: $relationParamCov\n");
                        
                        # [ lower_bounds ] : 1 / (median - maximum)
                        my $estremo_inferiore =FormatShort((-1)/($massimoCVariables{$covariate[$i]} - $medianCVariables{$covariate[$i]}));
                        
                        # [ initial_value ] : 0.001/(mediana - minimo)
                        #my $valore_iniziale = FormatShort($minimoCVariables{$covariate[$i]} *(10**(-5)));
                        my $valore_iniziale = FormatShort((-0.001)/($minimoCVariables{$covariate[$i]} - $medianCVariables{$covariate[$i]}));
                        #print("Valore iniziale: $valore_iniziale\n");
                        
                        # [ upper_bounds ] : 1 / (median - minimum)
                        my $estremo_superiore = FormatShort((-1)/($minimoCVariables{$covariate[$i]} - $medianCVariables{$covariate[$i]}));
                        
                        #print(" \n--- $covariate[$i] ---\nestremo inferiore $estremo_inferiore\nvalore iniziale $valore_iniziale\nestremo superiore $estremo_superiore\n");
                        
                        my $valHS = 1;
                        
                        given($relationParamCov)
                        {
                            # when ('1'): no relation!
                            
                            when('2') # Linear Model
                            {
                                #print("   > Relazione lineare\n");
                                print $FILE ("\$THETA ($estremo_inferiore,$valore_iniziale,$estremo_superiore); $parametri[$j]$covariate[$i]\n");
                            }
 
                            when('3') # Hockey - Stick or piece - wise linear
                            {
                                #print("   > Relazione hockey stick\n");
                                print $FILE ("\$THETA (-10,$valore_iniziale,$estremo_superiore); $parametri[$j]$covariate[$i]1\n");
                                my $valore_iniziale2 = FormatShort((-0.001)/($massimoCVariables{$covariate[$i]} - $medianCVariables{$covariate[$i]}));
                                if($valore_iniziale2 == 0)
                                {
                                	$valore_iniziale2 = 0.001;
                                }
                                print $FILE ("\$THETA ($estremo_inferiore,$valore_iniziale2,10); $parametri[$j]$covariate[$i]2\n");
                                
                            }
 
                            when('4') # Relazione esponenziale
                            {
                                #print("   > Relazione esponenziale\n");
                                
                                    #print("   > MI TROVO NEL MATCH CON CL\n");
                                    print $FILE ("\$THETA (-10,0.001,10); $parametri[$j]$covariate[$i]\n");
                            }
                            
                            when('5') # Power Model
                            {
                                #print("   > Relazione power\n");
                                if($parametri[$j] =~ m/CL/) # Per la cleareance
                                {
                                    #print("   > MI TROVO NEL MATCH CON CL\n");
                                    print $FILE ("\$THETA (-4,0.75,4); $parametri[$j]$covariate[$i]\n");
                                }
                                elsif($parametri[$j] =~ m/V/) # Per il volume
                                {
                                    #print("   > MI TROVO NEL MATCH CON V \n");
                                    print $FILE ("\$THETA (-4,1,4); $parametri[$j]$covariate[$i]\n");
                                }
                                else
                                {
                                    print $FILE ("\$THETA (-10,0.001,10); $parametri[$j]$covariate[$i]\n");
                                }
                            }# close when '5'

                        } #close given
                    } #close if any continuous
                } #close if name
            } # close for parameters
        } # close for covariates
    } # close for Vect cont
} # close subrutine



# ******************************** count_x_Position ########################################

sub count_x_Position
{
    my $e = $_[0];
    my $numberofStates = $_[1];
    my $CromosomaRef = $_[2]; my @cromosoma = @{$CromosomaRef};
    my $WholeGenesRef = $_[3]; my @WholeGenes = @{$WholeGenesRef};
    my $CountRef = $_[4]; my %count = %{$CountRef};
    
    my %countResults;
    
    unless($e!=0) # Solo alla prima iterazione inizializzo (e conta il numero di cromosomi)
    {
        for(my $index=0; $index<=$numberofStates+1; $index++)
        {
            for(my $jndex=0; $jndex<=$#cromosoma; $jndex++)
            {
                $count{$WholeGenes[$index]}[$jndex] = 0; # inizializzo il vettore count1 come un vettore di 0
            }
        }
    }
    #print Dumper \%count;
    #print("Numero geni x cromosomi $#cromosoma\n\n\n");#print Dumper \%count;
    
    # Calcolo le frequenze dei diversi geni per ogni cromosoma
    
    for(my $js= 0 ; $js<=$#cromosoma; $js++)
    {
        my $countResultsRef = countGenes($cromosoma[$js], \@WholeGenes, \%count, $js);
        %countResults = %{$countResultsRef};
        %count = %countResults;
    }
    return \%countResults;
}

# ******************************** COUNT GENES ########################################
sub countGenes
{
    my $geni = $_[0];
    my $AmGenesRef = $_[1]; my @AmGenes = @{$AmGenesRef};
    my $CountRef = $_[2]; my %Count = %{$CountRef};
    my $j = $_[3];
    my $i;
    
    for($i=0; $i<=$#AmGenes; $i++)
    {
        if($geni eq $AmGenes[$i])
        {
            $Count{$AmGenes[$i]}[$j] = $Count{$AmGenes[$i]}[$j]+1;
        }
        
    }
    
    #print Dumper \%Count;
    my $results =  \%Count;
    return $results;
    
}

# ******************************** GeneMostFrequent_x_Position ********************************
sub GenesMostFrequent_x_Position
{
    my $WholeGenesRef = $_[0]; my @WholeGenes = @{$WholeGenesRef};
    my $countRef = $_[1]; my %count = %{$countRef};
    
    my @geneMax;
    
    # Count max contiene le frequenze del gene 00 in tutte le posizioni del cromosma
    # Es. @countMax = (1 2 1 0 0) -> il gene 00 in posizione 0 si ripete una volta. In posizione
    # 1 si ripete 2 volte, in posizione 2 si ripete una volta. Nelle altre posizioni il gene
    # non compare mai
    # La lunghezza di @countMax = numero di geni nel cromosoma
    my @countMax = @{$count{$WholeGenes[0]}};
    
    # Frequenza del primo gene (es. 00)
    for(my $i=0; $i<=$#countMax; $i++)
    {
        $geneMax[$i]= $WholeGenes[0];
    }
    
    
    for(my $i=1; $i<=$#WholeGenes; $i++) # iterazioni sui geni che ho tranne il primo (es. 01, 10, 11)
    {
        for(my $j = 0; $j<=$#countMax; $j++)   # iterazioni sulle colonne (numero di geni nel cromosoma)
        {
            # Il conteggio non deve essere maggiore o uguale, ma è necessario anche che
            # la fitness associata a quel cromosoma sia maggiore
=a
            given($count{$WholeGenes[$i]}[$j])
            {
                when($count{$WholeGenes[$i]}[$j] > $countMax[$j])
                {
                    $countMax[$j] = $count{$WholeGenes[$i]}[$j];
                    $geneMax[$j] = $WholeGenes[$i];
                }
                when($count{$WholeGenes[$i]}[$j] < $countMax[$j])
                {
                    $countMax[$j] = $countMax[$j];
                    $geneMax[$j] = $geneMax[$j];
                }
                when($count{$WholeGenes[$i]}[$j] == $countMax[$j])
                {
                    if(any ($_ eq $WholeGenes[$i]) @WholeGenes)
 {
 }
                }
                
            }
=cut
            if($count{$WholeGenes[$i]}[$j] >= $countMax[$j])
            {
            
                $countMax[$j] = $count{$WholeGenes[$i]}[$j];
                $geneMax[$j] = $WholeGenes[$i];
                
            }
            else
            {
                $countMax[$j] = $countMax[$j];
                $geneMax[$j] = $geneMax[$j];
            }
            
        }
    }
    return @geneMax;
}

# ******************************** RUN NONMEM ********************************
sub runNONMEM
{
    my $e = $_[0];
    my @args = ("execute", "Cromosoma$e.mod" , "-dir=resultDirectory$e");
    system(@args)==0 or die "Error occurred in function runNONMEM for chromosome $e $_\n"; # prova con break (passa a run successivo)
    #or die "Error occurred in function runNONMEM for chromosome $e $_\n";
}

# ******************************** RouletteWheel ********************************

# La selezione dei genitori è il compito di allocare opportunità riproduttive a ciascun individuo
# Roulette Wheel : i genitori sono selezionati in base alla loro fitness: i migliori cromosomi hanno maggiore probabilità di essere selzionati
#La funzione deve prendere in ingresso il vettore delle fitness e deve restituire l indice del cromosoma scelto

sub RouletteWheel
{
    my $ref1 = $_[0]; my %FitnessHash = %{$ref1};
    my $ref2 = $_[1]; my %popolazione = %{$ref2};
    
    # Vettore delle fitness
    my @FitnessVect = values %FitnessHash;
    
    # Calcolo la somma delle fitness di tutti gli individui nella popolazione
    my $FitnessSum = sum(@FitnessVect);
    
    # Calcolo l'ampiezza di ogni slice della Roulette Wheel
    # Immagazzino questa informazione in un vettore perché devo avere un valore di
    # ampiezza per ogni individuo nella popolazione
    my @sliceWidth;
    for(my $i=0; $i<=$#FitnessVect; $i++)
    {
        $sliceWidth[$i] = $FitnessVect[$i] / $FitnessSum;
    }
    
    # Bisgona posizionare poi ogni cromosoma sulla roulette wheel
    # vogliamo calcolare qual è il valore degli UPPER BOUND per fetta (slice)
    # Anche in questo caso devo costruire un vettore: ad ogni cromosoma devo associare un bound
    my @WheelBounds;
    for(my $i=0; $i<=$#FitnessVect; $i++)
    {
        $WheelBounds[$i] = sum(@sliceWidth[0..$i]);
    }
    $WheelBounds[-1] = 1; # Ultimo bound
    
    # Genero numero casuale
    my $nrand = rand(1);
    
    # L'individuo selezionato è il primo che ha il valore di fitness cumulata (somma delle
    # fitness di tutti i cromosomi fino a lui) maggiore del numero casuale generato
    my $ParentIndex;
    for(my $i=0; $i<=$#WheelBounds; $i++)
    {
        if($WheelBounds[$i]>$nrand)
        {
            $ParentIndex = $i;
            last;
        }
    }
    
    #print("Indice del genitore $ParentIndex\n");
    my @ParentCromosom = @{$popolazione{"$ParentIndex"}};
    @ParentCromosom = split('', join('', @ParentCromosom));
    
    #print("Genitore: @ParentCromosom\n");
    return \@ParentCromosom;
}
#****************************** Tournament Selection ********************************

sub TournamentSelection
{
    my $popSize = $_[0];
    my $ref1 = $_[1]; my %FitnessHash = %$ref1;
    my $ref2 = $_[2]; my %popolazione = %$ref2;
    
    # Parametro del Tournament Selection
    my $parametroTS = 0.75;
    my $IndexGenitore;
    
    # Selezione di due indici casuali nella popolazione
    my $indice1 = irand($popSize);
    my $indice2 = irand($popSize);
    
    my $FitnessBest = max($FitnessHash{$indice1}, $FitnessHash{$indice2});
    
    # Scelgo numero casuale da 0 a 1
    my $randomValue = rand(1);
    
    given($randomValue)
    {
        when($randomValue < $parametroTS)
        {
            # Scelta del genitore con fitness maggiore
            if($FitnessHash{$indice1}>=$FitnessHash{$indice2})
            {
                $IndexGenitore = $indice1;
            }
            else
            {
                $IndexGenitore = $indice2;
            }
        }
        
        when($randomValue >= $parametroTS)
        {
            # Scelta del genitore con fitness minore
            if($FitnessHash{$indice1}>=$FitnessHash{$indice2})
            {
                $IndexGenitore = $indice2;
            }
            else
            {
                $IndexGenitore = $indice1;
            }
        }
    }
    my @ParentCromosom = @{$popolazione{$IndexGenitore}};
    @ParentCromosom = split('', join('', @ParentCromosom));
    return(\@ParentCromosom);
}


# ******************************** onePointCrossover ********************************

sub onePoinCrossOver
{
    my $ref1 = $_[0]; my @ParentCromosom1 = @{$ref1};
    my $ref2 = $_[1]; my @ParentCromosom2 = @{$ref2};
    
    # Seleziono casualmente l'indice di crossover
    # NB Bisogna tener conto che ci sono meno punti per il crossover rispetto al numero di elementi
    my $index = irand($#ParentCromosom1);
    #print("\nIndice di crossover: $index\n");
    
    # PROLE
    
    my @offspring1 = (@ParentCromosom1[0..$index] , @ParentCromosom2[$index+1..$#ParentCromosom2]);
    #print("\n[PROLE]\n");
    #print("Primo Figlio: @offspring1\n");
    my @offspring2 = (@ParentCromosom2[0..$index] , @ParentCromosom1[$index+1..$#ParentCromosom1]);
    #print("Secondo Figlio: @offspring2\n");
    return(\@offspring1, \@offspring2);
}


# ******************************** splitNewpop ********************************

sub splitNewpop
{
    my $nrelazioniContinue = $_[0];
    my $countsOffspring = $_[1];
    my $nGene = $_[2];
    my $ref1 = $_[3]; my %NewPopolazione = %{$ref1};
    my $ref2 =$_[4]; my @MostFrequent = @$ref2;
    my $ref3 = $_[5]; my @AmGenes = @$ref3;
    my $numberofBits = $_[6];
    
    
    my %NuovaPop;
    my $np2 = tie %NuovaPop, 'Tie::IxHash';
    
    for (my $i=0; $i< $countsOffspring; $i++) #ciclo sul numero di individui nella popolazione
    {
        my $ac=0;
        my $cromosoma = $NewPopolazione{"$i"};
        
        # Sostituzione con gene più frequente tra i cromosomi con miglior fitness
        for(my $j=0; $j<=$nrelazioniContinue; $j++) # ciclo sul numero di relazioni continue
        {
            $NuovaPop{"$i"}[$j] = substr $cromosoma, $ac, $numberofBits;
            
            if  (! any {$_ eq "$NuovaPop{$i}[$j]"} @AmGenes)
            {
                $NuovaPop{"$i"}[$j] = $MostFrequent[$j];# sostituisci con quello più frequente
            }
            
            $ac = $ac + $numberofBits;
        }
        
        for(my $j=$nrelazioniContinue+1; $j<=$nGene; $j++ )
        {
            $NuovaPop{"$i"}[$j] = substr $cromosoma, $ac, 1; #Per categoriche
            $ac++;
            
        }
    }
    
    return \%NuovaPop;
    
}



# ******************************** findFitness ********************************
sub findFitness
{
    my $e = $_[0];
    my $maxTheta = $_[1];
    my $MAXTHETA = $_[2];
    my $ref2 = $_[3]; my %FitnessHash= %{$ref2};
    
    $CWD ="resultDirectory$e";    # entra nella cartella dei risultati
    
    my $csvname = "raw_results_Cromosoma$e.csv";
    open FILE ,'<' ,$csvname , or die "Error occurred: the file raw_results_Cromosoma$e.csv does not exists \n";
    my @rawResultFile = <FILE>;
    close FILE;
    #print("@rawResultFile\n");
    # Nome delle variabili
    my $Resultsnames = $rawResultFile[0];
    my @splitResultNames = split(/,\D/, $Resultsnames);
    
    # Valore dei risultati
    my $Resultvalues = $rawResultFile[1];
    my @splitResultValues = split(/,/, $Resultvalues);
    
    # Creo un array associativo: associo ad ogni nome il suo valore
    my $var;
    my %Rawresults;
    
    for(my $i=0; $i<=$#splitResultNames; $i++)
    {
        $splitResultNames[$i] =~ s/[^a-zA-Z0-9]//g;
        $Rawresults{$splitResultNames[$i]} = $splitResultValues[$i];
    };
   
    # hash *indice cromosoma* -> objective function - numero di tetas stimati (considero tutti per il momento)
    
    # *******************************************************************************
    # ******************************** Fitness Value ********************************
    # *******************************************************************************
    
    # Per valutare la fitness di un cromosoma non valuto solamente il valore della fitness ma anche il numero di parametri utilizzati dal modello
    # In questo caso utilizzo quindi l'AIC (Akaike's information criterion) che tiene conto sia della bontà dell'adattamento  sia della complessità del modello)
    # Si preferiscono sempre modelli con AIC più basso
    my $numeroParametri = $maxTheta - $MAXTHETA;
    my $FitValue = $Rawresults{ofv};
    $FitnessHash{"$e"} =  (-1)*($FitValue+2*$numeroParametri);
    #print("\nObjective Function Cromosoma$e:\t\t$Rawresults{ofv}\n");
    my $result = $FitnessHash{"$e"};
    $CWD = '..';                  # Torna indietro nella cartella di partenza
    
    return ($result, $FitValue);
}

# ******************************** trim   ********************************************

#Elimina caratteri di spaziatura (spazio, tabulazione, new line ecc) da una stringa passata come input
sub trim
{
    my @input= @_;
    my $nelementi = $#input;
    
    given($nelementi)
    {
        # Trim di uno scalare
        when($nelementi eq 0)
        {
            my $string = $_[0];
    
            $string =~ s/\s+//;
            $string =~ s/\s*//;
            $string =~ s/\n//;
    
            return $string;
        }
        # Trim di un array
        when($nelementi > 0)
        {
            my @elemento;
            for(my $i=0; $i<=$nelementi; $i++)
            {
                $elemento[$i] = $_[$i];
                $elemento[$i] =~ s/\s+//;
                $elemento[$i] =~ s/\s*//;
                $elemento[$i] =~ s/\n//;
            }
            return \@elemento;
        }
    }
    
}


# ************************************* ZEROS *****************************************
# Subrutines per inizializzare i vettori con valori nulli
# Equivalente di zeros in Matlab
sub zeros{
    
    my @var = @_;
    if($#var == 0)
    {
        my $n = $_[0];
        my $i;
        my @vect;
        
        for($i=0; $i<$n; $i++)
        {
            $vect[$i] = 0 ;
        }
        
        return @vect;
    }
    
    else
    {
        
        die "Error occurred: number of input when call zeros subrutines must be = 1\n";
    }
}


# ******************************** LOGARITMO IN BASE 2 ***********************************
sub log2
{
    my $n = shift;
    if($n<=0){die "Error occurred in function log2: \$argument must be > 0 \n"};
    return log($n) / log(2);
}
# ******************************** CEIL: arrotonda ad intero superiore *********************
sub CEIL
{
    my @input = @_;
    if($#input > 0){die "Error Occurred in function CEIL: number of inputs must be eqaul to 1\n"}
    
    my $number = $_[0];
    
    if (($number)> int($number))
    {
        $number = int($number+1);
    }
    return $number;
}

# ***************************************** FormatShort ***********************************
# Prendo le prime 4 cifre decimali di un numero con la virgola
sub FormatShort
{
    my $numeroDecimale = $_[0];
    if($#_ >= 1)
    {  die "Error occurred in function FormatShort: too much input values!\n"  }
    
    my $fattoreMolt=0;
    
    if($numeroDecimale =~ m/\d+e\W?\d\d/) # Numero riportato in notazione scientifica (es. e-01)
    {
        $fattoreMolt = substr($numeroDecimale, index($numeroDecimale, 'e'), length($numeroDecimale)- index($numeroDecimale,'e'));
        $numeroDecimale = substr($numeroDecimale, 0, index($numeroDecimale, 'e')-1);   
     }

    #if($numeroDecimale =~ m/[a-zA-Z]/)
    #{  die "Error occurred in function FormatShort: input values are not valid!\n" }
    
    my $indiceVirgola = index($numeroDecimale, '.');
    my @numeroSplit = split(//, $numeroDecimale);
    my @parteDecimale = @numeroSplit[$indiceVirgola+1..$#numeroSplit];
    
    if ($indiceVirgola != -1 && $#parteDecimale >= 3)
    {
        my $newNumero = join('',@numeroSplit[0..$indiceVirgola+4]); 
        $newNumero = join('',($newNumero, $fattoreMolt));
        return($newNumero);
    }
    else
    {  return($numeroDecimale);  }
}
# **********************************  media ***********************************
# Calcolo della media aritmetica degli elementi di un vettore
sub media {
    my @input = @_;
    for(my $i=0; $i<$#input; $i++)
    {
        if ($input[$i] =~ m/[a-zA-Z]/)
        {
            my $n = $i+1;
            die "Error occurred in function media: function value n. $n is not valid!\n"
        }
    }
    return sum(@_)/@_;
}

# ********************************** var ***********************************
# Calcolo la varianza degli elementi di un vettore
sub varianza
{
    my @Vector = @_;
    
    for(my $i=0; $i<$#Vector; $i++)
    {
        if ($Vector[$i] =~ m/[a-zA-Z]/)
        {
            my $n = $i+1;
            die "Error occurred in function media: function value n. $n is not valid!\n"
        }
    }
    
    my $SUM = 0;
    
    for (my $i=0; $i<= $#Vector; $i++)
    {
        $SUM = $SUM + ( ($Vector[$i] - media(@Vector))**2);
    }
    
    my $VARIANCE = $SUM / ($#Vector+1);
    
    return $VARIANCE;
}



# ********************************** mirror ***********************************
# Specchia un vettore
sub mirror
{
    my @SubVect = @_;
    my @Mirror;
    my $j = $#SubVect;
    my $i = 0;
    
    while($i <= $#SubVect)
    {
        while($j >= 0)
        {
            $Mirror[$i] = $SubVect[$j];
            $j= $j-1;
            $i++;
        }
    }
    return @Mirror;
}

# **********************************  ascendHash  ***********************************
# Sort the keys of the hash according to their values
sub ascendHash
{
    my $HashRef = $_[0];
    my %Hash = %$HashRef;
    
    if (!ref(%Hash) eq "HASH")
    {
        die "Error occured in function ascendHash: subrutine argument is not valid\n";
    }
    
    my %HashOrdinato;
    my @IndiciOrdinati;
    my @ValoriOrdinati;
    my $i=0;
    
    foreach my $element(sort { $Hash{$a} <=> $Hash{$b}} keys %Hash)
    {
        $IndiciOrdinati[$i] = $element;
        $ValoriOrdinati[$i] = $Hash{$element};
        $i++;
    }
    
    return \@IndiciOrdinati, \@ValoriOrdinati;
    
}

# **********************************  descendHash  ***********************************

# Sort the keys of the hash according to their values
sub descendHash
{
    my $HashRef = $_[0];
    my %Hash = %$HashRef;
    
    if (!ref(%Hash) eq "HASH")
    {
        die "Error occured in function descendHash: subrutine argument is not valid\n";
    }
    
    my %HashOrdinato;
    my @IndiciOrdinati;
    my @ValoriOrdinati;
    my $i=0;
    
    foreach my $element(sort { $Hash{$b} <=> $Hash{$a}} keys %Hash)
    {
        $IndiciOrdinati[$i] = $element;
        $ValoriOrdinati[$i] = $Hash{$element};
        $i++;
    }
    
    return \@IndiciOrdinati, \@ValoriOrdinati;
    
}
# **********************************  NegazioneBit  ***********************************

# Nega bit: 1 -> 0 e 0 -> 1

sub NegazioneBit
{
    my $bit = $_[0];
    if($bit != 0 and $bit != 1)
    {   die "Error occurred in function NegazioneBit: input value is not valid: boolean data type is required  \n"}
    
    my $BitNegato;
    
    if($bit eq 0)
    {   $BitNegato= 1;  }
    else
    {   $BitNegato = 0; }
    
    return $BitNegato;
}

# **********************************  FindHessian  ***********************************
sub FindHessian
	{
		# POSSO ENTRARE NELLA CARTELLA resultDirectory -> NM_run1 e aprire il file psn.lst
		# Nel caso in cui fosse presente l'errore:
		# 0HESSIAN OF POSTERIOR DENSITY IS NON-POSITIVE-DEFINITE DURING SEARCH
		# Posso impostare il valore della fitness come valore molto basso: 99999999 ad esempio
		# in questo modo so che quel valore nel processo di ricerca verrà pescato con una probabilità molto bassa
	
        my $e = $_[0];
		my $ObjFunction = $_[1];
		my $FitnessValue = $_[2];
		my $TermMinimization ;
		 
        $CWD ="resultDirectory$e";
        $CWD="NM_run1";
        open my $FILE_psnlst ,'<' ,"psn.lst" , or die $!;
        my @psnlst = <$FILE_psnlst>;
        
        my $indice = 0;
        my $indexTerm;
    
        foreach my $riga_lst(@psnlst)
        {
            if($riga_lst =~ m/TERM:/)
			{
                $indexTerm = $indice;
				last;
			}
			$indice++;
		 }	
    
        $TermMinimization = $psnlst[$indexTerm+1];$TermMinimization=~s/\s+//;$TermMinimization=~s/\s+//;
        
        foreach my $riga(@psnlst)
        {
             if ($riga =~m/0HESSIAN OF POSTERIOR DENSITY IS NON-POSITIVE-DEFINITE DURING SEARCH/)
             {
                 $TermMinimization = $psnlst[$indexTerm-2];$TermMinimization=~s/\s+//;$TermMinimization=~s/\n+//;
                 $ObjFunction = 0;
                 $FitnessValue = 0;
             }
        }
        close $FILE_psnlst;
		
        $CWD="..";
		$CWD="..";
        
        return($ObjFunction, $FitnessValue, $TermMinimization);
	
}