/* LEXFIT.C */
/*-
VERSION OF XXXXFIT, for:
  
lsq for kinetic data on pH dependence of the breakdown of
	wt LexA and mutants;
fit to model described below;
weighted for equal error in log kobs; 
  
 */
/* MODEL */
/*-
 * simulation of breakdown of LexA:
 *
 * model:
 *
 *                Kestar                  kstar
 *      EstarH+    ---->        Estar     --->    products
 *
 *        ^                       ^
 *        |                       |
 *        |                       | Kconf
 *        |                       |
 *
 *       EH+       ---->          E
 *                  Ke
 *
 *
 * v = kstar . (Estar)
 *
 *   = k . (Etotal)
 *
 *
 * k = observed first order rate constant for appearance of products
 *
 * k = kstar / [ 1/Kconf + (H+)/{Kconf.Ke} + 1 + (H+)/Kestar ]
 *
 * apparent pK controlling k:
 * 	pK_app = - log( [ 1/Kconf + 1 ] / [ 1/{Kconf.Ke} + 1/Kestar ] )
 *
 * high pH limit of k:
 *	logk_limit = logkstar - log( 1/Kconf + 1 );
 */

/*-
 *
 * modification:
 *
 * assume that Ke and/or kstar vary with Kconf;
 *
 * let pKe and logkstar depend linearly on log(Kconf),
 * i.e., a linear free energy relationship;
 *
 * this is equivalent to:
 *
 * kstar = kstar,ref . (Kconf)^alpha
 * Ke = Ke,ref . (Kconf)^beta
 *
 */

/*-
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 LOGKSTAR_REF	0	/* defines for elements of struct pstruct */
#define PKESTAR		1	/* (mnemonic) */
#define PKE_REF		2
#define ALPHA		3
#define BETA		4
#define LOGKCONF_START	5	/* first of a set of logKconf values */

#define LOGKCONF(x)	(x + LOGKCONF_START)	/* index into set of logKconf
						 * 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          Kestar_1;
	double          Ke_1;
	double          Kconf_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 logKconf 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;
		}
		if (index != (nparm - LOGKCONF_START - 1)) {
			fprintf(stderr,
				"func: pass one; error in Kconf count\n");
			exit(1);/* ERROR in Kconf count */
		}
	}
	/*-	 reminder -- from math.h
#define	M_LOG10E	0.43429448190325182765
#define	M_LN10		2.30258509299404568402
	*/

	pnam->val = 0;
	/* 1/Kestar from pKestar input */
	Kestar_1 = exp(M_LN10 * pnam->parm[PKESTAR]);
	last_index = -1;
	for (n = 0; n < ndata; n++) {
		/* 1/Kconf from set of Kconf values in parms */
		/* and 1/Ke from pKe,ref input and Kconf */
		if ((index = (int) data[n].datval[KCONF_INDEX]) !=
		    last_index) {
			last_index = index;
			Kconf_1 = exp(-(M_LN10 *
					pnam->parm[LOGKCONF(index)]));
			Ke_1 = exp(M_LN10 *
				   (pnam->parm[PKE_REF] -
				    pnam->parm[BETA] *
				    pnam->parm[LOGKCONF(index)]));
		}
		/* 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] = pnam->parm[LOGKSTAR_REF] +
			pnam->parm[ALPHA] * pnam->parm[LOGKCONF(index)] -
			(1 / M_LN10) *
			log(Kconf_1 +
			    Kconf_1 * Ke_1 * H +
			    1 +
			    Kestar_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
*/

#define NUM_KCONF_MAX	10

char           *parm_name[] = {
			       "logk*,ref",
			       "pKe*",
			       "pKe,ref",
			       "alpha",
			       "beta",
			       "logKconf"
};

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();

	double          logKconf[NUM_KCONF_MAX];	/* array of logKconf
							 * values (results of
							 * fitting) */
	double          pKapp[NUM_KCONF_MAX];	/* array of pK_app values
						 * (calculated from results
						 * of fitting) */
	double          logklim[NUM_KCONF_MAX];	/* array of logk_limit values
						 * (high pH limiting rate)
						 * (calculated from results
						 * of fitting) */
	int             num_Kconf;	/* number of data sets with different
					 * logKconf == different mutants */
	/*
	 * temporaries 
	 */
	int             index;
	double          Kestar_1;
	double          Ke_1;
	double          Kconf_1;

	num_Kconf = nparm - LOGKCONF_START;

	/*
	 * fill logKconf[], logklim[], and pKapp[] arrays 
	 */
	Kestar_1 = exp(M_LN10 * pcent.parm[PKESTAR]);
	for (index = 0; index < num_Kconf; index++) {
		Kconf_1 = exp(-(M_LN10 * pcent.parm[LOGKCONF(index)]));
		Ke_1 = exp(M_LN10 *
			   (pcent.parm[PKE_REF] -
			    pcent.parm[BETA] * pcent.parm[LOGKCONF(index)]));
		logKconf[index] = pcent.parm[LOGKCONF(index)];
		logklim[index] = pcent.parm[LOGKSTAR_REF] +
			pcent.parm[ALPHA] * pcent.parm[LOGKCONF(index)] -
			(1 / M_LN10) * log(Kconf_1 + 1);
		pKapp[index] = -(1 / M_LN10) * (log(Kconf_1 + 1) -
					    log(Kconf_1 * Ke_1 + Kestar_1));
	}
	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 < LOGKCONF_START; index++)
		fprintf(fptr, "parm %d\t%-15.15s\t%8.5f  %s\n",
			index, parm_name[index], pcent.parm[index],
			std_dev_string(index)
			);
	for (index = 0; index < num_Kconf; index++)
		fprintf(fptr,
		"* %d\tpKapp %8.5f    logklim %8.5f    logKconf %8.5f  %s\n",
			LOGKCONF(index),
			pKapp[index], logklim[index], logKconf[index],
			std_dev_string(index + LOGKCONF_START)
			);

	fprintf(fptr, "\n\n");
	/*
	 * a label for use in plotting, or whatever 
	 */
	fprintf(fptr,
		"label logk*=%-4.2f pKe*=%-4.2f pKe=%-4.2f A=%-4.2f B=%-4.2f logKconf=",
		pcent.parm[LOGKSTAR_REF],
		pcent.parm[PKESTAR],
		pcent.parm[PKE_REF],
		pcent.parm[ALPHA],
		pcent.parm[BETA]);
	for (index = 0; index < num_Kconf; index++) {
		fprintf(fptr, "%-4.2f ", logKconf[index]);
	}
	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			 */