/*
 *  The program takes as input files "s_in.txt" and "s_out.txt", that are the two NxM stoichiometric matrices 
 *  of input and output, where N=59 and M=39 are respectively the reactions and metabolites numbers of the 
 *	Red Blood Cell metabolic network.
 *  The program produces as output the files "h.txt" which contains 10000 fluxes solutions of the Von Neumann problem. 
 *  
 *  Created by Daniele Granata on 10/06/08.
 */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>





typedef unsigned long long int RANDn;                 /* definizione del tipo unsigned long long int */

double gen();
void normalize();
void Consumo(int,int, double);
int Minover(int,int,double);
int trovamin (int);
void caricafile();


RANDn In;

int a[59][39]={0}, b[59][63]={0};
double s[59]={0},c[39]={0};
double s_pr[59]={0},c_pr[39]={0};
double x;



int main()  {

	int N=59,M=39,j,m,t,k,i=0,v=0,passi=0;

	caricafile();

	FILE *f;
	f=fopen("h.txt","w");
	FILE *pippo;
	pippo=fopen("c.txt","w");

	
	double somma=0.,p=0.;
	srand(time(0)); In=rand();            
	In=rand();


	for(k=0;k<10000;k++){		
		p=0.5;
		passi=0;
		//x=1.;
	
		for(m=0;m<N;m++)  {
			s[m]=gen();            
			}
	
		normalize();
		Consumo(N,M,p);
		t=trovamin(M);
	
	
	
		while(passi<1000000&&p<.99981){
		
			if (p<.999) p=p+0.001;
			else p=p+.0001;
		
			passi=Minover(N,M,p);
			
			}
	
		for (j=0;j<M;j++)  fprintf(pippo,"%f ", c_pr[j]); fprintf(pippo,"\n");
		for (i=0;i<N;i++)  fprintf(f,"%f ",s_pr[i]); fprintf(f,"\n\n");
		
		}

		
	fclose(f);
	fclose(pippo);

	}



void caricafile() {                  // input files load

	FILE *s_in;
	FILE *s_out;

	int i,j,N=59,M=39;

	s_in=fopen("s_in.txt","r");
	s_out=fopen("s_out.txt","r");

	for (i=0; i<N; i++) {
		for (j=0; j<M; j++) {
			fscanf(s_in,"%d",&a[i][j]);  
			fscanf(s_out,"%d",&b[i][j]);  
			}

				
			  }
	fclose(s_in);
	fclose(s_out);

	}



double gen() {                              //random numbers generator

	RANDn g = 16807;                           
	RANDn m = 2147483647;
	RANDn numero = 0;
	double x;

	numero= (g*In) % m;
	In = numero;
  
  
	x =  ((double)numero/(double)(m+1));                 
	return x;
}

		

int Minover(int N, int M, double p){         //  central algorithm

	int i,j,m,passi=0;
	
	Consumo(N,M,p);	

	m=trovamin(M);                                 
	
	while (c[m]<=0&&passi<1000000){                    
		
		x=5./(passi+1.);
		
		for (j=0; j<N; j++) {
			if (s[j]+(b[j][m]-p*a[j][m])*x>0) s[j]=s[j]+(b[j][m]-p*a[j][m])*x; 
			}
		
		normalize();			
		Consumo(N,M,p);				

		m=trovamin(M);
		passi++;
		
		}


	if(passi!=1000000) {            
		
		for(i=0;i<N;i++)	s_pr[i]=s[i];
		for(i=0;i<M;i++)	c_pr[i]=c[i];
		
		}

	
	return passi;

	}



	
void normalize() {

	double sum=0.,a=0;
	int j, N=59;
	
	for (j=0; j<N; j++) {
		sum+=s[j];
		}
	
	for (j=0; j<N; j++) {
		s[j]=s[j]/sum*59.;
		}
}


void Consumo(int n, int m, double p) {     
	int i,k;
    double sumout,sumin,j;
	
    for (i=0; i<m; i++) {
	    sumin = 0.0;
		sumout = 0.0;
		
		for (k=0; k<n; k++){
			sumout += b[k][i]*s[k] ;
			sumin += a[k][i]*s[k] ;
			}

		c[i] = sumout-p*sumin;

		}

	}


	

int trovamin (int M) {

	double c_min=c[0];
	int j=0, min=-1;
 
		
	for (j=0; j<M; j++){ 
		if (c[j]<=c_min) {
			c_min=c[j];
			min=j;
			}
		}
	

	return min; 

	}