////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///********* This code is devceloped by Ms. Mouli Das of Machine Intelligence Unit ****************///
///********* Indian Statistical Institute, Kolkata. The program was implemented in C ****************///
///********* in SunOS Release 5.9 Generic\_117171-07.  The system configurations ****************///
///********* are: Sun System Name: Sun Blade 2000; Platform Name: SUNW, ****************///
///********* Sun-Blade-1000; sparc; Platform Group: sun4u, Physical Memory (RAM) = 2048 MB. ****************///
///********* The code can not be comercialized for bussiness purpose. The code can be modified and use for research only. ****************///
///********* This code is completely documented.****************///
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
#include<stdio.h>
#include <stdlib.h>
#include <math.h>

#define noSample  100  //Total no of samples presented
#define lambda  0.1  //Regularization parameter
#define maxIteration  100  //Total number of iterations to be executed
#define noMetabolite  20  //Total number of metabolites of the stoichiometric matrix
#define noFlux  52  //  Total number of fluxes of the stoichiometric matrix
	

/*Variable declaration*/
int i,j,s;
int **stoichioMatrix; //denotes a two dimensional array 
double **fluxVector; //denotes a two dimensional array 
double constraint,z; //denotes the constraint lambda.(stoichiometric matrix.(weight factor.fluxvector)) and the objective function z
double newc[noFlux], steadyState[noMetabolite], del_c[noFlux], del2_c[noFlux];

//newc[noFlux] denotes the weight factor c 
//steadyState[noMetabolite] denotes the constraint stoichiometric matrix.(c.fluxvector)
//del_c[noFlux] denotes the first derivative of the reformulted objective function y w.r.t c
//del_c[noFlux] denotes the second derivative of the reformulted objective function y w.r.t c

float array_lambda[noMetabolite]; //denotes the storage of the regularization parameter lambda into the array
FILE *fpr,*fpr1,*fpr2;


/*Function declaration*/
void mem_alloc();
void read_data_file();
void initial_c();
void set_lambda();
void compute_constraint();
void calculate_delc();
void calculate_del2c();
double delz_delci(int);



/*Memory allocation for the stoichiometric matrix named "stoichio.txt and flux vector matrix named "fluxvector.txt"*/
void mem_alloc(){
	stoichioMatrix=(int**)malloc(noMetabolite*sizeof(int*));
 	for(i=0;i<noMetabolite;i++){
  		stoichioMatrix[i]=(int*)malloc(noFlux*sizeof(int));
  		if(stoichioMatrix[i]==NULL)
   		printf("Memory allocation failed for stoichioMatrix");
 	}        
	fluxVector=(double**)malloc(noSample*sizeof(double*));
 	for(s=0;s<noSample;s++){
  		fluxVector[s]=(double*)malloc(noFlux*sizeof(double));
  		if(fluxVector[s]==NULL)
   		printf("Memory allocation failed for fluxVector");
 	}
   
}

 
/*Reading and storing of the data file "stoichio.txt" and "fluxvector.txt" into the array*/
void read_data_file(){
	mem_alloc();
	for(i=0;i<noMetabolite;i++){ 
		for(j=0;j<noFlux;j++)
  			fscanf(fpr,"%d",&stoichioMatrix[i][j]);
  		fscanf(fpr,"\n");
  	}
  	for(s=0;s<noSample;s++){
  		for(j=0;j<noFlux;j++)
  			fscanf(fpr1,"%lf",&fluxVector[s][j]);
  		fscanf(fpr1,"\n");
  	}
  	return;
}


/*Initialization of the weight factor c-values*/
void initial_c(){
	double c;
	for(i=0;i<noFlux;i++){
		//c=((double)rand()/RAND_MAX)*0.02;
		//c=((double)rand()/RAND_MAX)*0.04 - 0.02;
            c=((double)rand()/RAND_MAX);
		newc[i] = c;
	}
}


/*storage of the regularization parameter lambda into the array*/
void set_lambda(){
    	for(i=0;i<noMetabolite;i++)
    		array_lambda[i]=lambda;
}


/*Formulation of the constraint lambda.(stoichiometric matrix.(c.fluxVector))*/
void compute_constraint(){
    	for(j=0;j<noMetabolite;j++){
    		steadyState[j] = 0.0;
    		for(i=0;i<noFlux;i++)
    			steadyState[j] += stoichioMatrix[j][i]*newc[i]*fluxVector[s][i];
    		constraint += array_lambda[j]*steadyState[j];
    	}
}


/*Calculation of the derivative z w.r.t c*/
double delz_delci(int ii){
	double result;	
	if((ii==32) || (ii==33))
		result = fluxVector[s][ii];
	else if(ii==34)
		result = 0.0 - fluxVector[s][ii];
	else
		result = 0.0;		
	return result;
}


/*Calculation of the first derivative of y w.r.t c*/
void calculate_delc(){
	for(i=0;i<noFlux;i++){
	 	double constraint_ci = 0.0;
		for(j=0;j<noMetabolite;j++)
			constraint_ci += array_lambda[j]*stoichioMatrix[j][i];
		del_c[i] += -1.0/(z*z)*delz_delci(i) + constraint_ci*fluxVector[s][i];
	}
}


/*Calculation of the second derivative of y w.r.t c*/
void calculate_del2c(){	
	for(i=0;i<noFlux;i++){
		for(j=0;j<noMetabolite;j++)
			del2_c[i] += (2.0/(pow(z,3)))*(pow(delz_delci(i),2));
	}
}
	

/*Main program for estimation of the y-value and c-value*/
int main(){
    	int iteration;
    	double y,total_flux,total_z;
      int seed;

      //y denotes the reformulted objective function
      //total_flux denotes the inverse of the objective function
    	//total_z denotes the total summation of the z value

	printf("\nEnter the initial seed:");
	scanf("%d",&seed);
	printf("\n%d",seed);
	srand(seed);
	
    	fpr = fopen("stoichio.txt","r");
    	if(!fpr)
       	{
       		printf("Error:Cannot open read data file : stoichio.txt\n");
       		exit(0);
       	}
       	
       	fpr1 = fopen("fluxvector.txt","r");
       	if(!fpr1)
       	{
       		printf("Error:Cannot open read data file : fluxvector.txt\n");
       		exit(0);
       	}
       	
       	fpr2 = fopen("Newc.txt","w");
       	if(!fpr2)
       	{
       	printf("Error:Cannot open write data file : Newc.txt\n");
       	exit(0);
       	}
       	fprintf(fpr2,"lambda = %.1f\n", lambda);       	
       	
    	read_data_file();
    	initial_c();      
    	set_lambda();
    	iteration = 0;
    	while(iteration<=maxIteration){
    		total_flux = 0.0;
            total_z = 0.0;
    	    	constraint = 0.0;
    	    	for(s=0;s<noSample;s++){
    			z = (double)((newc[32]*fluxVector[s][32]) + (newc[33]*fluxVector[s][33]) - (newc[34]*fluxVector[s][34]));
			total_flux+=1/z;
	            total_z+=z;
    			compute_constraint();
    		}
    		y = total_flux + constraint;   		
    		fprintf(fpr2,"\niteration = %d \t y = %lf\t z = %lf",iteration,y,total_z/noSample);    		
    		for(i=0;i<noFlux;i++){
    				del_c[i] = 0.0;
    				del2_c[i] = 0.0;
    		}
    		for(s=0;s<noSample;s++){
    			z = (double)((newc[32]*fluxVector[s][32]) + (newc[33]*fluxVector[s][33]) - (newc[34]*fluxVector[s][34]));
    			calculate_delc();
    			calculate_del2c();
    		}
    		for(i=0;i<noFlux;i++){
    			newc[i] += (-(del_c[i]/fabs(del2_c[i])))/noSample; 
    			if(newc[i]<0.0)
    				newc[i] = 0.0;
    			if(newc[i]>1.0)
    				newc[i] = 1.0;   		   		
       			fprintf(fpr2,"\nvalue of newc[%d] = %lf",i+1,newc[i]);
       		}
       		//printf("\nStaying at iteration %d", iteration);
       		iteration++;       		        			
	}
	fclose(fpr);
	fclose(fpr1);
	fclose(fpr2);
	return (0);
}