/*************************************************************************
* Report program results and save computations                           *
*   Commissioned by: R. Clifton Bailey, HCFA                             *
*                    Henry Krakauer, HCFA HSQB                           *
*   Based on code written by Jim Summe, NIH                              *
*   Author, Bill Rogers: Rand, Computing Resource Center                 *
*************************************************************************/

#include "bailey.h"
#include "ds.h"
#include <string.h>
#include <math.h>

char outbuf[130];

static DOUBLE last_ll_change;
static DOUBLE lastabsch;

/*************************************************************************
* void iter_report(void);                                                *
*                                                                        *
*  Report the results of one iteration.  This corresponds to the STATA   *
*  option "trace".                                                       *
*************************************************************************/

void iter_report() {

	DOUBLE ll_change;
	int goodsolution=0;
	int i, j;

/*************************************************************************
*  Output the parameter likelihood                                       *
*    Note, the output to "stderr" is designed to keep the user aware     *
*    of progress under conditions where other signs of activity are      *
*    a minimum.                                                          *
*************************************************************************/

	sprintf (outbuf,"Iteration %3d: Log Likelihood = %16.4f (%d)\n", iter, llr0, reasoniter);
	spout21o();
	if (*logfile) fprintf (stderr,outbuf);

/*************************************************************************
*                                                                        *
*    Look critically at convergence.  If quadratic convergence is        *
*    being achieved, then the likelihood should look quadratic about     *
*    the minimum.  This means that the LL function changes by a          *
*    reducing factor on each iteration (e.g >= 2). The corresponding     *
*    stepsizes^2 should be reducing by a similar factor.  In a bad       *
*    solution, the likelihood goes down linearly, but the stepsize       *
*    doesn't go down at all.                                             *
*                                                                        *
*    However, Clift reports that sometimes the criterion is fooled,      *
*    so this needs more work.                                            *
*************************************************************************/

	if (iter>0) {

/*************************************************************************
*    First, report some information on step length, ability to invert,   *
*    and relative change in likelihood.                                  *
*************************************************************************/

		if (trace) {
			dupblnk(15);
			sprintf(outbuf,"Norm Step = %12g\n", sqrt(absch));
			spout21o();
		}
		if (trace & singularity>=0) {
			dupblnk(15);
			sprintf(outbuf,"2nd partial Singularity found for parameter %s.%s\n",
				ident[parm_number[singularity%vars]],
				spname[singularity/vars]);
			spout21o();
		}
		ll_change = (llr0 - llrlast);
		if (trace) {
			dupblnk(15);
			sprintf(outbuf,"Gradient  = %12g\n", ll_change/absch);
			spout21o();
		}

/*************************************************************************
*    Second, evaluate iteration constants to try to get a picture of     *
*    the kind of progress being made.                                    *
*************************************************************************/

		if (iter>1 & ll_change<0.01*nobs/200.) {
			goodsolution = ll_change/last_ll_change < 0.5;
			goodsolution += (ll_change/last_ll_change < 0.1);
			goodsolution += (absch/lastabsch < 0.5);
			goodsolution += (absch/lastabsch < 0.1);
			goodsolution += (fabs(absch-lastabsch)/lastabsch
			  /(fabs(ll_change - last_ll_change)/nobs+1e-16) > 1e3);
						if (goodsolution>=3) {
				if (trace) {
					sprintf (outbuf,"               Predicting interior max(%d), step type %d\n",
						 goodsolution,found_better_answer);
					spout21o();
				}
			} else {
				if (trace) {
					sprintf (outbuf,"               Predicting infinite max(%d), step type %d\n",
					 goodsolution,found_better_answer);
					spout21o();
				}

/*************************************************************************
*    Third, fix (stop) a coefficient as appropriate.  Note:  we elect to *
*    stop a coefficient if the value appears to be proceeding to         *
*    infinity unabated.                                                  *
*************************************************************************/

				if (ll_change<=autofix*nobs/200. && id_largest_change>=0) {
					sprintf(outbuf,"               (parameter %s.%s stopped)\n",
						ident[parm_number[actprm[id_largest_change]%vars]],
						spname[actprm[id_largest_change]/vars]);
					spout21o();
					fix(id_largest_change);
				}
				else if (id_largest_change>=0) {
					sprintf (outbuf,"               Recommend fixing %s.%s\n",
					  ident[parm_number[actprm[id_largest_change]%vars]],
					  spname[actprm[id_largest_change]/vars]);
					spout21o();
				}
				for (i=0; i<actparm; ++i) {
					j = actprm[i];
					if ((pch[i]/(lpch[i]+1e-30))>0.5) {
						sprintf(outbuf,"                 Likely participant: %s.%s\n",
							ident[parm_number[j%vars]],
							spname[j/vars]);
						spout21o();
					}
				}
			}
		}

/*************************************************************************
*    Finally, record some important constants and actually print the     *
*    coefficient table.                                                  *
*************************************************************************/

		lastabsch = absch;
		last_ll_change = ll_change;
	}
	for (i=0; i<actparm; ++i) lpch[i] = pch[i];
	llrlast = llr0;
	if (trace) spout21("\n");
	if (trace) print_parms();
	if (trace) spout21("\n---------------------------------------------------------------------------\n");
	if (insight) outparm("_checkpt.spm");
}

/*************************************************************************
*    print_parms() prints the table of coefficients.                     *
*       The names of the structural parameters are printed across the    *
*       top, and the names of the independent variables are printed      *
*       down the side.  Variables that are fixed are followed with an    *
*       asterisk.                                                        *
*************************************************************************/

void print_parms() {
	int i, j;
	char *spn, *fixp;

	dupblnk(15);
	spout21("Variable |");
	for (i=0; i<NSP; ++i) {
		sprintf(outbuf,"%15s ",spname[i]);
		spout21o();
	}
	spout21("\n");
	dupblnk(15);
	spout21("---------+");
	dupch('-',16*NSP); spout21("\n");
	for (i = (in_constant_pass?(vars-1):0); i<vars; ++i) {
		dupblnk(15);
		sprintf(outbuf,"%8s |",ident[parm_number[i]]);
		spout21o();
		for (j=0; j<NSP; ++j) {
			fixp = fixed[j*vars+i]?"*":" ";
			sprintf(outbuf,"%15g%1s",theta[j*vars+i],fixp);
			spout21o();
		}
		spout21("\n");
	}
}

/*************************************************************************
*    mxerrors(f,i,s) is responsible for printing terminal error messages *
*                                                                        *
*       f is the file pointer                                            *
*       i is the error number                                            *
*       s is an additional information string                            *
*                                                                        *
*    mxerror(i,s) is the routine responsible for routing the messages    *
*    to the log and output files, as appropriate.                        *
*************************************************************************/

void mxerrors(f,i,s) FILE *f; int i; char *s; {
	fprintf(f,"Error %d: ", i);
	switch(i) {
		case 1: fprintf(f,"%s.\n%s\n",
		   "The Marquardt parameter has exceeded 1e17",
		   "Computationally a stationary point has been reached");
		   break;
		case 2: fprintf(f,"%s.\n%s %s.\n%s,\n%s.\n",
		   "Numerical singularity exists in the matrix of 2nd partials",
		   "The singularity is associated with ", s,
		   "If its concomitant variable is independent",
		   "fix or rescale the parameter.");
		   break;
		case 3: fprintf(f,"Unable to open data file: %s\n", s);
		   break;
		case 4: fprintf(f,"The probability of failing in the interval is zero\n");
		   break;
		case 7: fprintf(f,"Problem with STATA data file %s\n", s);
		   break;
		case 14: fprintf(f,"Out of range data values.\n");
		   break;
		case 17: fprintf(f,"Illegal use of keyword %s\n", s);
		   break;
		case 18: fprintf(f,"Independent variable <%s> missing from dataset\n", s);
		   break;
		case 19: fprintf(f,"Weighting variable <%s> missing from dataset\n", s);
		   break;
		case 20: fprintf(f,"Dependent variable <%s> missing from dataset\n", s);
		   break;
		case 21: fprintf(f,"Bad logfile %s.\n", s);

		default: fprintf(f,"Unknown error %d, text: %s\n", i, s);
		   break;
	}
}

void mxerror(i,s) int i; char *s; {
	mxerrors(logf,i,s);
	mxerrors(stderr,i,s);
	fclose(logf);
	error_flag = i;
	outparm(outdsn);
	exit(i);
}

/*************************************************************************
*    Utility routines for output to the active log file, or the screen   *
*    if there is no active log file.                                     *
*************************************************************************/

void spout21(s) char *s; {
	fprintf(logf,s);
}

void spout21o() {
	fprintf(logf,outbuf);
}

void dupch(c,i) char c; int i; {
	char s[2];
	s[1] = 0; s[0] = c;
	while(i--) spout21(s);
}

void dupblnk(i) int i; {
	dupch(' ',i);
}




/*************************************************************************
*    output(a,fa) prints the final report.                               *
*       This routine is only concerned with the appearance of the output *
*************************************************************************/

static double chip, chi2;
static int chidf;

void output(a, fa) DOUBLE *a, *fa; {
	int i, j, k, m, ii, iv, ivv;
	int isfix=0, isstop=0;
	LDOUBL t;

	for (i=0; i<actparm; ++i) adiag[i] = a[(i*(i+3))/2];
	j = matinvert (a, fa, ntheta, actparm);
	if (j>=0) {
		sprintf(outbuf,"Numerical singularity exists, parameter %d\n", j);
		spout21o();
	}
	if (debug) sprt(a,actparm,"inverse of second partials");
	sprintf(outbuf,"%s\n", model_name);
	spout21o();
	if (!onestep) {
		sprintf(outbuf,"Log Likelihood (C) =%14.3f", ll0);
		spout21o();
		chi2 = 2.0 * (llr0 - ll0);
		for (chidf=i=0; i<idcons; ++i) {
			for (j=0;j<NSP;++j) if (!fixed[j*vars+i]) ++chidf;
		}
		chip = chiprob(chi2,(double)chidf);
	} else dupblnk(38);
	sprintf(outbuf,"   Number of observations = %13.0f\n", nobs); spout21o();
	if (!onestep) {
		sprintf(outbuf, "Chi2(%3d)      = %17.3f", chidf, chi2); spout21o();
	}
	sprintf(outbuf,"   Log Likelihood         =%14.3f\n", llr0); spout21o();
	if (!onestep) {
		sprintf(outbuf, "Prob>chi2      = %17.4f\n", chiprob(chi2,(double)chidf));
			spout21o();
	}
	spout21("\n");
	spout21("Structural        |\n");
	spout21("parameter    Var. |        Coef.    Std. Err.      t      Sig.          Mean\n");
	spout21("------------------+----------------------------------------------------------\n");
	for (i=0; i<parms; ++i) {
		k = i/NSP; iv=parm_number[k];
		m = i%NSP;
		ivv = vars*m + k;
		for (ii=0; ii<actparm; ++ii) { if (ivv==actprm[ii]) break; }
		sprintf(outbuf,"%8s%9s | %12g",
			spname[m], ident[iv], theta[ivv]);
		spout21o();
		if(ii>=actparm){
			if (fixed[ivv]==1) { isfix = 1; spout21("           *"); }
			else if (fixed[ivv]==2) { isstop = 1; spout21("          **"); }
			dupblnk(20);
			sprintf(outbuf,"%12g\n",means[k]); spout21o();
			continue;
		} else {
			t = sqrt(a[(ii*(ii+3))/2]);
			sprintf(outbuf,"%12g  %8.3f  %6.4f  %12g\n",
				t, theta[ivv]/t, tprob(theta[ivv]/t,nobs-(LDOUBL)actparm), means[k]);
			spout21o();
		}
	}
	spout21("------------------+"); dupch('-',58); spout21("\n");
	if (isfix) spout21("                                       (*) parameter fixed, SE not estimated.\n");
	if (isstop) spout21("                           (**) paramter going to infinity, SE not estimated.\n");
	return;
}



/*************************************************************************
* outhdr(f) puts a Mata header file for the data to the file "f"         *
*************************************************************************/

int outhdr(f) FILE *f; {
	struct mtadf header;
	header.s1 = MTASIG1; header.s2 = MTASIG2;
	header.s3 = MTASIG3; header.s4 = MTASIG4;
	header.rlse = MTARLSE; header.byteorder = BYTORD;
	header.intsize=4;
	fwrite(&header,sizeof(struct mtadf),1,f);
}

/*************************************************************************
* outval(f,name,value,rows,cols) outputs a matrix of 8-byte IEEE numbers *
*   in Mata format.  See Mata documentation for details.  The idea is to *
*   make the output readable in Mata.                                    *
*                                                                        *
* Currently, this approach is NOT used.                                  *
*************************************************************************/

int outval(f,name,value,rows,cols)
FILE *f;
char *name;
DOUBLE *value;
Int4 rows, cols;
{
	struct dmmkey entry;
	entry.len = (8L*rows)*cols;
	entry.lvl = 0L;
	entry.user1 = rows;
	entry.user2 = cols;
	entry.suser2 = 0L;      /* real, not string */
	strncpy(entry.name,name,8);
	entry.typ = 1; /* data, not function */
	fwrite(&entry,sizeof(struct dmmkey),1,f);
	while (entry.len>32768L) {
		fwrite(&value,1,32768,f);
		entry.len -= 32768L;
		value += 4096;
	}
	fwrite(&entry,1,entry.len,f);
}

/*************************************************************************
* outparm(outdsn) outputs a checkpoint/result file to the named dataset  *
*   outdsn.                                                              *
*                                                                        *
* We must be careful to write constants in such a way that they can be   *
*   transferred across systems!  This means 4 byte ints and 8 byte       *
*   doubles.  Note that byteorder is specified in the header and should  *
*   be respected.                                                        *
*                                                                        *
* If we cannot write the parameters, at least we write a short file with *
*   an error code that STATA could read.                                 *
*************************************************************************/

int outparm(outdsn) char *outdsn; {
	/* Output in meta format */
	FILE *f;
	Int2 i, j;
	f = fopen(outdsn,"wb");
	if (!f) return 1;

	outhdr(f);
	fwrite(&error_flag,sizeof(Int2),1,f);
	fwrite(&modelid,sizeof(Int2),1,f);
	fwrite(&actparm,sizeof(Int2),1,f);
	fwrite(&nobs,sizeof(DOUBLE),1,f);
	i = NSP; fwrite(&i,sizeof(Int2),1,f);
	i = ndv; fwrite(&i,sizeof(Int2),1,f);
	fwrite(&vars,sizeof(Int2),1,f);
	fwrite(&n_extra_vals,sizeof(Int2),1,f);
	i=(ridp-ridbuf); fwrite(&i,sizeof(Int2),1,f);
	fwrite(&labheads,sizeof(Int2),1,f);
	/* End of constants that determine allocation */
	for (i=0; i<NSP; ++i) fwrite(spname[i],9,1,f);
	for (i=0; i<ndv; ++i) fwrite(ident[dvlist[i]],9,1,f);
	for (i=0; i<vars; ++i) fwrite(ident[parm_number[i]],9,1,f);
	fwrite(weightvar?ident[weightvar]:"",9,1,f);
	fwrite(extra_vals,sizeof(DOUBLE),n_extra_vals,f);
	fwrite(means,sizeof(DOUBLE),vars,f);
	fwrite(actprm,sizeof(Int2),actparm,f);
	fwrite(fixed,1,parms,f);
	fwrite(theta,sizeof(DOUBLE),parms,f);
	fwrite(q,sizeof(DOUBLE),actparm,f);
	fwrite(a,sizeof(DOUBLE),(actparm*(actparm+1))/2,f);
	fwrite(adiag,sizeof(DOUBLE),actparm,f);
	fwrite(&llr0,sizeof(DOUBLE),1,f);
	i = input_type; fwrite(&i,sizeof(Int2),1,f);
	fwrite(dsn,64,1,f);
	fwrite(format,255,1,f);
	fwrite(ridbuf,1,ridp-ridbuf,f);
	fclose(f);
	return 0;
}

int out_err(outdsn)
char *outdsn;
{
	FILE *f;
	int i, j;
	f = fopen(outdsn,"wb");
	if (!f) return 1;
	outhdr(f);
	fwrite(&error_flag,sizeof(int),1,f);
	fclose(f);
	return 0;
}