/* SIMPFIT.C */
/*-
simplex minimization by use of nelder-mead algorithm,
for model defined in  <func()> .

J.A. Rupley, Tucson, Arizona
rupley!local@cs.arizona.edu
*/
/*-
SIMPLEX METHOD OF FUNCTION MINIMIZATION

The strategy is that of Nelder and Mead (Computer Journal, vol 7, p
308, 1965).

A simplex is constructed in parameter space, with vertices that are
arbitrary except they describe a volume that includes the minimum
point.

Each vertex is associated with a set of parameter values, for which one
calculates the function value for that vertex.  The simplex has (M + 1)
vertices, where M is the number of free parameters.

In successive iterations the vertex with the highest function value is
moved to where it has a smaller function value.  The
movement is with respect to the center of the reduced simplex,
which is the simplex less the highest vertex.  The allowed movements
are:

     reflection;
     expansion beyond the reflected position;
     reflection after failed expansion;
     contraction toward the center of the reduced simplex from
          the original position;
     contraction from the reflected position;
     if none of these operations gives a lower function value,
          then the entire simplex is shrunk about the lowest
          vertex.

Exit from the iterative minimization follows if the fractional RMS
deviation of the function values at the vertices is less than a test
value and if the function value at the centroid of the simplex is
within two RMS deviations of the mean of the values at the vertices.

At exit from the minimization, the centroid of the simplex gives
an arbitrarily close approximation of the parameter values at the
function minimum.

The routine that calculates values for the function to be
optimized must be coded by the user.

The fitting of a model with variable parameters to a set of data is
handled by minimizing an error function, such as the weighted sum of
squares of the differences between observed values and values
calculated according to the model.
*/

#include "simpdefs.h"

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

int
simpfit()
{
	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     pcopy();

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

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

		/*
		 * 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 new centroid */
		centroid(&pcent, p_p, nvert);
		func(&pcent);
		mse = pcent.val / (ndata - nfree);
		rms_data = sqrt(mse);

		/*
		 * calculate 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);
		if (ratio_exit_test && mean_func > 0)
			test = rms_func / mean_func;
		else
			test = rms_func;

		/*
		 * exit test
		 */
		if (iter == maxiter) {
			break;
		}
		if (rms_func == 0.0) {
			prog_done = TRUE;
			break;
		}
		if (test < exit_test) {
			if ((pcent.val - mean_func) < (2 * rms_func)) {
				prog_done = TRUE;
				break;
			}
		}
	}
	/* END OF MAIN LOOP OF MINIMIZATION */

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