/* SIMPFIT.C */
/*-
simplex fitting routine			= simpfit()
  
J.A. Rupley, Tucson, Arizona
rupley!local@megaron.arizona.edu
*/

#define SIMPFIT

#include "simpdefs.h"

/*-
SIMPFIT
	SIMPLEX MINIMIZATION BY USE OF NELDER-MEAD ALGORITHM,
	FOR MODEL DEFINED IN  <FUNC()> .
*/

#define ILLEGAL		0
#define REFLECT		1
#define STARCONTRACT	2
#define HIGHCONTRACT	3
#define SHRINK		4
#define FAILREFLECT	5
#define STAREXPAND  	6

int
simpfit(fptr)
	FILE           *fptr;
{
	register int    j, i;

	int             opcode;
	int             nlow, nhigh;

	struct pstruct *p_best, pbar, pstar, p2star;

	extern void     qsort();
	extern double   sqrt();

	extern          pvalcomp();
	extern          centroid();
	extern          func();
	extern          ptrial();
	extern void     fpprint();
	extern void     pcopy();
	extern void     fspecial();

	/* expansion-contraction coefficients */
	double          alpha;
	double          beta;
	double          coef_gamma;

	/* expansion-contraction operations */
	char           *opname[7];
	opname[0] = "illegal";
	opname[1] = "reflect";
	opname[2] = "reflect & contract";
	opname[3] = "contract";
	opname[4] = "shrink on lowest vertex";
	opname[5] = "reflect after fail to expand";
	opname[6] = "reflect and expand";

	nlow = 0;
	nhigh = nvert - 1;

	/*
	 * initialization of coefficients and pointers 
	 */
	alpha = 1;
	beta = 0.5;
	coef_gamma = 2;
	for (j = 0; j < nvert; j++)
		p_p[j] = &p[j];

	/* MAIN LOOP OF MINIMIZATION */
	while (++iter) {

		/*
		 * ascending sort of pointers p_p to struct array p 
		 */
		qsort((char *) p_p, (unsigned) nvert,
		      sizeof(*p_p), pvalcomp);

		/*
		 * form reduced simplex, pbar, for (nvert - 1) vertices ie
		 * without high vertex 
		 */
		centroid(&pbar, p_p, (nvert - 1));

		/* form pstar, new reflection trial vertex */
		ptrial(&pstar, &pbar, p_p[nhigh], -alpha);
		func(&pstar);

		/*
		 * NELDER-MEAD ALGORITHM = test trial vertex vs old high and
		 * low vertices; 
		 *
		 * construct new trial vertices as appropriate; 
		 *
		 * set pointer to best new vertex 
		 */
		if (pstar.val < p_p[nlow]->val) {
			ptrial(&p2star, &pstar, &pbar, -coef_gamma);
			func(&p2star);
			if (p2star.val < p_p[nlow]->val) {
				p_best = &p2star;
				opcode = STAREXPAND;
			} else {
				p_best = &pstar;
				opcode = FAILREFLECT;
			}
		} else if (pstar.val <= p_p[nhigh - 1]->val) {
			p_best = &pstar;
			opcode = REFLECT;
		} else {
			if (pstar.val <= p_p[nhigh]->val) {
				pcopy(p_p[nhigh], &pstar);
				opcode = STARCONTRACT;
			} else {
				opcode = HIGHCONTRACT;
			}
			ptrial(&p2star, p_p[nhigh], &pbar, (1 - beta));
			func(&p2star);
			p_best = &p2star;
			if (p2star.val > p_p[nhigh]->val)
				opcode = SHRINK;
		}
		/* END OF NELDER-MEAD */

		/* print results */
		fprintf(fptr, "\n%5d %20.14e %20.14e %s\n",
			iter, p_p[nhigh]->val, p_best->val, opname[opcode]);
		fpprint(fptr, p_p[nhigh]);
		fpprint(fptr, p_best);

		/*
		 * adjust simplex 
		 *
		 * either loop over all vertices > lowest = [0] to shrink on
		 * lowest else copy best movement into high vertex 
		 */
		if (opcode == SHRINK)
			for (j = 1; j < nvert; j++) {
				for (i = 0; i < nparm; i++)
					p_p[j]->parm[i] = (p_p[nlow]->parm[i]
						     + p_p[j]->parm[i]) / 2;
				func(p_p[j]);
			}
		else
			pcopy(p_p[nhigh], p_best);

		/* calculate and print new centroid */
		centroid(&pcent, p_p, nvert);
		fpprint(fptr, &pcent);

		/*
		 * calculate and print mean and rms of function values and
		 * value of exit test 
		 */
		mean_func = 0;
		for (j = 0; j < nvert; j++)
			mean_func = mean_func + p[j].val;
		mean_func = mean_func / nvert;
		rms_func = 0;
		for (j = 0; j < nvert; j++)
			rms_func = rms_func
				+ (p[j].val - mean_func) * (p[j].val - mean_func);
		rms_func = sqrt(rms_func / nvert);
		test = rms_func / mean_func;
		rms_data = sqrt(p_p[nlow]->val / (ndata - nfree));
		fprintf(fptr, "mean=%20.14e rms=%20.14e test=%20.14e\n",
			mean_func, rms_func, test);
		fprintf(fptr,
			"root mean square weighted error of best fit to data = %20.14e\n\n",
			rms_data);
		fspecial(fptr);

		/*
		 * exit test calculate centroid function value if exit 
		 */
		if (iter == maxiter) {
			func(&pcent);
			break;
		}
		if (test < exit_test) {
			func(&pcent);
			if ((pcent.val - mean_func) < (2 * rms_func))
				break;
		}
	}
	/* END OF MAIN LOOP OF MINIMIZATION */

	return (OK);
}				/* END OF SIMPFIT      			 */