#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <iostream.h>
#include <fstream.h>
#include <iomanip.h>
#include <time.h>
#include <string.h>
#include <ctype.h>
#include <malloc.h>
#include "nrutil.h"

void fit(float x[], float y[], int ndata, float sig[], int mwt, float *a,
        float *b, float *siga, float *sigb, float *chi2, float *q);
void usage();

void usage()
{
	printf("Usage: fit2line inputfile\n");
	exit(0);
}

int main(int argc, char **argv)
{
	int i,j,k;
	if(argc<2) usage();
	FILE *input;
	input=fopen(argv[1],"r");
	if(!input) { printf("Can not open file: %s\n",argv[1]); exit(0); }
	float *data;
	float *x, *y, *sig,a,b,siga,sigb,chi2,q;
	int ndata,mwt=0;
	char line[1000];
	int index=0;
	while(fgets(line,1000,input)!=NULL) {
		if(strcmp(line,"")!=0 && (line[0]!='\n')) ++index;
	}
	rewind(input);
	//allocate memory to array x, y
	data=vector(1,index);
	x=vector(1,index);
	y=vector(1,index);
	sig=vector(1,index);

	index=0;
        while(fgets(line,1000,input)!=NULL) {
                if(strcmp(line,"")!=0 && (line[0]!='\n')) {
			++index;
			data[index]=atof(line);
			data[index]=-1.0*log(data[index]);
		}
	}
	fclose(input);
	
	float *variance, *var1, *var2, *A1, *B1, *A2, *B2;
	variance=vector(2,index-1);
	var1=vector(2,index-1);
	var2=vector(2,index-1);
	A1=vector(1,index-1);
	B1=vector(1,index-1);
	A2=vector(1,index-1);
	B2=vector(1,index-1);

	for(i=2;i<=index-1;i++) {
		//prepare to fit the first part
		for(j=1;j<=i;j++) {
			x[j]=j;
			y[j]=data[j];
		}
		ndata=i;
		fit(x,y,ndata,sig,mwt,&a,&b,&siga,&sigb,&chi2,&q);
		var1[i]=chi2; A1[i]=a; B1[i]=b;
		
		for(j=i;j<=index;j++) {
			x[j-i+1]=j;
			y[j-i+1]=data[j];
		}
		ndata=index-i+1;
		fit(x,y,ndata,sig,mwt,&a,&b,&siga,&sigb,&chi2,&q);
		var2[i]=chi2; A2[i]=a; B2[i]=b;
	}
	for(i=2;i<=index-1;i++) variance[i]=var1[i]+var2[i];

	float min=10000;
	int index_min;
	for(i=2;i<=index-1;i++) {
		if(variance[i]<min) { min=variance[i]; index_min=i; }
	}

	//output results
	printf("Data:\n");
	for(i=1;i<=index;i++) printf("%8.3f\n",data[i]);
	printf("  i    variance var1      var2     A1      B1     A2      B2\n");
	for(i=2;i<=index-1;i++)
		printf("%4d  %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n",i,variance[i],var1[i],var2[i],A1[i],B1[i],A2[i],B2[i]);
	printf("\n\nThe biggest variance is: %8.3f, locates at: %4d\n",min,index_min);
}
