/* LDHFIT.C */
/*-
VERSION OF XXXXFIT, for:
  
nonlinear least squares fit by simplex minimizaton (Nelder-Mead algorithm).
  
analysis of lactate dehydrogenase initial rate kinetic data:
	two substrate -- two product mechanism,
	with inhibition by first product,
	and with reactant and product abortive complexes;
  
see DATA_REDN.DOC for description of rate law.
 */
/*-
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@megaron.arizona.edu
*/

#define XXXXFIT

#include "simpdefs.h"		/* externals are read from <simpdefs> */

/* DEFINES SPECIAL TO FITTING FUNCTION */

#define Y	0		/* defines for elements of struct dat */
#define YC	1		/* (mnemonic) */
#define W	2
#define A	3
#define B	4
#define P	5

#define VMAX	0		/* defines for elements of struct pstruct */
#define KMA	1		/* (mnemonic) */
#define KMB	2
#define KMAB	3
#define KMQ_KMPQ 4
#define KIP_1	5
#define KIB_K3K4 6

/*-
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(g)
	struct pstruct *g;
{
	register int    n;

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

	/* bounds checks */
	for (n = 0; n < nparm; n++)
		if (g->parm[n] <= 0) {
			/* if bound violated, set function value HUGE */
			/* and return ERROR */
			g->val = 1.E38;
			fprintf(stderr, "function error: bounds violated\n");
			return (ERROR);
		}
	g->val = 0;
	for (n = 0; n < ndata; n++) {
		/* evaluation of dependent variable */
		data[n].datval[YC] =
			g->parm[VMAX] /
			(
			 1 + g->parm[KMQ_KMPQ] * data[n].datval[P]
			 + g->parm[KMA] / data[n].datval[A]
			 + g->parm[KMB] / data[n].datval[B]
			 * (1 + g->parm[KMQ_KMPQ] * data[n].datval[P])
			 * (1 + g->parm[KIP_1] * data[n].datval[P])
		     + g->parm[KMAB] / data[n].datval[A] / data[n].datval[B]
			 * (1 + g->parm[KMQ_KMPQ] * data[n].datval[P])
			 + g->parm[KIB_K3K4] * data[n].datval[B]
			);
		/* evaluation of weighted sum of residuals squared */
		g->val = g->val + (data[n].datval[Y] - data[n].datval[YC])
			* (data[n].datval[Y] - data[n].datval[YC])
			* 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;
{
	register int    j;

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

	extern void     fpointprint();

	fprintf(fptr, "\1\f\n%-s\n\niteration number %d\n\n", title, iter);
	fprintf(fptr,
		"        yobs   ycalc   del-y  weight       a         b         p\n");
	for (j = 0; j < ndata; j++)
		fprintf(fptr,
			"%4d %7.3f %7.3f %7.3f %7.3f %10.5f %10.5f %10.5f\n",
			(j + 1), data[j].datval[Y], data[j].datval[YC], (data[j].datval[Y] - data[j].datval[YC]),
			data[j].datval[W], data[j].datval[A], data[j].datval[B], data[j].datval[P]);

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

/*-
FPOINTPRINT
	PRINT POINTS FOR CONSTRUCTION OF PRIMARY AND SECONDARY PLOTS
	CODED ACCORDING TO MODEL AND DATA
*/

void
fpointprint(fptr)
	FILE           *fptr;
{

	extern struct pstruct pcent;
	extern char     title[];
	extern int      iter;
	double          v, k1, k2, k3, k4, k5, vint, v1, v2, s1, s2;

	v = pcent.parm[VMAX];
	k1 = pcent.parm[KMA];
	k2 = pcent.parm[KMB];
	k3 = pcent.parm[KMAB];
	k4 = pcent.parm[KMQ_KMPQ];
	k5 = pcent.parm[KIP_1];
	vint = (1 - k1 * k2 / k3) / v;
	v1 = (1 + k1 /.035) / v;
	v2 = v1 + k4 * 150 / v;
	s1 = (k2 + k3 /.035) / v;
	s2 = (s1 + k2 * k5 * 150 / v) * (1 + k4 * 150);

	fprintf(fptr, "\1\f\n%-s\n\niteration number %d\n\n", title, iter);
	fprintf(fptr,
		"For primary plot at fixed b = [Pyruvate]:\n");
	fprintf(fptr,
		"intersection (1/a = 1/[NADH], 1/V)      %10.2f%10.2f\n", -k2 / k3, vint);
	fprintf(fptr,
		"1/V intercepts at b = (.7,.35,.22,.15)      %8.2f%8.2f%8.2f%8.2f\n\n",
		(1 + k2 /.7) / v, (1 + k2 /.35) / v, (1 + k2 /.22) / v, (1 + k2 /.15) / v);

	fprintf(fptr,
		"For primary plot at fixed a = [NADH]:\n");
	fprintf(fptr,
		"intersection (1/b = 1/[Pyruvate], 1/V)  %10.2f%10.2f\n",
		-k1 / k3, vint);
	fprintf(fptr,
		"1/V intercepts at a = (.035,.022,.015,.011) %8.2f%8.2f%8.2f%8.2f\n\n",
		(1 + k1 /.035) / v, (1 + k1 /.022) / v, (1 + k1 /.015) / v, (1 + k1 /.011) / v);

	fprintf(fptr,
		"For secondary plot:\n");
	fprintf(fptr,
		"1/a = 1/[NADH] and 1/b = 1/[Pyruvate] intercepts %9.2f%10.2f\n", -1 / k1, -1 / k2);
	fprintf(fptr,
		"1/V intercept     \t\t\t\t%10.2f\n\n", 1 / v);

	fprintf(fptr, "For lactate inhibition, at a = [NADH] = 0.035\n");
	fprintf(fptr, "1/V values at 1/b = 1/[Pyruvate] = 0 and 5.0\n");
	fprintf(fptr, "\t\t\t- Lactate%7.2f%10.2f\n", v1, (v1 + 5 * s1));
	fprintf(fptr, "\t\t\t+ Lactate%7.2f%10.2f\n\n", v2, (v2 + 5 * s2));

	fprintf(fptr,
	  "R(+/-,intercept), R(+/-,slope)%10.2f%10.2f\n", v2 / v1, s2 / s1);

	fprintf(fptr,
		"\n\n\n\nValues of the kinetic constants\n");
	fprintf(fptr, "parm(0) = V        %17.11e\n", v);
	fprintf(fptr, "parm(1) = KmA      %17.11e\n", k1);
	fprintf(fptr, "parm(2) = KmB      %17.11e\n", k2);
	fprintf(fptr, "parm(3) = KmAB     %17.11e\n", k3);
	fprintf(fptr, "parm(4) = KmQ/KmPQ %17.11e\n", k4);
	fprintf(fptr, "parm(5) = 1/K(I,P) %17.11e\n", k5);

}				/* END OF FPOINTPRINT			 */

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

void
fspecial(fptr)
	FILE           *fptr;
{
	register        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, "    maxiter = %d  quad_test = %7.2e\n",
		maxiter, quad_test);

	return;
}				/* END OF FSPECIAL			 */