43 #ifndef O2SCL_MMIN_SIMP2_H 44 #define O2SCL_MMIN_SIMP2_H 51 #include <boost/numeric/ublas/vector.hpp> 53 #include <gsl/gsl_multimin.h> 55 #include <o2scl/mmin.h> 56 #include <o2scl/cblas.h> 57 #include <o2scl/vec_stats.h> 59 #ifndef DOXYGEN_NO_O2NS 130 #ifndef DOXYGEN_INTERNAL 164 for (j = 0; j <
dim; j++) {
166 for (i = 0; i < dim+1; i++) {
169 val /= ((double)dim)+1;
188 for(
size_t i=0;i<P;i++) {
189 for(
size_t j=0;j<
dim;j++) ws1[j]=x1[i][j];
197 return sqrt(ss/(
double)(P));
208 vec_t &xc, func_t &f,
size_t nvar,
217 double alpha=(1.0-coeff)*((
double)P)/((
double)
dim);
218 double beta=(((double)P)*coeff-1.0)/((double)
dim);
220 for(
size_t j=0;j<
dim;j++) {
221 xc[j]=center[j]*alpha;
222 xc[j]+=x1[corner][j]*beta;
233 const size_t P=
dim+1;
236 for(
size_t j=0;j<
dim;j++) {
242 for(
size_t j=0;j<
dim;j++) {
252 S2 += (2.0 / ((double)P)) * xmcd +
253 ((((
double)P) - 1.0) / ((double)P)) * (d * d / ((
double)P));
256 double alpha=1.0/((double)P);
257 for(
size_t j=0;j<
dim;j++) {
258 center[j]-=alpha*x1[i][j];
259 center[j]+=alpha*xx[j];
262 for(
size_t j=0;j<
dim;j++) {
294 for (i=0;i<
dim+1;i++) {
298 for (j=0;j<
dim;j++) {
299 x1[i][j]=0.5*(x1[i][j]+x1[best][j]);
302 if (!std::isfinite(y1[i])) {
303 std::string err=((std::string)
"Function not finite (returned ")+
304 dtos(y1[i])+
" in mmin_simp2::contract_by_best().";
380 template<
class vec2_t>
int set_step(
size_t nv, vec2_t &step) {
383 for(
size_t i=0;i<nv;i++) step_vec[i]=step[i];
410 virtual int mmin(
size_t nn, vec_t &xx,
double &fmin,
414 O2SCL_ERR2(
"Tried to min over zero variables ",
418 int ret=0,status,iter=0;
423 for (
size_t is=0;is<nn;is++) ss[is]=step_vec[is % step_vec.size()];
424 ret=
set(ufunc,nn,xx,ss);
441 status=gsl_multimin_test_size(size,this->
tol_abs);
443 }
while(status == GSL_CONTINUE && iter<this->
ntrial);
445 for (
size_t i=0;i<nn;i++) xx[i]=x[i];
450 if(iter>=this->ntrial) {
451 std::string str=
"Exceeded maximum number of iterations ("+
452 itos(this->ntrial)+
") in mmin_simp2::mmin().";
463 virtual int mmin_twovec(
size_t nn, vec_t &xx, vec_t &xx2,
double &fmin,
466 int ret=0,i,status,iter=0;
471 for (
size_t is=0;is<nn;is++) ss[is]=xx2[is]-xx[is];
472 ret=
set(ufunc,nn,xx,ss);
489 status=gsl_multimin_test_size(size,this->
tol_abs);
491 }
while(status == GSL_CONTINUE && iter<this->
ntrial);
493 for (i=0;i<((int)nn);i++) xx[i]=x[i];
498 if(iter>=this->ntrial) {
499 std::string str=
"Exceeded maximum number of iterations ("+
500 itos(this->ntrial)+
") in mmin_simp2::mmin_twovec().";
510 template<
class mat_t>
514 int ret=0,i,status,iter=0;
534 status=gsl_multimin_test_size(size,this->
tol_abs);
536 }
while(status == GSL_CONTINUE && iter<this->
ntrial);
538 for (i=0;i<((int)nn);i++) sx(0,i)=x[i];
543 if (iter>=this->ntrial) {
544 std::string str=
"Exceeded maximum number of iterations ("+
545 itos(this->ntrial)+
") in mmin_simp2::mmin_simplex().";
560 for(
size_t i=0;i<n+1;i++) x1[i].resize(n);
573 virtual int set(func_t &ufunc,
size_t n, vec_t &ax,
581 for (i=0;i<
dim;i++) x[i]=ax[i];
586 if (!std::isfinite(y1[0])) {
587 std::string err=((std::string)
"Function not finite (returned ")+
588 dtos(y1[0])+
" in mmin_simp2::set().";
591 for(i=0;i<
dim;i++) x1[0][i]=ax[i];
595 for (i=1;i<dim+1;i++) {
596 for(
size_t j=0;j<
dim;j++) x1[i][j]=x[j];
597 x1[i][i-1]=x1[i][i-1]+step_size[i-1];
598 y1[i]=ufunc(dim,x1[i]);
612 template<
class mat_t>
622 for(
size_t i=0;i<dim+1;i++) {
623 for(
size_t j=0;j<
dim;j++) {
626 y1[i]=ufunc(dim,x1[i]);
627 if (!std::isfinite(y1[i])) {
628 std::string err=((std::string)
"Function not finite (returned ")+
629 dtos(y1[i])+
" in mmin_simp2::set_simplex().";
650 double dhi,ds_hi,dlo;
666 }
else if (val > dhi) {
671 }
else if(val > ds_hi) {
681 if (avoid_nonzero && ret1!=0) {
682 std::cout <<
"Found problem move1: " << std::endl;
683 for (
size_t it=0;it<20 && ret1!=0;it++) {
686 std::cout <<
"it,ret: " << it <<
" " << ret1 << std::endl;
694 if (std::isfinite(val) && val<y1[lo]) {
700 if (avoid_nonzero && ret2!=0) {
701 std::cout <<
"Found problem move2: " << std::endl;
702 for (
size_t it=0;it<20 && ret2!=0;it++) {
705 std::cout <<
"it,ret: " << it <<
" " << ret2 << std::endl;
713 if (std::isfinite(val2) && val2<y1[lo]) {
719 }
else if (!std::isfinite(val) || val > y1[s_hi]) {
723 if (std::isfinite(val) && val <= y1[hi]) {
735 if (avoid_nonzero && ret3!=0) {
736 std::cout <<
"Found problem move3: " << std::endl;
737 for (
size_t it=0;it<20 && ret3!=0;it++) {
740 std::cout <<
"it,ret: " << it <<
" " << ret3 << std::endl;
748 if (std::isfinite(val2) && val2 <= y1[hi]) {
757 O2SCL_ERR(
"Function contract_by_best() failed in iterate().",
774 for(i=0;i<
dim;i++) x[i]=x1[lo][i];
800 double value,
double limit,
801 std::string comment) {
803 if (this->
verbose<=0)
return 0;
808 (*this->
outs) << comment <<
" Iteration: " << iter << std::endl;
809 (*this->
outs) <<
"x: ";
810 for(i=0;i<nv;i++) (*this->
outs) << xx[i] <<
" ";
811 (*this->
outs) << std::endl;
812 if (print_simplex>0) {
813 if (print_simplex==2) {
814 (*this->
outs) <<
"Simplex Values:" << std::endl;
815 for(i=0;i<nv+1;i++) {
816 (*this->
outs) << i <<
": ";
817 for(
size_t j=0;j<nv;j++) {
818 (*this->
outs) << simp[i][j] <<
" ";
820 (*this->
outs) <<
": " << y1[i] << std::endl;
823 (*this->
outs) <<
"Simplex Values:" << std::endl;
824 for(i=0;i<nv+1;i++) {
825 (*this->
outs) << y1[i] <<
" ";
827 (*this->
outs) << std::endl;
830 (*this->
outs) <<
"y: " << y <<
" Val: " << value <<
" Lim: " 831 << limit << std::endl;
833 (*this->
outs) <<
"Press a key and type enter to continue. ";
841 virtual const char *
type() {
return "mmin_simp2";}
843 #ifndef DOXYGEN_INTERNAL 850 (
const mmin_simp2<func_t,vec_t>&);
856 #ifndef DOXYGEN_NO_O2NS vec_t ws1
Workspace vector 1.
double S2
Squared simplex size.
std::function< double(size_t, const boost::numeric::ublas::vector< double > &)> multi_funct
Multi-dimensional function typedef.
virtual int update_point(size_t i, vec_t &xx, double val)
Update point i in the simplex with values xx.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
std::istream * ins
Stream for verbose input.
ubvector step_vec
Vector of step sizes.
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
std::ostream * outs
Stream for verbose output.
virtual int allocate(size_t n)
Allocate the memory.
invalid argument supplied by user
virtual int mmin_twovec(size_t nn, vec_t &xx, vec_t &xx2, double &fmin, func_t &ufunc)
Calculate the minimum min of func w.r.t the array x of size nvar, using xx and xx2 to specify the sim...
exceeded max number of iterations
Multidimensional minimization by the Simplex method (v2) (GSL)
int mmin_simplex(size_t nn, mat_t &sx, double &fmin, func_t &ufunc)
Calculate the minimum min of func w.r.t the array x of size nvar, given an initial simplex...
void vector_min(size_t n, const vec_t &data, size_t &index, data_t &val)
Compute the minimum of the first n elements of a vector.
int verbose
Output control.
double size
Size of current simplex computed by iterate()
double dnrm2(const size_t N, const vec_t &X)
Compute the norm of the vector X.
virtual int iterate()
Perform an iteration.
bool err_nonconv
If true, call the error handler if the routine does not "converge".
matrix, vector lengths are not conformant
Multidimensional minimization [abstract base].
vec_t xmc
Distance of vector from center.
bool avoid_nonzero
If true, try to automatically avoid regions where the function returns a non-zero value (default fals...
vec_t x
Present minimum vector computed by iterate()
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
size_t dim
Number of variables to be mind over.
virtual int try_corner_move(const double coeff, size_t corner, vec_t &xc, func_t &f, size_t nvar, double &newval)
Move a corner of a simplex.
double ddot(const size_t N, const vec_t &X, const vec2_t &Y)
Compute .
vec_t * x1
An array of n+1 vectors containing the simplex.
double compute_size()
Compute the size of the simplex.
virtual int mmin(size_t nn, vec_t &xx, double &fmin, func_t &ufunc)
Calculate the minimum min of func w.r.t the array x of size nvar.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
int last_ntrial
The number of iterations for in the most recent minimization.
double tol_abs
The independent variable tolerance.
int set_simplex(func_t &ufunc, mat_t &sx)
Set the function and initial simplex.
void daxpy(const double alpha, const size_t N, const vec_t &X, vec2_t &Y)
Compute .
int set_step(size_t nv, vec2_t &step)
Set the step sizes for each independent variable.
vec_t center
Center of simplex.
vec_t ws3
Workspace vector 3.
problem with user-supplied function
double fval
Function value at minimum computed by iterate()
ubvector y1
The n+1 function values at the simplex points.
virtual int contract_by_best(size_t best, func_t &f, size_t nvar)
Contract the simplex towards the best point.
vec_t ws2
Workspace vector 2.
std::string itos(int x)
Convert an integer to a string.
int print_simplex
Print simplex information in print_iter() (default 0)
int compute_center()
Compute the center of the simplex.
bool set_called
True if set() has been called.
virtual int print_iter(size_t nv, vec_t &xx, vec_t *simp, double y, int iter, double value, double limit, std::string comment)
Print out iteration information.
virtual const char * type()
Return string denoting type("mmin_simp2")
int ntrial
Maximum number of iterations.