#include <stdio.h>#include <math.h>#define WANT_STREAM#define WANT_MATH#include "newmatap.h"#include "newmatio.h"#define ANSI#define TINY 1.0e-10#define NMAX 5000//#define GET_PSUM \// for (j=1;j<=ndim;j++) {\// for (sum=0.0,i=1;i<=mpts;i++) sum += p[i][j];\// psum(j)=sum;}#define SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;}void freeColumnVector(ColumnVector rv, int dim){ for (int i=1; i<= dim; i++) { rv(i) = 0; }}/*Multidimensional minimization of the function funk(x) where x[1..ndim] is a vector in ndim dimensions, by the downhill simplex method of Nelder and Mead. The matrix p[1..ndim+1][1..ndim] is input. Its ndim+1 rows are ndim-dimensional vectors which are the vertices ofthe starting simplex. Also input is the vector y[1..ndim+1], whose components must be preinitialized to the values of funk evaluated at the ndim+1 vertices (rows) of p; and ftol the fractional convergence tolerance to be achieved in the function value (n.b.!). On output, p andy will have been reset to ndim+1 new points all within ftol of a minimum function value, and nfunk gives the number of function evaluations taken.*/void amoeba(Matrix p, ColumnVector y, int ndim, float ftol, float (*funk)(ColumnVector), int *nfunk){ float amotry(Matrix p, ColumnVector y, ColumnVector psum, int ndim, float (*funk)(ColumnVector), int ihi, float fac); int i,ihi,ilo,inhi,j,mpts=ndim+1; float rtol,sum,swap,ysave,ytry; ColumnVector psum(ndim); *nfunk=0; for (j=1;j<=ndim;j++) { psum(j)=p.Column(j).Sum(); } for (;;) { ilo=1; /*First we must determine which point is the highest (worst), next-highest, and lowest (best), by looping over the points in the simplex.*/ ihi = y(1)>y(2) ? (inhi=2,1) : (inhi=1,2); for (i=1;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; } rtol=2.0*fabs(y(ihi)-y(ilo))/(fabs(y(ihi))+fabs(y(ilo))+TINY); /*Compute the fractional range from highest to lowest and return if satisfactory.*/ if (rtol < ftol) { /*If returning, put best point and value in slot 1.*/ SWAP(y(1),y(ilo)) for (i=1;i<=ndim;i++) SWAP(p(1,i),p(ilo,i)) break; } if (*nfunk >= NMAX) return; *nfunk += 2; /*Begin a new iteration. First extrapolate by a factor −1 through the face of the simplex across from the high point, i.e., reflect the simplex from the high point.*/ ytry=amotry(p,y,psum,ndim,funk,ihi,-1.0); if (ytry <= y(ilo)) /*Gives a result better than the best point, so try an additional extrapolation by a factor 2.*/ ytry=amotry(p,y,psum,ndim,funk,ihi,2.0); else if (ytry >= y(inhi)) { /*The reflected point is worse than the second-highest, so look for an intermediate lower point, i.e., do a one-dimensional contraction.*/ ysave=y(ihi); ytry=amotry(p,y,psum,ndim,funk,ihi,0.5); if (ytry >= ysave) { /*Can't seem to get rid of that high point. Better contract around the lowest (best) point.*/ for (i=1;i<=mpts;i++) { if (i != ilo) { for (j=1;j<=ndim;j++) p(i,j)=psum(j)=0.5*(p(i,j)+p(ilo,j)); float a; a = (*funk)(psum); y(i) = a; //y(i)=(*funk)psum; } } *nfunk += ndim; /*Keep track of function evaluations.*/ /*Recompute psum.*/ for (j=1;j<=ndim;j++) { psum(j)=p.Column(j).Sum(); } } } else --(*nfunk); /*Correct the evaluation count.*/ } /*Go back for the test of doneness and the next iteration.*/ //ColumnVector psum(ndim); freeColumnVector(psum, ndim); //free_vector(psum,1,ndim); }/*Extrapolates by a factor fac through the face of the simplex across from the high point, tries it, and replaces the high point if the new point is better.*/float amotry(Matrix p, ColumnVector y, ColumnVector psum, int ndim, float (*funk)(ColumnVector), int ihi, float fac){ int j; float fac1,fac2,ytry; ColumnVector ptry(ndim);// ptry=vector(1,ndim); fac1=(1.0-fac)/ndim; fac2=fac1-fac; for (j=1;j<=ndim;j++) ptry(j)=psum(j)*fac1-p(ihi,j)*fac2; ytry=(*funk)(ptry); /* Evaluate the function at the trial point. If it's better than the highest, then replace the highest. */ if (ytry < y(ihi)) { y(ihi)=ytry; for (j=1;j<=ndim;j++) { psum(j) += ptry(j)-p(ihi,j); p(ihi,j)=ptry(j); } } //free_vector(ptry,1,ndim); freeColumnVector(ptry, ndim); return ytry;}float x_squared(ColumnVector x){ return x(1)*x(1);}void main(int argc, char* argv){ printf("bu"); int ndim = 1; Matrix p(2,1); p(1,1) = 30; p(2,1) = 40; ColumnVector y(ndim+1); y(1) = 30*30; y(2) = 40*40; int nfunk = 100; float ftol = 0.1; amoeba(p, y, ndim, ftol, x_squared, &nfunk); for (int i=1; i<=ndim+1; i++) printf("\ny %d: %.2f", i, y(i)); for (int i=1; i<=ndim+1; i++) for (int j=1; j<=ndim; j++) printf("\np(%d,%d): %.2f", i, j, p(i,j)); printf("\nnfunk: %d", nfunk); getchar();}