多変数関数の極小値求解プログラム２

/*		amoeba.c		多次元の滑降シンプレックス法	*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define		NMAX		5000		/* 関数評価回数の上限 */
#define		SWAP(a, b)		{swap = (a); (a) = (b); (b) = swap;}

void amoeba(double **p, double y[], int ndim, double ftol,
	double (*funk)(double []), int *nfunk);
double amotry(double **p, double y[], double psum[], int ndim,
	double (*funk)(double []), int ihi, double fac);
double *vector(long size);
void free_vector(double *v);
double **matrix(long rsize, long csize);
void free_matrix(double **m);
double f(double x[]);

/*
Nelder-Mead の滑降シンプレックス法による関数 funk(x) の多次元の最小化。
x[ ] は ndim 次元のベクトル。入力する行列 p[ ][ ] の ndim+1 個の行は、出発点
でのシンプレックスの頂点の座標を表す ndim 次元のベクトル。入力するベクトル
ftol は関数値（極小点の座標ではない！）の収束を判定するための相対許容値。
戻り時の p と y は、最終位置での ndim+1 個の頂点の座標と関数値（どの関数値も
最小の関数値との相対的な外れが ftol 以内）。nfunk は関数の評価回数。
*/
void amoeba(double **p, double y[], int ndim, double ftol,
	double (*funk)(double []), int *nfunk)
{
	int i, ihi, ilo, inhi, j, mpts = ndim + 1;
	double rtol, sum, swap, ysave, ytry, *psum;

	for(j = 0; j < mpts; j++)	y[j] = (*funk)(p[j]);
	psum = vector(ndim);
	*nfunk = 0;
	for(j = 0; j < ndim; j++)
	{
		sum = 0.0;
		for(i = 0; i < mpts; i++)	sum += p[i][j];
		psum[j] = sum;
	}
	for(;;)
	{
		ilo = 0;
/* まずシンプレックスの頂点についてループし、最悪、最悪の次、最良の点を調べる。*/
		if(y[0] > y[1])
		{
			inhi = 1;
			ihi = 0;
		}
		else
		{
			inhi = 0;
			ihi = 1;
		}
		for(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;
		}
		rtol = 2.0 * fabs(y[ihi] - y[ilo]) / (fabs(y[ihi]) + fabs(y[ilo]));

/* 関数値の最大の比を求め、これが十分小さいなら終了。*/
		if(rtol < ftol)				/* 戻る際には最良の点を頭に持って行く */
		{
			SWAP(y[0], y[ilo]);
			for(i = 0; i < ndim; i++)	SWAP(p[0][i], p[ilo][i]);
			break;
		}
		if(*nfunk >= NMAX)
		{
			printf("Process Stop : NMAX(=%d) exceeded.\n", NMAX);
			SWAP(y[0], y[ilo]);
			for(i = 0; i < ndim; i++)	SWAP(p[0][i], p[ilo][i]);
			break;
		}
		*nfunk += 2;
/* 新しい反復を始める。まず最悪の点を残りの点の重心の反対側に対称移動する。*/
		ytry = amotry(p, y, psum, ndim, funk, ihi, -1.0);
		if(ytry <= y[ilo])
		{
/* 現在の最良の点より良くなったので、さらに２倍だけ進んでみる。*/
			ytry = amotry(p, y, psum, ndim, funk, ihi, 2.0);
		}
		else if(ytry >= y[inhi])
		{
/* 対称移動した点は２番目に悪い点よりも悪いので１次元の収縮を試みる。*/
			ysave = y[ihi];
			ytry = amotry(p, y, psum, ndim, funk, ihi, 0.5);
			if(ytry >= ysave)
			{
/* それでも良くないなら全体を最良点に向かって収縮。*/
				for(i = 0; i < mpts; i++)
				{
					if(i != ilo)
					{
						for(j = 0; j < ndim; j++)
							p[i][j] = psum[j] = 0.5 * (p[i][j] + p[ilo][j]);
						y[i] = (*funk)(psum);
					}
				}
				*nfunk += ndim;			/* 関数の評価回数を数える。 */
				for(j = 0; j < ndim; j++)		/* psum を再計算。 */
				{
					sum = 0.0;
					for(i = 0; i < mpts; i++)	sum += p[i][j];
					psum[j] = sum;
				}
			}
		}
		else	--(*nfunk);				/* 評価回数を補正。 */
	}
	free_vector(psum);
}

/*
最悪の点を、それ以外の点の重心について fac 倍に伸ばした点に仮に移し、その点の
方が良ければ最悪の点を捨てて置き換える。
*/
double amotry(double **p, double y[], double psum[], int ndim,
	double (*funk)(double []), int ihi, double fac)
{
	int j;
	double fac1, fac2, ytry, *ptry;

	ptry = vector(ndim);
	fac1 = (1.0 - fac) / ndim;
	fac2 = fac1 - fac;
	for(j = 0; j < ndim; j++)	ptry[j] = psum[j] * fac1 - p[ihi][j] * fac2;
	ytry = (*funk)(ptry);				/* 新しい点で関数評価 */
	if(ytry < y[ihi])					/* 新しい点の方が良ければ元の最悪の */
	{									/* 捨てて置き換える。 */
		y[ihi] = ytry;
		for(j = 0; j < ndim; j++)
		{
			psum[j] += ptry[j] - p[ihi][j];
			p[ihi][j] = ptry[j];
		}
	}
	free_vector(ptry);
	return ytry;
}

/* allocate a double vector with subscript range v[nl..nh] */
double *vector(long size)
{
	double *v;

	v = (double *)malloc((size_t)((size + 1) * sizeof(double)));
	if(v == NULL)	fprintf(stderr, "Error : allocation failure in vector().\n");
	return v;
}

/* free a double vector allocated with vector() */
void free_vector(double *v)
{
	free((char *)v);
}

/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
double **matrix(long rsize, long csize)
{
	double **m, **p;
	long i;

	m = (double **)malloc((size_t)(rsize * sizeof(double *)));
	if(m == NULL)
		fprintf(stderr, "Error : allocation failure 1 in matrix().\n");
	*m = (double *)malloc((size_t)(rsize * csize * sizeof(double)));
	if(*m == NULL)
		fprintf(stderr, "Error : allocation failure 2 in matrix().\n");
	for(p = m, i = 1; i < rsize; i++, p++)	*(p + 1) = *p + csize;
	return m;
}

/* free a double matrix allocated by matrix() */
void free_matrix(double **m)
{
	free((char *)*m);
	free((char *)m);
}

double f(double x[])
{
	double a0, a1, a2;

	a0 = (x[0] - 1.0);
	a1 = (x[1] - 2.0);
	a2 = (x[2] - 3.0);
	return fabs(a0 * a0 * a0) + a1 * a1 + fabs(a2);
}

int main(void)
{
	double **p, *y, ftol;
	int i, j, ndim, mpts, nfunk;

	ndim = 3;
	mpts = ndim + 1;
	ftol = 1.0e-12;

	p = matrix(mpts, ndim);
	y = vector(mpts);
	srand(1);
	for(i = 0; i < mpts; i++)
		for(j = 0; j < ndim; j++)	p[i][j] = (double)rand() / (double)RAND_MAX;

	printf("f(x,y,z) = |(x-1)^3| + (y-2)^2 + |z-3|\n");
	printf("極小点：x=1, y=2, z=3\n\n");
	printf("[ input ]\n");
	for(i = 0; i < mpts; i++)
	{
		printf("f(%f", p[i][0]);
		for(j = 1; j < ndim; j++)	printf(",%f", p[i][j]);
		printf(") = %f\n", f(p[i]));
	}
	printf("tolerance = %e\n\n", ftol);

	amoeba(p, y, ndim, ftol, f, &nfunk);

	printf("[ output ]\n");
	printf("count=%d\n", nfunk);
	for(i = 0; i < mpts; i++)
	{
		printf("f(%f", p[i][0]);
		for(j = 1; j < ndim; j++)	printf(",%f", p[i][j]);
		printf(") = %f\n", y[i]);
	}

	free_matrix(p);
	free_vector(y);

	return 1;
}

f(x,y,z) = |(x-1)^3| + (y-2)^2 + |z-3|
極小点：x=1, y=2, z=3

[ input ]
f(0.010559,0.003967,0.335154) = 7.617647
f(0.033265,0.355724,0.217200) = 6.389932
f(0.536973,0.195776,0.700339) = 5.654155
f(0.949919,0.274789,0.444288) = 5.532191
tolerance = 1.000000e-12

Process Stop : NMAX(=5000) exceeded.
[ output ]
count=5000
f(0.999505,2.000010,3.000000) = 0.000000
f(0.999505,2.000010,3.000000) = 0.000000
f(0.999505,2.000010,3.000000) = 0.000000
f(0.999505,2.000010,3.000000) = 0.000000