#!/usr/bin/perl

#---------------------------------------------------------------------------------
# Daniel Segre'
# May 17, 2002
# dsegre@genetics.med.harvard.edu
#---------------------------------------------------------------------------------
#
print "\n\n";
print "******************************************************************\n";
print "Program doko_May17_2002.pl\n";
print "by Daniel Segre', dsegre\@genetics.med.harvard.edu\n";
print scalar localtime;
print "\n";
print "******************************************************************\n";

#-----------------------------------------------------------------------------------
# HOME DIRECTORY AND ENVIRONMENTAL VARIABLES:
#-----------------------------------------------------------------------------------
$homedir="/home/dsegre";
# ( On Llama: substitute /home1/dsegre/ )
print "Home directory: $homedir\n";
print "\n";

#print "NOTE: QP will run only after defining the following sys variables:\n";
#$syscall1=qq(export OSL_HOME=IBM_optimization/Linux_QP/solutions);
#$syscall2=qq(export PATH=\$PATH:\$OSL_HOME/bin);
#print "$syscall1\n";
#print "$syscall2\n";
#print "\n";

# NOTE: Initialization of system variables for QP:
# The Quadratic Programming part will not work unless the 
# following environmental variables are set:
#
# Current setup on lab laptop:
# export OSL_HOME=IBM_optimization/Linux_QP/solutions
# export PATH=$PATH:$OSL_HOME/bin
#
# On Cloister (for userid dsegre) the followng is currently done 
# automatically when logging in (in file .bash_profile)
# 
# export OSL_HOME="/home/dsegre/IBM_optimization/Linux_QP/solutions"
# export PATH="$PATH:$OSL_HOME/bin"
#
# On Llama, the path is /home1/... instead:
# export OSL_HOME="/home1/dsegre/IBM_optimization/Linux_QP/solutions"
# export PATH="$PATH:$OSL_HOME/bin"
#

#-----------------------------------------------------------------------------------
# READ PARAMETERS FILES:
#-----------------------------------------------------------------------------------
print "Reading parameter files...\n";
$S_matrix_filename="S.par";
$objvector_filename="objvector.par";
$constraints_filename="constraintsJ.par";
read_S_par($S_matrix_filename);
print "S matrix read from $S_matrix_filename \n";
read_constraints_par($constraints_filename);
print "Constraints read from $constraints_filename\n";
read_objvector_par($objvector_filename);
print "Objective function read from $objvector_filename\n";
print "Done.\n";
print "N. of fluxes : $N \n";
print "N. of metabolites : $M \n";

#-----------------------------------------------------------------------------------
# BASIC SIMULATION PARAMETERS:
#-----------------------------------------------------------------------------------

# OPTIONS:
$MODE=2; # 0=sensitivity loop ; 1=basic LP & QP loop ; 2=single mutant
print "MODE=$MODE  (0=sensitivity loop ; 1=basic LP & QP loop ; 2=single mutant) \n";

# SWITCH FOR DELETION OF ISOZYMES: (0=don't delete ; 1=delete)
$elimin_isozymes=1;
print "Isozymes deletion switch (0=don't delete ; 1=delete) = $elimin_isozymes\n";

# QP kind of computation: (NOTE: as of today Jan 29, Normalized doesn't really work...)
$normalized_QP=0;   # (0=normal QP ; 1=QP with fluxes normalized with respect to their wt)
if ($normalized_QP==1) {
  print "Running QP in NORMALIZED MODE\n";
} else {
  print "Running QP in standard mode\n";
}

# Loop parameters: (between 1 and $N)
$first_mut=1;
$last_mut=$N;
#  $first_mut=1;
#  $last_mut=11;
if ($MODE!=2) {
  print "Loop will be performed from flux n.$first_mut to flux n.$last_mut. \n";
}

# Mutant being analyzed in sensitivity analysis
if ($MODE==0) {
  $n_mut_for_sensitivity=7; # order number (1 to N)
  $index_mut_sensitivity=$n_mut_for_sensitivity-1;  # array position (0 to N-1)
  print "MUTANT BEING ANALYZED FOR SENSITIVITY =n.$n_mut_for_sensitivity ($enzyme[$index_mut_sensitivity])\n";
}

# Fraction of flux remaining in partial k.o.:
if ($MODE==0) {
  $activity_factor=0.9;
  print "Activity Factor = $activity_factor\n";
}

#$Vgr_position=970; # position of growth flux for complete Ecoli (with f and b separate)
$Vgr_position=602; # position of growth flux for complete Ecoli
print "Postion of Vgro for QP = $Vgr_position\n";

# Boundary to be used in (quasi) complete knock-out:
$ko_threshold_value=0;

$write_fluxes_flag=1;


#----------------------------------------------------------------------------------
# ASSIGN mps FILENAMES 
#----------------------------------------------------------------------------------
$LP_mps_wt_filename="Ecoli_complete_wt.mps";

$QP_mps_filename="ecoli_QP_mut.mps";
$LP_mps_filename="ecoli_LP_mut.mps";


#-----------------------------------------------------------------------------------
# INITIALIZE OUTPUT FILES
#-----------------------------------------------------------------------------------
print "File initialization...";
open(FID, ">QP_LP_Vgrowth.dat") or die "Can't open file: $!\n";
close(FID);
open(FID, ">QP_fluxes.dat") or die "Can't open file: $!\n";
close(FID);
open(FID, ">LP_fluxes.dat") or die "Can't open file: $!\n";
close(FID);
open(TMPID,">tmp")  or die "Can't open file: $!\n";
close(TMPID);
open(FID, ">summary_table.dat") or die "Can't open file: $!\n";
close(FID);
print "Done. \n";


#-----------------------------------------------------------------------------------
#  Calculate wild type fluxes (LP):
#-----------------------------------------------------------------------------------
print "Preliminary part : Calculate wild type with LP \n";

# Copy original constraints into LP specific constraints
foreach $a (keys %constraints){
  $constraints_LP{$a}=$constraints{$a};
}
data2mps_LP($LP_mps_wt_filename);                # Generate LP mps file
do_LP($LP_mps_wt_filename);                      # Performs optimization
$LP_wt_obj=0;
($LP_wt_obj,@LP_wt_flux)=read_LP_output();       # Read output

print  "\n";
print  "========================================================================\n";
print  "WILD TYPE PREDICTION OUTPUT:\n";
print  "V_Growth = \t$LP_wt_obj\n";

# NOTE: following part is S.par specific for complete E.coli!
print  "Optimized UPTAKE fluxes:\n";
print  "\tGLCxtI (764) = \t$LP_wt_flux[763]  \n";
print  "\tO2xtI (924) = \t$LP_wt_flux[923]  \n";
print  "\tNH3xtI (898) = \t$LP_wt_flux[897] \n";
print  "========================================================================\n";

print  "\n";

# Writes wt fluxes to file:
open(FID,">wt_fluxes.dat");
foreach $oneflux (@LP_wt_flux) {
  print FID "$oneflux \n";
}


#-----------------------------------------------------------------------------------
# DO SINGLE MUTANT (FROM STDIN)
#-----------------------------------------------------------------------------------
# Old part: get mutant number and upper bound from stdin:
if ($MODE==2) {

  # Prepare for LP:
  foreach $a (keys %constraints){
    $constraints_LP{$a}=$constraints{$a};
  }


  # Prepare for QP:
  if ($normalized_QP==1) {  # prepare for normalized QP
    $c=0;
    foreach $flux(@enzyme) {
      if ($LP_wt_flux[$c]!=0) {  
        $objhash_QP{$flux}=1;   # for nonzero fluxes, normalized coord. is 1 
      } else {
        $objhash_QP{$flux}=0;   # for zero fluxes, stays zero
      } 
      $c++;
    }

    $c=0;
    foreach $a (keys %constraints){
      for ($i_const=0 ; $i_const<=1 ; $i_const++) {
        $checkandum=$constraints{$a}[$i_const];
        if (($checkandum != "Inf")&&($checkandum != "-Inf")&&($checkandum != "0")&&($LP_wt_flux[$c]!=0)) {         
          $constraints_QP{$a}[$i_const]=$constraints{$a}[$i_const]/$LP_wt_flux[$c]; 
        } else {
          $constraints_QP{$a}[$i_const]=$constraints{$a}[$i_const];
        }
      }
      $c++;
    }

  } else {  # ... or prepare for REGULAR QP (not normalized)

    $c=0;
    foreach $flux(@enzyme) {
      $objhash_QP{$flux}=$LP_wt_flux[$c];
      $c++;
    }

    foreach $a (keys %constraints){
      $constraints_QP{$a}=$constraints{$a};
    }
  }

  print "Number of genes to be deleted or constrained: ";
  $howmanymuts=<STDIN>;
  chomp($howmanymuts);

   # Following part reads from STDIN the genes knocked out,
   # as modified by Wayne.
   # The format of input is the gene NAMES, not the numbers, as 
   # in previous versions, or in loop MODEs

    for ($i_manymuts=1 ; $i_manymuts<=$howmanymuts ; $i_manymuts ++){
      $mutant_allowed=0;
      while ($mutant_allowed==0) {
        print  "Gene (n.$i_manymuts) = ";
        $name_of_mutant=<STDIN>;
        chomp($name_of_mutant);
        for $enz_check (@enzyme) {
	  if ($enz_check eq $name_of_mutant) {
	    $mutant_allowed = 1;  
          }
        }
        if ($mutant_allowed == 0) {
          for $enz_check (@enzyme) {
	    if (index($enz_check, $name_of_mutant) == 0 && $mutant_allowed == 0) { 
	      $name_of_mutant = $enz_check;
              print  "changing name of gene to $name_of_mutant\n";
	      $mutant_allowed=1;
            }
          }      
        }
        if ($mutant_allowed == 0) {
          $name_of_mutant = uc($name_of_mutant);
          for $enz_check (@enzyme) {
            if (index($enz_check, $name_of_mutant) == 0 && $mutant_allowed == 0) { 
	        $name_of_mutant = $enz_check;
                print  "changing name of gene to $name_of_mutant\n";
	        $mutant_allowed=1;
            }
          }
        }      
        if ($mutant_allowed == 0) {
              print  "Gene name not allowed. Try again.\n";
        }
      }
    
      print  "Name of gene: $name_of_mutant\n";  

  # OLD PART (receiving number of genes rather than names)
  #for ($i_manymuts=1 ; $i_manymuts<=$howmanymuts ; $i_manymuts ++){
  #  print "Mutant (n.$_manymuts) = ";
  #  $i_mut=<STDIN>;
  #  chomp($i_mut);#
  #   $i_mut_array=$i_mut-1;
  #   $name_of_mutant=$enzyme[$i_mut_array];
  #   print "Name of mutant: $name_of_mutant\n";  

    print "Lower bound for flux ($name_of_mutant) = ";
    $lowerbound=<STDIN>;
    chomp($lowerbound);

    print "Upper bound for flux ($name_of_mutant) = ";
    $upperbound=<STDIN>;
    chomp($upperbound);

    # SINGLE LP extra constraint:  
    $constraints_LP{"$name_of_mutant"}=["$lowerbound", "$upperbound"];

    # SINGLE QP:
    $constraints_QP{"$name_of_mutant"}=["$lowerbound", "$upperbound"];

  }

  data2mps_LP($LP_mps_filename);
  do_LP($LP_mps_filename);
  ($LP_obj,@LP_flux)=read_LP_output();


  data2mps_QP($QP_mps_filename);
  do_QP($QP_mps_filename);
  $QP_obj=0;
  ($QP_obj,@QP_flux)=read_QP_output();

#  if ($normalized_QP==1) {  # retransform fluxes in absoulte coordinates
#    for ($c=0 ; $c<N ; $c++) {
#      if ($QP_flux[$c] !=0){
#        $QP_flux[$c]=$QP_flux[$c]*$LP_wt_flux[$c];
#      }
#    }
#  }


  $LP_Vgr=$LP_flux[$Vgr_position-1];
  $QP_Vgr=$QP_flux[$Vgr_position-1];

  write_output_to_file();
  
  print  "\n";
  print  "========================================================================\n";
  print  "OPTIMAL MUTANT PREDICTION OUTPUT (LP):\n";
  print  "V_Growth = \t$LP_Vgr\n";
  # NOTE: following part is S.par specific for complete E.coli!
  print  "MUTANT (LP) UPTAKE fluxes:\n";    
  print  "\tGLCxtI (764) = \t$LP_flux[763]  \n";
  print  "\tO2xtI (924) = \t$LP_flux[923]  \n";
  print  "\tNH3xtI (898) = \t$LP_flux[897] \n";
  print  "========================================================================\n";


  
  print  "\n";
  print  "========================================================================\n";
  print  "SUBOPTMAL MUTANT PREDICTION OUTPUT (QP):\n";
  print  "V_Growth = \t$QP_Vgr\n";
  # NOTE: following part is S.par specific for complete E.coli!
  print  "MUTANT (QP) UPTAKE fluxes:\n";    
  print  "\tGLCxtI (764) = \t$QP_flux[763]  \n";
  print  "\tO2xtI (924) = \t$QP_flux[923]  \n";
  print  "\tNH3xtI (898) = \t$QP_flux[897] \n";
  print  "========================================================================\n";

  print "\n";

  
}


#-----------------------------------------------------------------------------------
# DO LOOP OF MUTANTS
#-----------------------------------------------------------------------------------

if ($MODE==1){

  #-----------------------------------------------------------------------------------
  #-----------------------------------------------------------------------------------
  # HERE START ACTUAL COMPUTATION OF ALL MUTANTS
  #-----------------------------------------------------------------------------------
  #-----------------------------------------------------------------------------------

  print "\nLoop on mutants start here ------------------------------------\n\n";
  for ($i_mut=$first_mut ; $i_mut<=$last_mut ; $i_mut++) {
    $i_mut_array=$i_mut-1;
    $name_of_mutant=$enzyme[$i_mut_array];
    print "Current mutant: enzyme n.$i_mut  ($name_of_mutant) \n";

    if ($elimin_isozymes==1) {
      # LOOKS FOR ISOENZYMES OF THE CURRENT ENZYME  ------------ --------------------------
      # reminder $S[metabolite][enzyme]
      undef %grouped_flux;
      print "111111 Searching for isoenzymes of the current enzyme\n";
      $num_alias=$i_mut_array;
        $grouping_counter=0;
        for ($flux_loop=0 ; $flux_loop<$N  ; $flux_loop++){
          $sum_squares=0;
          for ($metab_loop=0 ; $metab_loop<$M ; $metab_loop++){
            $sum_squares=$sum_squares + ( $S[$metab_loop][$flux_loop] != $S[$metab_loop][$num_alias] );
          }
          if ($sum_squares==0) {
            $grouped_flux{$grouping_counter}=$enzyme[$grouping_counter];
            print "111111 Isozyme $enzyme[$grouping_counter] added to hash \n";
          }
        $grouping_counter++;
        } # end loop on all fluxes, looking for isoenzymes of current alias
    
      print "111111 Summary of Isozymes:\n";
      for $num_grouped(keys %grouped_flux) {
        print "111111 Isoenzyme n.$num_grouped ($grouped_flux{$num_grouped}) added\n";
      }
    } elsif ($elimin_isozymes==0) {
      undef %grouped_flux;
      print "000000 NOT searching for isoenzymes...\n";
      $grouped_flux{$i_mut_array}=$name_of_mutant;
    }
    

    # FOR EACH ISOZYME SEARCH FOR ALL REACTIONS IT CATALYZES (BY NAME SIMILARITY)  ---------------------------     
    undef %delenda_fluxa;    # empty hash
    for $num_grouped(keys %grouped_flux) { # starts loop on all isozymes found before 
      $name_of_iso=$enzyme[$num_grouped];
      print "222222 Now looking for multiple functions of $name_of_iso...\n";
      $delenda_fluxa{$num_grouped}=$name_of_iso; # enters in list the isozymne itself
      $elab_mutant=$name_of_iso;
      $elab_mutant =~ s/R$//;    # eliminates final R, if present (reversible)
      if ($elab_mutant =~ /(\D+)(\d+)/) {
        $pure_name=$1;
        $gene_version=$2;
        #print "Pure gene name: $pure_name  (Current version: $gene_version)\n";
        # Now searches through all fluxes for those with same enzyme:
        $alias_counter=0;
        for $pot_al_enzyme (@enzyme) {  # loop on potential alias enzymes...
          $pot_al_enz_tmp= $pot_al_enzyme ;
          $pot_al_enz_tmp =~ s/R$//;
          if ($pot_al_enz_tmp =~ /(\D+)(\d+)/) {
            $alias_pure_name=$1;
            $alias_gene_version=$2;
            if ($alias_pure_name eq $pure_name) {              
              $delenda_fluxa{$alias_counter}=$pot_al_enzyme;
              print "222222 Flux $pot_al_enzyme added to hash of fluxes to set to zero\n";
            }
          }
          $alias_counter++;
        }
      } 
    }
     
    open (FID,">>summary_table.dat");
    print FID "<$i_mut> ";
    print "222222 Final report of fluxes being set to zero: \n";
    for $num_alias(sort keys %delenda_fluxa) {
      print FID "$delenda_fluxa{$num_alias} ";
      print "222222  Alias n.$num_alias ($delenda_fluxa{$num_alias})\n";
    }
    print FID "\n";
    close(FID);

    # PREPARE CONSTRAINTS DEFAULT FOR LP----------------------------------------------
    # Sets to default vaue the constraints 
    undef %constraints_LP;
    foreach $a (keys %constraints){
      $constraints_LP{$a}=$constraints{$a};
    }
  
    # PREPARE CONSTRAINTS DEFAULT FOR QP--------------------------------
    # Sets to default vaue the constraints
    undef %constraints_QP;
    foreach $a (keys %constraints){
      $constraints_QP{$a}=$constraints{$a};
    }

   
    # STARTS LOOP ON ALL FLUXES TO BE DELETED TOGETHER, FOR DEFINING CONSTRAINTS
    for $num_alias(keys %delenda_fluxa) {
      if ($LP_wt_flux[$num_alias]>0) {
        $ko_bound_min=0;
        $ko_bound_max=$ko_threshold_value;
      } elsif ($LP_wt_flux[$num_alias]<0) {
        $ko_bound_min=-$ko_threshold_value;
        $ko_bound_max=0;
      } else {
        $ko_bound_min=0;
        $ko_bound_max=0;
      }

      #-----------------------------------------------------------------------------------------
      # Generates Mutant mps file for LP:
      #-----------------------------------------------------------------------------------------


      $constraints_LP{$enzyme[$num_alias]}=[$ko_bound_min, $ko_bound_max];
      $constraints_QP{$enzyme[$num_alias]}=[$ko_bound_min, $ko_bound_max];
      print "Constraint set for LP and QP: constraint($enzyme[$num_alias])= $ko_bound_min  $ko_bound_max \n";

    } # end loop on all fluxes to be deleted together

    data2mps_LP($LP_mps_filename);
    do_LP($LP_mps_filename);
    ($LP_obj,@LP_flux)=read_LP_output();
    print "Immediate LP result: obj=$LP_obj\n"; 

    #-----------------------------------------------------------------------------------------
    # Generates Mutant mps file for QP:
    #----------------------------------------------------------------------------------------- 

    # PREPARE OBJECTIVE FOR QP--------------------------------
    undef %objhash_QP;

    $c=0;
    foreach $flux(@enzyme) {
      $objhash_QP{$flux}=$LP_wt_flux[$c];
      $c++;
    }
    

    data2mps_QP($QP_mps_filename);
    do_QP($QP_mps_filename);
    $QP_obj=0;
    ($QP_obj,@QP_flux)=read_QP_output();
 
    $LP_Vgr=$LP_flux[$Vgr_position-1];
    $QP_Vgr=$QP_flux[$Vgr_position-1];

    push @Vgrowth_LP, $LP_Vgr;
    push @Vgrowth_QP, $QP_Vgr;
    
    $LP_Vgr_norm=$LP_Vgr/$LP_wt_obj;   
    $QP_Vgr_norm=$QP_Vgr/$LP_wt_obj;
 
    print "V_Growth (LP) =$LP_Vgr  (Normalized: $LP_Vgr_norm)\n";
    print "V_Growth (QP) =$QP_Vgr  (Normalized: $QP_Vgr_norm)\n";
    
    open(FID,">>summary_table.dat");
    printf FID "V_Growth (LP) =%5.3g (Normalized: %5.3g)\n",$LP_Vgr,$LP_Vgr_norm;
    printf FID "V_Growth (QP) =%5.3g (Normalized: %5.3g)\n",$QP_Vgr,$QP_Vgr_norm;
    print FID "\n";
    close(FID);

    write_output_to_file();

    print "Mutation Done\n\n";
  }; # End of loop on mutants 

  print "@Vgrowth_LP\n";
  print "@Vgrowth_QP\n";


}; # end if $MODE =1



#-----------------------------------------------------------------------------------
# DO SENSITIVITY ANALYSIS  (MODE=0)
#-----------------------------------------------------------------------------------

if ($MODE==0){

  #-----------------------------------------------------------------------------------
  #-----------------------------------------------------------------------------------
  # HERE START ACTUAL COMPUTATION OF HOW A SPECIFIC MUTANT (QP) VARIES BY
  # MOVING AROUND THE ORIGINAL WILD-TYPE OPTIMUM
  #-----------------------------------------------------------------------------------
  #-----------------------------------------------------------------------------------

  # LOOKS FOR ISOENZYMES OF THE ENZYME AROUND WHICH WE WANT TO STUDY SENSITIVITY  --------------------------
  # reminder $S[metabolite][enzyme]
  undef %grouped_flux;
  print "111111 Searching for isoenzymes of the current enzyme\n";
  $num_alias=$index_mut_sensitivity;
  $grouping_counter=0;
  for ($flux_loop=0 ; $flux_loop<$N  ; $flux_loop++){
    $sum_squares=0;
    for ($metab_loop=0 ; $metab_loop<$M ; $metab_loop++){
      $sum_squares=$sum_squares + ( $S[$metab_loop][$flux_loop] != $S[$metab_loop][$num_alias] );
    }
    if ($sum_squares==0) {
      $grouped_flux{$grouping_counter}=$enzyme[$grouping_counter];
      print "111111 Isozyme $enzyme[$grouping_counter] added to hash \n";
    }
    $grouping_counter++;
  } # end loop on all fluxes, looking for isoenzymes of current alias
  print "111111 Summary of Isozymes:\n";
  for $num_grouped(keys %grouped_flux) {
    print "111111 Isoenzyme n.$num_grouped ($grouped_flux{$num_grouped}) added\n";
  }

  # FOR EACH ISOZYME SEARCH FOR ALL REACTIONS IT CATALYZES (BY NAME SIMILARITY)  ---------------------------
  undef %delenda_fluxa;    # empty hash
  for $num_grouped(keys %grouped_flux) { # starts loop on all isozymes found before
    $name_of_iso=$enzyme[$num_grouped];
    print "222222 Now looking for multiple functions of $name_of_iso...\n";
    $delenda_fluxa{$num_grouped}=$name_of_iso; # enters in list the isozymne itself
    $elab_mutant=$name_of_iso;
    $elab_mutant =~ s/R$//;    # eliminates final R, if present (reversible)
    if ($elab_mutant =~ /(\D+)(\d+)/) {
      $pure_name=$1;
      $gene_version=$2;
      #print "Pure gene name: $pure_name  (Current version: $gene_version)\n";
      # Now searches through all fluxes for those with same enzyme:
      $alias_counter=0;
      for $pot_al_enzyme (@enzyme) {  # loop on potential alias enzymes...
        $pot_al_enz_tmp= $pot_al_enzyme ;
        $pot_al_enz_tmp =~ s/R$//;
        if ($pot_al_enz_tmp =~ /(\D+)(\d+)/) {
          $alias_pure_name=$1;
          $alias_gene_version=$2;
          if ($alias_pure_name eq $pure_name) {
            $delenda_fluxa{$alias_counter}=$pot_al_enzyme;
            print "222222 Flux $pot_al_enzyme added to hash of fluxes to set to zero\n";
          }
        }
        $alias_counter++;
      }
    }
  }

  print "222222 Final report of fluxes being set to zero: \n";
  for $num_alias(sort keys %delenda_fluxa) {
    print "222222  Alias n.$num_alias ($delenda_fluxa{$num_alias})\n";
  }
  print "\n";



  undef %constraints_QP;
  foreach $a (keys %constraints){
    $constraints_QP{$a}=$constraints{$a};
  }
  # INSTEAD OF NEXT LINE, SET TO ZERO SET OF FLUXES, AS IN MODE=1 XXXXXXXXXXXXXXXXXXXXXX
  # $constraints_QP{"$mutant_for_sensitivity"}=["0", "0"];


  # STARTS LOOP ON ALL FLUXES TO BE DELETED TOGETHER, FOR DEFINING CONSTRAINTS
    for $num_alias(keys %delenda_fluxa) {
      $constraints_QP{$enzyme[$num_alias]}=[0,0];
      print "Constraint set for QP: constraint($enzyme[$num_alias])= 0 0\n";
    } # end loop on all fluxes to be deleted together




  print "\nLoop on all directions for sensitivity start here ------------------------------------\n\n";
  for ($i_mut=$first_mut ; $i_mut<=$last_mut ; $i_mut++) {
    $i_mut_array=$i_mut-1;
    $name_of_mutant=$enzyme[$i_mut_array];
    
    # Checks whether flux for displacement is being contracted or expanded
    if ($activity_factor<1) {
      $upperbound=$activity_factor*$LP_wt_flux[$i_mut_array];
    } elsif ($activity_factor>1) {
      if ($LP_wt_flux[$i_mut_array]==0) {  
        $upperbound="0.1";
      } else {
        $upperbound=$activity_factor*$LP_wt_flux[$i_mut_array];
      } 
    }
   
    print "Current flux being bounded: enzyme n.$i_mut  ($name_of_mutant) ; ";

    print "FBA flux = $LP_wt_flux[$i_mut_array]  ;  ";
    print "Current upper bound = $upperbound \n";
    
    if ($upperbound!=0) {

      #-----------------------------------------------------------------------------------------
      # Generates Mutant mps file for LP:
      #-----------------------------------------------------------------------------------------

      undef %constraints_LP;
      foreach $a (keys %constraints){
        $constraints_LP{$a}=$constraints{$a};
      }

      # Actually sets constraint for flux along which the displacement is being done
      if ($activity_factor<1) {
        if ($upperbound>0) {
          if ($name_of_mutant =~ /R$/) {  
            print "******** Reversible Flux \n";
            $constraints_LP{"$name_of_mutant"}=["-Inf", "$upperbound"];
          } else {
            $constraints_LP{"$name_of_mutant"}=["0", "$upperbound"];
          }
        } elsif ($upperbound<0) {  # must be a reversible flux
          $constraints_LP{"$name_of_mutant"}=["$upperbound", "Inf"];
        }
      } elsif ($activity_factor>1) {
        $constraints_LP{"$name_of_mutant"}=["$upperbound", "10000"]; 
      }

      data2mps_LP($LP_mps_filename);
      do_LP($LP_mps_filename);
      ($LP_obj,@LP_flux)=read_LP_output();
      print "Immediate LP result: obj=$LP_obj\n"  ;

      #-----------------------------------------------------------------------------------------
      # Generates Mutant mps file for QP:
      #generate_mps_mut_QP();
      #-----------------------------------------------------------------------------------------

      undef %objhash_QP;
      $c=0;
      foreach $flux(@enzyme) {
        $objhash_QP{$flux}=$LP_flux[$c];
        #print "QPobj> c=$c ; flux=$flux ; LP_flux[$c]=$LP_flux[$c]\n";
        open(TMPID,">>tmp");
        printf TMPID "%5.3e\t",$objhash_QP{$flux};
        $c++;
      }
      print TMPID "\n";
      close(TMPID);

      data2mps_QP($QP_mps_filename);
      do_QP($QP_mps_filename);
      $QP_obj=0;
      ($QP_obj,@QP_flux)=read_QP_output();
 
      $LP_Vgr=$LP_flux[$Vgr_position-1];
      $QP_Vgr=$QP_flux[$Vgr_position-1];

    } else {
      print "Trivial case, skipping calculations...\n";
      $LP_Vgr=0;
      $QP_Vgr=0;
      @LP_flux=@LP_wt_flux;
      @QP_flux=@LP_wt_flux;

    }

    push @Vgrowth_LP, $LP_Vgr;
    push @Vgrowth_QP, $QP_Vgr;

    print "V_Growth (LP) =$LP_Vgr\n";
    print "V_Growth (QP) =$QP_Vgr\n";

    write_output_to_file();

    print "Displacement Done\n\n";
  }; # End of loop on mutants 

  print "@Vgrowth_LP\n";
  print "@Vgrowth_QP\n";


}; # end if $MODE =0






#----------------------------------------------------------------------------------------
#                            S U B R O U T I N E S
#----------------------------------------------------------------------------------------

sub do_LP {

print "Running Linear Programming...";
# EXECUTES Linear Programming using GNU Linear Programming Kit:
system("$homedir/gnusimplex/glpk-2.4.1/glpsol_dan.exe $_[0] > output_LP_trash");
system("cp dan_output LP_output.txt");
print "Done.\n";
}

#----------------------------------------------------------------------------------------

sub do_QP {

print "Running Quadratic Programming...";
# Removes previous QP output file"
if ( -e "QP_output.txt") {
	system("rm QP_output.txt");
}
# EXECUTES Quadratic Programming using IBM QP solutions:
# NOTE: this requires preliminary initialization with:
# export OSL_HOME=/home/dsegre/IBM_optimization/Linux_QP/solutions
# export PATH=$PATH:$OSL_HOME/bin
system("oslqslv < $_[0] -walltime=yes -imaxiter=20000 -dspace=1000000 -maxmin=min -output=QP_output.txt > output_QP_trash");
print "Done.\n";
}


#----------------------------------------------------------------------------------------

sub read_LP_output {
my (@LP_flux);
#print "Reading LP output file...";
open(LP_DATA, "LP_output.txt") or die "Can't open file: $!\n";

@lp_data = <LP_DATA>;
close(LP_DATA);
chomp(@lp_data);

# DETECTS Objective function value:
$counter=0;
$flag=0;
while ($flag==0) {
	$line=@lp_data[$counter];
	if ($line=~/Objective/) {
		@line_pieces=split /\s+/,$line;
		$LP_obj=-$line_pieces[3];
		#print "LP obj = $LP_obj\n";
		$flag=1;
	}
	$counter++;
}


# DETECTS Beginning of Fluxes values:
$flag=0;
while ($flag==0) {
        $line=@lp_data[$counter];
	if ($line=~/Column name/) {
		$flag=1;
	}
	$counter++;
}


# READS LP fluxes

$counter++;
while ($counter<$#lp_data) {
        $line=@lp_data[$counter];
        @line_pieces = split /\s+/, $line;
	$flux=$line_pieces[4];
	push @LP_flux, $flux;
	$counter++;
}

#print "Done.\n";
return ($LP_obj,@LP_flux);

}

#----------------------------------------------------------------------------------------

sub read_QP_output {
my (@QP_flux);
#print "Reading QP output file...";

open(QP_DATA, "QP_output.txt") or die "Can't open file: $!\n";

@qp_data = <QP_DATA>;
close(QP_DATA);
chomp(@qp_data);

# Check if a solution was found:
$counter=0;
$flag=0;
while ($flag==0) {
	$line=@qp_data[$counter];
	if ($line=~/1EKK0010I    Rows Section/) {
		$line_obj=@qp_data[$counter-1];
                @line_pieces=split /\s+/,$line_obj;
                #print "@line_pieces\n";
                $obj_piece=$line_pieces[7];
                $obj_piece =~ m/((\d|\.|-)+)--Optimal/;
                $QP_obj = $1;
		if ($line_pieces[7]=="infeasibilities") {
			print "*** INFEASIBLE ! ***\n";
		    print "$line_obj\n";
		};
                #print "QP Obj=$QP_obj\n";
		$flag=1;
	}
	$counter++;
}

if ($flag==0) {
print "PROBLEM with QP!!n";
print "supposedly QP_obj=$QP_obj\n";
}


# DETECTS Beginning of Fluxes values:
$flag=0;
while ($flag==0) {
        $line=@qp_data[$counter];
	if ($line=~/1EKK0011I/) {
		#print "$line\n";
		$flag=1;
	}
	$counter++;
}


# READS FLUXES:
$flag=0;
while ($counter <= $#qp_data) {
  $line=@qp_data[$counter];
  if ($line=~/EKK0064I/) {
    #$line=@qp_data[$counter];
    @line_pieces = split /\s+/, $line;
    $flux=$line_pieces[5];
    $number=$line_pieces[2];
    if ($flux eq '.') { $flux=0;};
    push @QP_flux, $flux;	
  }
  $counter++;
}
#print "$line\n";

#print "Done.\n";
return ($QP_obj,@QP_flux);

}

#----------------------------------------------------------------------------------------

sub write_output_to_file {

#print "Writing to file...";
open(FID, ">>QP_LP_Vgrowth.dat") or die "Can't open file: $!\n";
#print FID "$LP_Vgr\t$QP_Vgr\n";

#print FID "i_mut=$i_mut\t($name_of_mutant)\t\t";
#print FID "wt_flux=$LP_wt_flux[$i_mut_array]\t\t";
#print FID "flux_up_bnd=$upperbound \t\t";
#print FID "Vgr_LP=$LP_Vgr\t\tVgr_QP=$QP_Vgr\n";

print FID "$i_mut\t\t";
print FID "$LP_wt_flux[$i_mut_array]\t\t";
print FID "$upperbound \t\t";
print FID "$LP_Vgr\t\t$QP_Vgr\n";

close(FID);

if ($write_fluxes_flag==1){
  open(FID, ">>LP_fluxes.dat") or die "Can't open file: $!\n";
  for ($_=0; $_<$#LP_flux; $_++) {
    print FID "$LP_flux[$_]\t";
  }
  print FID "\n";
  close(FID);

  open(FID, ">>QP_fluxes.dat") or die "Can't open file: $!\n";
  for ($_=0; $_<($#QP_flux+1); $_++) {
    print FID "$QP_flux[$_]\t";
  }
  print FID "\n";
  close(FID);
};

#print "Writing done.\n";
}
#----------------------------------------------------------------------------------------

#****************************************************************************
# FOLLOWING SUBROUTINES ORIGINALLY FROM read_par_files_Jan12.pl
#****************************************************************************

sub read_S_par {

# READING THE STOICHIOMETRIC MATRIX PARAMETER FILE
#---------------------------------------------------------------------

# print "Reading file of S matrix...";
open(FID,"$_[0]");
@alllines=<FID>;
close(FID);
chomp(@alllines);

@enzyme = split /\s+/, $alllines[0];
$N=scalar(@enzyme);

shift(@alllines);
$M=scalar(@alllines);

$counter1=0;
foreach $line(@alllines) {
  ($currmetabolite, $rest)= split /\s+/, $line, 2;
  push @metabolite, $currmetabolite;
  @Sline= split /\s+/, $rest;
  for ($counter2=0 ; $counter2<$N ; $counter2++){
    $S[$counter1][$counter2]=$Sline[$counter2];
  }
  $counter1++;
}


# Generate hash of hashes with all fluxes:
for ($i=0 ; $i<$N ; $i++){
  #print "$enzyme[$i] ) ";
  for ($j=0 ; $j<$M ; $j++) {
    if ($S[$j][$i] != 0 ) {
      #print " $S[$j][$i] $metabolite[$j] ";
      $Shash{$enzyme[$i]}{$metabolite[$j]} = $S[$j][$i];
    }
  }
  #print "\n";	
}

# Prints list of enzymes into a matlab file:
open(FID,">enzyme.m");
print FID "Enzyme=\{\n";
for $k1 (@enzyme) {
  print FID "\'$k1\'\n";
}
print FID "\};\n";
close(FID);


# Prints list of reactions into a matlab file:
open(FID,">reaction.m");
print FID "Reaction=\{\n";
for $k1 (@enzyme) {
  print FID "\'$k1 ) ";
  for $k2 (keys %{$Shash{$k1}} ) {
   print FID "$Shash{$k1}{$k2} $k2 ";
  }
  print FID "\'\n";
}
print FID "\};\n";
close(FID);

# print "Done.\n";

}

#****************************************************************************

sub read_constraints_par {

# READING THE CONSTRAINTS PARAMETER FILE
#---------------------------------------------------------------------

# print "Reading file of constraints...";
open(FID,"$_[0]");
@alllines=<FID>;
close(FID);
chomp(@alllines);

foreach $line(@alllines) {
($number, $name, $min, $max)= split /\s+/, $line;
$constraints{$name} = [ $min , $max ];
}

#for $_ (keys %constraints) {
#  print "$_ : min=$constraints{$_}[0]  -  max=$constraints{$_}[1] \n";
#}

# print "Done.\n";

}

#****************************************************************************

sub read_objvector_par {


# READING THE OBJECTIVE FUNCTION PARAMETER FILE
#---------------------------------------------------------------------

# print "Reading file of objective function...";
open(FID,"$_[0]");
@alllines=<FID>;
close(FID);
chomp(@alllines);


foreach $line(@alllines) {
  ($name, $value)= split /\s+/, $line;
  $objhash{$name} = $value;
}


#foreach $_(@enzyme) {
#  print "$_ -> $objhash{$_}\n";
#}

# print "Done.\n";

}


#****************************************************************************

sub data2mps_LP {

# Generate mps file for LP

#print "Entered subroutine data2mps_LP \n";
#print "This will generate file $_[0] \n";

open(FID,">$_[0]");
close(FID);

open(FID,">$_[0]");

$mps_name='EColi';
$next_line="NAME          $mps_name";
print FID "$next_line\n";

# ROWS
$next_line="ROWS";
print FID "$next_line\n";

# Objective Function
$next_line="* Objective Function:";
print FID "$next_line\n";
$next_line=" N  COST";
print FID "$next_line\n";

for ($i=0 ; $i<$M ; $i++) {
  $next_line=" E  $metabolite[$i]";
  print FID "$next_line\n";
} 

$next_line="COLUMNS";
print FID "$next_line\n";

for $k1 (@enzyme) {
  $costcoeff=$objhash{$k1};
  if ($costcoeff != "0"){
    $next_line=sprintf "    %-8s  COST      %-5.3g",$k1,-$costcoeff;
    print FID "$next_line\n";
  }
  $ccc=0;
  for $k2 (keys %{$Shash{$k1}} ) {
    $coeff=$Shash{$k1}{$k2};
    if ($ccc==0){
      $next_line=sprintf "    %-8s  %-8s  %-5.3e", $k1,$k2,$coeff;
    } else {
      $next_line=sprintf "    %-8s  %-8s  %-5.3e", $k1,$k2,$coeff; 
    }    
    print FID "$next_line\n";
  }
}

$next_line="RHS";
print FID "$next_line\n";

$next_line="RANGES";
print FID "$next_line\n";

$next_line="BOUNDS";
print FID "$next_line\n";

for $c (sort keys %constraints_LP) {
  $min=$constraints_LP{$c}[0];
  $max=$constraints_LP{$c}[1];

  if ($min==$max){
    $next_line=sprintf " FX BND1      %-8s  %-5.3e", $c,$min; 
    print FID "$next_line\n";
  } elsif (($min=="-Inf")&&($max=="Inf")) {
    $next_line=sprintf " FR BND1      %-8s", $c; 
    print FID "$next_line\n";
  } elsif (($min=="-Inf")&&($max!="Inf")) {  # bound only from above
    $next_line=sprintf " MI BND1      %-8s", $c;
    print FID "$next_line\n";
    $next_line=sprintf " UP BND1      %-8s  %-5.3e", $c,$max;
    print FID "$next_line\n";
  } elsif (($min!="-Inf")&&($max=="Inf")) {  # bound only from below
    $next_line=sprintf " LO BND1      %-8s  %-5.3e", $c,$min;
    print FID "$next_line\n";
    $next_line=sprintf " PL BND1      %-8s", $c;
    print FID "$next_line\n";
  } else {
    $next_line=sprintf " LO BND1      %-8s  %-5.3e", $c,$min; 
    print FID "$next_line\n";
    $next_line=sprintf " UP BND1      %-8s  %-5.3e", $c,$max; 
    print FID "$next_line\n";
  }
}

$next_line="ENDATA";
print FID "$next_line\n";

close(FID);

#print "Subroutine data2mps_LP completed.\n";

}

#****************************************************************************

sub data2mps_QP {

# Generate mps file for LP

#print "Entered subroutine data2mps_QP \n";

open(FID,">$_[0]");
close(FID);

open(FID,">$_[0]");

$mps_name='EColiQP';
$next_line="NAME          $mps_name";
print FID "$next_line\n";

# ROWS
$next_line="ROWS";
print FID "$next_line\n";

# Objective Function
$next_line="* Objective Function:";
print FID "$next_line\n";
$next_line=" N  COST";
print FID "$next_line\n";

for ($i=0 ; $i<$M ; $i++) {
  $next_line=" E  $metabolite[$i]";
  print FID "$next_line\n";
} 

$next_line="COLUMNS";
print FID "$next_line\n";


for $k1 (@enzyme) {
  $costcoeff=$objhash_QP{$k1};
  if ($costcoeff != "0"){
    $next_line=sprintf "    %-8s  COST      %-5.3e",$k1,-$costcoeff;
    print FID "$next_line\n";
  }
  for $k2 (keys %{$Shash{$k1}} ) {
    $coeff=$Shash{$k1}{$k2};
    $next_line=sprintf "    %-8s  %-8s  %-5.3e", $k1,$k2,$coeff;
    print FID "$next_line\n";
  }
}

$next_line="RHS";
print FID "$next_line\n";

$next_line="RANGES";
print FID "$next_line\n";

$next_line="BOUNDS";
print FID "$next_line\n";

for $c (keys %constraints_QP) {
  $min=$constraints_QP{$c}[0];
  $max=$constraints_QP{$c}[1];

  if ($min==$max){
    $next_line=sprintf " FX BND1      %-8s  %-5.3e", $c,$min; 
    print FID "$next_line\n";
  } elsif (($min=="-Inf")&&($max=="Inf")) {
    $next_line=" FR BND1      $c";
    print FID "$next_line\n";
  } elsif (($min=="-Inf")&&($max!="Inf")) {  # bound only from above
    $next_line=sprintf " MI BND1      %-8s", $c;
    print FID "$next_line\n";
    $next_line=sprintf " UP BND1      %-8s  %-5.3e", $c,$max;
    print FID "$next_line\n";
  } elsif (($min!="-Inf")&&($max=="Inf")) {  # bound only from below
    $next_line=sprintf " LO BND1      %-8s  %-5.3e", $c,$min;
    print FID "$next_line\n";
    $next_line=sprintf " PL BND1      %-8s", $c;
    print FID "$next_line\n";
  } else {
    $next_line=sprintf " LO BND1      %-8s  %-5.3e", $c,$min; 
    print FID "$next_line\n";
    $next_line=sprintf " UP BND1      %-8s  %-5.3e", $c,$max; 
    print FID "$next_line\n";
  }
}

$next_line="ENDATA";
print FID "$next_line\n";

# Now add quadratic part of objective function:
$next_line="NAME          $mps_name";
print FID "$next_line\n";
$next_line="QSECTION";
print FID "$next_line\n";

for $k1 (@enzyme) {
  $next_line="                        1.0";
  substr($next_line,4,length($k1))=$k1;
  substr($next_line,14,length($k1))=$k1;
  print FID "$next_line\n";
}

$next_line="ENDATA";
print FID "$next_line\n";

close(FID);

}