#include "iostream"
#include "stdlib.h"
#include "math.h"
#include "time.h"
#include "fstream"
#include "vector"
#include "string"
using namespace std;
                         //INPUT/OUTPUT
#define Stoichiometric_matrix_filename  "matrix.dat"     // stoichiometric matrix filename
#define Nomifile   "metabolites.txt"    // list of metabolites filename
#define N 1668                                         // number of metabolites
#define M 2381                                           // number of reactions

                     // ALGORITHM PARAMETERS
#define max 10               // max number of montecarlo search before starting relaxation                     
#define initT 1;             // initial temperature of the Montecarlo
#define coolrate 1e-3        // cooling rate of the simulated annealing 
#define lambda 1.9           // lenght of the relaxation step  
#define minimum 1e-9         
#define maximum 1e9
#define relaxationmax  1e6  // maximum number of relaxation steps 

                         // VARIABLES DECLARATION
int kerneldim;
int intkerneldim;
std::vector<int> Nsolutions[N];
std::vector<double> Nsolutions2[N];
std::vector<int> matched;
std::vector<int> intmatched;
std::vector<int> J[N];
std::vector<double> J2[N];
double fields[N];
                            // FUNCTIONS
double casual(){
      int m = random();
      double r = double(m)/(RAND_MAX+1.0);                      // random number generator between 0 and 1
      return r;
      }
void quicksort(int k, int km, int orders[], int pivots[]){                               // sorting algorithm (quicksort)
                            int centre;
                            if(k-km>=1){
                               int pivot = km+int((k-km)/2);
                               int l,p;
                               l=0;
                               p=k-km-1;
                               double   neworders[k-km];
                               for(int i=km;i<=k-1;i++){
                                                   if(i != pivot){
                                                                      if(pivots[orders[i]]<pivots[orders[pivot]]){
                                                                               neworders[l]=orders[i];
                                                                               l++;
                                                                              }
                                                                      else{
                                                                               neworders[p]=orders[i];
                                                                               p--;
                                                                           }
                                                              }
                                                   }
                               
                               neworders[p]=orders[pivot];
                               for(int i=km;i<=k-1;i++) orders[i]=neworders[i-km];
                               centre = p+km;
                               quicksort(k,centre+1,orders,pivots);
                               quicksort(centre,km, orders,pivots);
                               }   
                           }
         
void kernel(){                                                                        // kernel calculation by gaussian elimination
        fstream input_file;
        input_file.open(Stoichiometric_matrix_filename, ios::in);                    // read the stoichiometric matrix   rows: reaction   columns: columns 
        int i1=0;
        int j1=0;        
        double pepe;
        std::vector<int> matrix[N];
        std::vector<double> matrix2[N];
        while(input_file >> pepe){
           if(pepe!=0){
                       matrix[j1].push_back(i1);
                       matrix2[j1].push_back(pepe);
                      }
           j1++;
           if(j1==N){
                    j1=0;
                    i1++;
                    }
           }        
        input_file.close();  
        for(int i=0;i<=N-1;i++){
                 matrix[i].push_back(M+i);
                 matrix2[i].push_back(1);
                 }

       int ok=0;
       int order[N];
       for(int i=0;i<=N-1;i++) order[i]=i;
       int pivots[N];
       for(int i=0;i<=N-1;i++){
                            if(matrix[i].size()>0) pivots[i]=matrix[i][0];
                            else pivots[i]=maximum;
                           }
      
        while(ok==0){
        
         quicksort(N,0,order,pivots);
         for(int j=0;j<N-1;j++){
          if(pivots[order[j+1]]==pivots[order[j]] && pivots[order[j]]!=maximum){                                                 // scaled pivoting
             double  min1=100000000;
             if(matrix[order[j]].size()>1) for(int i=0;i<=matrix[order[j]].size()-1;i++) if(fabs(matrix2[order[j]][0]/matrix2[order[j]][i])<min1) min1 = fabs(matrix2[order[j]][0]/matrix2[order[j]][i]);
             double  min2=100000000;
             if(matrix[order[j+1]].size()>1) for(int i=0;i<=matrix[order[j+1]].size()-1;i++)if(fabs(matrix2[order[j+1]][0]/matrix2[order[j+1]][i])<min2) min2 = fabs(matrix2[order[j+1]][0]/matrix2[order[j+1]][i]); 
             if(min2>min1){
                 int k2 = order[j+1];
                 order[j+1]=order[j];
                 order[j] = k2;              
                 }
             }
            }
          ok=1;
          for(int j=0;j<N-1;j++){

                      if(pivots[order[j+1]]==pivots[order[j]] && pivots[order[j]]!=maximum){
                                        int k1 = order[j+1];
                                        int k2 = order[j];
                                        double   colonna[N+M];
                                        for(int i=0;i<=N+M-1;i++) colonna[i]=0;
                                        double    g = matrix2[k2][0]/matrix2[k1][0];
                                        for(int i=1;i<=matrix[k1].size()-1;i++) colonna[matrix[k1][i]]  = matrix2[k1][i]*g;
                                        for(int i=1;i<=matrix[k2].size()-1;i++) colonna[matrix[k2][i]] -= matrix2[k2][i];
                                        matrix[k1].clear();
                                        matrix2[k1].clear();
                                        for(int i=0;i<=N+M-1;i++){
                                                  if(fabs(colonna[i])>minimum){
                                                        matrix[k1].push_back(i);
                                                        matrix2[k1].push_back(colonna[i]);
                                                        }
                                                  }           
                                        ok=0;
                                        if(matrix[order[j+1]].size()>0) pivots[order[j+1]]=matrix[order[j+1]][0];
                                        else pivots[order[j+1]]=maximum;
                                        }
                       }
                  };
   
   std::vector<int> Rsolutions[N];
   std::vector<double> Rsolutions2[N];
   kerneldim=0;
   for(int i=0;i<=N-1;i++){
                   int ok=1;
                   if(matrix[i].size()>0)for(int j=0;j<=matrix[i].size()-1;j++) if(matrix[i][j]<M) ok=0;
                   if(ok==1 && matrix[i].size()>0){
                             for(int j=0;j<=matrix[i].size()-1;j++){
                                                              Rsolutions[kerneldim].push_back(matrix[i][j]-M);
                                                              Rsolutions2[kerneldim].push_back(matrix2[i][j]);
                                                              }
                              kerneldim++;          
                          }
                   }
    
    for(int i=0;i<=N-1;i++){
                       matrix[i].clear();
                       matrix2[i].clear();
                         }
    int i2=0;  
    for(int i=0;i<=kerneldim-1;i++){
                             int ok2=1;
                             if(Rsolutions[i].size()>0) for(int j=0;j<=Rsolutions[i].size()-1;j++){
                                                                          if(Rsolutions2[i][j]*Rsolutions2[i][0]<0) ok2=0;
                                                                          if(matched.size()==0) matched.push_back(Rsolutions[i][j]);
                                                                          else{
                                                                              int ok3=1;  
                                                                              for(int k=0;k<=matched.size()-1;k++) if(matched[k]==Rsolutions[i][j]) ok3=0; 
                                                                              if(ok3==1)  matched.push_back(Rsolutions[i][j]);
                                                                              }
                                                                          
                                                                           }
                             if(ok2==1 && Rsolutions[i].size()>0){
                                        double min=maximum;
                                        for(int j=0;j<=Rsolutions[i].size()-1;j++){
                                                              Nsolutions[i2].push_back(Rsolutions[i][j]);
                                                              Nsolutions2[i2].push_back(fabs(Rsolutions2[i][j]));
                                                              if(min>fabs(Rsolutions2[i][j])) min = fabs(Rsolutions2[i][j]);
                                                              if(intmatched.size()==0) intmatched.push_back(Nsolutions[i2][j]);
                                                              else{
                                                                              int ok3=1;  
                                                                              for(int k=0;k<=intmatched.size()-1;k++) if(intmatched[k]==Nsolutions[i2][j]) ok3=0; 
                                                                              if(ok3==1)  intmatched.push_back(Nsolutions[i2][j]);
                                                                              }
                                                              }
                                        
                                        for(int j=0;j<=Nsolutions[i2].size()-1;j++) Nsolutions2[i2][j]/=min;
                                        i2++;
                                        }
                               }   
    intkerneldim=i2; 
    cout << "Kernel found" << endl;
} 

void fill(){                                                  // interaction matrix construction
         
        fstream input_file;
        input_file.open(Stoichiometric_matrix_filename, ios::in);                    // read the stoichiometric matrix   rows: reactions   columns: metabolites 
        int i1=0;
        int j1=0;    
        int dim=matched.size();     
        double pepe;
        std::vector<int> matrix[dim];
        std::vector<double> matrix2[dim];
        while(input_file >> pepe){
           if(pepe!=0){
                       int prendo=dim;
                       if(dim>0) for(int i=0;i<=dim-1;i++) if(j1==matched[i]) prendo=i;
                       if(prendo<dim){
                            matrix[prendo].push_back(i1);
                            matrix2[prendo].push_back(pepe);
                            }
                      }
           j1++;
           if(j1==N){
                    j1=0;
                    i1++;
                    }
           }        
        input_file.close(); 
       
        for(int i=0;i<=dim-1;i++){                                                             
                               for(int j=i;j<=dim-1;j++){
                                                    double interactions=0;
                                                    if(matrix[i].size()>0)for(int po=0;po<=matrix[i].size()-1;po++){ 
                                                    if(matrix[j].size()>0)     for(int pu=0;pu<=matrix[j].size()-1;pu++) if(matrix[i][po]==matrix[j][pu]) interactions += matrix2[i][po]*matrix2[j][pu];
                                                        }
                                                    if(j==i) fields[i]=interactions;
                                                    else{
                                                             if(fabs(interactions)>minimum){
                                                                   J[i].push_back(j);
                                                                   J2[i].push_back(interactions);
                                                                   J[j].push_back(i);
                                                                   J2[j].push_back(interactions);
                                                                   }
                                                    }
                               
                               }
 }

         
}


int LinearDependence(int vettore[]){                                              // check if the solution found with montecarlo is linearly independent with respect to the previous ones
        int K=intkerneldim+1;
        std::vector<int> matrix[K];
        std::vector<double> matrix2[K];
        
        for(int i=0;i<=K-2;i++){
                   for(int j=0;j<=Nsolutions[i].size()-1;j++){
                             matrix[i].push_back(Nsolutions[i][j]);
                             matrix2[i].push_back(Nsolutions2[i][j]);
                             }
                   }
        int order2[matched.size()];
        int pivots2[matched.size()];
         for(int i=0;i<=matched.size()-1;i++){
                                       order2[i]=i;
                                       pivots2[i]=matched[i];
                                       }
        quicksort(matched.size(),0,order2,pivots2);
        for(int i=0;i<=matched.size()-1;i++) if(vettore[order2[i]]>minimum){
                                                                                     matrix[K-1].push_back(matched[order2[i]]);
                                                                                     matrix2[K-1].push_back(double(vettore[order2[i]]));
                                                             }
                                                   
          
        int ok=0;
        int order[K];
        for(int i=0;i<=K-1;i++) order[i]=i;
        int pivots[K];
        for(int i=0;i<=K-1;i++){
                            if(matrix[i].size()>0) pivots[i]=matrix[i][0];
                            else pivots[i]=maximum;
                           }
      
        while(ok==0){
        
         quicksort(K,0,order,pivots);
         for(int j=0;j<K-1;j++){
          if(pivots[order[j+1]]==pivots[order[j]] && pivots[order[j]]!=maximum){ 
             double  min1=maximum;
             if(matrix[order[j]].size()>1) for(int i=0;i<=matrix[order[j]].size()-1;i++) if(fabs(matrix2[order[j]][0]/matrix2[order[j]][i])<min1) min1 = fabs(matrix2[order[j]][0]/matrix2[order[j]][i]);
             double  min2=maximum;
             if(matrix[order[j+1]].size()>1) for(int i=0;i<=matrix[order[j+1]].size()-1;i++)if(fabs(matrix2[order[j+1]][0]/matrix2[order[j+1]][i])<min2) min2 = fabs(matrix2[order[j+1]][0]/matrix2[order[j+1]][i]); 
             if(min2>min1){
                 int k2 = order[j+1];
                 order[j+1]=order[j];
                 order[j] = k2;              
                 }
             }
            }
          ok=1;
          for(int j=0;j<K-1;j++){

                      if(pivots[order[j+1]]==pivots[order[j]] && pivots[order[j]]!=maximum){
                                        int k1 = order[j+1];
                                        int k2 = order[j];
                                        double   colonna[N];
                                        for(int i=0;i<=N-1;i++) colonna[i]=0;
                                        double    g = matrix2[k2][0]/matrix2[k1][0];
                                        for(int i=1;i<=matrix[k1].size()-1;i++) colonna[matrix[k1][i]]  = matrix2[k1][i]*g;
                                        for(int i=1;i<=matrix[k2].size()-1;i++) colonna[matrix[k2][i]] -= matrix2[k2][i];
                                        matrix[k1].clear();
                                        matrix2[k1].clear();
                                        for(int i=0;i<=N-1;i++){
                                                  if(fabs(colonna[i])>minimum){
                                                        matrix[k1].push_back(i);
                                                        matrix2[k1].push_back(colonna[i]);
                                                        }
                                                  }           
                                        ok=0;
                                        if(matrix[k1].size()>0) pivots[k1]=matrix[k1][0];
                                        else pivots[k1]=maximum;
                                        }
                       }
                  };
            int K1=0;
            for(int i=0;i<=K-1;i++) if(matrix[i].size()>0) K1++;
            int yes=0;
            if(K==K1) yes=1;
            return yes;

}

void Reduce(){                                      // reducing the solution found by montecarlo
        int K = intkerneldim;
        int ok=0;
        int order[K];
        for(int i=0;i<=K-1;i++) order[i]=i;
        int pivots[K];
        for(int i=0;i<=K-1;i++) pivots[i]=-Nsolutions[i].size();
        do{
          quicksort(K,0,order,pivots);
          ok=1;
          for(int i=0;i<K-1;i++){
                   for(int j=i+1;j<=K-1;j++){
                                        int k1 = order[i];
                                        int k2 = order[j];
                                        long double   colonna[N];
                                        for(int i=0;i<=N-1;i++) colonna[i]=0;
                                        int ok1=1;
                                        for(int i=0;i<=Nsolutions[k1].size()-1;i++) colonna[Nsolutions[k1][i]]  = Nsolutions2[k1][i];
                                        for(int i=0;i<=Nsolutions[k2].size()-1;i++){
                                                             colonna[Nsolutions[k2][i]] -= Nsolutions2[k2][i];
                                                             if(colonna[Nsolutions[k2][i]]<-minimum) ok1=0;
                                                            }
                                        if(ok1==1){
                                          ok=0;
                                          Nsolutions[k1].clear();
                                          Nsolutions2[k1].clear();
                                          for(int i=0;i<=N-1;i++){
                                                  if(fabs(colonna[i])>minimum){
                                                        Nsolutions[k1].push_back(i);
                                                        Nsolutions2[k1].push_back(colonna[i]);
                                                        }
                                                  }
                                          pivots[k1]=-Nsolutions[k1].size();
                                          }           
                         }               
                                        
                   }
            }while(ok==0);
}


int Montecarlo(){                                                                               // Montecarlo simulated annealing 
                                       int dim = matched.size();
                                       int num[dim];                    // variables
                                       int numtot=0;                       
                                       for(int i=0;i<=dim-1;i++){
                                                             if(J[i].size()>0)  num[i]=int(2*casual());           // random initialization
                                                             else num[i]=0; 
                                                             numtot += num[i];
                                                             }
                                      
                                       double H=0;
                                       for(int i=0;i<=dim-1;i++){
                                                 H += fields[i]*num[i]*num[i];
                                                 if(J[i].size()>0) for(int j=0;j<=J[i].size()-1;j++) H += J2[i][j]*num[i]*num[J[i][j]];            
                                                 }     
                                         
                                       int stop=0;
                                       int count=0; 
                                       double T1 = initT;
                                       int howmany=0;
                                       double e = exp(-1/T1);     
                                       cout << "Montecarlo start " << endl;   
                                       do{  

                                          int en = int(casual()*dim);
                                          while(J[en].size()==0) en = int(casual()*dim);
                                          int p=1;
                                          if(num[en]>0 && casual()<0.5) p=-1;
                                          double delta = fields[en]*num[en]; 
                                          for(int i=0;i<=J[en].size()-1;i++) delta += J2[en][i]*num[J[en][i]];
                                          delta = 2*p*delta + fields[en];
                                          if(delta<0 || casual()<pow(e,delta)){                         // metropolis rule
                                                       num[en] += p;
                                                       numtot += p;
                                                       H += delta;
                                                    }
                                          count++;
                                          if(count%int(dim)==0){
                                                          T1 -= coolrate;
                                                          if(T1<=0) T1=coolrate;
                                                          e = exp(-1/T1);
                                                          }
                                          if(count==int(double(dim)/coolrate)){ 
                                                               T1=initT;
                                                               e = exp(-1/T1);   
                                                               count=0;
                                                               numtot=0;
                                                               for(int i=0;i<=dim-1;i++) num[i]=0;
                                                               int en = int(casual()*dim);
                                                               while(J[en].size()==0) en = int(casual()*dim);
                                                               num[en]=1;
                                                               numtot=1;
                                                               H=0;
                                                               for(int i=0;i<=dim-1;i++){
                                                                           H += fields[i]*num[i]*num[i];
                                                                           if(J[i].size()>0) for(int j=0;j<=J[i].size()-1;j++) H += J2[i][j]*num[i]*num[J[i][j]];
                                                                           }   
                                                               
                                      
                                                                howmany++;
                                                               }
                                          
                                          if((H<minimum && numtot>0) || howmany==10*max) stop=1;
                                          
                                        }while(stop==0);                                               
                                        int yes;
                                       
                                        if(howmany<10*max){
                                           if(intmatched.size()>0) yes = LinearDependence(num);
                                           else yes=1;  
                                           if(yes==1){
                                                    int order2[matched.size()];
                                                    int pivots2[matched.size()];
                                                    for(int i=0;i<=matched.size()-1;i++){
                                                                   order2[i]=i;
                                                                   pivots2[i]=matched[i];
                                                                     }
                                                    quicksort(matched.size(),0,order2,pivots2);
                                                    for(int i=0;i<=matched.size()-1;i++) if(num[order2[i]]>0){
                                                                                        Nsolutions[intkerneldim].push_back(matched[order2[i]]);
                                                                                        Nsolutions2[intkerneldim].push_back(num[order2[i]]);
                                                                                        }
                                                    intkerneldim++;
                                                    int yes2 = LinearDependence(num);
                                                    Reduce();
                                                    double min=1000;
                                                    for(int i=0;i<=Nsolutions[intkerneldim-1].size()-1;i++){
                                                                                        if(intmatched.size()==0) intmatched.push_back(Nsolutions[intkerneldim-1][i]);
                                                                                        else{
                                                                                              int ok3=1;  
                                                                                              for(int k=0;k<=intmatched.size()-1;k++) if(intmatched[k]==Nsolutions[intkerneldim-1][i]) ok3=0; 
                                                                                              if(ok3==1)  intmatched.push_back(Nsolutions[intkerneldim-1][i]);
                                                                                              
                                                                                               }
                                                                                        if(Nsolutions2[intkerneldim-1][i]<min) min=Nsolutions2[intkerneldim-1][i];
                                                                                        }
                                                   for(int i=0;i<=Nsolutions[intkerneldim-1].size()-1;i++) Nsolutions2[intkerneldim-1][i]/=min;
                                                   
                                                   cout << " I have found a linearly indipendent moiety! Now they are " << intkerneldim  << " engaging " << intmatched.size() << " metabolites" << endl;
                                                   }
                                              else cout << " I have found a moiety but it is linearly dependent!" << endl;
                                        }
                                        else yes=0;
                                        return yes;
}

int Relaxation(){                                                                          // relaxation algorithm: it tests if all the metabolites not classified in conserved moieties are producible
        fstream input_file;
        input_file.open(Stoichiometric_matrix_filename, ios::in);                    // read the stoichiometric matrix   rows: metabolites   columns: reactions 
        int i1=0;
        int j1=0;   
        double pepe;
        int K = intmatched.size();
        std::vector<int> matrix[K];
        std::vector<double> matrix2[K];
        while(input_file >> pepe){
           if(pepe!=0){
                       int prendo=K;
                       if(K>0) for(int i=0;i<=K-1;i++) if(j1==intmatched[i]) prendo=i;                              // reading the stoichiometry of metabolites in conserved moieties: conservation induces flux dependencies             
                       if(prendo<K){
                            matrix[prendo].push_back(i1);
                            matrix2[prendo].push_back(pepe);
                            }
                      }
           j1++;
           if(j1==N){
                    j1=0;
                    i1++;
                    }
           }        
        input_file.close();
        int ok=0;
        int order[K];
        for(int i=0;i<=K-1;i++) order[i]=i;
        int pivots[K];
        for(int i=0;i<=K-1;i++){
                 if(matrix[i].size()>0) pivots[i]=matrix[i][0];
                 else pivots[i]=maximum;
                 }
        while(ok==0){                                                                     // reducing the stoichiometric matrix of conserved moieties to row echelon form by gaussian elimination
            quicksort(K,0,order,pivots);
            for(int j=0;j<K-1;j++){
               if(pivots[order[j+1]]==pivots[order[j]] && pivots[order[j]]!=maximum){ 
               double  min1=maximum;
               if(matrix[order[j]].size()>1) for(int i=0;i<=matrix[order[j]].size()-1;i++) if(fabs(matrix2[order[j]][0]/matrix2[order[j]][i])<min1) min1 = fabs(matrix2[order[j]][0]/matrix2[order[j]][i]);
               double  min2=maximum;
               if(matrix[order[j+1]].size()>1) for(int i=0;i<=matrix[order[j+1]].size()-1;i++)if(fabs(matrix2[order[j+1]][0]/matrix2[order[j+1]][i])<min2) min2 = fabs(matrix2[order[j+1]][0]/matrix2[order[j+1]][i]); 
               if(min2>min1){
                 int k2 = order[j+1];
                 order[j+1]=order[j];
                 order[j] = k2;              
                 }
               }
             }
           ok=1;
           int j=0;
           for(int j=0;j<K-1;j++){
                       if(pivots[order[j+1]]==pivots[order[j]] && pivots[order[j]]!=maximum){
                                        int k1 = order[j+1];
                                        int k2 = order[j];
                                        double   colonna[M];
                                        for(int i=0;i<=M-1;i++) colonna[i]=0;
                                        double    g = matrix2[k2][0]/matrix2[k1][0];
                                        for(int i=1;i<=matrix[k1].size()-1;i++) colonna[matrix[k1][i]]  = matrix2[k1][i]*g;
                                        for(int i=1;i<=matrix[k2].size()-1;i++) colonna[matrix[k2][i]] -= matrix2[k2][i];
                                        matrix[k1].clear();
                                        matrix2[k1].clear();
                                        for(int i=0;i<=M-1;i++){
                                                  if(fabs(colonna[i])>minimum){
                                                        matrix[k1].push_back(i);
                                                        matrix2[k1].push_back(colonna[i]);
                                                        }
                                                  }            
                                        ok=0;
                                        if(matrix[order[j+1]].size()>0) pivots[order[j+1]]=matrix[order[j+1]][0];
                                        else pivots[order[j+1]]=maximum;
                                        }
                       }
           };
       for(int i=0;i<=K-1;i++) {
                            if(matrix[i].size()>0){
                                      double   norm = matrix2[i][0];
                                      for(int j=0;j<=matrix[i].size()-1;j++) matrix2[i][j]/=norm;
                                      }
                          } 
       for(int k1=K-2;k1>=0;k1--){
                           int k = order[k1];
                           if(matrix[k].size()>1){
                              for(int i=1;i<=matrix[k].size()-1;i++){
                                               for(int j1=k1+1;j1<=K-1;j1++){
                                                                int j = order[j1];
                                                                if(matrix[j].size()>0){
                                                                    if(matrix[j][0]==matrix[k][i]){
                                                                                              double   rigak[M];
                                                                                              for(int a=0;a<=M-1;a++) rigak[a]=0;
                                                                                              for(int a=0;a<=matrix[k].size()-1;a++) rigak[matrix[k][a]] = matrix2[k][a];        
                                                                                              for(int a=0;a<=matrix[j].size()-1;a++) rigak[matrix[j][a]]-= matrix2[j][a]*matrix2[k][i];
                                                                                              matrix[k].clear();
                                                                                              matrix2[k].clear(); 
                                                                                              for(int a=0;a<=M-1;a++){
                                                                                                            if(rigak[a]!=0){
                                                                                                                   matrix[k].push_back(a); 
                                                                                                                   matrix2[k].push_back(rigak[a]);
                                                                                                                   }
                                                                                                            }
                                                                                               }
                                                                                       }
                                                                       }           
                                                                   }
                                                 }
                                 }
       
       int indip[M];
       for(int i=0;i<=M-1;i++) indip[i]=K+1;
       for(int i=0;i<=K-1;i++) if(matrix[i].size()>0) indip[matrix[i][0]]=i;
       int M1=0;
       for(int i=0;i<=M-1;i++) if(indip[i]==K+1){ 
                                                  indip[i]=K+M1;
                                                  M1++;
                                                 } 
       std::vector<int> matrixAus[M1];
       std::vector<double> matrixAus2[M1];
       i1=0;
       for(int i=0;i<=M-1;i++){
               if(indip[i]>=K){
                        matrixAus[i1].push_back(i);
                        matrixAus2[i1].push_back(1);
                        i1++;
                        }
               else{
                       int t = indip[i];
                       if(matrix[t].size()>1){
                                     for(int k=1;k<=matrix[t].size()-1;k++){
                                         int quelo = indip[matrix[t][k]]-K;
                                        
                                         matrixAus[quelo].push_back(i);
                                         matrixAus2[quelo].push_back(-matrix2[t][k]);
                                        }
                                    }
                        }
               }            
       for(int i=0;i<=K-1;i++){
                      matrix[i].clear();
                      matrix[i].clear();
                      }  
       int N1=N-K;
       std::vector<int> matrixaus[N1];
       std::vector<double> matrixaus2[N1];
       input_file.open(Stoichiometric_matrix_filename, ios::in);
       int k1=0;                 
       i1=j1=k1=0;
       while(input_file >> pepe){
            int prendo=1;
            if(intmatched.size()>0) for(int i=0;i<=intmatched.size()-1;i++) if(j1==intmatched[i]) prendo--;               // reading the stoichiometry of metabolites NOT in conserved moieties
            if(pepe!=0){
                    if(prendo==1){
                         matrixaus[k1].push_back(i1);
                         matrixaus2[k1].push_back(pepe);
                         }
                    }
           j1++;
           k1+=prendo;
           if(j1==N){
                    j1=k1=0;
                    i1++;
                    }
           }        
       input_file.close();
       std::vector<int> matrixb[N1];
       std::vector<double> matrixb2[N1];
       for(int i=0;i<=M1-1;i++){                                                                 // matrix multiplication
            for(int j=0;j<=N1-1;j++){
                 double prod=0;
                 if(matrixaus[j].size()*matrixAus[i].size()>0){
                                       for(int ib=0;ib<=matrixAus[i].size()-1;ib++) for(int jb=0;jb<=matrixaus[j].size()-1;jb++) if(matrixAus[i][ib]==matrixaus[j][jb]) prod += matrixAus2[i][ib]*matrixaus2[j][jb];
                                       if(fabs(prod)>minimum ){
                                                         matrixb[j].push_back(i);
                                                         matrixb2[j].push_back(prod);
                                                         }
                                       }
                   }
             }
       cout << " Relaxation start " << endl;
       for(int i=0;i<=M1-1;i++){
                      matrixAus[i].clear();
                      matrixAus2[i].clear();
                      }  
       for(int i=0;i<=N1-1;i++){
                      matrixaus[i].clear();
                      matrixaus2[i].clear();
                      }
       double var[M1];
       for(int i=0;i<=M1-1;i++) var[i]=minimum;
       ok=0;
       int tiempo=0;
       int min;
       do{                                                                                   // relaxation algorithm
          double cmin=1000;
          for(int j=0;j<=N1-1;j++){
                         double constr=0;
                         if(matrixb[j].size()>0) for(int i=0;i<=matrixb[j].size()-1;i++) constr += matrixb2[j][i]*var[matrixb[j][i]];
                         if(constr<cmin){  
                                         min = j;
                                         cmin = constr;
                                        }
                                
                         }
          tiempo++;
          if(cmin>=0) ok=1;    // if TRUE constraints satisfied                 
          else{
               double alpha=-1.9*cmin;                                // Motzkin relaxation
               double fact=0;
               for(int j=0;j<=matrixb[min].size()-1;j++) fact += matrixb2[min][j]*matrixb2[min][j];                
               alpha /= fact;
               if(alpha<1e-9*minimum) alpha=1e-9*minimum;
               for(int j=0;j<=matrixb[min].size()-1;j++) var[matrixb[min][j]] += alpha*matrixb2[min][j];              // update
               }
           }while(ok==0 && tiempo<relaxationmax);
       int yes=0;
       if(ok==1) yes=1;
       return yes;
}






int main(){
   srand(time(0));                             // seed for random numbers
   kernel();
   cout << "There are " << kerneldim << " conservation laws, engaging " << matched.size() << " metabolites " << endl;
   cout << intkerneldim << " are integers (conserved moieties), engaging   " << intmatched.size() << " metabolites " << endl;
   fill();
   int finish=0;
   if(intkerneldim==kerneldim) finish=1;
   int timer=0;
   while(finish==0){
                  int yes;
                  yes = Montecarlo();
                  if(intkerneldim==kerneldim) finish=1;
                  if(yes==0) timer++;
                  else timer = 0;
                  if(timer==max){
                                 cout << "I have found too many linearly dependent moieties, let's check if they are all with Relaxation" << endl;
                                 finish = Relaxation();
                                 if(finish==1) cout << "they are all" << endl;
                                 else cout << "probably not, back to Montecarlo" << endl;
                                 timer=0;
                              }
                  }
   Reduce();
   fstream file2;
   file2.open(Nomifile,ios::in);
   string nomi[N];
   for(int i=0;i<=N-1;i++) file2 >> nomi[i];
   file2.close();
   cout << "moiety basis found" << endl;
   cout << "There are " << intkerneldim << " linearly independent conserved moieties, engaging   " << intmatched.size() << " metabolites " << endl;
   if(intkerneldim==kerneldim) cout << " They generate all the conservation laws " << endl;
   else cout << " They don't generate all the conservation laws, " << kerneldim -intkerneldim << " of them are not reducible to moieties"  << endl;
   for(int i=0;i<=intkerneldim-1;i++){
           cout << "moiety number " << i+1 << " engages "  << Nsolutions[i].size() << " metabolites: " << endl;
           if(Nsolutions[i].size()>0) for(int j=0;j<=Nsolutions[i].size()-1;j++) cout << nomi[Nsolutions[i][j]] << "  "<< Nsolutions[i][j] << "  " << Nsolutions2[i][j] << endl;
           }
   cout << "machine time: "<< double(clock())/double(CLOCKS_PER_SEC) << "s" << endl;           // print on standard output the machine time (s)                  

      return 0;
 }