/* demo-energy-outline.c
 *
 * Rouben Rostamian
 * March 2015
 *
 * This is a sketch of the contents of the file demo-energy.c
 * from Project 18.2 on page 207.
 *
 * The book's description of that file is somewhat terse.  This
 * file provides a little more of the details to help you along
 * with completing the project.
 */

#include <stdio.h>
#include <math.h>
#include "nelder-mead.h"

#define SQ(x)	((x)*(x))

/* The elastic stored energy function.

   Note: Our stored energy function blows up at x=0, but the
   Nelder-Mead algorithm is not designed to handle singularities.
   To avoid the singularity, we replace W(x) near x=0 by a sufficiently
   large constant.

   The cut-off criterion x^2 < 10^-6 flattens W near
   the origin at the W(10^-3) = 83333.21... level.
*/
static double W(double x)
{
	double eps = 1.0e-3;
	if (x < eps)
		x = eps;
	... calculate and return W(x) according to formula (18.8)
}

/* Energy of a two-link truss with links (0,0)-(1,1) and (1,0)-(1,1)
   and applied force (F1,F2) at (1,1).
*/
static double energy(double *X, int n, void *params)
{
	double x = X[0], y = X[1];
	double *F = params;
	double lambda1 = sqrt(SQ(x)+SQ(y))/sqrt(2);	// stretch of link 1
	double lambda2 = sqrt(SQ(x-1) + SQ(y));		// stretch of link 2
	... calculate and return E(x,y) according to formula (18.9)
	... note that a = F[0], b = F[1].
}

// Case study 1: Unique minimum at D1 = ((0.8319642234, -1.2505278265)
void CaseStudy1(void)
{
	double p[] = { -1.0, 1.0 };
	double F[] = { 0.0, -0.3 };
	struct nelder_mead NM = { 	// C99-style initialization
		.f	= energy,
		.n	= 2,
		.s	= NULL,
		.x	= p,
		.h	= 1.0,
		.tol	= 1.0e-3,
		.maxevals = 1000,
		.params	= F
	};
	int evalcount;

	printf("Case study 1\n");
	printf("Expected answer:   ");
	printf("min = -0.6395281605 at (0.8319642234, -1.2505278265)\n");
	evalcount = nelder_mead(&NM);

	if (evalcount > NM.maxevals)
		printf("No convergence after %d function evaluation\n",
			evalcount);
	else {
		printf("Computed solution: min = %.10f at (%.10f, %.10f)\n",
				NM.minval, p[0], p[1]);
		printf("Converged after %d function evaluations\n", evalcount);
	}
	putchar('\n');
}

// Case study 2: minumum at D1
// min = -0.2047379473 at (0.9208762303, -1.0926128103)
void CaseStudy2a(void)
{
	...
}

// Case study 2: minumum at D2
// min = -0.0055676481 at (1.1554672673,  0.8776603419)
void CaseStudy2b(void)
{
	...
}

int main(void)
{
	CaseStudy1();
	CaseStudy2a();
	CaseStudy2b();
	return 0;
}