/* BIPHASIC.C */
/*-
VERSION OF XXXXFIT, for:
  
lsq for kinetic data on the possibly biphasic breakdown of
	wt LexA and mutants;
fit to model described below;
weighted for equal error in concentration (percent) of uncleaved protein;
  
 */
/* MODEL */
/*-
 * 
 * 	    k2
 *         E2 ---> P
 *         ^       ^
 *     k12 |       |
 *         |       |
 *         E1 -----|
 * 	    k1
 * 
 * e1(t)  =  e1(t=0) . exp(-(k1 + k12)t)
 * 
 * e2(t)  =  exp(-k2t) ( e2(t=0) - 
 * 
 * 	e1(t=0) . k12 / (k2 - k1 - k12) . ( 1 - exp((k2 - k1 - k12)t) ) )
 * 
 * etotal(t)   =  e1(t) + e2(t)
 * 
 * p(t)   =  etotal(t=0) - e1(t) - e2(t)
 * 
 */

/*-
contents of the XXXXFIT module  =
declarations and routines for simplex fitting special to function to be fit
(no other functions/files need to be rewritten or adapted for a new model):
  
	mnemonic defines for members of data and parameter structures
  
	function for calculation of dependent variable and
		weighted sum of residuals squared		= func()
  
	print of <data> records					= fdatprint()
  
	additional display called by fdatprint(),
		specially written for this model 		= fpointprint()
  
	customizable display called by <simpfit()>,
		following the output tracking each cycle
		of minimization					= fspecial()
  
J.A. Rupley, Tucson, Arizona
rupley!local@cs.arizona.edu
*/

#define XXXXFIT

#include "simpdefs.h"		/* externals */

/* DEFINES SPECIAL TO FITTING FUNCTION */

#define E_TOTAL_OBS	0	/* defines for elements of struct dat */
#define E_TOTAL_CALC	1	/* (mnemonic) */
#define W		2
#define T		3

#define K1		0	/* defines for elements of struct pstruct */
#define K2		1	/* (mnemonic) */
#define K12		2
#define FRACT_E1_0	3
#define ETOTAL_0	4

/*-
FUNC
	CALCULATION OF LEAST SQUARES FUNCTION
	CODED ACCORDING TO MODEL BEING FIT
	IT SHOULD BE EFFICIENT
	DURING THE FIT, TIME IS MOSTLY SPENT HERE
*/

int
func(pnam)
	struct pstruct *pnam;
{
	static int      pass_one = 0;

	int		n;

	extern char     title[];
	extern struct dat data[];
	extern int      nparm;
	extern int      ndata;

	double		e1, e2, k1, k2, k12, e1_0, e2_0, etotal_0;
	double		c1, c2, c3, t;

	/*
	 * here set/test bounds on parms... if applicable 
	 */
	for (n = 0; n < nparm; n++)
		if (pnam->parm[n] < 0) {
			/* if bound violated, set function value HUGE */
			/* and return ERROR */
/* 			pnam->val = HUGE; */
			pnam->val = 1.E38;
			fprintf(stderr, "function error: parm[%d] = %g\n",
				n, pnam->parm[n]);
			return (ERROR);
		}
		if (pnam->parm[FRACT_E1_0] > 1.0) {
			pnam->val = 1.E38;
			fprintf(stderr, "function error: FRACT_E1_0 = %g\n",
				n, pnam->parm[n]);
			return (ERROR);
		}
		if (pnam->parm[ETOTAL_0] > 1.0) {
			pnam->val = 1.E38;
			fprintf(stderr, "function error: ETOTAL_0 = %g\n",
				n, pnam->parm[n]);
			return (ERROR);
		}

	/*
	 * set weights and other stuff on first pass through func(),
	 * if applicable
	 */

	/*-	 reminder -- from math.h
#define	M_LOG10E	0.43429448190325182765
#define	M_LN10		2.30258509299404568402
	*/

	pnam->val = 0;
	etotal_0 = pnam->parm[ETOTAL_0];
	e1_0 = etotal_0 * pnam->parm[FRACT_E1_0];
	e2_0 = etotal_0 - e1_0;
	k1 = pnam->parm[K1];
	k2 = pnam->parm[K2];
	k12 = pnam->parm[K12];
	c1 = - (k1 + k12);
	c2 = k2 + c1;
	c3 = e1_0 * k12 / c2;
	for (n = 0; n < ndata; n++) {
		/*
		 * calculation of dependent variable
		 * based on model given above
		 */
		t = data[n].datval[T];
		e1 = e1_0 * exp(c1 * t);
		e2 = exp(-k2 * t) * (e2_0 - c3 * (1 - exp(c2 * t)));
		data[n].datval[E_TOTAL_CALC] = e1 + e2;
		/* evaluation of weighted sum of residuals squared */
		pnam->val = pnam->val
			+ (data[n].datval[E_TOTAL_OBS] -
			   data[n].datval[E_TOTAL_CALC])
			* (data[n].datval[E_TOTAL_OBS] -
			   data[n].datval[E_TOTAL_CALC])
			* data[n].datval[W]
			* data[n].datval[W];
	}
	return (OK);
}				/* END OF FUNC        			 */

/*-
FDATPRINT
	PRINT DATA AND COMPARE WITH CALCULATED VALUES
	CODED ACCORDING TO MODEL AND DATA
*/

void
fdatprint(fptr)
	FILE           *fptr;
{
	int             j;

	extern void     fpointprint();

	extern int      iter, ndata, maxiter;
	extern struct dat data[];
	extern char     title[];

	fprintf(fptr, "\1\f\n%-s\n\niteration number %d\n\n", title, iter);
	fprintf(fptr,
	   "     Etotobs Etotcal    diff  weight       time\n");
	/* "1111 2222222 3333333 4444444 5555555 6666666666"); */
	for (j = 0; j < ndata; j++)
		fprintf(fptr,
			"%4d %7.3f %7.3f %7.3f %7.3g %10.3f\n",
			(j + 1),
			data[j].datval[E_TOTAL_OBS],
			data[j].datval[E_TOTAL_CALC],
			(data[j].datval[E_TOTAL_OBS] -
			 data[j].datval[E_TOTAL_CALC]),
			data[j].datval[W],
			data[j].datval[T]);

	if (maxiter != 0 || nvert == 1)
		fpointprint(fptr);
}				/* END OF FDATPRINT			 */

/*-
FPOINTPRINT
	CODE FOR SPECIAL OUTPUT, IF ANY
*/

char           *parm_name[] = {
				"k1",
				"k2",
				"k12",
				"fract_e1_0",
				"etotal_0"
};

void
fpointprint(fptr)
	FILE           *fptr;
{
	extern struct pstruct pcent;
	extern char     title[];
	extern int      iter;
	extern double   rms_data;
	extern int      nfree, ndata;

	extern char    *parm_name[];
	extern char    *std_dev_string();

	/*
	 * temporaries 
	 */
	int             index;

	fprintf(fptr, "\n\f\nSTART SPECIAL OUTPUT\n\n");
	fprintf(fptr, "%s\n\n", title);
	fprintf(fptr, "iteration number %d\n\n", iter);
	fprintf(fptr,
		"number of data points %d\t\trms weighted error %-15.5e\n\n",
		ndata, rms_data);
	for (index = 0; index < nparm; index++)
		fprintf(fptr, "parm %d\t%-15.15s\t%8.5f  %s\n",
			index, parm_name[index], pcent.parm[index],
			std_dev_string(index)
			);

	fprintf(fptr, "\n\n");
	/*
	 * label(s) for use in plotting, or whatever 
	 */
	fprintf(fptr,
		"label k1=%-4.2f k2=%-4.2f k12=%-4.2f f_e1_0=%-4.2f etot_0=%-4.2f",
		pcent.parm[K1],
		pcent.parm[K2],
		pcent.parm[K12],
		pcent.parm[FRACT_E1_0],
		pcent.parm[ETOTAL_0]);
	fprintf(fptr, "\n\n");
	fprintf(fptr, "label1 %s\n", title);
	fprintf(fptr, "\nEND SPECIAL OUTPUT\n\n");
	return;
}				/* END OF FPOINTPRINT			 */


char           *
std_dev_string(i_parm)
	int             i_parm;
{
	extern int      nfree;
	extern struct qstruct q;

	static char     string_store[BUFSIZ];

	int             i;

	for (i = 0; i < nfree; i++) {
		if (i_parm == q.parmndx[i]) {
			sprintf(string_store,
				"+- %8.5f", q.std_dev[i]);
			return string_store;
		}
	}
	return "  fixed";
}


/*-
FSPECIAL
	DISPLAY ADDITIONAL INFORMATION DURING
	TRACKING OF MINIMIZATION, IN SIMPFIT()
*/

void
fspecial(fptr)
	FILE           *fptr;
{
	int             j;

	extern int      maxiter;
	extern double   ypmin, yzero, quad_test;

	if (ypmin > 0.99E38) {
		fprintf(fptr, "y-pmin error");
		for (j = 0; j < nvert; j++)
			if (pmin.parm[j] < 0)
				fprintf(fptr, " (parm(%d) < 0)", j);
	} else if (ypmin > yzero)
		fprintf(fptr, "y-pmin > yzero");

	fprintf(fptr, "    next_prt= %d  next_quad= %7.2e\n",
		maxiter, quad_test);

	return;
}				/* END OF FSPECIAL			 */