// 08/19/02 add -G parameter to correct estimated alpha calcualation
#include "mpi.h"
#include "ancestral.h"
#include <sys/ddi.h>

double dbrent1(double ax, double bx, double cx, double (*f)(int,int,double), double tol,
        double *xmin,int site,int root);
void dpowell(double p[], double **xi, int n, double ftol, int *iter, double *fret,
	double (*func)(double []));
void worker();
void manager(int numprocs);
int main_function(int argc, char **argv);
double func_oneLength(double *x);


char aa[NP+1]={'-','A','R','N','D','C','Q','E','G','H','I','L','K','M','F','P','S','T','W','Y','V'};
char **listArray;
double *frequency;
int listArrayIndex;
FILE *resultFile=stdout;
int numprocs1;

void usage()
{
	printf("Usage: mpirun -np N /home1/UTSW/wcai/likelihood/likely2_optimize_pei -L listFile -F freqFile\n");
	exit(1);
}
	

int main (int argc, char **argv)
{
	int numprocs, myrank;
	MPI_Init(&argc, &argv);
	MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
	MPI_Comm_rank(MPI_COMM_WORLD, &myrank);

	main_function(argc, argv);

	if(myrank!=0) worker();
	else manager(numprocs);

	MPI_Finalize();
	return 0;
}

void manager(int numprocs)
{
//	int i,j,k;
	int sender, numsent=0;
	int send;
	MPI_Status status;
	int namelen; char processor_name[MPI_MAX_PROCESSOR_NAME];
	int myrank;
	MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
	MPI_Get_processor_name(processor_name, &namelen);
	numprocs1=numprocs;
	printf("numprocs: %d\n",numprocs);

	double ran1(long *idum);
	int i,iter,j,ndim=NP;
	double fret,**xi;
	double *p;
	p=dvector(1,ndim);
	xi=dmatrix(1,ndim,1,ndim);
	for(i=1;i<=NP-1;i++) p[i]=frequency[i]; p[NP]=0.5;
//	for(i=1;i<=ndim;i++) p[i]=0.05;
	for(i=1;i<=ndim;i++) fprintf(resultFile,"%8.4f ",p[i]);
	fprintf(resultFile,"\n");
//	fprintf(resultFile," Val: %12.6f\n",func_oneLength(p));
	fflush(stdout);
	for(i=1;i<=ndim;i++) for(j=1;j<=ndim;j++) xi[i][j]=(i==j? 1.0:0.0);
	dpowell(p,xi,ndim,FTOL,&iter,&fret,func_oneLength);
	fprintf(resultFile,"Iterations: %3d\n\n",iter);
	fprintf(resultFile,"Minimum found at:  ");
	for(i=1;i<=ndim;i++) fprintf(resultFile,"%12.6f",p[i]);
	fprintf(resultFile,"  MinfuncValue = %12.6f ",fret);
//	fprintf(resultFile,"  Check: %12.6f\n\n",func_oneLength(p));
	double *p1, optimizedDist;
	p1=dvector(1,NP);
	p1[NP]=1.0;
	for(i=1;i<NP;i++) { p1[i]=p[i]; p1[NP]-=p1[i];}
	optimizedDist=p[NP];
	//output result
	fprintf(resultFile, "Result from optimization:\n");
	fprintf(resultFile, "Amino Acid    Original    Optimized\n");
	for(i=1;i<=NP;i++) fprintf(resultFile,"    %c    %12.6f    %12.6f\n",aa[i],frequency[i],p1[i]);
	fprintf(resultFile,"optimized distanceRatio: %12.6f\n",optimizedDist);
	fprintf(resultFile,"logL after optimization: %12.6f\n",fret);
	fflush(stdout);
	free_dmatrix(xi,1,NP,1,NP);

	//terminate all worker processes
	for(i=1;i<numprocs;i++)
		MPI_Send(MPI_BOTTOM,0,MPI_INT,i,i,MPI_COMM_WORLD);
}

double func_oneLength(double *x)
{
	int i,j,k;
	int sender, numsent=0;
	int send;
	MPI_Status status;
	int namelen; char processor_name[MPI_MAX_PROCESSOR_NAME];
	int myrank;
	MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
	MPI_Get_processor_name(processor_name, &namelen);

	FILE *output;
	double normalPI=0.0,totalLikely=0.0;
	for(i=1;i<NP;i++) {
		if(x[i]<0.0) return 1.0e100;
		normalPI+=x[i];
	}
	if(normalPI>1.00000001) return 1.0e100;
	if(x[NP]<0.0 || x[NP]>1.0) return 1.0e100;
	if(normalPI>1.0) for(i=1;i<NP;i++) x[i]/=normalPI;
	//write out frequency and ratio to file frequency.txt
	output=fopen("frequency.txt","w");
	if(!output) { printf(" can not open frequency.txt file to write\n"); exit(1); }
	for(i=1;i<=NP;i++) fprintf(output,"%-16.10f\n",x[i]);
	fclose(output);

	//use mpi to calculate logL for all alignment in the list file
	for(i=1;i<min(listArrayIndex+1,numprocs1);i++) {
		send=i;
		MPI_Send(&send,1,MPI_INT,i,i,MPI_COMM_WORLD);
	}
	numsent=min(listArrayIndex+1,numprocs1);

	double receive_logL;
	for(i=1;i<=listArrayIndex;i++){ 
		MPI_Recv(&receive_logL,1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
		sender=status.MPI_SOURCE;
//		int site=status.MPI_TAG;
		totalLikely+=receive_logL;

		if(numsent<=listArrayIndex) {
			send=numsent;
			MPI_Send(&send,1,MPI_INT,sender,numsent,MPI_COMM_WORLD);
			numsent++;
		}
	}
	for(i=1;i<=NP;i++) fprintf(resultFile,"%16.10f ",x[i]);
	fprintf(resultFile," totalLikely: %20.10f\n",totalLikely);
	fflush(resultFile);
	return totalLikely;
}


void worker()
{
	char command[1000], logL_file[100],line[1000];
	FILE *input;
	double send_logL;

//	int i,j,k;
	int myrank;
	MPI_Status status;
	int send;
	int namelen; char processor_name[MPI_MAX_PROCESSOR_NAME];
	MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
	MPI_Get_processor_name(processor_name, &namelen);
	printf("processor: %s rank: %d is running\n",processor_name,myrank);

	if(myrank<=listArrayIndex) {
		MPI_Recv(&send,1,MPI_INT,0,MPI_ANY_TAG,MPI_COMM_WORLD, &status);
		while(status.MPI_TAG>0){
			//use system call to calculate the logL and return to manager
			strcpy(command,"/home1/UTSW/wcai/likelihood/likely2_pei_morpho -L -i ");
			strcat(command,listArray[send]);
			strcat(command,".aln -t ");
			strcat(command,listArray[send]);
			strcat(command,".tre -o junk");
//			printf("%s\n",command);
			system(command);
			//read logL file
			strcpy(logL_file,listArray[send]);
			strcat(logL_file,".aln.logL");
			input=fopen(logL_file,"r");
			if(!input) { printf(" can not open logL file: %s\n",logL_file); exit (1); }
			fgets(line,1000,input);
			send_logL=atof(line);
			fclose(input);
			//delete logL file after read it
			strcpy(command, "rm -f ");
			strcat(command, logL_file);
			system(command);

			int row=send;
			MPI_Send(&send_logL,1,MPI_DOUBLE,0,row,MPI_COMM_WORLD);
			MPI_Recv(&send,1,MPI_INT,0,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
		}
	}

}


int main_function(int argc, char **argv)
{
	int i,j;
	char listFile[100], freqFile[100];
	bool listFlag=false, freqFlag=false;
	FILE *input;
	char line[2000];

	if(argc<=1) usage();
	for(i=1;i<argc;++i){
		char option,optidx;
		if(argv[i][0]=='-') option=argv[i][1];
		else usage();
		switch(option){
			case 'L':
				listFlag=true;
				++i;
				optidx=0;
				if(i<argc) {
					strcpy(listFile,argv[i]);
					input=fopen(listFile,"r");
					if(!input) { printf("can not open file:%s\n",listFile); exit(1); }
					fclose(input);
				}
				else usage();
				break;
			case 'F':
				freqFlag=true;
				++i;
				optidx=0;
				if(i<argc) {
					strcpy(freqFile,argv[i]);
					input=fopen(freqFile,"r");
					if(!input) { printf("can not open file:%s\n",freqFile); exit(1); }
                                        fclose(input);
                                }
                                else usage();
                                break;
			case 'o':
				++i;
				optidx=0;
				if(i<argc) {
					resultFile=fopen(argv[i],"w");
					if(!resultFile) {
						printf("can not open file:%s\n",argv[i]);
						exit(1);
					}
				}
				else usage();
				break;
			default:
				usage();
		}
	}
	if(!listFlag || !freqFlag) usage();

	//read list file
	input=fopen(listFile,"r");
	if(!input) { printf("can not open file: %s\n",listFile); exit(1); }
	listArrayIndex=0;
	while(fgets(line,2000,input)!=NULL){
		bool flag=false;
		for(i=0;i<strlen(line);i++) if(line[i]!=' ') { flag=true; break; }
		if(strlen(line)>1 && flag) {
			++listArrayIndex;
		}
	}
	rewind(input);
	//define listArray
	listArray = new char * [listArrayIndex+1];
	for(i=0;i<=listArrayIndex;i++) listArray[i]=new char [100];
	
	listArrayIndex=0;
	while(fgets(line,2000,input)!=NULL){
		bool flag=false;
		for(i=0;i<strlen(line);i++) if(line[i]!=' ') { flag=true; break; }
		if(strlen(line)>1 && flag) {
			strcpy(listArray[++listArrayIndex],line);
		}
	}
	fclose(input);

	for(i=1;i<=listArrayIndex;i++) for(j=0;j<strlen(listArray[i]);j++) if(listArray[i][j]=='\n') listArray[i][j]='\0';
	for(i=1;i<=listArrayIndex;i++) printf("%s\n",listArray[i]);

	//define frequency vector
	frequency=dvector(1,NP);

	//read frequency
	input=fopen(freqFile,"r");
	if(!input) { printf("can not open file: %s\n",freqFile); exit(1); }
	double totalFreq=0.0;
	for(i=1; i<=NP; i++) {
		fgets(line,2000,input);
		frequency[i]=atof(line);	
		totalFreq+=frequency[i];
	}
//	for(i=1;i<=NP;i++) printf("%f\n",frequency[i]); 

	if(totalFreq>1.001 || totalFreq<0.999) { printf("totalFreq :%f does not equal to 1, please check file: %s\n", totalFreq,freqFile); exit(1); }
	fclose(input);

} 

