all the .h files sophie hopefully pulled off adding

ugh

adapt.h

@@ -0,0 +1,86 @@
+struct Adapt {
+	Doub TOL,toler;
+	static const Doub alpha,beta,x1,x2,x3,x[12];
+	bool terminate,out_of_tolerance;
+	Adapt(Doub tol);
+	template <class T>
+	Doub integrate(T &func, const Doub a, const Doub b);
+	template <class T>
+	Doub adaptlob(T &func, const Doub a, const Doub b, const Doub fa,
+		const Doub fb, const Doub is);
+};
+
+Adapt::Adapt(Doub tol) : TOL(tol),terminate(true),out_of_tolerance(false)
+{
+	const Doub EPS=numeric_limits<Doub>::epsilon();
+	if (TOL < 10.0*EPS)
+		TOL=10.0*EPS;
+}
+
+template <class T>
+Doub Adapt::integrate(T &func, const Doub a, const Doub b)
+{
+	Doub m,h,fa,fb,i1,i2,is,erri1,erri2,r,y[13];
+	m=0.5*(a+b);
+	h=0.5*(b-a);
+	fa=y[0]=func(a);
+	fb=y[12]=func(b);
+	for (Int i=1;i<12;i++)
+		y[i]=func(m+x[i]*h);
+	i2=(h/6.0)*(y[0]+y[12]+5.0*(y[4]+y[8]));
+	i1=(h/1470.0)*(77.0*(y[0]+y[12])+432.0*(y[2]+y[10])+
+		625.0*(y[4]+y[8])+672.0*y[6]);
+	is=h*(0.0158271919734802*(y[0]+y[12])+0.0942738402188500*
+		(y[1]+y[11])+0.155071987336585*(y[2]+y[10])+
+		0.188821573960182*(y[3]+y[9])+0.199773405226859*
+		(y[4]+y[8])+0.224926465333340*(y[5]+y[7])+
+		0.242611071901408*y[6]);
+	erri1=abs(i1-is);
+	erri2=abs(i2-is);
+	r=(erri2 != 0.0) ? erri1/erri2 : 1.0;
+	toler=(r > 0.0 && r < 1.0) ? TOL/r : TOL;
+	if (is == 0.0)
+		is=b-a;
+	is=abs(is);
+	return adaptlob(func,a,b,fa,fb,is);
+}
+
+template <class T>
+Doub Adapt::adaptlob(T &func, const Doub a, const Doub b, const Doub fa,
+		const Doub fb, const Doub is)
+{
+	Doub m,h,mll,ml,mr,mrr,fmll,fml,fm,fmrr,fmr,i1,i2;
+	m=0.5*(a+b);
+	h=0.5*(b-a);
+	mll=m-alpha*h;
+	ml=m-beta*h;
+	mr=m+beta*h;
+	mrr=m+alpha*h;
+	fmll=func(mll);
+	fml=func(ml);
+	fm=func(m);
+	fmr=func(mr);
+	fmrr=func(mrr);
+	i2=h/6.0*(fa+fb+5.0*(fml+fmr));
+	i1=h/1470.0*(77.0*(fa+fb)+432.0*(fmll+fmrr)+625.0*(fml+fmr)+672.0*fm);
+	if (abs(i1-i2) <= toler*is || mll <= a || b <= mrr) {
+		if ((mll <= a || b <= mrr) && terminate) {
+			out_of_tolerance=true;
+			terminate=false;
+		}
+		return i1;
+	}
+	else
+		return adaptlob(func,a,mll,fa,fmll,is)+
+			adaptlob(func,mll,ml,fmll,fml,is)+
+			adaptlob(func,ml,m,fml,fm,is)+
+			adaptlob(func,m,mr,fm,fmr,is)+
+			adaptlob(func,mr,mrr,fmr,fmrr,is)+
+			adaptlob(func,mrr,b,fmrr,fb,is);
+}
+const Doub Adapt::alpha=sqrt(2.0/3.0);
+const Doub Adapt::beta=1.0/sqrt(5.0);
+const Doub Adapt::x1=0.942882415695480;
+const Doub Adapt::x2=0.641853342345781;
+const Doub Adapt::x3=0.236383199662150;
+const Doub Adapt::x[12]={0,-x1,-alpha,-x2,-beta,-x3,0.0,x3,beta,x2,alpha,x1};

amebsa.h

@@ -0,0 +1,147 @@
+template <class T>
+struct Amebsa {
+	T &funk;
+	const Doub ftol;
+	Ranq1 ran;
+	Doub yb;
+	Int ndim;
+	VecDoub pb;
+	Int mpts;
+	VecDoub y;
+	MatDoub p;
+	Doub tt;
+	Amebsa(VecDoub_I &point, const Doub del, T &funkk, const Doub ftoll) :
+		funk(funkk), ftol(ftoll), ran(1234),
+		yb(numeric_limits<Doub>::max()), ndim(point.size()), pb(ndim),
+		mpts(ndim+1), y(mpts), p(mpts,ndim)
+	{
+		for (Int i=0;i<mpts;i++) {
+			for (Int j=0;j<ndim;j++)
+				p[i][j]=point[j];
+			if (i != 0) p[i][i-1] += del;
+		}
+		inity();
+	}
+	Amebsa(VecDoub_I &point, VecDoub_I &dels, T &funkk, const Doub ftoll) :
+		funk(funkk), ftol(ftoll), ran(1234),
+		yb(numeric_limits<Doub>::max()), ndim(point.size()), pb(ndim),
+		mpts(ndim+1), y(mpts), p(mpts,ndim)
+	{
+		for (Int i=0;i<mpts;i++) {
+			for (Int j=0;j<ndim;j++)
+				p[i][j]=point[j];
+			if (i != 0) p[i][i-1] += dels[i-1];
+		}
+		inity();
+	}
+	Amebsa(MatDoub_I &pp, T &funkk, const Doub ftoll) : funk(funkk),
+		ftol(ftoll), ran(1234), yb(numeric_limits<Doub>::max()),
+		ndim(pp.ncols()), pb(ndim), mpts(pp.nrows()), y(mpts), p(pp)
+	{ inity(); }
+
+	void inity() {
+		VecDoub x(ndim);
+		for (Int i=0;i<mpts;i++) {
+			for (Int j=0;j<ndim;j++)
+				x[j]=p[i][j];
+			y[i]=funk(x);
+		}
+	}
+	Bool anneal(Int &iter, const Doub temperature)
+	{
+		VecDoub psum(ndim);
+		tt = -temperature;
+		get_psum(p,psum);
+		for (;;) {
+			Int ilo=0;
+			Int ihi=1;
+			Doub ylo=y[0]+tt*log(ran.doub());
+			Doub ynhi=ylo;
+			Doub yhi=y[1]+tt*log(ran.doub());
+			if (ylo > yhi) {
+				ihi=0;
+				ilo=1;
+				ynhi=yhi;
+				yhi=ylo;
+				ylo=ynhi;
+			}
+			for (Int i=3;i<=mpts;i++) {
+				Doub yt=y[i-1]+tt*log(ran.doub());
+				if (yt <= ylo) {
+					ilo=i-1;
+					ylo=yt;
+				}
+				if (yt > yhi) {
+					ynhi=yhi;
+					ihi=i-1;
+					yhi=yt;
+				} else if (yt > ynhi) {
+					ynhi=yt;
+				}
+			}
+			Doub rtol=2.0*abs(yhi-ylo)/(abs(yhi)+abs(ylo));
+			if (rtol < ftol || iter < 0) {
+				SWAP(y[0],y[ilo]);
+				for (Int n=0;n<ndim;n++)
+					SWAP(p[0][n],p[ilo][n]);
+				if (rtol < ftol)
+					return true;
+				else
+					return false;
+			}
+			iter -= 2;
+			Doub ytry=amotsa(p,y,psum,ihi,yhi,-1.0);
+			if (ytry <= ylo) {
+				ytry=amotsa(p,y,psum,ihi,yhi,2.0);
+			} else if (ytry >= ynhi) {
+				Doub ysave=yhi;
+				ytry=amotsa(p,y,psum,ihi,yhi,0.5);
+				if (ytry >= ysave) {
+					for (Int i=0;i<mpts;i++) {
+						if (i != ilo) {
+							for (Int j=0;j<ndim;j++) {
+								psum[j]=0.5*(p[i][j]+p[ilo][j]);
+								p[i][j]=psum[j];
+							}
+							y[i]=funk(psum);
+						}
+					}
+					iter -= ndim;
+					get_psum(p,psum);
+				}
+			} else ++iter;
+		}
+	}
+	inline void get_psum(MatDoub_I &p, VecDoub_O &psum)
+	{
+		for (Int n=0;n<ndim;n++) {
+			Doub sum=0.0;
+			for (Int m=0;m<mpts;m++) sum += p[m][n];
+			psum[n]=sum;
+		}
+	}
+	Doub amotsa(MatDoub_IO &p, VecDoub_O &y, VecDoub_IO &psum,
+		const Int ihi, Doub &yhi, const Doub fac)
+	{
+		VecDoub ptry(ndim);
+		Doub fac1=(1.0-fac)/ndim;
+		Doub fac2=fac1-fac;
+		for (Int j=0;j<ndim;j++)
+			ptry[j]=psum[j]*fac1-p[ihi][j]*fac2;
+		Doub ytry=funk(ptry);
+		if (ytry <= yb) {
+			for (Int j=0;j<ndim;j++) pb[j]=ptry[j];
+			yb=ytry;
+		}
+		Doub yflu=ytry-tt*log(ran.doub());
+		if (yflu < yhi) {
+			y[ihi]=ytry;
+			yhi=yflu;
+			for (Int j=0;j<ndim;j++) {
+				psum[j] += ptry[j]-p[ihi][j];
+				p[ihi][j]=ptry[j];
+			}
+		}
+		return yflu;
+	}
+};

amoeba.h

@@ -0,0 +1,116 @@
+struct Amoeba {
+	const Doub ftol;
+	Int nfunc;
+	Int mpts;
+	Int ndim;
+	Doub fmin;
+	VecDoub y;
+	MatDoub p;
+	Amoeba(const Doub ftoll) : ftol(ftoll) {}
+	template <class T>
+	VecDoub minimize(VecDoub_I &point, const Doub del, T &func)
+	{
+		VecDoub dels(point.size(),del);
+		return minimize(point,dels,func);
+	}
+	template <class T>
+	VecDoub minimize(VecDoub_I &point, VecDoub_I &dels, T &func)
+	{
+		Int ndim=point.size();
+		MatDoub pp(ndim+1,ndim);
+		for (Int i=0;i<ndim+1;i++) {
+			for (Int j=0;j<ndim;j++)
+				pp[i][j]=point[j];
+			if (i !=0 ) pp[i][i-1] += dels[i-1];
+		}
+		return minimize(pp,func);
+	}
+	template <class T>
+	VecDoub minimize(MatDoub_I &pp, T &func)
+	{
+		const Int NMAX=5000;
+		const Doub TINY=1.0e-10;
+		Int ihi,ilo,inhi;
+		mpts=pp.nrows();
+		ndim=pp.ncols();
+		VecDoub psum(ndim),pmin(ndim),x(ndim);
+		p=pp;
+		y.resize(mpts);
+		for (Int i=0;i<mpts;i++) {
+			for (Int j=0;j<ndim;j++)
+				x[j]=p[i][j];
+			y[i]=func(x);
+		}
+		nfunc=0;
+		get_psum(p,psum);
+		for (;;) {
+			ilo=0;
+			ihi = y[0]>y[1] ? (inhi=1,0) : (inhi=0,1);
+			for (Int i=0;i<mpts;i++) {
+				if (y[i] <= y[ilo]) ilo=i;
+				if (y[i] > y[ihi]) {
+					inhi=ihi;
+					ihi=i;
+				} else if (y[i] > y[inhi] && i != ihi) inhi=i;
+			}
+			Doub rtol=2.0*abs(y[ihi]-y[ilo])/(abs(y[ihi])+abs(y[ilo])+TINY);
+			if (rtol < ftol) {
+				SWAP(y[0],y[ilo]);
+				for (Int i=0;i<ndim;i++) {
+					SWAP(p[0][i],p[ilo][i]);
+					pmin[i]=p[0][i];
+				}
+				fmin=y[0];
+				return pmin;
+			}
+			if (nfunc >= NMAX) throw("NMAX exceeded");
+			nfunc += 2;
+			Doub ytry=amotry(p,y,psum,ihi,-1.0,func);
+			if (ytry <= y[ilo])
+				ytry=amotry(p,y,psum,ihi,2.0,func);
+			else if (ytry >= y[inhi]) {
+				Doub ysave=y[ihi];
+				ytry=amotry(p,y,psum,ihi,0.5,func);
+				if (ytry >= ysave) {
+					for (Int i=0;i<mpts;i++) {
+						if (i != ilo) {
+							for (Int j=0;j<ndim;j++)
+								p[i][j]=psum[j]=0.5*(p[i][j]+p[ilo][j]);
+							y[i]=func(psum);
+						}
+					}
+					nfunc += ndim;
+					get_psum(p,psum);
+				}
+			} else --nfunc;
+		}
+	}
+	inline void get_psum(MatDoub_I &p, VecDoub_O &psum)
+	{
+		for (Int j=0;j<ndim;j++) {
+			Doub sum=0.0;
+			for (Int i=0;i<mpts;i++)
+				sum += p[i][j];
+			psum[j]=sum;
+		}
+	}
+	template <class T>
+	Doub amotry(MatDoub_IO &p, VecDoub_O &y, VecDoub_IO &psum,
+		const Int ihi, const Doub fac, T &func)
+	{
+		VecDoub ptry(ndim);
+		Doub fac1=(1.0-fac)/ndim;
+		Doub fac2=fac1-fac;
+		for (Int j=0;j<ndim;j++)
+			ptry[j]=psum[j]*fac1-p[ihi][j]*fac2;
+		Doub ytry=func(ptry);
+		if (ytry < y[ihi]) {
+			y[ihi]=ytry;
+			for (Int j=0;j<ndim;j++) {
+				psum[j] += ptry[j]-p[ihi][j];
+				p[ihi][j]=ptry[j];
+			}
+		}
+		return ytry;
+	}
+};