/****************************************************************
* Utility program to read in matrix of prediction information   *
* used by the Bailey-Makeham estimation routines.  
* step.                                                         *
****************************************************************/

#include "ds.h"
#include <stdio.h>
#include "inmat.h"

/****************************************************************
* The following globals are used to store information from the  *
* dataset that contains the information from the maximization   *
* step.                                                         *
****************************************************************/

DOUBLE nobs;
int modelid;
Int2 actparm;
Int2 nsp;
Int2 ndv;
Int2 vars;
Int2 parms;
Int2 n_extra_vals;
char *spnames;
char *dvnames;
char *weightname;
char *varnames;
DOUBLE *means;
DOUBLE *extra_vals;
Int2 *actprm;
unsigned char *fixed;
DOUBLE *theta;
DOUBLE *q;
DOUBLE *a;
DOUBLE *adiag;
int error_flag;
int result;
DOUBLE weight;
DOUBLE llr0;
Int2 intype;
char dsn1[64];
char format1[255];
int idcons = (-1);
Int2 lreadv;
Int2 labheads1;
int nread=0;

FILE *f;

/************************************************************
 * Input routines:  Compatible with MATA (c) CRC            *
 ************************************************************/

struct mtadf header;

int inhdr(f) FILE *f; {
	fread(&header,sizeof(struct mtadf),1,f);
	if (header.s1!=MTASIG1 || header.s2!=MTASIG2 || header.s3!=MTASIG3)
		mxerror(23,"Not a valid output file.");
}

#include <stdlib.h>

myalloc() {
	means = (DOUBLE *) malloc(sizeof(DOUBLE)*vars);
	actprm = (Int2 *) malloc(sizeof(Int2)*actparm);
	fixed = (unsigned char *) malloc(parms);
	theta = (DOUBLE *) malloc(sizeof(DOUBLE)*parms);
	spnames = (char *) malloc(nsp*9);
	varnames = (char *) malloc(vars*9);
	weightname = (char *) malloc(9);
	dvnames = (char *) malloc(ndv*9);
	extra_vals = (DOUBLE *) malloc(sizeof(DOUBLE)*n_extra_vals);
	adiag = (DOUBLE *) malloc(sizeof(DOUBLE)*actparm);
	ridbuf = (char *) malloc(lreadv);
	q = (DOUBLE *) malloc(sizeof(DOUBLE)*actparm);
	a = (DOUBLE *) malloc(sizeof(DOUBLE)*(actparm*(actparm+1))/2);
	if (!means || !actprm || !q || !a) mxerror(10,"Can't malloc");
}

void ridpsetup() 
{
	int i; 
	char *p;
	if (nread>0) return; /* a read has already been set up */
	p = ridbuf;
	while (p-ridbuf<lreadv) {
		rvec[nread++] = p;
		/* printf("Setting up %s, %d\n", p, nread-1); */
		p += strlen(p)+1;
	}
}

void byte_order_check(buf,size)
char *buf;
int size;
{
	int i; char c;
	if (header.byteorder!=BYTORD) {
		for (i=0;i<size/2;++i) {
			c = buf[i];
			buf[i]=buf[size-1-i];
			buf[size-1-i] = c;
		}
	}
}

int cread(buf, dsiz, num, f)
void *buf;
unsigned dsiz,num;
FILE *f;
{
	int i, j;
	j = fread(buf,dsiz,num,f);
	for (i=0; i<num; ++i) byte_order_check((char *)buf+i*dsiz,dsiz);
	return j;
}

int inparmf(f)
FILE *f;
{
	inhdr(f);
	cread(&error_flag, sizeof(Int2), 1, f);
	cread(&modelid,sizeof(Int2),1,f);
	cread(&actparm,sizeof(Int2),1,f);
	cread(&nobs,sizeof(DOUBLE),1,f);
	cread(&nsp,sizeof(Int2),1,f);
	cread(&ndv,sizeof(Int2),1,f);
	cread(&vars,sizeof(Int2),1,f);
	cread(&n_extra_vals,sizeof(Int2),1,f);
	cread(&lreadv,sizeof(Int2),1,f);
	cread(&labheads1,sizeof(Int2),1,f);
	parms = vars*nsp;
	myalloc();
	fread(spnames,9,nsp,f);
	fread(dvnames,9,ndv,f);
	fread(varnames,9,vars,f);
	fread(weightname,9,1,f);
	cread(extra_vals,sizeof(DOUBLE),n_extra_vals,f);
	cread(means,sizeof(DOUBLE),vars,f);
	cread(actprm,sizeof(Int2),actparm,f);
	/* for (i=0; i<actparm; ++i) printf("Parameter %d = %d\n", i, actprm[i]); */
	cread(fixed,1,parms,f);
	cread(theta,sizeof(DOUBLE),parms,f);
	cread(q,sizeof(DOUBLE),actparm,f);
	cread(a,sizeof(DOUBLE),(actparm*(actparm+1))/2,f);
	cread(adiag,sizeof(DOUBLE),actparm,f);
	cread(&llr0,sizeof(DOUBLE),1,f);
	cread(&intype,sizeof(Int2),1,f);
	fread(dsn1,64,1,f);
	fread(format1,255,1,f);
	fread(ridbuf,1,lreadv,f);
	if (!labheads1) ridpsetup();
	fclose(f);
	return 0;
}

int inparm(indsn) char *indsn; {
	FILE *f;
	int i, j;
	f = fopen(indsn,"rb");
	if (!f) return 1;
	return inparmf(f);
}

/**********************************************************************/
void indv_setup(need_dv)
int need_dv;
{
	int k, j;
	char *p;
	for (k=j=0; !k && j<vars; ++j) {
		p = &varnames[9*j];
		if (!strcmp("_cons",p)) {
			idcons = j;
			continue;
		}
		k = set_variable(p,&d[j]);
	}
	if (k) mxerror(18,p);
	if (!need_dv) return;
	for (k=j=0; !k && j<ndv; ++j) {
		p = &dvnames[9*j];
		k = set_variable(p,&d[vars+j]);
	}
	if (k) mxerror(19,p);
	if (*weightname) {
		k = set_variable(weightname,&d[vars+ndv]);
		if (k) mxerror(20,weightname);
	}
}