#include "inmat.h"
#include "predict.h"
#include "ds.h"
#include <math.h>

char outbuf[130];
char dsn[64];
char outdsn[64];
char logfile[64];
char parmdsn[64];

DOUBLE ltab[1000];
char *ident[255];
char format[130];
int debug=0;
int input_type = -1;
int outascii=0;
char quotechar = '\"';
int labheads=0;

short prediction_id[20];
DOUBLE prediction_times[20];
float prediction_sum[30];
int ntimes=0;
unsigned char prediction_p_type[10];
int prediction_p_id[10];
float prediction_p_sum[10];
int npred=0;

static FILE *f;

int main(argc,argv) int argc; char **argv; {
	int i;
	if (!lexopen(argv[1],ltab,(DOUBLE *)&ltab[LTAB])) exit(1);

	KeywordInstall();
	result = yyparse();
	if (result) {
		printf("Parse was rejected\n");
		return 1;
	}
	if (yydebug) printf("Parse was accepted\n");
	if (debug) {
		printf("The dataset is %s\n", dsn);
		printf("The output dataset is %s\n", outdsn);
	}
	if (inparm(parmdsn)) mxerror(21,parmdsn);
	if (input_type==STATAIN) {
		stata_setup(dsn);
	} else if(labheads){
		result=raw_setup(dsn,format);
		if (result) mxerror(10,dsn);
		headread();
	}
	indv_setup(0);
	if (outascii) {
		f=fopen(outdsn,"w");
		if (!f) mxerror(12,outdsn);
		fprintf(f,"   Obs_Num");
		for (i=0; i<ntimes; ++i) fprintf(f,"%12s",ident[prediction_id[i]]);
		for (i=0; i<npred; ++i) fprintf(f,"%12s",ident[prediction_p_id[i]]);
		fprintf(f,"\n");
	} else stata_openout(outdsn);
	/* Times to predict are in prediction_times[0,...,ntimes-1]
	   Names of variables are pointed to by prediction_id[]
	*/
	predict_run();
}

struct stata_header headout;

int stata_openout(dsn) char *dsn; {
	int i;
	short j=0;
	char *pn;

	f = fopen(dsn,"wb");
	if (!f) mxerror(11,dsn);

	headout.release = 103;       /* backward compatible */
	headout.byteorder = BYTORD;
	headout.filetype = 1;        /* raw data */
	headout.nobs = 0L;           /* will need correct value */
	headout.nvar = ntimes+npred;
	strcpy(headout.data_label,"Bailey-Makeham Survival Est.");
	fwrite(&headout,sizeof(struct stata_header),1,f);
	for (i=0; i<ntimes+npred; ++i) fwrite("f",1,1,f);
	for (i=0; i<ntimes; ++i) fwrite(ident[prediction_id[i]],9,1,f);
	for (i=0; i<npred; ++i) fwrite(ident[prediction_p_id[i]],9,1,f);
	for (i=0; i<=ntimes+npred; ++i) fwrite(&j,2,1,f);
	for (i=0; i<ntimes+npred; ++i) fwrite("%9.0g",7,1,f);
	for (i=0; i<ntimes+npred; ++i) fwrite("",9,1,f);
	for (i=0; i<ntimes; ++i) {
		sprintf(outbuf,"Prediction at t==%-12g", prediction_times[i]);
		fwrite(outbuf,32,1,f);
	}
	for (i=0; i<npred; ++i) {
		j = prediction_p_type[i];
		switch(j) {
			case 0: pn="Structural log(alpha)";         break;
			case 1: pn="Structural log(gamma)";         break;
			case 2: pn="Structural log(delta)";         break;
			case 3: pn="Var(log(alpha))";               break;
			case 4: pn="Cov(log(alpha),log(gamma))";    break;
			case 5: pn="Var(log(gamma))";               break;
			case 6: pn="Cov(log(alpha),log(delta))";    break;
			case 7: pn="Cov(log(gamma),log(delta))";    break;
			case 8: pn="Var(log(delta))";               break;
		}
		fwrite(pn,32,1,f);
	}
}

/**********************************************************************/

predict_run() {
	int i, j, k;
	int ii, kk, ik, jj, ij;
	int iii, kkk;
	long nobs=0L;
	long atend;
	DOUBLE sp[NSP], esp[NSP], cov[NSP*(NSP+1)/2], t, p;
	while (1) {
		i = dsread();
		j = missings();
		if (j) {
			if (!outascii) {
				for (k=0; k<ntimes+npred; ++k) prediction_sum[k]=fltmiss;
				while(j--) fwrite(prediction_sum,sizeof(float),ntimes+npred,f);
				++nobs;
			} else {
				while(j--) {
					++nobs ;
					fprintf(f,"%10ld ",nobs);
					for (k=0; k<ntimes+npred; ++k) fprintf(f, "%11s ",misssym);
					fprintf(f,"\n");
				}
			}
		}
		if (i) break;
		++nobs;
		if (idcons>=0) d[idcons] = 1.0;
		for (j=0; j<NSP; ++j) {
			sp[j] = 0.0;
			for (i=0; i<vars; ++i) {
				if (debug & yydebug) printf("nobs = %ld Data = %15.8g, parm = %15.8g\n",
					nobs, d[i], theta[j*vars+i]);
				sp[j] += theta[j*vars+i]*d[i];
			}
		}
		for (j=0; j<NSP*(NSP+1)/2; ++j) cov[j] = 0.0;
		for (i=0; i<actparm; ++i) for (k=0; k<actparm; ++k) {
			ii = actprm[i]%vars;
			iii = actprm[i]/vars;
			kk = actprm[k]%vars;
			kkk = actprm[k]/vars;
			j = (kkk>iii) ? kkk : iii;
			jj = kkk + iii - j;
			ik = (k>i)? k:i;
			/* printf("Updating %d %d %d %d %22.7e %22.7e\n", j*(j+1)/2+jj,
				ij=(ik*(ik-1))/2+i+k, iii, kkk, d[ii], d[kk], a[ij]); */
			cov[j*(j+1)/2+jj] += a[(ik*(ik-1))/2+i+k]*d[ii]*d[kk];
		}
		/* For Bailey-Makeham */
		for (j=0; j<NSP; ++j) {
			if (debug) printf("For nobs=%ld, j=%d, sp[j]=%15g\n", (long)nobs, j, sp[j]);
			esp[j] = exp(sp[j]);
		}
		for (j=0; j<ntimes; ++j) {
			t = prediction_times[j];
			p = esp[1]*t;
			if (p>300.0)p=300.0;
			p = exp(-esp[0]/esp[1]*(1.0-exp(-p)) - esp[2]*t);
			prediction_sum[j] = 1.0-p;
		}
		if (outascii) {
			fprintf(f,"%10ld ",nobs);
			for (j=0; j<ntimes; ++j) fprintf(f,"%11.9g ",prediction_sum[j]);
		} else fwrite(prediction_sum,sizeof(float),ntimes,f);
		for (j=0; j<npred; ++j) {
			i = prediction_p_type[j];
			if (i<3) prediction_p_sum[j] = sp[i];
			else prediction_p_sum[j] = cov[i-3];
		}
		if (outascii) {
			for (j=0; j<npred; ++j) fprintf(f,"%11.9g ", prediction_p_sum[j]);
			fprintf(f,"\n");
		} else fwrite(prediction_p_sum,sizeof(float),npred,f);
	}
	if (!outascii) {
		fseek(f,0L,0);
		headout.nobs = nobs;
		fwrite(&headout,sizeof(struct stata_header),1,f);
	}
	fclose(f);
}

void mxerrors(f,i,s) FILE *f; int i; char *s; {
	fprintf(f,"Error %d: ", i);
	switch(i) {
		case 21: fprintf(f,"Unable to open parameter file %s\n", s);
			break;
		
		default: fprintf(f,"Unknown error %d, text: %s\n", i, s);
		   break;
	}
}

void mxerror(i,s) int i; char *s; {
	mxerrors(stderr,i,s);
	/* mxerrors(logf,i,s); */
	/* fclose(logf); */
	exit(i);
}