#!/usr/bin/perl
use 5.014;
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;
use Time::HiRes qw( time );
use IPC::Shareable qw(:lock);

srand(123);

# CONSTANTS
my $BAD_FITNESS = -10000;
my $BAD_OBJFUN = -10000;

# PARAMETERS
my $timer = 10000 ; #initial timer value
my $prctile = 95;
my $popfile = '../popfile.csv'; # use '' when popfile.csv is missing for random initialization

# GA PARAMS
my $sizepop=20;
my $ElitismPercentage = 0.1;
my $NumOfGen = 100;
my $pc = 0.7;           # crossing over probability
my $pm = 0.025;         # mutation probability

# PARALLELIZATION PARAMS 
my $MAXCP=8; # NUMBER OF CPUs in which the GA will run

# SHAREABLES
my $glue = "glue".time(); #creation of an almost unique glue, in this way it is not possible that different instances of this script can share the variables (that would be a mess)
tie my @ObjFun, 'IPC::Shareable', { key => $glue.'ObjFun', create => 1, destroy => 1};
tie my @FitnessValue, 'IPC::Shareable', { key => $glue.'FitnessValue', create => 1, destroy => 1};
tie my @normalizedUPVs, 'IPC::Shareable', { key => $glue.'normalizedUPVs', create => 1, destroy => 1};
tie my @TermMinimization, 'IPC::Shareable', { key => $glue.'TermMinimization', create => 1, destroy => 1};
tie my @chrTime, 'IPC::Shareable', {key => $glue.'chrTime', create => 1, destroy => 1};
tie my $activecount, 'IPC::Shareable', {key => $glue.'activecount', create => 1, destroy => 1};
tie my @timed_out_flag, 'IPC::Shareable', {key => $glue.'timed_out_flag', create => 1, destroy => 1}; 
#****************************************************************************************
#                                    FILE TO WRITE                                      #
#****************************************************************************************

my @argsNew = ("touch", "SumUpFile.csv");
system(@argsNew); # CHECK DOESN'T WORK DUE TO IPC SHAREABLE == 0 or die "Error occurred: is impossible to create a new file for model informations\n";
open my $SUMUPFILE, '>', "SumUpFile.csv", or die $!;

my @argsNews = ("touch", "Worst_SumUpFile.csv");
system(@argsNews); # CHECK DOESN'T WORK DUE TO IPC SHAREABLE  == 0 or die "Error occurred: is impossible to create a new file for model informations\n";
open my $WSUMUPFILE, '>', "Worst_SumUpFile.csv", or die $!;

my @args = ("touch", "InformazioniFitness.csv");
system(@args); # CHECK DOESN'T WORK DUE TO IPC SHAREABLE  == 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("Generation, FitnessMassima,FitnessMimina,Fitnessmean,varianceFitness,DeviazioneStandard,Time\n");

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

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

my $MAXTheta = $maxTheta;

#****************************************************************************************
#                                 READ CONFIGURATION FILE                               #
#****************************************************************************************

open my $FILE ,'<' ,'ConfigurationFile.scm' , or die 'configuration file not found';
my @configFile = <$FILE>;
close FILE;

# Delete comments from configuration file
dropComment(\@configFile,'ConfigurationFileWithoutComments.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 covariate -> parameters                            #
#****************************************************************************************

my $refHash = hashCovariateParameter(\@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]);
    }
}

my $toprint = join(',', @to_print);
my $statusCromosom = 'created';


print $SUMUPFILE ("#Chromosom, #Generation, FitnessValue, ObjectiveFunction, TermOptimization, $toprint, StatusCromosom, Time, indexCorr, NormalizedUPVs\n");
print $WSUMUPFILE ("#Chromosom, #Generation, FitnessValue, ObjectiveFunction, TermOptimization, $toprint\n");

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

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

#****************************************************************************************
#                                      READ data.csv                                    #
#****************************************************************************************

open FILE ,'<' ,$FileDataName , or die 'file.csv not found';
my @wholeCsv=<FILE>;
close FILE;

my $names = $wholeCsv[0];
my @splitNames = split(/,|;|\s/, $names); # data can be separated by ; or , or tab

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

(my $indexContinuousRef , my $indexCategoricalRef, my $indexID) = subInfoDataCsv(\@splitNames,\@continuousVariables,\@categoricalVariables);
my %indexContinuous = %{$indexContinuousRef};
my %indexCategorical = %{$indexCategoricalRef};

#  ------------------------------- DATA MATRIX -------------------------------------
my @splitValues;
my @ID;
my %matrix;

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];   # first row: variable names
    }
    
    $ID[$i-1] =$matrix{$i}{$indexID};       # column ID
}


#  ---------------------------- meanNS of CONTINUOUS DATA  -------------------------------

(my $meannCVariablesRef, my $minRef, my $maxRef) = meannContinuousData(\@continuousVariables, \@ID,\%matrix, \%indexContinuous);
(my $sdCVariablesRef) = sdContinuousData(\@continuousVariables, \@ID,\%matrix, \%indexContinuous);
my %meannCVariables = %$meannCVariablesRef;   
my %minimumCVariables = %$minRef;                
my %maximumCVariables = %$maxRef;               
my %sdCVariables = %$sdCVariablesRef;               

#  --------------------------- CORRELATIONS of CONTINUOUS DATA  -------------------------------

my $CorrMatrixRef = correlazioneContinuousData(\@continuousVariables, \@ID,\%matrix, \%indexContinuous);
my %CorrMatrix = %{$CorrMatrixRef};

#  --------------------------- FREQUENCES of CATEGORICAL DATA -----------------------------

(my $FreqTotRef, my $FreqCovRef, my $refUnique)=FrequenceCategoricalData(\@categoricalVariables, \@ID, \%matrix, \%indexCategorical);
my %FreqTot = %{$FreqTotRef};
my %FreqCov = %{$FreqCovRef};
my %VarUnique = %{$refUnique};

#****************************************************************************************
#                           ESTIMATE OF MODEL WITHOUT COVARIATES                        #
#****************************************************************************************

print "\n******* [BASE MODEL ESTIMATION] *******\n";
my @arg1 = ("mkdir", "BaseModel");system(@arg1);
my @arg2 = ("cp", "MODEL.mod", "BaseModel");  system(@arg2);
my @arg3 = ("cp", "DataFile.csv", "BaseModel"); system(@arg3);
$CWD ="BaseModel";

@args = ("execute", "MODEL.mod", "-dir=resultDirectory", "-silent");
system(@args); # CHECK DOESN'T WORK DUE TO IPC SHAREABLE ==0 or die "Error occurred in function runNONMEM for MODEL.mod $_\n";
    
$CWD ="resultDirectory";
    
my $csvname = "raw_results_MODEL.csv";
open FILE ,'<' ,$csvname , or die "Error occurred: the file raw_results_MODEL.csv does not exists \n";
my @rawResultFile = <FILE>;
close FILE;

my $Resultsnames = $rawResultFile[0];
my @splitResultNames = split(/,\D/, $Resultsnames);
    
my $Resultvalues = $rawResultFile[1];
my @splitResultValues = split(/,/, $Resultvalues);
    
my $theta_index;
my $first_omega_index;
my $last_omega_index;

for (my $i=0; $i<=$#splitResultNames; $i++)
{
    $splitResultNames[$i] =~ s/[^a-zA-Z0-9]//g; 
    if ($splitResultNames[$i] eq 'ofv'){
        $theta_index = $i+1;
        next;
    }
    if ($splitResultNames[$i] eq 'OMEGA11'){
        $first_omega_index=$i;
        next;
    }
    if ($splitResultNames[$i] eq 'SIGMA11'){
        $last_omega_index=$i-1;
    }
}

#UPV
my $num_theta = (-1+sqrt(1+8*($last_omega_index-$first_omega_index+1)))/2; 
my @UPVbase;

for (my $i=0; $i<$num_theta; $i++){
    $UPVbase[$i] =  (sqrt(exp($splitResultValues[$first_omega_index+$i*$num_theta])-1)*$splitResultValues[$theta_index+$i])**2;
}
    
$CWD = '../..';

#****************************************************************************************
#                               DEFINITION OF THE GENOMA                                #
#****************************************************************************************

my $numberofStates = $#ValidStatesContinuous+1;
my $numberofBits = CEIL(log2($numberofStates));
my @AmGeni = genome($numberofBits);
my @AmGeniCat = ('0', '1');
my @WholeGenes = (@AmGeni, @AmGeniCat);
my $nAmGeni = $#AmGeni+1;


#****************************************************************************************
#                            INITIAL POPULATION - Generation 0                          #
#****************************************************************************************

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

my %FitnessHash;
my %population;

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

my $e; # Counts the individuals
my $nContinuousRelations;
my $NumChromElitism = int($ElitismPercentage*$sizepop);
my $TermMinimization;
my $infoElitism;
my $FitnessValue;

tie my @Chromosom, 'IPC::Shareable', {key => $glue.'Chromosom', create => 1, destroy => 1}; 
my @statusCromosom = ();
my @indexCorr = ();
my @vecttoprintComplete;
my @ObjFunVect = ();
my @FitnessVect = ();
my @TermMinimizationVect =();
my $toprintVect;

my $ObjFun;

@arg1 = ("mkdir", "Generation0");system(@arg1);
@arg2 = ("cp", "MODEL.mod", "Generation0");  system(@arg2);
@arg3 = ("cp", "DataFile.csv", "Generation0"); system(@arg3);
$CWD ="Generation0";

my $G=0;

#****************************************************************************************
#                        0. IMPORTING CHROMOSOMES FROME POPFILE                         #
#****************************************************************************************

my @chr=();

if ($popfile ne ''){
    open FILE ,'<' ,$popfile , or die "Error occurred: the initial population file $popfile does not exist \n";
    while (my $line = <FILE>)
    {
        chomp($line);
        $line =~ s/,//g;
        @chr=(@chr, $line);
    };
    close FILE;
}
#**************************************************************************************** 

    
my @maxTheta=();
my $begin_time = time();
my $usefulTime = 0;
my $usefulChromosomes = 0;
$activecount = 0;
my @usefulTimes = ();

for($e=0; $e<$sizepop; $e++)
{
    my %relationstart;
    my %paramCovList;
	my @gene = (); 				
    my @Vect_to_print = ();
    @Chromosom = ();
    my $rawchr = shift @chr;
	$maxTheta[$e] = $MAXTheta;
    
    my @args = ("touch", "Chromosom$e.mod");
    system(@args); # CHECK DOESN'T WORK DUE TO IPC SHAREABLE  == 0 or die "system @args failed: $?";
    
    open my $FILE ,'>' ,"Chromosom$e.mod" , or die $!;
    
    $line[$indexPROBLEM] = substr($line[$indexPROBLEM], index($line[$indexPROBLEM], '$PROBLEM'));
    print $FILE("@line[$indexPROBLEM..$indexPK]\n");
    

    #****************************************************************************************
    #                             1. WRITE NEW NONMEM FILE                                  #
    #****************************************************************************************

    (my $RelationStartRef, my $ChromosomRef, $maxTheta[$e], $nContinuousRelations, my $refhash, my $paramCovListRef) =  write_NONMEMfile(\@param, \@AmGeni, \@AmGeniCat, \%relationstart, $maxTheta[$e], \%meannCVariables, $FILE, \@ValidStatesContinuous, \@Chromosom, \%FreqCov, \%FreqTot, \@gene, $nAmGeni, $indexPK, $indexOmega, \@line, $G, \@Covariate,\%hash, \@continuousVariables, \@categoricalVariables, \%VarUnique, \%minimumCVariables, \%maximumCVariables, $popfile, $rawchr,\%sdCVariables);

    %relationstart = %$RelationStartRef;
    @Chromosom = @$ChromosomRef;
    for(my $K=0; $K<=$#Chromosom;$K++)
    {
        $population{$e}[$K] = $Chromosom[$K];
    }
	%paramCovList = %{$paramCovListRef};
    my %FenotypeHash = %{$refhash};
    
    for (my $b=0; $b<= $#to_print; $b++)
    {
        @Vect_to_print = (@Vect_to_print, $FenotypeHash{$to_print[$b]});
    }
    my $vecttoprint = join(',', @Vect_to_print);
    $vecttoprintComplete[$e] = join(',', @Vect_to_print);
    
    close $FILE;
    
    $indexCorr[$e] = compute_correlation_index(\%paramCovList, \%CorrMatrix);
    
    #****************************************************************************************
    #                                  2. RUN NONMEM                                        #
    #****************************************************************************************
    
    #FORK
    my $pid = fork;
    die "failed to fork: $!" unless defined $pid;
    #MAIN PROCESS
    if ($pid){
        tied($activecount)->lock;
        $activecount++;
        tied($activecount)->unlock;
        while ($activecount>=$MAXCP){
                sleep 0.5;
        }
        next;
    }; #the parent goes on forking in the next iteration of the loop
    
    #CHILD PROCESS
    print("BEGIN Chromosome $e of Generation $G\n");
    my $begin_c_time = time();
    tie my @timed_out_flag, 'IPC::Shareable', {key => $glue.'timed_out_flag'}; 
    $timed_out_flag[$e] = runNONMEM($e, $timer);
    
    tie my @chrTime, 'IPC::Shareable', {key => $glue.'chrTime'}; 
    tied(@chrTime)->lock;
    $chrTime[$e] = time() - $begin_c_time;
    tied(@chrTime)->unlock;
    
    tie my $activecount, 'IPC::Shareable', {key => $glue.'activecount'}; 
    tied($activecount)->lock;
    $activecount--;
    tied($activecount)->unlock;
	print("E N D Chromosome $e of Generation $G\t($chrTime[$e] s)\n");
    exit; #CHILD exits the loop
    
} #close parallelized for sizePop

#PARENT WAITS FOR CHILDREN 
my $kid;
do {
$kid = waitpid -1, 0;
} while($kid > 0);

for(my $e=0; $e<$sizepop; $e++)
{

    #****************************************************************************************
    #                     3.  FIND FITNESS in ResultDirectory                               #
    #****************************************************************************************

    my $result, my $ObjFunction, my $normalizedUPVs;
    if ($timed_out_flag[$e] == 0){
	push @usefulTimes, $chrTime[$e];
        ($result, $ObjFunction, $normalizedUPVs) = findFitness($e, $maxTheta[$e], $MAXTheta, \%FitnessHash, \@UPVbase,$indexCorr[$e]);
        $FitnessHash{$e} = $result;
        $FitnessValue = $result;
        $usefulChromosomes++;
        $usefulTime+=$chrTime[$e];
        # Find if Hessian
       ($ObjFunction, $FitnessValue, $TermMinimization) = FindHessian($e, $ObjFunction, $FitnessValue);
    }else{
        $result = $BAD_FITNESS;
        $FitnessHash{$e} = $result;
        $FitnessValue = $result;
        $ObjFunction = $BAD_OBJFUN;
        $normalizedUPVs = '';
        $TermMinimization = 'TIMED OUT';
    }
    
    print $SUMUPFILE ("$e, $G, $FitnessValue, $ObjFunction, $TermMinimization,  $vecttoprintComplete[$e], $statusCromosom, $chrTime[$e], $indexCorr[$e], $normalizedUPVs\n");
    
    @TermMinimizationVect = (@TermMinimizationVect, $TermMinimization);
    @FitnessVect = (@FitnessVect, $FitnessValue);
    @ObjFunVect = (@ObjFunVect, $ObjFunction);
    

} #close for sizePop

@usefulTimes = sort { $a <=> $b } @usefulTimes;
$timer = $usefulTimes[ceil($prctile/100*$#usefulTimes)];
print "\n Timer: --> $timer \n";

$CWD ="..";

my @nGenesVect= @{$population{0}};
my $nGene = $#nGenesVect;

#****************************************************************************************
#                     Select best chromosomes (higher fitness)                          #
#****************************************************************************************

my %FitnessOrder;

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

my @Index = keys %FitnessOrder;
my @BestIndex = @Index[0..2]; 

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

my @Fitnesssss = values %FitnessHash;
my @TrueFitness = ();
for (my $i=0; $i<=$#FitnessVect; $i++)
{
	
	if($FitnessVect[$i] != 0)
	{
		@TrueFitness = (@TrueFitness, $FitnessVect[$i]);
	}
}
my $variance1 = computevariance(@TrueFitness);
my $minFit = min(@TrueFitness);
my $maxFit = max(@TrueFitness);
my $meanFit = computemean(@TrueFitness);
my $devFit = sqrt($variance1);

my $genTime = time() - $begin_time;

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

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

my @crom;
my %countResults;
my %count;

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

#****************************************************************************************
#                         MOST FREQUENT GENES IN EACH POSITION                          #
#****************************************************************************************

my @geneMax = GenesMostFrequent_x_Position(\@WholeGenes, \%count);

#****************************************************************************************
#                               Selection + Crossing Over                               #
#****************************************************************************************

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

my @Table= @vecttoprintComplete; 
my @ObjTable = @ObjFunVect;
my @FitnessTab = @FitnessVect;
my @TermMinimizationTab = @TermMinimizationVect;

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

my $countsOffspring;

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# 					              		FURTHER GENERATIONS											 %
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for (my $G=1; $G<=$NumOfGen; $G++)
{
    my $begin_time = time();
    #****************************************************************************************
    #                                 CREO - COPIO FILE                                     #
    #****************************************************************************************

    print("\n\n*********** [Generation $G] *********** \n");
    
    my @arg1 = ("mkdir", "Generation$G");              system(@arg1);
    my @arg2 = ("cp", "MODEL.mod", "Generation$G");  system(@arg2);
    my @arg3 = ("cp", "DataFile.csv", "Generation$G"); system(@arg3);
    $CWD ="Generation$G";
    
    $countsOffspring = 0;
    @Vect_to_printBest = ();
    
    #****************************************************************************************
    #                                 ELITISM - STEP 1                                     #
    #****************************************************************************************

    (my $OrderedIndxBestRef,my $OrderedFitBestRef) = descendHash(\%FitnessHash);
    my @OrderedIndxBest = @$OrderedIndxBestRef;
    my @OrderedFitBest = @$OrderedFitBestRef;
    
    my @BestIdxElitism = @OrderedIndxBest[0..$NumChromElitism];
    
    my @FitnessBestElitismo = @OrderedFitBest[0..$NumChromElitism];
    
    for (my $obix =0; $obix<=$#BestIdxElitism; $obix++)
    {
        $ObjectiveBest[$obix] = $ObjFunVect[$BestIdxElitism[$obix]];
        $FitBest[$obix] = $FitnessVect[$BestIdxElitism[$obix]];
		$TermBest[$obix] = $TermMinimizationVect[$BestIdxElitism[$obix]];
        @Vect_to_printBest = (@Vect_to_printBest, $vecttoprintComplete[$BestIdxElitism[$obix]]);
    }
    

    for(my $i=0; $i< $sizepop/2; $i++)
    {
        #******************************************************************
        #          PARENT SELECTION - Tournament selection                #
        #******************************************************************

        my $ParentCromosomRef1 = TournamentSelection($sizepop, \%FitnessHash, \%population);
        my @ParentCromosom1 = @{$ParentCromosomRef1};
        
        #my $ParentCromosomRef2 = RouletteWheel(\%FitnessHash, \%population);
        my $ParentCromosomRef2 = TournamentSelection($sizepop, \%FitnessHash, \%population);
        my @ParentCromosom2 = @{$ParentCromosomRef2};
	
        #******************************************************************
        #         CROSSING OVER - OnePoint Crossover Recombination        #
        #******************************************************************

        if(rand(1)<= $pc)
        {
            (my $offspring1Ref, my $offspring2Ref) = onePoinCrossOver(\@ParentCromosom1, \@ParentCromosom2);
            $Newpopulation{"$countsOffspring"} = join('',@$offspring1Ref);
            $countsOffspring++;
            @offspring2 = @$offspring2Ref;
            $Newpopulation{"$countsOffspring"} = join('', @$offspring2Ref);
            $countsOffspring++;
        }
        else
        {
            $Newpopulation{"$countsOffspring"} = join('', @ParentCromosom1);
            $countsOffspring++;
            $Newpopulation{"$countsOffspring"} = join('', @ParentCromosom2);
            $countsOffspring++;
        }
        
    }
    my @MutatedIndividual;

    #******************************************************************
    #                            MUTATION                            #
    #******************************************************************

    for(my $i=0; $i<$sizepop; $i++)
    {
        my @individual = split('', $Newpopulation{"$i"});
        
        for(my $j=0; $j<=$#individual; $j++)
        {
            my $randomicNumber = rand(1);
            
            if($randomicNumber <= $pm)
            {
                $MutatedIndividual[$j] = NegazioneBit($individual[$j]);
            }
            else
            {
                $MutatedIndividual[$j] = $individual[$j];
            }
        }
        
        $Newpopulation{"$i"} = join('', @MutatedIndividual);
    }
    
    my $NewPopRef = splitNewpop($nContinuousRelations,$countsOffspring,$nGene,\%Newpopulation, \@geneMax, \@AmGeni, $numberofBits);
    %NewPop = %{$NewPopRef};
    
    @ObjFunVect = ();
    @FitnessVect = ();
	@TermMinimizationVect = ();
	@vecttoprintComplete = ();
	my @indexCorrVect = ();
	my @normalizedUPVsVect = ();
    my @statusCromosoms = ();
	#my @chrTime = ();
	my @indexCorr = ();  
    for(my $e=0; $e<$sizepop; $e++)
    {
        
        my @gene=();
        @Chromosom = ();
        my @Vect_to_print= ();
        my $nAmGeni = $#AmGeni+1;
        my %relationstart;
        my $b;
        my @ret;
        
        my %paramCovList;
        
        my @args = ("touch", "Chromosom$e.mod");
        system(@args); # CHECK DOESN'T WORK DUE TO IPC SHAREABLE  == 0 or die "system @args failed: $?";
        
        open my $FILE ,'>' ,"Chromosom$e.mod" , or die $!;
        print $FILE("@line[$indexPROBLEM..$indexPK]\n");
        
        
        #****************************************************************************************
        #                             1. WRITE NEW NONMEM FILE                                  #
        #****************************************************************************************
 
        my @ChromosomIN = @{$NewPop{"$e"}};
        
       (my $RelationStartRef, my $ChromosomRef, $maxTheta[$e], my $nnn, my $refhash, my $paramCovListRef) =  write_NONMEMfile(\@param, \@AmGeni, \@AmGeniCat, \%relationstart, $MAXTheta, \%meannCVariables, $FILE, \@ValidStatesContinuous, \@ChromosomIN, \%FreqCov, \%FreqTot, \@gene, $nAmGeni, $indexPK, $indexOmega, \@line, $G,\@Covariate,\%hash, \@continuousVariables, \@categoricalVariables, \%VarUnique, \%minimumCVariables, \%maximumCVariables, $popfile, \@chr,\%sdCVariables);
        %relationstart = %$RelationStartRef;
        %paramCovList = %{$paramCovListRef};
        my %FenotypeHash = %{$refhash};
 
        for (my $b=0; $b<= $#to_print; $b++)
        {
            @Vect_to_print = (@Vect_to_print, $FenotypeHash{$to_print[$b]});
        }
        $vecttoprintComplete[$e] = join(',', @Vect_to_print);
        @Chromosom = @$ChromosomRef;
        
        close $FILE;
        
        $indexCorr[$e] = compute_correlation_index(\%paramCovList, \%CorrMatrix);
        #****************************************************************************************
        #                                  2. RUN NONMEM                                        #
        #****************************************************************************************
        
        $chrTime[$e]=0; #duplicated and elited chromosomes don't need estimation
        
        $normalizedUPVs[$e] = '';
        if( any {$_ ~~ $vecttoprintComplete[$e]}@Table)
        {
			$statusCromosoms[$e] = 'duplicated';
            
			my $IndiceTrue;
        
            for(my $indiceTable=0; $indiceTable<=$#Table; $indiceTable++)
            {
                if($vecttoprintComplete[$e] eq $Table[$indiceTable])
                {
                    $IndiceTrue = $indiceTable;
                    last;
                }
            }
            
            $FitnessValue[$e] = $FitnessTab[$IndiceTrue];
            $ObjFun[$e] = $ObjTable[$IndiceTrue];
            $vecttoprintComplete[$e] = $Table[$IndiceTrue];
			#$toprintVect = $Table[$IndiceTrue]; # fenotipo
			$TermMinimization[$e] = $TermMinimizationTab[$IndiceTrue];
			
			$FitnessHash{$e} = $FitnessTab[$IndiceTrue];
        }
        else
        {
				$statusCromosoms[$e] = 'created';
			
			    
                #FORK
                my $pid = fork;
                die "failed to fork: $!" unless defined $pid;
                #MAIN PROCESS
                if ($pid){
                    tied($activecount)->lock;
                    $activecount++;
                    tied($activecount)->unlock;
                    while ($activecount>=$MAXCP){
                            sleep 0.5;
                    }
                    next;
                }; #the parent goes on forking in the next iteration of the loop
                
                #CHILD PROCESS
                
                #SHAREABLES
                tie my @ObjFun, 'IPC::Shareable', { key => $glue.'ObjFun'};
                tie my @FitnessValue, 'IPC::Shareable', { key => $glue.'FitnessValue'};
                tie my @normalizedUPVs, 'IPC::Shareable', { key => $glue.'normalizedUPVs'};
                tie my @TermMinimization, 'IPC::Shareable', { key => $glue.'TermMinimization'};
                #tie my %FitnessHash, 'IPC::Shareable', {key => $glue.'FitnessHash'};
                tie my @chrTime, 'IPC::Shareable', {key => $glue.'chrTime'};
                tie my $activecount, 'IPC::Shareable', {key => $glue.'activecount'};
                
                print("BEGIN Chromosome $e of Generation $G\n");
                my $begin_c_time = time();
                
                my $timed_out_flag= runNONMEM($e, $timer);
              
                tied(@chrTime)->lock;
                $chrTime[$e] = time() - $begin_c_time;
                tied(@chrTime)->unlock;

                tied($activecount)->lock;
                $activecount--;
                tied($activecount)->unlock;  
			
		#****************************************************************************************
            #                     3.  FIND FITNESS in ResultDirectory                               #
            #****************************************************************************************
                my $result;
                if ($timed_out_flag == 0){

                    ($result, $ObjFun[$e], $normalizedUPVs[$e]) = findFitness($e, $maxTheta[$e], $MAXTheta, \%FitnessHash, \@UPVbase, $indexCorr[$e]);
                    
                    $FitnessValue[$e] = $result;
			        ($ObjFun[$e], $FitnessValue[$e], $TermMinimization[$e]) = FindHessian($e, $ObjFun[$e], $FitnessValue[$e]);
			    }else{
                    $FitnessValue[$e] = $BAD_FITNESS;
                    $ObjFun[$e] = $BAD_OBJFUN;
                    $normalizedUPVs[$e] = '';
                    $TermMinimization[$e] = 'TIMED OUT';
                }
                
			    print("E N D Chromosome $e of Generation $G\t($chrTime[$e] s)\n");
			    exit; #CHILD exits the loop
		}
		
        
    } #close for sizePop (over $e)
    
    #PARENT WAITS FOR CHILDREN 
    my $kid;
    do {
    $kid = waitpid -1, 0;
    } while($kid > 0);
    
    for(my $e=0; $e<$sizepop; $e++)
    {
        $FitnessHash{$e}=$FitnessValue[$e];
        @ObjFunVect = (@ObjFunVect, $ObjFun[$e]);
        @FitnessVect = (@FitnessVect, $FitnessValue[$e]);
	    @TermMinimizationVect = (@TermMinimizationVect, $TermMinimization[$e]);
		@normalizedUPVsVect = (@normalizedUPVsVect, $normalizedUPVs[$e]);
		@indexCorrVect = (@indexCorrVect, $indexCorr[$e]);
		if ($timed_out_flag[$e] == 0){
		    push @usefulTimes, $chrTime[$e];
		}
	}
    
   @usefulTimes = sort { $a <=> $b } @usefulTimes;
   $timer = $usefulTimes[$prctile/100*$#usefulTimes];
   print "\n Timer: --> $timer \n";
    
    #****************************************************************************************
    #                                  ELITISM - STEP 2                                    #
    #****************************************************************************************

    # Sort in acending order for ELITISM
    (my $refIndex, my $refValues)  = ascendHash(\%FitnessHash);
    my @OrderedIndexWorst = @$refIndex;
    my @FitnessOrderedWorst = @$refValues;
    
    # Migliori 10% degli indici
    my @WorstIndexElitism = @OrderedIndexWorst[0..$NumChromElitism];
    
    # Miglior 10% delle fitness
    my @FitnessWorstElitism = @FitnessOrderedWorst[0..$NumChromElitism];
    
    #****************************************************************************************
    #                                  ELITISM - STEP 3                                    #
    #****************************************************************************************
    
    for(my $l=0; $l <= $#BestIdxElitism; $l++)
    {
		
		$statusCromosoms[$WorstIndexElitism[$l]] = 'elited';
		
		print $WSUMUPFILE("$WorstIndexElitism[$l], $G, $FitnessVect[$WorstIndexElitism[$l]],$ObjFunVect[$WorstIndexElitism[$l]],$TermMinimizationVect[$WorstIndexElitism[$l]],$vecttoprintComplete[$WorstIndexElitism[$l]]\n");
        $NewPop{$WorstIndexElitism[$l]} = $population{$BestIdxElitism[$l]};
        $FitnessHash{$WorstIndexElitism[$l]} = $FitnessBestElitismo[$l];
        
        $vecttoprintComplete[$WorstIndexElitism[$l]] = $Vect_to_printBest[$l];
		$ObjFunVect[$WorstIndexElitism[$l]] = $ObjectiveBest[$l]; 
        $FitnessVect[$WorstIndexElitism[$l]] = $FitBest[$l];
		$TermMinimizationVect[$WorstIndexElitism[$l]] = $TermBest[$l];
    }
    
    @ObjTable = (@ObjTable, @ObjFunVect);
    @FitnessTab = (@FitnessTab, @FitnessVect);
	@TermMinimizationTab = (@TermMinimizationTab, @TermMinimizationVect);
	@Table = (@Table, @vecttoprintComplete);
	
    # population finale
    %population = %NewPop;
    $CWD ="..";
    
	my @TrueFitness = ();
	for (my $iq=0; $iq<=$#FitnessVect; $iq++)
	{
		if($FitnessVect[$iq] != 0)
		{
			@TrueFitness = (@TrueFitness, $FitnessVect[$iq]);
		}
	}	
    my $MaxValueFit = max(@TrueFitness);
    my $MinValueFit = min(@TrueFitness);
    my $MeanValueFit = computemean(@TrueFitness);
    my $VarianceValueFit = computevariance(@TrueFitness);
	my $DeviazioneFit = sqrt($VarianceValueFit);
    my $genTime = time() - $begin_time; 
    print $FILEFITNESS("$G, $MaxValueFit,$MinValueFit,$MeanValueFit,$VarianceValueFit, $DeviazioneFit, $genTime\n");
    
	for(my $individualE=0; $individualE<$sizepop; $individualE++)
    {
        print $SUMUPFILE ("$individualE, $G, $FitnessVect[$individualE], $ObjFunVect[$individualE] ,$TermMinimizationVect[$individualE] ,$vecttoprintComplete[$individualE], $statusCromosoms[$individualE], $chrTime[$individualE], $indexCorrVect[$individualE], $normalizedUPVsVect[$individualE]\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 $row(@FileNonmem)
    {
        $line[$k] = $row;
        $line[$k] =~ s/\\//;
        
        if($row  =~ m/\$PROBLEM\s*/) {
            $indexPROBLEM = $k;
        }
        
        if($row  =~ m/\$DATA\s*/) {
            $indexDATA = $k;
        }
        
        if($row  =~ m/\$PK\s*/) {
            $indexPK = $k;
        }
        
        if($row  =~ m/\$OMEGA\s*/)
        {
            $indexOmega = $k;
        }
        
        if($row  =~ m/\$THETA\s*/)
        {
            $indexTheta = $k;
        }
        
        $k++;
        
        if($row =~ m/TV.+THETA/)
        {
            $flag = 1;
            
            $param[$i] = trim(substr $row, 2, (index($row, '=')-2));
            chomp($param[$i]);
            $param[$i]=~s/TV//;
            
            $i++;
        }
        
        my $position = 0;
        my $info;
        
        if($row =~ m/.+THETA/)
        {
            for (my $indice = 0; $indice< length($row); $indice++)
            {
                $info = index($row, 'THETA(', $position+1);
                $position = $info;
                if ($position != -1)
                {
                    $tetas[$j] = substr $row, $position+6, 1;
                    $j++;
                }
            }
        }
    }
    
    my $maxTheta = max(@tetas);
    
    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  ***********************************

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 $row(@FileVect)
    {
        my $toCopy = substr($row, 0, index($row, ';'));
        my $toDrop = substr($row, index($row, ';'));
        $toDrop = "\n";
        $FileVect[$q] = $toCopy.$toDrop;
        print $FILE $FileVect[$q];
        $q++;
    }
    
    close $FILE;
}


# ********************************   subConfigFile ********************************

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 $row(@configFile)
        {
            if($row =~ m/\[$search[$j]\]/)
            {
                $ix{$search[$j]} = $i;
            }
            
            if($row=~ m/categorical_covariates/)
            {
                $ct = substr $row, (index($row, '=')+1), -1;
                $ct =~ s/\\//;
                $flagCategorical = 0;
            }
            
            if($row=~ m/continuous_covariates/)
            {
                $cc = substr $row, (index($row, '=')+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);
}



# ********************************* hashCovariateParameter  ***********************************
sub hashCovariateParameter
{
    # 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;
    
    for(my $i=0; $i<= $#continuousVariables; $i++)
    {
        my $k = 0;
        my $j = 0;
        
        my $covariata = trim($continuousVariables[$i]);
        
        foreach my $row(@configFile)
        {
            if($j> $v[$indTR] && $j < $v[$indVS] && $row =~ m/\w\s*\=\w\W*/)
            {
                
                my $parametro = substr($row, 0, index($row, '='));
                my $covariate = trim(substr($row, index($row, '=')+1, -1));
                $covariate =~ s/\\//;
                
                my @covariateVect = split(/,/ , $covariate);
                
                if(any {$_ eq $covariata} @covariateVect)
                {
                    $hash{"$covariata"}[$k] = $parametro;
                    $k++;
                }
            }
            $j++;
        }
    }
    
    for(my $i=0; $i<= $#categoricalVariables; $i++)
    {
        my $j = 0;
        my $k = 0;
        my $covariata = $categoricalVariables[$i];
        
        foreach my $row(@configFile)
        {
            if($j> $v[$indTR] && $j < $v[$indVS] && $row =~ m/\w\s*\=\w\W*/)
            {
                my $parametro = substr($row, 0, index($row, '='));
                my $covariate = trim(substr($row, index($row, '=')+1, -1));
                $covariate =~ s/\\//;
                
                my @covariateVect = split(/,/ , $covariate);

                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 $row(@configFile)
    {
        if($j > $v[$indVS] && $row =~ m/continuous\=(\d+\W*)/)
        {
            my $relazioni = substr($row, index($row, '=')+1, -1);
            $relazioni =~ s/\\//;
            @ValidStatesContinuous = split(/,/,$relazioni);
        }
        $j++;
    }
    for (my $element=0; $element<=$#ValidStatesContinuous; $element++)
    {
        $ValidStatesContinuous[$element]= trim($ValidStatesContinuous[$element]);
    }
    return \@ValidStatesContinuous;
}


# **********************************  subInfoDataCsv **********************************
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++)
    {
        $splitNames[$i]=~ s/\W//;
        $splitNames[$i]=trim($splitNames[$i]);
        
        for(my $j=0; $j<= $#continuousVariables; $j++)
        {
            my $a = ($continuousVariables[$j]);
            $a=trim($a);
            
            if($splitNames[$i] eq $a)
            {
                $indexContinuous{$a}= $i;
            }
        }
        
        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);
}


# ******************************** correlazioneContinuousData ********************************
sub correlazioneContinuousData
{
    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 %CorrMatrix = ();
    
    my @ValoreCovariata;
    my @ValoreCovariata2; 
    for(my $j=0; $j<=$#continuousVariables; $j++){
        my $a = $continuousVariables[$j];$a= trim($a);
        for(my $k=$j+1; $k<=$#continuousVariables; $k++){
        
            @ValoreCovariata= ();
            @ValoreCovariata2= ();
            
            my $b = $continuousVariables[$k];$b= trim($b);
            
            for(my $i=1; $i<= $#ID; $i++)
            {
                    $ValoreCovariata[$i-1] = $matrix{$i}{$indexContinuous{$a}};
                    $ValoreCovariata2[$i-1] = $matrix{$i}{$indexContinuous{$b}};
                    
            }
            $CorrMatrix{$a}{$b} = correlation(vector(@ValoreCovariata),vector(@ValoreCovariata2));
            $CorrMatrix{$b}{$a} = correlation(vector(@ValoreCovariata),vector(@ValoreCovariata2));
        }
        #print "\n";
    }
    
    return (\%CorrMatrix);
}

# ******************************** meannContinuousData ********************************
sub meannContinuousData
{
    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 %meannCVariables;  # ---> hash: covariatacontinua -> valore meanno
    my %minimumCVariables;  # ---> hash: covariatacontinua -> valore minimo
    my %maximumCVariables; # ---> 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++;
            }
        }
        $meannCVariables{$a} = computemean(@ValoreCovariata);
        $minimumCVariables{$a} = min(@ValoreCovariata);
        $maximumCVariables{$a} = max(@ValoreCovariata);
    }

    return (\%meannCVariables, \%minimumCVariables, \%maximumCVariables);
}

# ******************************** sdContinuousData ********************************
sub sdContinuousData
{
    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 %sdCVariables;  # ---> hash: covariatacontinua -> sd

    
    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++;
            }
        }
        $sdCVariables{$a} = stddev(@ValoreCovariata);
    }

    return (\%sdCVariables);
}

# ******************************** FrequenceCategoricalData ********************************
sub FrequenceCategoricalData
{
    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++){
        
        @ColonnaCategorical= ();
        my $a = $categoricalVariables[$j];$a=~ s/\n//;$a=~ s/\s+//;
        my $o=0;
        
        for(my $i=1; $i< $#ID; $i++)
        {
            if($ID[$i-1] != $ID[$i])
            {
                $ColonnaCategorical[$o] = $matrix{$i}{$indexCategorical{$a}};
                $o++;
                
            }
        }
        @varfjoin = grep {$_ =~ /[\d*]/} @ColonnaCategorical;
        my $varfjoin = join("", @varfjoin);trim($varfjoin);
        
        my @unico =(uniq(@ColonnaCategorical));
        $ValUnique{$categoricalVariables[$j]} = $#unico;
        
        
        for(my $q=0; $q<= $#unico; $q++)
        {
            $FreqTot{$a}[$q] = $unico[$q];
            
            my $count0=0;
            
            for(my $s=0; $s<= $#ColonnaCategorical; $s++)
            {
                
                if($ColonnaCategorical[$s] == $unico[$q])
                {
                    $count0++;         
                }
            }
            $counts[$q] = $count0;    
        }
        
        my $m = max(@counts);          
        
        for(my $q=0; $q<= $#unico; $q++)
        {
            if($counts[$q] == $m)
            {
                $indexMax = $q;                    
            }
        }
        $FreqCov{$a} = $unico[$indexMax];
    }
    return(\%FreqTot, \%FreqCov, \%ValUnique);
}


# ******************************** DS_continuous ********************************

sub DS_continuous
{
    my $gene = $_[0];
    my $covariata = $_[1];
    my $maxTheta = $_[2];
    my $meanna = $_[3]; $meanna =~ s/,/./;
    my $param = $_[4];
    my $FILE = $_[5];
    my $relationsRef = $_[6]; my @relations = @$relationsRef;
    my $AmGeneRef = $_[7]; my @AmGene =  @$AmGeneRef;
    my $paramCovListRef = $_[8]; my %paramCovList = %{$paramCovListRef};
    my $sd = $_[9]; $sd =~ s/,/./;
    my %relazioni;
    my $VectDS_continuous;
    my $r = tie %relazioni, 'Tie::IxHash';
    
    for(my $i=0; $i<=$#relations; $i++)
    {
        $relations[$i] = trim($relations[$i]);
        $relazioni{$AmGene[$i]} =  $relations[$i];
    }

    given($gene)
    {
        
        when($relazioni{$gene} ne 1)
        {
            push @{$paramCovList{$param}}, $covariata;
        }
    }
    
    given($gene)
    {      
        
        #    **** 2: LINEAR RELATIONSHIP ****
        when($relazioni{$gene} eq 2)
        {
            
            $maxTheta = $maxTheta +1;
            print $FILE(";;; $param$covariata-DEFINITION START\n$param$covariata = (1+THETA($maxTheta)*(($covariata-$meanna)/$sd))\n;;;$param$covariata-DEFINITION END\n\n");
            $VectDS_continuous = "$param$covariata";
            return ($maxTheta, \%relazioni, $VectDS_continuous,  \%paramCovList);
        }
        
        #    **** 3: PIECE-WISE LINEAR RELATIONSHIP  ****
        when($relazioni{$gene} eq 3)
        {
            
            $maxTheta = $maxTheta +1;
            
            print $FILE(";;; $param$covariata-DEFINITION START\nIF($covariata.LE.$meanna) $param$covariata=(1+THETA($maxTheta)*(($covariata-$meanna)/$sd))\n");
            $maxTheta = $maxTheta +1;
            print $FILE("IF($covariata.GT.$meanna) $param$covariata=(1+THETA($maxTheta)*(($covariata-$meanna)/$sd))\n;;;$param$covariata-DEFINITION END\n\n");
            $VectDS_continuous = "$param$covariata";
            return ($maxTheta, \%relazioni, $VectDS_continuous, \%paramCovList);
        }
        
        
        #     **** 4: EXPONENTIAL RELATIONSHIP ****
        when($relazioni{$gene} eq 4)
        {         
            $maxTheta = $maxTheta +1;
            print $FILE(";;; $param$covariata-DEFINITION START\n$param$covariata =EXP(THETA($maxTheta)*(($covariata - $meanna)/$sd))\n;;;$param$covariata-DEFINITION END\n\n");
            $VectDS_continuous = "$param$covariata";
            return ($maxTheta, \%relazioni, $VectDS_continuous, \%paramCovList);
        }
        
        
        #     **** 5: POWER RELATIONSHIP ****
        when($relazioni{$gene} eq 5)
        {
            
            $maxTheta = $maxTheta +1;
            print $FILE(";;; $param$covariata-DEFINITION START\n$param$covariata = (($covariata/$meanna)**(THETA($maxTheta)/$sd))\n;;;$param$covariata-DEFINITION END\n\n");
            $VectDS_continuous = "$param$covariata";
            return ($maxTheta, \%relazioni, $VectDS_continuous, \%paramCovList);
        }
    }
}

# ********************************* DS_categorical **********************************

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}};
    my $VectDS_categorical;
    
    #    **** 1: LINEAR RELATIONSHIP ****
    
    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})
        {
            $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: NO RELATIONSHIP ****
    
    return($maxTheta, $VectDS_categorical);
}

# ****************************************** genome **************************************
sub genome
{
    my $numberofBits = $_[0];
    
    my @AmGeni;
    
    given($numberofBits)
    {
        when($numberofBits==1)
        {
            @AmGeni = ('0', '1');
        }
        when($numberofBits==2)
        {
            @AmGeni = ('00', '01', '10', '11');
        }
        when($numberofBits==3)
        {
            @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 %meannCVariables = %$ref6;
    my $FILE = $_[6];
    my $ref7 =$_[7]; my @ValidStatesContinuous = @$ref7;
    my $ref8 = $_[8]; my @ChromosomIN = @$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 %minimumCVariables = %$ref18;
    my $ref19 = $_[23]; my %maximumCVariables = %$ref19;
    my $popfile = $_[24];
    my $rawchr = $_[25];
    my $ref20 = $_[26]; my %sdCVariables = %$ref20;
    
    my $nContinuousRelations;
    my $k=0;
    my $nelementi;
    
    my %ParamCov_Gene;
    my %relazioni;
    my @VectDSContinuous;
    my $RelDSContinuous;
    my %paramCov_Relazione;
    
    my %paramCovList;
    
    # Inizializza Chromosom
    my @Chromosom =();
    my $nAmGeniCat = $#AmGeniCat+1;
    
    my $pos = 0;
    my $nbits = length($AmGeni[0]);
    #****************************************************************************************
    #                         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;
            my $paramCovListRef;

            $nelementi = $#valore1+1;
            
            #@VectDSContinuous =();
            for(my $i=0; $i<$nelementi; $i++) # iterazioni sui parametri
            {
                my $parametro = $hash{$covariata}[$i];
                
                #****************************************************************************************
                #                                    find gene                                          #
                #****************************************************************************************
                if ($G==0)
                {
                    if ($popfile ne ''){ 
                        $gene[$i] = substr($rawchr,$pos,$nbits);
                        $pos = $pos + $nbits;
                    }
                    else{
                        my $gene_casuale = irand($nAmGeni);
                        $gene[$i] = $AmGeni[$gene_casuale];
                    }
                }
                else
                {
                    $gene[$i] = $ChromosomIN[$k];
                    $k++;
                }
                my $name = "$parametro$covariata";
                $ParamCov_Gene{$name} = $gene[$i];
                
                #****************************************************************************************
                #                     find relations for continuous covariates                          #
                #****************************************************************************************
                unless ($gene[$i] eq $AmGeni[0])
                {
                    if(exists $relationstart{$parametro})
                    {
                        my @covariatePerParametro= @{$relationstart{$parametro}};
                        $indiceCovariata = $#covariatePerParametro;
                        $indiceCovariata++;
                    }
                    else
                    {
                        $indiceCovariata = 0;
                    }
                    
                    $relationstart{$parametro}[$indiceCovariata] = $covariata;
                    
                    ($maxTheta, $relazioniRef, $RelDSContinuous, $paramCovListRef) = DS_continuous($gene[$i],$covariata,$maxTheta,$meannCVariables{$covariata} ,$parametro, $FILE, \@ValidStatesContinuous, \@AmGeni, \%paramCovList, $sdCVariables{$covariata});
                    %relazioni = %{$relazioniRef};
                    %paramCovList= %{$paramCovListRef};
                    @VectDSContinuous = (@VectDSContinuous, $RelDSContinuous);
                    $paramCov_Relazione{$covariata.$parametro} = $relazioni{$gene[$i]};
                    
                }
                if ($gene[$i] eq $AmGeni[0])
                {
                    $paramCov_Relazione{$covariata.$parametro} = 0;
                }
                
                @Chromosom = (@Chromosom, $gene[$i]);
            }
        }
        $nContinuousRelations = $#Chromosom;
        
        #****************************************************************************************
        #                         DEFINITION START CATEGORICAL COVARIATE                        #
        #****************************************************************************************

        if (exists $hash{$covariata} && any{$_ eq $covariata} @categoricalVariables)
        {
            my $refCov = $hash{$covariata};
            my @valore1 = @{$refCov};
            my $relazioniRef;
            
            my $nelementi = $#valore1+1;
            
            for(my $i=0; $i<$nelementi; $i++)
            {
                my $parametro = $hash{$covariata}[$i];
                my $ixx = 0;
                #****************************************************************************************
                #                                    find gene                                          #
                #****************************************************************************************
                if($G==0)
                {
                    if ($popfile ne ''){
                        $gene[$i] = substr($rawchr,$pos,1);
                        $pos = $pos + 1;
                    }
                    else{
                    my $gene_casuale = irand($nAmGeniCat);
                    $gene[$i] = $AmGeniCat[$gene_casuale];
                    }
                }
                else
                {
                    $gene[$i] = $ChromosomIN[$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];
    
                @Chromosom = (@Chromosom, $gene[$i]);
            }
        }
    }
   
    my @c= @Chromosom;

    #****************************************************************************************
    #                              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, \%minimumCVariables, \%maximumCVariables, \%meannCVariables, \@VectDSContinuous, \%sdCVariables);
    
    print $FILE("@line[$indexOmega..$#line]\n");
    return (\%relationstart, \@c, $maxTheta, $nContinuousRelations, \%paramCov_Relazione, \%paramCovList);
}

# ************************************* 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);
    
    for(my $i=0; $i<= $#valhash; $i++)
    {
        my $parametro = $valhash[$i];
        
        for (my $row = 0; $row<=$#line; $row++)
        {
            if($line[$row] =~ m/TV$parametro=THETA\(\d+\)/)
            {
                # Rewrite relation
                my @rewrite = ("$line[$row]\n","TV$parametro=$parametro","COV*TV$parametro\n");
                $line[$row] = join("", @rewrite);
            }
        }
    }
    
    # stampa il modello
    print $FILE("@line[$indexPK+1..$indexOmega-1]");
}

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

sub printInitialTheta
{
    
    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 %minimumCVariables = %$minref;
    my $maxref = $_[8]; my %maximumCVariables = %$maxref;
    my $meannref = $_[9]; my %meannCVariables = %$meannref;
    my $dsref= $_[10]; my @VectDSCont= @$dsref;
    my $sdref = $_[11]; my %sdVariables = %$sdref;
    
    my $k=0;
 
    my @covariate = (@categoricalVariables, @continuousVariables);
    
    for(my $ix=0; $ix<=$#VectDSCont; $ix++)
    {
        for(my $i=0; $i<= $#covariate; $i++)
        {
            my $covariata = $covariate[$i];
            for(my $j=0; $j<=$#parametri; $j++)
            {
                my $name = "$parametri[$j]$covariate[$i]";
                
                if($name eq $VectDSCont[$ix])
                {

                    if(any {$_ eq $covariate[$i]} @categoricalVariables)
                    {

                        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");
                        }
                    }
 

                    if(any {$_ eq $covariate[$i]} @continuousVariables)
                    {

                        my $gene = $ParamCov_Gene{$name};
 
                        my $relationParamCov = $Relazioni{$gene};

                        my $estremo_inferiore =FormatShort((-1)/($maximumCVariables{$covariate[$i]} - $meannCVariables{$covariate[$i]}));

                        my $valore_iniziale = FormatShort((-0.001)/($minimumCVariables{$covariate[$i]} - $meannCVariables{$covariate[$i]}));

                        my $estremo_superiore = FormatShort((-1)/($minimumCVariables{$covariate[$i]} - $meannCVariables{$covariate[$i]}));          
                        
                        my $valHS = 1;
                        
                        given($relationParamCov)
                        {   
                            when('2') # Linear Model
                            {
                                print $FILE ("\$THETA ($estremo_inferiore,$valore_iniziale,$estremo_superiore); $parametri[$j]$covariate[$i]\n");
                            }
 
                            when('3') # Hockey - Stick or piece - wise linear
                            {
                                print $FILE ("\$THETA (-10,$valore_iniziale,$estremo_superiore); $parametri[$j]$covariate[$i]1\n");
                                my $valore_iniziale2 = FormatShort((-0.001)/($maximumCVariables{$covariate[$i]} - $meannCVariables{$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')
                            {
                                    print $FILE ("\$THETA (-10,0.001,10); $parametri[$j]$covariate[$i]\n");
                            }
                            
                            when('5')
                            {
                                if($parametri[$j] =~ m/CL/)
                                {
                                    print $FILE ("\$THETA (-4,0.75,4); $parametri[$j]$covariate[$i]\n");
                                }
                                elsif($parametri[$j] =~ m/V/) # Per il volume
                                {
                                    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 $ChromosomRef = $_[2]; my @Chromosom = @{$ChromosomRef};
    my $WholeGenesRef = $_[3]; my @WholeGenes = @{$WholeGenesRef};
    my $CountRef = $_[4]; my %count = %{$CountRef};
    
    my %countResults;
    
    unless($e!=0)
    {
        for(my $index=0; $index<=$numberofStates+1; $index++)
        {
            for(my $jndex=0; $jndex<=$#Chromosom; $jndex++)
            {
                $count{$WholeGenes[$index]}[$jndex] = 0; # inizializzo il vettore count1 come un vettore di 0
            }
        }
    }
    
    for(my $js= 0 ; $js<=$#Chromosom; $js++)
    {
        my $countResultsRef = countGenes($Chromosom[$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;
        }
        
    }
    
    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;
    
    my @countMax = @{$count{$WholeGenes[0]}};
    
    for(my $i=0; $i<=$#countMax; $i++)
    {
        $geneMax[$i]= $WholeGenes[0];
    }
    
    
    for(my $i=1; $i<=$#WholeGenes; $i++)
    {
        for(my $j = 0; $j<=$#countMax; $j++)
        {
            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 $timer = $_[1];
    my @args = ("timeout","-k","5","$timer","execute", "Chromosom$e.mod" , "-dir=resultDirectory$e", "-silent");
    system(@args);
    
        my @out = qx(ls ./resultDirectory$e/raw_results_Chromosom$e.csv);
        chomp @out;
        if (defined $out[0]){
            return 0;
        }
        else{
            return 1;
        }

}

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

sub RouletteWheel
{
    my $ref1 = $_[0]; my %FitnessHash = %{$ref1};
    my $ref2 = $_[1]; my %population = %{$ref2};
    
    my @FitnessVect = values %FitnessHash;
    
    my $FitnessSum = sum(@FitnessVect);
    
    my @sliceWidth;
    for(my $i=0; $i<=$#FitnessVect; $i++)
    {
        $sliceWidth[$i] = $FitnessVect[$i] / $FitnessSum;
    }
    
    my @WheelBounds;
    for(my $i=0; $i<=$#FitnessVect; $i++)
    {
        $WheelBounds[$i] = sum(@sliceWidth[0..$i]);
    }
    $WheelBounds[-1] = 1;
    
    my $nrand = rand(1);
    
    my $ParentIndex;
    for(my $i=0; $i<=$#WheelBounds; $i++)
    {
        if($WheelBounds[$i]>$nrand)
        {
            $ParentIndex = $i;
            last;
        }
    }
    
    my @ParentCromosom = @{$population{"$ParentIndex"}};
    @ParentCromosom = split('', join('', @ParentCromosom));
    
    return \@ParentCromosom;
}
#****************************** Tournament Selection ********************************

sub TournamentSelection
{
    my $popSize = $_[0];
    my $ref1 = $_[1]; my %FitnessHash = %$ref1;
    my $ref2 = $_[2]; my %population = %$ref2;
  
    my $parametroTS = 0.75;
    my $IndexGenitore;
    
    my $indice1 = irand($popSize);
    my $indice2 = irand($popSize);

    
    my $FitnessBest = max($FitnessHash{$indice1}, $FitnessHash{$indice2});
    
    my $randomValue = rand(1);
    
    given($randomValue)
    {
        when($randomValue < $parametroTS)
        {

            if($FitnessHash{$indice1}>=$FitnessHash{$indice2})
            {
                $IndexGenitore = $indice1;
            }
            else
            {
                $IndexGenitore = $indice2;
            }
        }
        
        when($randomValue >= $parametroTS)
        {
            if($FitnessHash{$indice1}>=$FitnessHash{$indice2})
            {
                $IndexGenitore = $indice2;
            }
            else
            {
                $IndexGenitore = $indice1;
            }
        }
    }
    my @ParentCromosom = @{$population{$IndexGenitore}};
    @ParentCromosom = split('', join('', @ParentCromosom));
    return(\@ParentCromosom);
}


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

sub onePoinCrossOver
{
    my $ref1 = $_[0]; my @ParentCromosom1 = @{$ref1};
    my $ref2 = $_[1]; my @ParentCromosom2 = @{$ref2};

    my $index = irand($#ParentCromosom1);

    my @offspring1 = (@ParentCromosom1[0..$index] , @ParentCromosom2[$index+1..$#ParentCromosom2]);

    my @offspring2 = (@ParentCromosom2[0..$index] , @ParentCromosom1[$index+1..$#ParentCromosom1]);

    return(\@offspring1, \@offspring2);
}


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

sub splitNewpop
{
    my $nContinuousRelations = $_[0];
    my $countsOffspring = $_[1];
    my $nGene = $_[2];
    my $ref1 = $_[3]; my %Newpopulation = %{$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++)
    {
        my $ac=0;
        my $Chromosom = $Newpopulation{"$i"};
        
        for(my $j=0; $j<=$nContinuousRelations; $j++)
        {
            $NuovaPop{"$i"}[$j] = substr $Chromosom, $ac, $numberofBits;
            
            if  (! any {$_ eq "$NuovaPop{$i}[$j]"} @AmGenes)
            {
                $NuovaPop{"$i"}[$j] = $MostFrequent[$j];
            }
            
            $ac = $ac + $numberofBits;
        }
        
        for(my $j=$nContinuousRelations+1; $j<=$nGene; $j++ )
        {
            $NuovaPop{"$i"}[$j] = substr $Chromosom, $ac, 1;
            $ac++;
            
        }
    }
    
    return \%NuovaPop;
    
}

# ******************************** findFitness ********************************
sub findFitness
{
    my $e = $_[0];
    my $maxTheta = $_[1];
    my $MAXTHETA = $_[2];
    my $ref2 = $_[3]; my %FitnessHash= %{$ref2};
    my $ref = $_[4];    my @UPVbase = @{$ref};
    my $indexCorr = $_[5];
    $CWD ="resultDirectory$e";
    
    my $csvname = "raw_results_Chromosom$e.csv";
    open FILE ,'<' ,$csvname , or die "Error occurred: the file raw_results_Chromosom$e.csv does not exists \n";
    my @rawResultFile = <FILE>;
    close FILE;

    my $Resultsnames = $rawResultFile[0];
    my @splitResultNames = split(/,\D/, $Resultsnames);
    
    my $Resultvalues = $rawResultFile[1];
    my @splitResultValues = split(/,/, $Resultvalues);
     
    my %Rawresults;
    
    my $theta_index;
    my $first_omega_index;
    my $last_omega_index;

    for (my $i=0; $i<=$#splitResultNames; $i++)
    {
        $splitResultNames[$i] =~ s/[^a-zA-Z0-9]//g;
        $splitResultNames[$i] =~ s/[^a-zA-Z0-9]//g;
        $Rawresults{$splitResultNames[$i]} = $splitResultValues[$i]; 
        if ($splitResultNames[$i] eq 'ofv'){
            $theta_index = $i+1;
            next;
        }
        if ($splitResultNames[$i] eq 'OMEGA11'){
            $first_omega_index=$i;
            next;
        }
        if ($splitResultNames[$i] eq 'SIGMA11'){
            $last_omega_index=$i-1;
        }
    }

    my $num_theta = (-1+sqrt(1+8*($last_omega_index-$first_omega_index+1)))/2; 
    my @UPV;
	my @modifiedUPV;
    my @normalUPV;
	my @modifiednormalUPV;
    my $NewFitValue=0;
    
    for (my $i=0; $i<$num_theta; $i++){
        $UPV[$i] =  (sqrt(exp($splitResultValues[$first_omega_index+$i*$num_theta])-1)*$splitResultValues[$theta_index+$i])**2;
        $modifiedUPV[$i] =(sqrt(exp($splitResultValues[$first_omega_index+$i*$num_theta])-1))**2;
		$normalUPV[$i] = 1 - $UPV[$i]/$UPVbase[$i];
		$modifiednormalUPV[$i] = 1 - $modifiedUPV[$i]/$UPVbase[$i];
        $NewFitValue += $normalUPV[$i];
    }
    
    # *******************************************************************************
    # ******************************** Fitness Value ********************************
    # *******************************************************************************

    my $numeroParametri = $maxTheta - $MAXTHETA;
    my $objfun = $Rawresults{ofv};
    my $FitValue =  (-1)*($objfun+3.84*($numeroParametri)); #
    
    $CWD = '..';
    
	push @normalUPV, @modifiednormalUPV;
    my $normalizedUPVs = join (',' , @normalUPV);
    
    return ($FitValue, $objfun, $normalizedUPVs);
}

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

sub trim
{
    my @input= @_;
    my $nelementi = $#input;
    
    given($nelementi)
    {
        when($nelementi eq 0)
        {
            my $string = $_[0];
    
            $string =~ s/\s+//;
            $string =~ s/\s*//;
            $string =~ s/\n//;
    
            return $string;
        }

        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 *****************************************

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";
    }
}


# ******************************** LOG2 ***********************************
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 *********************
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 ***********************************
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/)
    {
        $fattoreMolt = substr($numeroDecimale, index($numeroDecimale, 'e'), length($numeroDecimale)- index($numeroDecimale,'e'));
        $numeroDecimale = substr($numeroDecimale, 0, index($numeroDecimale, 'e')-1);   
     }
    
    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);  }
}
# **********************************  mean ***********************************

sub computemean {
    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 mean: function value n. $n is not valid!\n"
        }
    }
    return sum(@_)/@_;
}

# ********************************** var ***********************************

sub computevariance
{
    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 mean: function value n. $n is not valid!\n"
        }
    }
    
    my $SUM = 0;
    
    for (my $i=0; $i<= $#Vector; $i++)
    {
        $SUM = $SUM + ( ($Vector[$i] - computemean(@Vector))**2);
    }
    
    my $VARIANCE = $SUM / ($#Vector+1);
    
    return $VARIANCE;
}

# ********************************** mirror ***********************************

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  ***********************************

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  ***********************************

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  ***********************************

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
	{

        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 $row_lst(@psnlst)
        {
            if($row_lst =~ m/TERM:/)
			{
                $indexTerm = $indice;
				last;
			}
			$indice++;
		 }	
    
        $TermMinimization = $psnlst[$indexTerm+1];$TermMinimization=~s/\s+//;$TermMinimization=~s/\s+//;
        
        foreach my $row(@psnlst)
        {
             if ($row =~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);
	
}
# ****************************  compute_correlation_index  *****************************
sub compute_correlation_index {
    my $ref0 = $_[0]; my %paramCovList = %{$ref0};
    my $ref1 = $_[1]; my %CorrMatrix = %{$ref1};
    
    my $sum=0;
  
    foreach my $param (keys %paramCovList){
        my @covs=@{$paramCovList{$param}};
        for (my $i=0; $i<= $#covs; $i++){
        #print $covs[$i];
            for (my $j=$i+1; $j<= $#covs; $j++){
                $sum = $sum + abs($CorrMatrix{$covs[$i]}{$covs[$j]});
            }
        }
    }

    return $sum;

}