/* LEXFIT.C */
/*-
VERSION OF XXXXFIT, for:
  
lsq for kinetic data on pH dependence of the breakdown of
	wt LexA and mutants;
fit to simple titration dependence of k observed;
weighted for equal error in log kobs; 
  
 */
/* MODEL */
/*-
 * simulation of breakdown of LexA:
 *
 * model:
 *
 *                 Ke                    klim
 *      EH+       ---->        E         --->    products
 *
 *
 * v = klim  . (E)
 *
 *   = k . (Etotal)
 *
 *
 * k = observed first order rate constant for appearance of products
 *
 * k = klim / [ 1 + (H+)/{Ke} ]
 *
 */

/*-
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 LOGK_OBS	0	/* defines for elements of struct dat */
#define LOGK_CALC	1	/* (mnemonic) */
#define W		2
#define PH		3
#define KCONF_INDEX	4

#define LOGKLIM(x)	(2 * x)		/* index into set of logklim values */
#define PK(x)		(2 * x + 1)	/* index into set of pK values */

/*-
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;

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

	int             n;
	int             index;
	int             last_index;

	double          H;
	double		logk;
	double		pk;
	double		Ke_1;

	/*
	 * here set/test bounds on parms... if applicable 
	 */
#if 0				/* not needed if fit log constants */
	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\n");
			return (ERROR);
		}
#endif

	/*
	 * set weights and other stuff on first pass through func() 
	 */
	if (pass_one == 0) {
		pass_one = pass_one + 1;
		/*
		 * reset the index values so they run from zero up
		 * by unit increments 
		 */
		last_index = -1;
		index = -1;
		for (n = 0; n < ndata; n++) {
			if ((int) data[n].datval[KCONF_INDEX] != last_index) {
				last_index = (int) data[n].datval[KCONF_INDEX];
				index++;
			}
			data[n].datval[KCONF_INDEX] = (double) index;
		}
	}
	/*-	 reminder -- from math.h
#define	M_LOG10E	0.43429448190325182765
#define	M_LN10		2.30258509299404568402
	*/

	pnam->val = 0;
	last_index = -1;
	for (n = 0; n < ndata; n++) {

		if ((index = (int) data[n].datval[KCONF_INDEX]) !=
		    last_index) {
			last_index = index;
			logk = pnam->parm[LOGKLIM(index)];
			pk = pnam->parm[PK(index)];
			Ke_1 = exp(M_LN10 * pk);
		}
		/* H+ concn from pH */
		H = exp(-(M_LN10 * data[n].datval[PH]));
		/*
		 * calculation of log observed rate constant, logk; based on
		 * model given above; someone should check this! 
		 */
		data[n].datval[LOGK_CALC] = logk -
			(1 / M_LN10) * log(1 + Ke_1 * H);
		/* evaluation of weighted sum of residuals squared */
		pnam->val = pnam->val
			+ (data[n].datval[LOGK_OBS] -
			   data[n].datval[LOGK_CALC])
			* (data[n].datval[LOGK_OBS] -
			   data[n].datval[LOGK_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,
	   "     logkobs logkcal    diff  weight         pH Kconf_index\n");
	/* "1111 2222222 3333333 4444444 5555555 6666666666 7777777"); */
	for (j = 0; j < ndata; j++)
		fprintf(fptr,
			"%4d %7.3f %7.3f %7.3f %7.3g %10.3f %7d\n",
			(j + 1),
			data[j].datval[LOGK_OBS],
			data[j].datval[LOGK_CALC],
			(data[j].datval[LOGK_OBS] -
			 data[j].datval[LOGK_CALC]),
			data[j].datval[W],
			data[j].datval[PH],
			(int) data[j].datval[KCONF_INDEX]);

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

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

void
fpointprint(fptr)
	FILE           *fptr;
{
	return;
}				/* END OF FPOINTPRINT			 */

/*-
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			 */