diff --git a/.gitignore b/.gitignore index 3161cea79..1f3a89a9f 100644 --- a/.gitignore +++ b/.gitignore @@ -7,5 +7,4 @@ tegra/ .sw[a-z] .*.swp tags -build/ Thumbs.db diff --git a/modules/optim/doc/downhill_simplex_method.rst b/modules/optim/doc/downhill_simplex_method.rst new file mode 100644 index 000000000..94d084c23 --- /dev/null +++ b/modules/optim/doc/downhill_simplex_method.rst @@ -0,0 +1,161 @@ +Downhill Simplex Method +======================= + +.. highlight:: cpp + +optim::DownhillSolver +--------------------------------- + +.. ocv:class:: optim::DownhillSolver + +This class is used to perform the non-linear non-constrained *minimization* of a function, given on an *n*-dimensional Euclidean space, +using the **Nelder-Mead method**, also known as **downhill simplex method**. The basic idea about the method can be obtained from +(`http://en.wikipedia.org/wiki/Nelder-Mead\_method `_). It should be noted, that +this method, although deterministic, is rather a heuristic and therefore may converge to a local minima, not necessary a global one. +It is iterative optimization technique, which at each step uses an information about the values of a function evaluated only at +*n+1* points, arranged as a *simplex* in *n*-dimensional space (hence the second name of the method). At each step new point is +chosen to evaluate function at, obtained value is compared with previous ones and based on this information simplex changes it's shape +, slowly moving to the local minimum. + +Algorithm stops when the number of function evaluations done exceeds ``termcrit.maxCount``, when the function values at the +vertices of simplex are within ``termcrit.epsilon`` range or simplex becomes so small that it +can enclosed in a box with ``termcrit.epsilon`` sides, whatever comes first, for some defined by user +positive integer ``termcrit.maxCount`` and positive non-integer ``termcrit.epsilon``. + +:: + + class CV_EXPORTS Solver : public Algorithm + { + public: + class CV_EXPORTS Function + { + public: + virtual ~Function() {} + //! ndim - dimensionality + virtual double calc(const double* x) const = 0; + }; + + virtual Ptr getFunction() const = 0; + virtual void setFunction(const Ptr& f) = 0; + + virtual TermCriteria getTermCriteria() const = 0; + virtual void setTermCriteria(const TermCriteria& termcrit) = 0; + + // x contain the initial point before the call and the minima position (if algorithm converged) after. x is assumed to be (something that + // after getMat() will return) row-vector or column-vector. *It's size and should + // be consisted with previous dimensionality data given, if any (otherwise, it determines dimensionality)* + virtual double minimize(InputOutputArray x) = 0; + }; + + class CV_EXPORTS DownhillSolver : public Solver + { + public: + //! returns row-vector, even if the column-vector was given + virtual void getInitStep(OutputArray step) const=0; + //!This should be called at least once before the first call to minimize() and step is assumed to be (something that + //! after getMat() will return) row-vector or column-vector. *It's dimensionality determines the dimensionality of a problem.* + virtual void setInitStep(InputArray step)=0; + }; + +It should be noted, that ``optim::DownhillSolver`` is a derivative of the abstract interface ``optim::Solver``, which in +turn is derived from the ``Algorithm`` interface and is used to encapsulate the functionality, common to all non-linear optimization +algorithms in the ``optim`` module. + +optim::DownhillSolver::getFunction +-------------------------------------------- + +Getter for the optimized function. The optimized function is represented by ``Solver::Function`` interface, which requires +derivatives to implement the sole method ``calc(double*)`` to evaluate the function. + +.. ocv:function:: Ptr optim::DownhillSolver::getFunction() + + :return: Smart-pointer to an object that implements ``Solver::Function`` interface - it represents the function that is being optimized. It can be empty, if no function was given so far. + +optim::DownhillSolver::setFunction +----------------------------------------------- + +Setter for the optimized function. *It should be called at least once before the call to* ``DownhillSolver::minimize()``, as +default value is not usable. + +.. ocv:function:: void optim::DownhillSolver::setFunction(const Ptr& f) + + :param f: The new function to optimize. + +optim::DownhillSolver::getTermCriteria +---------------------------------------------------- + +Getter for the previously set terminal criteria for this algorithm. + +.. ocv:function:: TermCriteria optim::DownhillSolver::getTermCriteria() + + :return: Deep copy of the terminal criteria used at the moment. + +optim::DownhillSolver::setTermCriteria +------------------------------------------ + +Set terminal criteria for downhill simplex method. Two things should be noted. First, this method *is not necessary* to be called +before the first call to ``DownhillSolver::minimize()``, as the default value is sensible. Second, the method will raise an error +if ``termcrit.type!=(TermCriteria::MAX_ITER+TermCriteria::EPS)``, ``termcrit.epsilon<=0`` or ``termcrit.maxCount<=0``. That is, +both ``epsilon`` and ``maxCount`` should be set to positive values (non-integer and integer respectively) and they represent +tolerance and maximal number of function evaluations that is allowed. + +Algorithm stops when the number of function evaluations done exceeds ``termcrit.maxCount``, when the function values at the +vertices of simplex are within ``termcrit.epsilon`` range or simplex becomes so small that it +can enclosed in a box with ``termcrit.epsilon`` sides, whatever comes first. + +.. ocv:function:: void optim::DownhillSolver::setTermCriteria(const TermCriteria& termcrit) + + :param termcrit: Terminal criteria to be used, represented as ``TermCriteria`` structure (defined elsewhere in openCV). Mind you, that it should meet ``(termcrit.type==(TermCriteria::MAX_ITER+TermCriteria::EPS) && termcrit.epsilon>0 && termcrit.maxCount>0)``, otherwise the error will be raised. + +optim::DownhillSolver::getInitStep +----------------------------------- + +Returns the initial step that will be used in downhill simplex algorithm. See the description +of corresponding setter (follows next) for the meaning of this parameter. + +.. ocv:function:: void optim::getInitStep(OutputArray step) + + :param step: Initial step that will be used in algorithm. Note, that although corresponding setter accepts column-vectors as well as row-vectors, this method will return a row-vector. + +optim::DownhillSolver::setInitStep +---------------------------------- + +Sets the initial step that will be used in downhill simplex algorithm. Step, together with initial point (givin in ``DownhillSolver::minimize``) +are two *n*-dimensional vectors that are used to determine the shape of initial simplex. Roughly said, initial point determines the position +of a simplex (it will become simplex's centroid), while step determines the spread (size in each dimension) of a simplex. To be more precise, +if :math:`s,x_0\in\mathbb{R}^n` are the initial step and initial point respectively, the vertices of a simplex will be: :math:`v_0:=x_0-\frac{1}{2} +s` and :math:`v_i:=x_0+s_i` for :math:`i=1,2,\dots,n` where :math:`s_i` denotes projections of the initial step of *n*-th coordinate (the result +of projection is treated to be vector given by :math:`s_i:=e_i\cdot\left`, where :math:`e_i` form canonical basis) + +.. ocv:function:: void optim::setInitStep(InputArray step) + + :param step: Initial step that will be used in algorithm. Roughly said, it determines the spread (size in each dimension) of an initial simplex. + +optim::DownhillSolver::minimize +----------------------------------- + +The main method of the ``DownhillSolver``. It actually runs the algorithm and performs the minimization. The sole input parameter determines the +centroid of the starting simplex (roughly, it tells where to start), all the others (terminal criteria, initial step, function to be minimized) +are supposed to be set via the setters before the call to this method or the default values (not always sensible) will be used. + +.. ocv:function:: double optim::DownhillSolver::minimize(InputOutputArray x) + + :param x: The initial point, that will become a centroid of an initial simplex. After the algorithm will terminate, it will be setted to the point where the algorithm stops, the point of possible minimum. + + :return: The value of a function at the point found. + +optim::createDownhillSolver +------------------------------------ + +This function returns the reference to the ready-to-use ``DownhillSolver`` object. All the parameters are optional, so this procedure can be called +even without parameters at all. In this case, the default values will be used. As default value for terminal criteria are the only sensible ones, +``DownhillSolver::setFunction()`` and ``DownhillSolver::setInitStep()`` should be called upon the obtained object, if the respective parameters +were not given to ``createDownhillSolver()``. Otherwise, the two ways (give parameters to ``createDownhillSolver()`` or miss the out and call the +``DownhillSolver::setFunction()`` and ``DownhillSolver::setInitStep()``) are absolutely equivalent (and will drop the same errors in the same way, +should invalid input be detected). + +.. ocv:function:: Ptr optim::createDownhillSolver(const Ptr& f,InputArray initStep, TermCriteria termcrit) + + :param f: Pointer to the function that will be minimized, similarly to the one you submit via ``DownhillSolver::setFunction``. + :param step: Initial step, that will be used to construct the initial simplex, similarly to the one you submit via ``DownhillSolver::setInitStep``. + :param termcrit: Terminal criteria to the algorithm, similarly to the one you submit via ``DownhillSolver::setTermCriteria``. diff --git a/modules/optim/doc/optim.rst b/modules/optim/doc/optim.rst index bead2122a..cdbaeac25 100644 --- a/modules/optim/doc/optim.rst +++ b/modules/optim/doc/optim.rst @@ -8,3 +8,4 @@ optim. Generic numerical optimization :maxdepth: 2 linear_programming + downhill_simplex_method diff --git a/modules/optim/include/opencv2/optim.hpp b/modules/optim/include/opencv2/optim.hpp index 0bbedad38..9acd95d02 100644 --- a/modules/optim/include/opencv2/optim.hpp +++ b/modules/optim/include/opencv2/optim.hpp @@ -47,6 +47,45 @@ namespace cv{namespace optim { +class CV_EXPORTS Solver : public Algorithm +{ +public: + class CV_EXPORTS Function + { + public: + virtual ~Function() {} + //! ndim - dimensionality + virtual double calc(const double* x) const = 0; + }; + + virtual Ptr getFunction() const = 0; + virtual void setFunction(const Ptr& f) = 0; + + virtual TermCriteria getTermCriteria() const = 0; + virtual void setTermCriteria(const TermCriteria& termcrit) = 0; + + // x contain the initial point before the call and the minima position (if algorithm converged) after. x is assumed to be (something that + // after getMat() will return) row-vector or column-vector. *It's size and should + // be consisted with previous dimensionality data given, if any (otherwise, it determines dimensionality)* + virtual double minimize(InputOutputArray x) = 0; +}; + +//! downhill simplex class +class CV_EXPORTS DownhillSolver : public Solver +{ +public: + //! returns row-vector, even if the column-vector was given + virtual void getInitStep(OutputArray step) const=0; + //!This should be called at least once before the first call to minimize() and step is assumed to be (something that + //! after getMat() will return) row-vector or column-vector. *It's dimensionality determines the dimensionality of a problem.* + virtual void setInitStep(InputArray step)=0; +}; + +// both minRange & minError are specified by termcrit.epsilon; In addition, user may specify the number of iterations that the algorithm does. +CV_EXPORTS_W Ptr createDownhillSolver(const Ptr& f=Ptr(), + InputArray initStep=Mat_(1,1,0.0), + TermCriteria termcrit=TermCriteria(TermCriteria::MAX_ITER+TermCriteria::EPS,5000,0.000001)); + //!the return codes for solveLP() function enum { diff --git a/modules/optim/src/debug.hpp b/modules/optim/src/debug.hpp new file mode 100644 index 000000000..fe5d00e87 --- /dev/null +++ b/modules/optim/src/debug.hpp @@ -0,0 +1,18 @@ +namespace cv{namespace optim{ +#ifdef ALEX_DEBUG +#define dprintf(x) printf x +static void print_matrix(const Mat& x){ + printf("\ttype:%d vs %d,\tsize: %d-on-%d\n",x.type(),CV_64FC1,x.rows,x.cols); + for(int i=0;i(i,j)); + } + printf("]\n"); + } +} +#else +#define dprintf(x) +#define print_matrix(x) +#endif +}} diff --git a/modules/optim/src/lpsolver.cpp b/modules/optim/src/lpsolver.cpp index 924b756a9..a046ddae1 100644 --- a/modules/optim/src/lpsolver.cpp +++ b/modules/optim/src/lpsolver.cpp @@ -2,16 +2,12 @@ #include #include #include +#include namespace cv{namespace optim{ using std::vector; #ifdef ALEX_DEBUG -#define dprintf(x) printf x -static void print_matrix(const Mat& x){ - print(x); - printf("\n"); -} static void print_simplex_state(const Mat& c,const Mat& b,double v,const std::vector N,const std::vector B){ printf("\tprint simplex state\n"); @@ -32,8 +28,6 @@ static void print_simplex_state(const Mat& c,const Mat& b,double v,const std::ve printf("\n"); } #else -#define dprintf(x) -#define print_matrix(x) #define print_simplex_state(c,b,v,N,B) #endif diff --git a/modules/optim/src/simplex.cpp b/modules/optim/src/simplex.cpp new file mode 100644 index 000000000..f45d0ce0b --- /dev/null +++ b/modules/optim/src/simplex.cpp @@ -0,0 +1,273 @@ +#include "precomp.hpp" +#include "debug.hpp" +#include "opencv2/core/core_c.h" + +namespace cv{namespace optim{ + + class DownhillSolverImpl : public DownhillSolver + { + public: + void getInitStep(OutputArray step) const; + void setInitStep(InputArray step); + Ptr getFunction() const; + void setFunction(const Ptr& f); + TermCriteria getTermCriteria() const; + DownhillSolverImpl(); + void setTermCriteria(const TermCriteria& termcrit); + double minimize(InputOutputArray x); + protected: + Ptr _Function; + TermCriteria _termcrit; + Mat _step; + private: + inline void createInitialSimplex(Mat_& simplex,Mat& step); + inline double innerDownhillSimplex(cv::Mat_& p,double MinRange,double MinError,int& nfunk, + const Ptr& f,int nmax); + inline double tryNewPoint(Mat_& p,Mat_& y,Mat_& coord_sum,const Ptr& f,int ihi, + double fac,Mat_& ptry); + }; + + double DownhillSolverImpl::tryNewPoint( + Mat_& p, + Mat_& y, + Mat_& coord_sum, + const Ptr& f, + int ihi, + double fac, + Mat_& ptry + ) + { + int ndim=p.cols; + int j; + double fac1,fac2,ytry; + + fac1=(1.0-fac)/ndim; + fac2=fac1-fac; + for (j=0;jcalc((double*)ptry.data); + if (ytry < y(ihi)) + { + y(ihi)=ytry; + for (j=0;j& p, + double MinRange, + double MinError, + int& nfunk, + const Ptr& f, + int nmax + ) + { + int ndim=p.cols; + double res; + int i,ihi,ilo,inhi,j,mpts=ndim+1; + double error, range,ysave,ytry; + Mat_ coord_sum(1,ndim,0.0),buf(1,ndim,0.0),y(1,ndim,0.0); + + nfunk = 0; + + for(i=0;icalc(p[i]); + } + + nfunk = ndim+1; + + reduce(p,coord_sum,0,CV_REDUCE_SUM); + + for (;;) + { + ilo=0; + /* find highest (worst), next-to-worst, and lowest + (best) points by going through all of them. */ + ihi = y(0)>y(1) ? (inhi=1,0) : (inhi=0,1); + for (i=0;i y(ihi)) + { + inhi=ihi; + ihi=i; + } + else if (y(i) > y(inhi) && i != ihi) + inhi=i; + } + + /* check stop criterion */ + error=fabs(y(ihi)-y(ilo)); + range=0; + for(i=0;i p(j,i) ) min = p(j,i); + if( max < p(j,i) ) max = p(j,i); + } + d = fabs(max-min); + if(range < d) range = d; + } + + if(range <= MinRange || error <= MinError) + { /* Put best point and value in first slot. */ + std::swap(y(0),y(ilo)); + for (i=0;i= nmax){ + dprintf(("nmax exceeded\n")); + return y(ilo); + } + nfunk += 2; + /*Begin a new iteration. First, reflect the worst point about the centroid of others */ + ytry = tryNewPoint(p,y,coord_sum,f,ihi,-1.0,buf); + if (ytry <= y(ilo)) + { /*If that's better than the best point, go twice as far in that direction*/ + ytry = tryNewPoint(p,y,coord_sum,f,ihi,2.0,buf); + } + else if (ytry >= y(inhi)) + { /* The new point is worse than the second-highest, but better + than the worst so do not go so far in that direction */ + ysave = y(ihi); + ytry = tryNewPoint(p,y,coord_sum,f,ihi,0.5,buf); + if (ytry >= ysave) + { /* Can't seem to improve things. Contract the simplex to good point + in hope to find a simplex landscape. */ + for (i=0;icalc((double*)coord_sum.data); + } + } + nfunk += ndim; + reduce(p,coord_sum,0,CV_REDUCE_SUM); + } + } else --(nfunk); /* correct nfunk */ + dprintf(("this is simplex on iteration %d\n",nfunk)); + print_matrix(p); + } /* go to next iteration. */ + res = y(0); + + return res; + } + + void DownhillSolverImpl::createInitialSimplex(Mat_& simplex,Mat& step){ + for(int i=1;i<=step.cols;++i) + { + simplex.row(0).copyTo(simplex.row(i)); + simplex(i,i-1)+= 0.5*step.at(0,i-1); + } + simplex.row(0) -= 0.5*step; + + dprintf(("this is simplex\n")); + print_matrix(simplex); + } + + double DownhillSolverImpl::minimize(InputOutputArray x){ + dprintf(("hi from minimize\n")); + CV_Assert(_Function.empty()==false); + dprintf(("termcrit:\n\ttype: %d\n\tmaxCount: %d\n\tEPS: %g\n",_termcrit.type,_termcrit.maxCount,_termcrit.epsilon)); + dprintf(("step\n")); + print_matrix(_step); + + Mat x_mat=x.getMat(); + CV_Assert(MIN(x_mat.rows,x_mat.cols)==1); + CV_Assert(MAX(x_mat.rows,x_mat.cols)==_step.cols); + CV_Assert(x_mat.type()==CV_64FC1); + + Mat_ proxy_x; + + if(x_mat.rows>1){ + proxy_x=x_mat.t(); + }else{ + proxy_x=x_mat; + } + + int count=0; + int ndim=_step.cols; + Mat_ simplex=Mat_(ndim+1,ndim,0.0); + + simplex.row(0).copyTo(proxy_x); + createInitialSimplex(simplex,_step); + double res = innerDownhillSimplex( + simplex,_termcrit.epsilon, _termcrit.epsilon, count,_Function,_termcrit.maxCount); + simplex.row(0).copyTo(proxy_x); + + dprintf(("%d iterations done\n",count)); + + if(x_mat.rows>1){ + Mat(x_mat.rows, 1, CV_64F, (double*)proxy_x.data).copyTo(x); + } + return res; + } + DownhillSolverImpl::DownhillSolverImpl(){ + _Function=Ptr(); + _step=Mat_(); + } + Ptr DownhillSolverImpl::getFunction()const{ + return _Function; + } + void DownhillSolverImpl::setFunction(const Ptr& f){ + _Function=f; + } + TermCriteria DownhillSolverImpl::getTermCriteria()const{ + return _termcrit; + } + void DownhillSolverImpl::setTermCriteria(const TermCriteria& termcrit){ + CV_Assert(termcrit.type==(TermCriteria::MAX_ITER+TermCriteria::EPS) && termcrit.epsilon>0 && termcrit.maxCount>0); + _termcrit=termcrit; + } + // both minRange & minError are specified by termcrit.epsilon; In addition, user may specify the number of iterations that the algorithm does. + Ptr createDownhillSolver(const Ptr& f, InputArray initStep, TermCriteria termcrit){ + DownhillSolver *DS=new DownhillSolverImpl(); + DS->setFunction(f); + DS->setInitStep(initStep); + DS->setTermCriteria(termcrit); + return Ptr(DS); + } + void DownhillSolverImpl::getInitStep(OutputArray step)const{ + _step.copyTo(step); + } + void DownhillSolverImpl::setInitStep(InputArray step){ + //set dimensionality and make a deep copy of step + Mat m=step.getMat(); + dprintf(("m.cols=%d\nm.rows=%d\n",m.cols,m.rows)); + CV_Assert(MIN(m.cols,m.rows)==1 && m.type()==CV_64FC1); + if(m.rows==1){ + m.copyTo(_step); + }else{ + transpose(m,_step); + } + } +}} diff --git a/modules/optim/test/test_downhill_simplex.cpp b/modules/optim/test/test_downhill_simplex.cpp new file mode 100644 index 000000000..95b2c6e9e --- /dev/null +++ b/modules/optim/test/test_downhill_simplex.cpp @@ -0,0 +1,63 @@ +#include "test_precomp.hpp" +#include +#include +#include + +static void mytest(cv::Ptr solver,cv::Ptr ptr_F,cv::Mat& x,cv::Mat& step, + cv::Mat& etalon_x,double etalon_res){ + solver->setFunction(ptr_F); + int ndim=MAX(step.cols,step.rows); + solver->setInitStep(step); + cv::Mat settedStep; + solver->getInitStep(settedStep); + ASSERT_TRUE(settedStep.rows==1 && settedStep.cols==ndim); + ASSERT_TRUE(std::equal(step.begin(),step.end(),settedStep.begin())); + std::cout<<"step setted:\n\t"<minimize(x); + std::cout<<"res:\n\t"<getTermCriteria().epsilon; + ASSERT_TRUE(std::abs(res-etalon_res)::iterator it1=x.begin(),it2=etalon_x.begin();it1!=x.end();it1++,it2++){ + ASSERT_TRUE(std::abs((*it1)-(*it2)) solver=cv::optim::createDownhillSolver(); +#if 1 + { + cv::Ptr ptr_F(new SphereF()); + cv::Mat x=(cv::Mat_(1,2)<<1.0,1.0), + step=(cv::Mat_(2,1)<<-0.5,-0.5), + etalon_x=(cv::Mat_(1,2)<<-0.0,0.0); + double etalon_res=0.0; + mytest(solver,ptr_F,x,step,etalon_x,etalon_res); + } +#endif +#if 1 + { + cv::Ptr ptr_F(new RosenbrockF()); + cv::Mat x=(cv::Mat_(2,1)<<0.0,0.0), + step=(cv::Mat_(2,1)<<0.5,+0.5), + etalon_x=(cv::Mat_(2,1)<<1.0,1.0); + double etalon_res=0.0; + mytest(solver,ptr_F,x,step,etalon_x,etalon_res); + } +#endif +}