public class amoeba
{

public delegate double function (vector x);

public int d,hi,lo;
public vector[] p;
vector p_ce,zero;
public vector values;
public function f;
public double size;

void update() {
	hi=0; for(int i=1;i<d+1;i++) if(values[i]>values[hi]) hi=i;
	lo=0; for(int i=1;i<d+1;i++) if(values[i]<values[hi]) lo=i;
	p_ce = zero;
	for(int i=0;i<d+1;i++)if(i!=hi) p_ce += p[i];
	p_ce /= d;
	calculate_size();
	System.Console.Write("top = {0:e5} \tsize={1:e5}\n",values[hi],size);
	}

public amoeba(function ofun, vector x, vector ss) {
	f = ofun; d = ss.size;
	values = new vector(d+1);
	p_ce = new vector(d);
	p = new vector[d+1];
	for(int i=0;i<d+1;i++) p[i] = x.copy();
	for(int i=0;i<d;i++) p[i][i]+= ss[i];
	for(int i=0;i<d+1;i++) values[i] = f(p[i]);
	zero = new vector(d); for(int i=0;i<d;i++) zero[i]=0;
	update();
	}

void calculate_size() {
	size=0;
	for(int i=0;i<p.Length;i++)if(i!=lo){
		double n = (p[i]-p[lo]).norm();
		if(n>size) size=n;}
	}

public void downhill(double simplex_size_goal) {
	while(size>simplex_size_goal) {
		vector p_re = p_ce+(p_ce-p[hi]); // try reflection
		double f_re = f(p_re);
		if(f_re<values[lo]){ // try expansion
			vector p_ex = p_ce+2.0*(p_ce-p[hi]);
			double f_ex = f(p_ex);
			if(f_ex<f_re){ // accept expansion
				values[hi]=f_ex; p[hi]=p_ex;
				update(); continue;
				}
			}
		if(f_re<values[hi]){
			values[hi]=f_re; p[hi]=p_re; // accept reflection
			update(); continue;
			}
		vector p_co = p_ce+0.5*(p[hi]-p_ce); // try reflection
		double f_co = f(p_co);
		if(f_co<values[hi]){// accept contraction
				values[hi]=f_co; p[hi]=p_co;
				update(); continue;}
		for(int i=0;i<d+1;i++)if(i!=lo){
			p[i]=0.5*(p[i]+p[lo]);
			values[i]=f(p[i]);}
		update();continue;
		}// end while
	}// end downhill

}