/*
File name: simplex.c
Created by: Ljubomir Buturovic
Created: 07/02/2004
Purpose: implementation of Nelder-Mead simplex algorithm for
function optimization without derivatives.
*/
/*
Copyright 2004 Ljubomir J. Buturovic
Permission is hereby granted, free of charge, to any person
obtaining a copy of this software and associated documentation files
(the "Software"), to deal in the Software without restriction,
including without limitation the rights to use, copy, modify, merge,
publish, distribute, sublicense, and/or sell copies of the Software,
and to permit persons to whom the Software is furnished to do so,
subject to the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/
#include <math.h>
#include <stdlib.h>
#include <errno.h>
#include <float.h>
#include "simplex.h"
#define SMALL 1E-10
static int simplex_alloc(float **centroid, float **xr, float **xe, float **xc,
float **xcc, int len)
{
int status;
status = -1;
*centroid = malloc(len*sizeof(float));
if (*centroid)
{
*xr = malloc(len*sizeof(len));
if (*xr)
{
*xe = malloc(len*sizeof(len));
if (*xe)
{
*xc = malloc(len*sizeof(len));
if (*xc)
{
*xcc = malloc(len*sizeof(len));
if (*xcc)
status = 0;
}
}
}
}
else
status = -1;
return status;
}
static void simplex_free(float *centroid, float *xr, float *xe, float *xc,
float *xcc)
{
free(centroid);
free(xr);
free(xe);
free(xc);
free(xcc);
}
/*
Replace the worst point (imax) in simplex 'smx' with 'xr'.
*/
static void accept_point(float **smx, int len, float *xr, int imax, float *fval,
int iteration, crit_func *func, void *parameters, int *errc)
{
int j;
for (j = 0; j < len; j++)
smx[imax][j] = xr[j];
fval[imax] = func(xr, len, iteration, parameters, errc);
}
/*
Implementation of the shrink step in simplex algorithm of
Nelder-Mead.
In case of success, set *errc to 0, otherwise set it to error code
set by func(). The only possible errors are function evaluation
errors set by func().
*/
static void shrink(float **smx, int len, float *fval, float *xc, int imin,
float fmin, float sigma, crit_func *func, void *parameters,
int iteration, int *errc)
{
int i;
int j;
if (errc)
{
*errc = 0;
/*
Save the min.
*/
for (j = 0; j < len; j++)
xc[j] = smx[imin][j];
for (i = 1; (i < len+1) && (*errc == 0); i++)
{
for (j = 0; j < len; j++)
smx[i][j] = xc[j]+sigma*(smx[i][j]-xc[j]);
fval[i] = func(smx[i], len, iteration, parameters, errc);
}
if (*errc == 0)
{
for (j = 0; j < len; j++)
smx[0][j] = xc[j];
fval[0] = fmin;
}
}
}
/*
simplex() returns minimum value of function 'func' of 'len'
variables, using simplex method of Nelder-Mead. On input, 'smx' is
the (len+1)*len matrix (array of len+1 pointers to len-long arrays)
representing the initial simplex; 'fval' are the corresponding len+1
values of 'func'; 'itmax' is the maximum allowed number of
iterations; 'parameters' is a pointer to a structure used to pass
additional arguments to the criterion function. On output, 'smx'
contains the simplex at time of termination. If the algorithm
converged (as opposed to reaching 'itmax'), all points in 'smx' will
have function values within 'ftol' of each other. 'fval' has the
corresponding function values and 'iter' the number of iterations it
took to terminate.
The signature of func() must be:
float func(float *x, int n, int iteration, void *parameters, int *errc);
Here, 'x' is the n-dimensional point at which the function is
evaluated and 'n' is length of the vector. 'parameters' is a
pointer to the structure used to pass additional arguments to
func(). 'iteration' is the current iteration, set by simplex(), and
'errc' contains the error code set by func(). In case of successful
evaluation, func() must set *errc to 0, and non-zero otherwise.
On success, simplex() returns the minimum value, sets 'iter' to the
number of iterations, and sets 'errc' to 0; on error, it returns
MAX_FLOAT and sets error code in 'errc'. The error codes are EINVAL
for bad arguments, malloc() error codes, or function evaluation
error code as set by func().
The implementation follows: Margaret H. Wright, Optimization Methods
For Base Station Placement In Wireless Applications, Proceedings of
1998 Vehicular Technology Conference, Vol. 89, pp. 11513--11517,
1998.
*/
float simplex(float **smx, int len, float *fval, crit_func *func,
void *parameters, int itmax, float ftol, int *iter,
int *errc)
{
int i;
int j;
int done;
int iteration;
int imin;
int imax;
int in;
int status;
float rdiff;
float val;
float fmin;
float fmax;
float fn;
float fr;
float fe;
float fc;
float fcc;
float rho;
float kappa;
float gamma;
float sigma;
float *centroid;
float *xr;
float *xe;
float *xc;
float *xcc;
val = FLT_MAX;
iteration = 1;
imin = 0;
if (smx && (len > 0) && fval && func && (itmax > 0) && errc)
{
done = 0;
rho = 1.0;
kappa = 2.0;
gamma = 0.5;
sigma = 0.5;
status = simplex_alloc(¢roid, &xr, &xe, &xc, &xcc, len);
if (!status)
{
while (!done)
{
/*
Find fmin, fmax, and fn (next-worst point after fmax).
*/
fmin = fval[0];
fmax = fval[0];
fn = fval[0];
imin = 0;
imax = 0;
in = 0;
for (i = 0; i < len+1; i++)
{
if (fval[i] < fmin)
{
fmin = fval[i];
imin = i;
}
else if (fval[i] > fmax)
{
fn = fmax;
in = imax;
fmax = fval[i];
imax = i;
}
else if (fval[i] > fn)
{
fn = fval[i];
in = i;
}
}
/*
Compute relative difference between the extreme vertices of the
simplex. If it is less than ftol, we are done.
*/
rdiff = 2.0*fabsf(fmax-fmin)/(fabsf(fmax)+fabsf(fmin)+SMALL);
if (rdiff < ftol)
done = 1;
else if (iteration >= itmax+1)
done = 1;
else
{
/*
Compute the reflection point xr and the corresponding
function value fr.
*/
for (i = 0; i < len; i++)
centroid[i] = 0.0;
for (i = 0; i < len+1; i++)
{
if (i != imax)
for (j = 0; j < len; j++)
centroid[j] += smx[i][j];
}
for (j = 0; j < len; j++)
centroid[j] = centroid[j]/len;
for (j = 0; j < len; j++)
xr[j] = centroid[j]+rho*(centroid[j]-smx[imax][j]);
fr = func(xr, len, iteration, parameters, errc);
if (*errc)
done = 1;
else
{
/*
If fmin <= fr < fn, replace the worst point with xr.
*/
if ((fmin <= fr) && (fr < fn))
accept_point(smx, len, xr, imax, fval, iteration, func,
parameters, errc);
else if (fr < fmin)
{
/*
Expansion step.
*/
for (j = 0; j < len; j++)
xe[j] = centroid[j]+rho*kappa*(centroid[j]-smx[imax][j]);
fe = func(xe, len, iteration, parameters, errc);
if (*errc)
done = 1;
else
{
if (fe < fr)
accept_point(smx, len, xe, imax, fval, iteration,
func, parameters, errc);
else
accept_point(smx, len, xr, imax, fval, iteration,
func, parameters, errc);
}
}
else
{
/*
Contraction step.
*/
if (fr < fmax)
{
/*
Outside contraction.
*/
for (j = 0; j < len; j++)
xc[j] = centroid[j]+rho*gamma*(centroid[j]-smx[imax][j]);
fc = func(xc, len, iteration, parameters, errc);
if (*errc)
done = 1;
else
{
if (fc <= fr)
accept_point(smx, len, xc, imax, fval, iteration,
func, parameters, errc);
else
/*
Shrink step.
*/
shrink(smx, len, fval, xc, imin, fmin, sigma,
func, parameters, iteration, errc);
}
}
else
{
/*
Inside contraction.
*/
for (j = 0; j < len; j++)
xcc[j] = centroid[j]-gamma*(centroid[j]-smx[imax][j]);
fcc = func(xcc, len, iteration, parameters, errc);
if (*errc)
done = 1;
else
{
if (fcc < fmax)
accept_point(smx, len, xcc, imax, fval, iteration,
func, parameters, errc);
else
/*
Shrink step.
*/
shrink(smx, len, fval, xc, imin, fmin, sigma,
func, parameters, iteration, errc);
}
}
}
}
iteration++;
}
}
simplex_free(centroid, xr, xe, xc, xcc);
val = fval[imin];
}
else if (errc)
*errc = errno;
if (errc && *errc)
val = FLT_MAX;
}
else if (errc)
*errc = EINVAL;
if (iter)
*iter = iteration-1;
return val;
}
syntax highlighted by Code2HTML, v. 0.9.1