/* LDHFIT.C */
/*-
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;

J.A. Rupley, Tucson, Arizona
rupley!local@cs.arizona.edu
*/

#define STANDALONE

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

int
func(g)
	struct pstruct *g;
{
	int             n;

	extern struct dat data[];
	extern double   bounds[2][NPARM];
	extern int      nparm;
	extern int      ndata;

	/* bounds checks */
	for (n = 0; n < nparm; n++)
		if (g->parm[n] <= bounds[0][n] || g->parm[n] >= bounds[1][n]) {
			/* if bound violated, set function value HUGE */
			/* and return ERROR */
			g->val = 1.E38;
			if (verbose >= 3)
				fprintf(stderr, "function error\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;
{
	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
	CALLED FROM FDATPRINT
*/

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