/****************************************************************
 * Programs for conducting Bailey's Surface Analysis            *
 *                                                              *
 *  This program is placed in the public domain                 *
 ********************************* Bill Rogers *** 04/23/91 ****/

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

char parmdsn[64];
char dsn[64];
char format[256]="";
FILE *f_parmdsn;

DOUBLE sp[20], dv[20];
DOUBLE interval;
int has_interval=0;
int labheads=0;

int debug=0;
FILE *macopen(long *iddir, char *s);

int main(argc,argv) int argc; char **argv; {
	int i;
	long iddir;

#if defined(THINK_C)
	f_parmdsn = macopen(&iddir,"rb");
#else
	if (argv[1]) strcpy(parmdsn,argv[1]);
	f_parmdsn = fopen(parmdsn,"rb");
#endif
	if (!f_parmdsn) {
		printf("Unable to open dataset %s\n", parmdsn);
		exit(1);
	}
	inparmf(f_parmdsn);
	strcpy(dsn,dsn1);
	if (intype==RAW) {
		strcpy(format,format1);
		labheads = labheads1;
		i=raw_setup(dsn,format);
		if (labheads) headread();
	}
	else i=stata_setup(dsn);
	if (i) mxerror(i,dsn);
	indv_setup(1);
	surface_run();
}

/**********************************************************************
* Evaluate likelihood function at +- 1 SE and +-2 SE                  *
*                                                                     *
*   +/- 2 SE is the confidence range for the parameter, allowing      *
*   other parameters to vary.  The likelihood needs to be compared    *
*   with what we expect based on the quadratic approximation          *
**********************************************************************/

#include <stdlib.h>

surface_run() {
	int i, j, k;
	DOUBLE *LL, *se, *lldiff;
	DOUBLE spsave;
	DOUBLE llsave;

	LL = (DOUBLE *) malloc(sizeof(DOUBLE)*parms*4);
	se = (DOUBLE *) malloc(sizeof(DOUBLE)*parms);
	lldiff = (DOUBLE *) malloc(sizeof(DOUBLE)*parms);

	for (j=0; j<4*parms; ++j) LL[j]=0.0;

	/* calculate the standard error distances */
	for (j=0; j<actparm; ++j) {
		i = actprm[j];
		se[i] = sqrt(a[(j*(j+3))/2]);
		lldiff[i] = 0.5*adiag[j]*se[i]*se[i];
	}

	while (1) {

		/* Read in the observations and skip over the missing values */

		i = dsread();
		if (i) break;
		j = missings();
		if (idcons>=0) d[idcons] = 1.0;

		/* Compute the structural parameters */

		for (j=0; j<NSP; ++j) {
			sp[j] = 0.0;
			for (i=0; i<vars; ++i) {
				/* printf("Setting %d %d %g %g\n", i, j, theta[j*vars+i], d[i]); */
				sp[j] += theta[j*vars+i]*d[i];
			}
		}

		/* Set up the dependent variable list and the weight */

		for (j=0; j<ndv; ++j) {
			dv[j] = d[vars+j];
		}
		for (j=0; j<n_extra_vals; ++j) dv[ndv+j] = extra_vals[j];
		if (*weightname) weight = d[vars+ndv]; else weight = 1.0;

		/* Debugging code
		printf("Evaluating with sp: %10.3f %10.3f %10.3f\n",
			sp[0], sp[1], sp[2]);
		printf(" and with dv: %10.3f %10.3f %10.3f\n",
			dv[0], dv[1], dv[2]);
		printf(" and with weight: %10.3f\n", weight);
		printf("extra_vals[0] = %10.3f\n", extra_vals[0]);
		*/

		for (j=0; j<NSP; ++j) {
			for (i=0; i<vars; ++i) {
				if (fixed[j*vars+i]) continue;
				spsave = sp[j];
				sp[j] += se[j*vars+i]*d[i];
				LL[(j*vars+i)*4+2] += (llsave=LL_compute(dv,sp,(DOUBLE *)NULL,
					(DOUBLE *)NULL,(DOUBLE *)NULL, &weight,0));
				if (debug>=2) printf("llsave,2 = %10.4f\n", llsave);
				sp[j] += se[j*vars+i]*d[i];
				LL[(j*vars+i)*4+3] += (llsave=LL_compute(dv,sp,(DOUBLE *)NULL,
					(DOUBLE *)NULL,(DOUBLE *)NULL, &weight,0));
				if (debug>=2) printf("llsave,3 = %10.4f\n", llsave);
				sp[j] = spsave - se[j*vars+i]*d[i];
				LL[(j*vars+i)*4+1] += LL_compute(dv,sp,(DOUBLE *)NULL,
					(DOUBLE *)NULL,(DOUBLE *)NULL, &weight,0);
				sp[j] -= se[j*vars+i]*d[i];
				LL[(j*vars+i)*4+0] += LL_compute(dv,sp,(DOUBLE *)NULL,
					(DOUBLE *)NULL,(DOUBLE *)NULL, &weight,0);
				sp[j] = spsave;
			}
		}
	}
	printf("\n");
	printf("Structural        |\n");
	printf("parameter    Var. |          -2 SE         -1 SE         +1 SE         +2 SE\n");
	printf("------------------+----------------------------------------------------------\n");
	for (i=0; i<vars; ++i) {
		for (j=0; j<NSP; ++j) {
			printf("%8s%9s | ",&spnames[j*9],&varnames[i*9]);
			if (!fixed[j*vars+i]) {
				for (k=0;k<4;++k) {
					printf("%8.2f",llsave=llr0-LL[(j*vars+i)*4+k]);
					if (k==0 || k==3) llsave = llsave/4.0;
					printf("(%4.1f)",llsave/lldiff[j*vars+i]);
				}
			}
			else {
				if (fixed[j*vars+i]==1) printf("   fixed");
				if (fixed[j*vars+i]==2) printf(" stopped");
			}
		printf("\n");
		}
	}
	printf("------------------+"); dupch('-',58); printf("\n");
}

void mxerrors(f,i,s) FILE *f; int i; char *s; {
	fprintf(f,"Error %d: ", i);
	switch(i) {
		case 3: fprintf(f,"Unable to open data file: %s\n", s);
		   break;
		case 18: fprintf(f,"Independent variable <%s> missing from dataset\n", s);
		   break;
		case 19: fprintf(f,"Dependent variable <%s> missing from dataset\n", s);
		   break;
		case 20: fprintf(f,"Weighting variable <%s> missing from dataset\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);
}

int dupch(c,n)
char c; int n;
{
	while(n--) putchar(c);
}

char outbuf[132];

int spout21o() {
	char *p;
	for (p=outbuf;*p;++p) putchar(*p);
}