/*
  This program is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

/* A Scalable Algorithm to Explore the Gibbs Energy Landscape of
   Genoma-Scale Metabolic Networks
   D. De Martino, M. Figliuzzi, A. De Martino and E. Marinari,
   PLoS Computational Biology, 
   code released May 2012.
*/

// This code corrects thermodynamic unfeasibilities among fluxes 
// and chemical potentials in  constraint-based models of metabolic networks

#include "iostream"
#include "math.h"
#include "stdlib.h"
#include "time.h"
#include "fstream"
#include "vector"
#include "string"
#include "sstream"

using namespace std;

int sign(double v){
  int s = ( v > 0 ) ? 1:-1;
  return s;
}

void usage(){
  cout << " Usage:  -m <stochiometric matrix filename> -p <chemical potentials filename> -f <fluxes filename>" << endl;

}

int main (int argc,char** argv){
  int M = 0;   
  int N = 0; 
  int ok = 0;
  // reading the input files
  fstream file;
  char *matrice = NULL;
  char *potenziali = NULL;
  char *flussi = NULL;
  int c;
  while((c = getopt(argc, argv, "m:p:f:"))!=-1) {
    switch(c) {			   
    case 'm':
      matrice = optarg;
      break;
    case 'p':
      potenziali = optarg;
      break;
    case 'f':
      flussi = optarg;
      break;  
    default:
      cout << "error parsing command line" << endl;
      usage();
      exit(1);
    }
  }

  file.open(matrice, ios::in);
  if(file.is_open()==0){
    cout << "the file " << matrice <<  " doesn't exist" << endl;
    exit(1);
  }

  do{
    string papa;
    while(file.eof()!=1){
      getline(file,papa);   // computing the size of the stochiometric matrix
      N++;
    }
    N--;  
    file.close();
    file.open(matrice, ios::in);
    double a;
    while(file >> a) M++;
    M/=N;
    file.close();
  }while(ok==1);

  std::vector<int> react[N];
  std::vector<int> metab[M];
  std::vector<double> react2[N];                      // definition of the stochiometrix matrix, flux and chemical potential vector
  std::vector<double> metab2[M];
  double flux[N];
  double chempot0[M], chempot[M];

  do{
    file.open(matrice, ios::in);       // reading the stochiometric matrix
    int i1 = 0;
    int j1 = 0;
    int y=0;        
    double pepe;
    while(file >> pepe){    
      if( pepe != 0 ){
	react[i1].push_back(j1);
	react2[i1].push_back(pepe);
	metab[j1].push_back(i1);
	metab2[j1].push_back(pepe);
      }
      j1++;
      if( j1 == M ){
                
	j1 = 0;
	i1++;  
      }
    }       
    file.close();     

    file.open(flussi,ios::in); // reading the flux vector
    if(file.is_open()==0){
      cout << "the file " << flussi <<  " doesn't exist" << endl; 
      exit(1);
    }
    for(int i = 0; i <= N-1; i++) file >> flux[i];
    file.close();
    file.open(potenziali,ios::in); // reading the chemical potential vector
    if(file.is_open()==0){
      cout << "the file " << potenziali <<  " doesn't exist" << endl;
      exit(1);
    }
    for(int i = 0; i <= M-1; i++)  file >> chempot0[i];
    file.close();                           
    for(int i = 0; i <= M-1; i++) chempot[i] = chempot0[i];
  } while(ok==1);

  int count = 0;
  int ourMax = 30 * N; 
  // ourMax is the time when the algorithm starts searching for loops
  int unsat[ourMax];
  do{  
    if( count == ourMax ){
      // for exceedingly large times search for loops among the unsat reactions
      cout << "Exceedingly large times: searching for loops" << endl;
      count = 0;
      std::vector <int> distinti;
      distinti.push_back(unsat[0]);
      for(int i = 1; i <= ourMax-1; i++){                             
	int oko = 0;
	for(int j = 0; j <= i-1; j++) if(unsat[i] == unsat[j])oko++;
	if(oko == 0) distinti.push_back(unsat[i]);
      }                   
      std::vector<int> matched;
      for(int i = 0; i <= distinti.size()-1; i++){
	int r = distinti[i]; 
	// search for distinct and  connected reactions among the unsat
	int match = 1;
	for(int j = 0; j <= react[r].size()-1; j++){
	  int k = 0;
	  int okei = 0;
	  do{
	    if(k != i){
	      int r1 = distinti[k];
	      for(int j1 = 0; j1 <= react[r1].size()-1; j1++){
		if(react[r][j] == react[r1][j1]) okei = 1;
	      }
	    }
	    k++;
	  } while( (k<distinti.size()) && (okei == 0) );
	  if(okei == 0) match = 0;
	}
	if(match == 1) matched.push_back(r);
      }
      distinti.clear();
      if( (matched.size() <= 26) && (matched.size()>0) ){
	cout << "searching among " << matched.size() << " reactions" << endl;
	std::vector<int> cycle;
	int check = 0;
	int lenght = 2; // exaustive search (max lenght 26)
	int index[matched.size()];
	for(int i = 0; i <= matched.size()-1; i++) index[i] = i;
	cout << "checking for loops large " << lenght << endl;
	do{
	  int nodes[lenght];
	  for(int i = 0; i <= lenght-1; i++) nodes[i] = matched[index[i]];
      	  int metab = 0;
	  int bell = 0;
	  do{
	    for(int i = 0; i <= lenght-1; i++){
	      int r = nodes[i];
	      int s = 1; 
	      if(flux[r]<0) s = -1;
	      for(int j = 0; j <= react[r].size()-1; j++){
		if(react[r][j] == metab) bell += s * react2[r][j];
	      }
	    }        
	    metab++;
	  } while( (metab<M) && (bell == 0) );
	  if(bell == 0){
	    check = 1;
	    for(int i = 0; i <= lenght-1; i++) cycle.push_back(nodes[i]);
	  }
	  index[lenght-1]++;
	  for(int i = 1; i <= lenght-1; i++){
	    if(index[lenght-i]>matched.size()-i){
	      index[lenght-i-1]++;
	      for(int k = 0; k <= i-1; k++){
		index[lenght-i+k] = index[lenght-i-1+k]+1;
	      }
	    }                     
	  }        
	  if(index[0]>matched.size()-lenght){                 
	    lenght++;
	    if( (lenght <= matched.size()) && (check == 0) ){
	      cout << "checking for loops large " << lenght << endl;
	    }
	    for(int i = 0; i <= matched.size()-1; i++) index[i] = i;        
	  }
	} while( (check == 0) && (lenght<matched.size()+1) );
	if(check == 1){ 
	  cout << "loop of size " << cycle.size() << " among" << endl; 
	  int corr;
	  double minimum = 1000;
	  for(int i = 0; i <= cycle.size()-1; i++){
	    if(fabs(flux[cycle[i]])<minimum){
	      minimum = fabs(flux[cycle[i]]);
	      corr = cycle[i];
	    }
	    cout << "reaction " << cycle[i] << " v_" << cycle[i] << "=" << flux[cycle[i]] << endl;
	  }                                              
	  cout << "correct them by " << minimum << endl;
	  for(int i = 0; i <= cycle.size()-1; i++){
	    int s = 1;
	    if( flux[cycle[i]]<0 ) s=-1;                     
	    flux[cycle[i]] -= s * 1.001 * minimum; 
	    // inverting the flux closest to equilibrium 
	    cout << "v_" << cycle[i] << "=" << flux[cycle[i]] << endl;
	  }
	}
	for(int i = 0; i <= M-1; i++) chempot[i] = chempot0[i];
	// once the loop is eliminated we restore the initial values 
	// of chemical potential
      }
      else cout << "I need to continue..." << endl;
    }

    // Minover

    int min;
    double xmin = 100000;
    for(int i = 0; i <= N-1; i++){        
      double free = 0;
      for(int j = 0; j <= react[i].size()-1; j++){
	free += react2[i][j] * chempot[react[i][j]];
      }
      double x = -free * sign(flux[i]);
      if(x<xmin){ // calculation of the most unsat constraint
	min = i;
	xmin = x;
      }
    }
    unsat[count] = min; // recording the least unsat
    count++;
    if(xmin >= 0) ok = 1;                 
    else{
      cout << "reaction "<< min << " has flux v_" << min << "=" << flux[min] << " but the free energy is G_" << min << "=" << -xmin * sign(flux[min]) << endl;
      double alpha = 0.1; // alpha is the lenght of the update step 
      for(int j = 0; j <= react[min].size()-1; j++){
	int cip = react[min][j];  
	cout << "from g_" << cip+1 << "=" << chempot[cip] ;
	if(flux[min]>0) chempot[cip] -= alpha * react2[min][j]; // update
	else chempot[cip] += alpha * react2[min][j];       
	cout << " to g_" << cip+1 << "=" << chempot[cip] << endl;
      }
    }
  } while(ok == 0);
  cout << "Chemical potentials: "<< endl;
  for(int j = 0; j <= M-1; j++) cout << chempot[j] << "   "; 
  // writing a thermodynamically feasible chemical potential and flux vector
  cout << endl;
  cout << "Fluxes: "<< endl;
  for(int j = 0; j <= N-1; j++) cout << flux[j] << "   ";
  cout << endl;

  return 0 ;
}