optim.hpp 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302
  1. /*M///////////////////////////////////////////////////////////////////////////////////////
  2. //
  3. // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
  4. //
  5. // By downloading, copying, installing or using the software you agree to this license.
  6. // If you do not agree to this license, do not download, install,
  7. // copy or use the software.
  8. //
  9. //
  10. // License Agreement
  11. // For Open Source Computer Vision Library
  12. //
  13. // Copyright (C) 2013, OpenCV Foundation, all rights reserved.
  14. // Third party copyrights are property of their respective owners.
  15. //
  16. // Redistribution and use in source and binary forms, with or without modification,
  17. // are permitted provided that the following conditions are met:
  18. //
  19. // * Redistribution's of source code must retain the above copyright notice,
  20. // this list of conditions and the following disclaimer.
  21. //
  22. // * Redistribution's in binary form must reproduce the above copyright notice,
  23. // this list of conditions and the following disclaimer in the documentation
  24. // and/or other materials provided with the distribution.
  25. //
  26. // * The name of the copyright holders may not be used to endorse or promote products
  27. // derived from this software without specific prior written permission.
  28. //
  29. // This software is provided by the copyright holders and contributors "as is" and
  30. // any express or implied warranties, including, but not limited to, the implied
  31. // warranties of merchantability and fitness for a particular purpose are disclaimed.
  32. // In no event shall the OpenCV Foundation or contributors be liable for any direct,
  33. // indirect, incidental, special, exemplary, or consequential damages
  34. // (including, but not limited to, procurement of substitute goods or services;
  35. // loss of use, data, or profits; or business interruption) however caused
  36. // and on any theory of liability, whether in contract, strict liability,
  37. // or tort (including negligence or otherwise) arising in any way out of
  38. // the use of this software, even if advised of the possibility of such damage.
  39. //
  40. //M*/
  41. #ifndef OPENCV_OPTIM_HPP
  42. #define OPENCV_OPTIM_HPP
  43. #include "opencv2/core.hpp"
  44. namespace cv
  45. {
  46. /** @addtogroup core_optim
  47. The algorithms in this section minimize or maximize function value within specified constraints or
  48. without any constraints.
  49. @{
  50. */
  51. /** @brief Basic interface for all solvers
  52. */
  53. class CV_EXPORTS MinProblemSolver : public Algorithm
  54. {
  55. public:
  56. /** @brief Represents function being optimized
  57. */
  58. class CV_EXPORTS Function
  59. {
  60. public:
  61. virtual ~Function() {}
  62. virtual int getDims() const = 0;
  63. virtual double getGradientEps() const;
  64. virtual double calc(const double* x) const = 0;
  65. virtual void getGradient(const double* x,double* grad);
  66. };
  67. /** @brief Getter for the optimized function.
  68. The optimized function is represented by Function interface, which requires derivatives to
  69. implement the calc(double*) and getDim() methods to evaluate the function.
  70. @return Smart-pointer to an object that implements Function interface - it represents the
  71. function that is being optimized. It can be empty, if no function was given so far.
  72. */
  73. virtual Ptr<Function> getFunction() const = 0;
  74. /** @brief Setter for the optimized function.
  75. *It should be called at least once before the call to* minimize(), as default value is not usable.
  76. @param f The new function to optimize.
  77. */
  78. virtual void setFunction(const Ptr<Function>& f) = 0;
  79. /** @brief Getter for the previously set terminal criteria for this algorithm.
  80. @return Deep copy of the terminal criteria used at the moment.
  81. */
  82. virtual TermCriteria getTermCriteria() const = 0;
  83. /** @brief Set terminal criteria for solver.
  84. This method *is not necessary* to be called before the first call to minimize(), as the default
  85. value is sensible.
  86. Algorithm stops when the number of function evaluations done exceeds termcrit.maxCount, when
  87. the function values at the vertices of simplex are within termcrit.epsilon range or simplex
  88. becomes so small that it can enclosed in a box with termcrit.epsilon sides, whatever comes
  89. first.
  90. @param termcrit Terminal criteria to be used, represented as cv::TermCriteria structure.
  91. */
  92. virtual void setTermCriteria(const TermCriteria& termcrit) = 0;
  93. /** @brief actually runs the algorithm and performs the minimization.
  94. The sole input parameter determines the centroid of the starting simplex (roughly, it tells
  95. where to start), all the others (terminal criteria, initial step, function to be minimized) are
  96. supposed to be set via the setters before the call to this method or the default values (not
  97. always sensible) will be used.
  98. @param x The initial point, that will become a centroid of an initial simplex. After the algorithm
  99. will terminate, it will be set to the point where the algorithm stops, the point of possible
  100. minimum.
  101. @return The value of a function at the point found.
  102. */
  103. virtual double minimize(InputOutputArray x) = 0;
  104. };
  105. /** @brief This class is used to perform the non-linear non-constrained minimization of a function,
  106. defined on an `n`-dimensional Euclidean space, using the **Nelder-Mead method**, also known as
  107. **downhill simplex method**. The basic idea about the method can be obtained from
  108. <http://en.wikipedia.org/wiki/Nelder-Mead_method>.
  109. It should be noted, that this method, although deterministic, is rather a heuristic and therefore
  110. may converge to a local minima, not necessary a global one. It is iterative optimization technique,
  111. which at each step uses an information about the values of a function evaluated only at `n+1`
  112. points, arranged as a *simplex* in `n`-dimensional space (hence the second name of the method). At
  113. each step new point is chosen to evaluate function at, obtained value is compared with previous
  114. ones and based on this information simplex changes it's shape , slowly moving to the local minimum.
  115. Thus this method is using *only* function values to make decision, on contrary to, say, Nonlinear
  116. Conjugate Gradient method (which is also implemented in optim).
  117. Algorithm stops when the number of function evaluations done exceeds termcrit.maxCount, when the
  118. function values at the vertices of simplex are within termcrit.epsilon range or simplex becomes so
  119. small that it can enclosed in a box with termcrit.epsilon sides, whatever comes first, for some
  120. defined by user positive integer termcrit.maxCount and positive non-integer termcrit.epsilon.
  121. @note DownhillSolver is a derivative of the abstract interface
  122. cv::MinProblemSolver, which in turn is derived from the Algorithm interface and is used to
  123. encapsulate the functionality, common to all non-linear optimization algorithms in the optim
  124. module.
  125. @note term criteria should meet following condition:
  126. @code
  127. termcrit.type == (TermCriteria::MAX_ITER + TermCriteria::EPS) && termcrit.epsilon > 0 && termcrit.maxCount > 0
  128. @endcode
  129. */
  130. class CV_EXPORTS DownhillSolver : public MinProblemSolver
  131. {
  132. public:
  133. /** @brief Returns the initial step that will be used in downhill simplex algorithm.
  134. @param step Initial step that will be used in algorithm. Note, that although corresponding setter
  135. accepts column-vectors as well as row-vectors, this method will return a row-vector.
  136. @see DownhillSolver::setInitStep
  137. */
  138. virtual void getInitStep(OutputArray step) const=0;
  139. /** @brief Sets the initial step that will be used in downhill simplex algorithm.
  140. Step, together with initial point (givin in DownhillSolver::minimize) are two `n`-dimensional
  141. vectors that are used to determine the shape of initial simplex. Roughly said, initial point
  142. determines the position of a simplex (it will become simplex's centroid), while step determines the
  143. spread (size in each dimension) of a simplex. To be more precise, if \f$s,x_0\in\mathbb{R}^n\f$ are
  144. the initial step and initial point respectively, the vertices of a simplex will be:
  145. \f$v_0:=x_0-\frac{1}{2} s\f$ and \f$v_i:=x_0+s_i\f$ for \f$i=1,2,\dots,n\f$ where \f$s_i\f$ denotes
  146. projections of the initial step of *n*-th coordinate (the result of projection is treated to be
  147. vector given by \f$s_i:=e_i\cdot\left<e_i\cdot s\right>\f$, where \f$e_i\f$ form canonical basis)
  148. @param step Initial step that will be used in algorithm. Roughly said, it determines the spread
  149. (size in each dimension) of an initial simplex.
  150. */
  151. virtual void setInitStep(InputArray step)=0;
  152. /** @brief This function returns the reference to the ready-to-use DownhillSolver object.
  153. All the parameters are optional, so this procedure can be called even without parameters at
  154. all. In this case, the default values will be used. As default value for terminal criteria are
  155. the only sensible ones, MinProblemSolver::setFunction() and DownhillSolver::setInitStep()
  156. should be called upon the obtained object, if the respective parameters were not given to
  157. create(). Otherwise, the two ways (give parameters to createDownhillSolver() or miss them out
  158. and call the MinProblemSolver::setFunction() and DownhillSolver::setInitStep()) are absolutely
  159. equivalent (and will drop the same errors in the same way, should invalid input be detected).
  160. @param f Pointer to the function that will be minimized, similarly to the one you submit via
  161. MinProblemSolver::setFunction.
  162. @param initStep Initial step, that will be used to construct the initial simplex, similarly to the one
  163. you submit via MinProblemSolver::setInitStep.
  164. @param termcrit Terminal criteria to the algorithm, similarly to the one you submit via
  165. MinProblemSolver::setTermCriteria.
  166. */
  167. static Ptr<DownhillSolver> create(const Ptr<MinProblemSolver::Function>& f=Ptr<MinProblemSolver::Function>(),
  168. InputArray initStep=Mat_<double>(1,1,0.0),
  169. TermCriteria termcrit=TermCriteria(TermCriteria::MAX_ITER+TermCriteria::EPS,5000,0.000001));
  170. };
  171. /** @brief This class is used to perform the non-linear non-constrained minimization of a function
  172. with known gradient,
  173. defined on an *n*-dimensional Euclidean space, using the **Nonlinear Conjugate Gradient method**.
  174. The implementation was done based on the beautifully clear explanatory article [An Introduction to
  175. the Conjugate Gradient Method Without the Agonizing
  176. Pain](http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf) by Jonathan Richard
  177. Shewchuk. The method can be seen as an adaptation of a standard Conjugate Gradient method (see, for
  178. example <http://en.wikipedia.org/wiki/Conjugate_gradient_method>) for numerically solving the
  179. systems of linear equations.
  180. It should be noted, that this method, although deterministic, is rather a heuristic method and
  181. therefore may converge to a local minima, not necessary a global one. What is even more disastrous,
  182. most of its behaviour is ruled by gradient, therefore it essentially cannot distinguish between
  183. local minima and maxima. Therefore, if it starts sufficiently near to the local maximum, it may
  184. converge to it. Another obvious restriction is that it should be possible to compute the gradient of
  185. a function at any point, thus it is preferable to have analytic expression for gradient and
  186. computational burden should be born by the user.
  187. The latter responsibility is accomplished via the getGradient method of a
  188. MinProblemSolver::Function interface (which represents function being optimized). This method takes
  189. point a point in *n*-dimensional space (first argument represents the array of coordinates of that
  190. point) and compute its gradient (it should be stored in the second argument as an array).
  191. @note class ConjGradSolver thus does not add any new methods to the basic MinProblemSolver interface.
  192. @note term criteria should meet following condition:
  193. @code
  194. termcrit.type == (TermCriteria::MAX_ITER + TermCriteria::EPS) && termcrit.epsilon > 0 && termcrit.maxCount > 0
  195. // or
  196. termcrit.type == TermCriteria::MAX_ITER) && termcrit.maxCount > 0
  197. @endcode
  198. */
  199. class CV_EXPORTS ConjGradSolver : public MinProblemSolver
  200. {
  201. public:
  202. /** @brief This function returns the reference to the ready-to-use ConjGradSolver object.
  203. All the parameters are optional, so this procedure can be called even without parameters at
  204. all. In this case, the default values will be used. As default value for terminal criteria are
  205. the only sensible ones, MinProblemSolver::setFunction() should be called upon the obtained
  206. object, if the function was not given to create(). Otherwise, the two ways (submit it to
  207. create() or miss it out and call the MinProblemSolver::setFunction()) are absolutely equivalent
  208. (and will drop the same errors in the same way, should invalid input be detected).
  209. @param f Pointer to the function that will be minimized, similarly to the one you submit via
  210. MinProblemSolver::setFunction.
  211. @param termcrit Terminal criteria to the algorithm, similarly to the one you submit via
  212. MinProblemSolver::setTermCriteria.
  213. */
  214. static Ptr<ConjGradSolver> create(const Ptr<MinProblemSolver::Function>& f=Ptr<ConjGradSolver::Function>(),
  215. TermCriteria termcrit=TermCriteria(TermCriteria::MAX_ITER+TermCriteria::EPS,5000,0.000001));
  216. };
  217. //! return codes for cv::solveLP() function
  218. enum SolveLPResult
  219. {
  220. SOLVELP_UNBOUNDED = -2, //!< problem is unbounded (target function can achieve arbitrary high values)
  221. SOLVELP_UNFEASIBLE = -1, //!< problem is unfeasible (there are no points that satisfy all the constraints imposed)
  222. SOLVELP_SINGLE = 0, //!< there is only one maximum for target function
  223. SOLVELP_MULTI = 1 //!< there are multiple maxima for target function - the arbitrary one is returned
  224. };
  225. /** @brief Solve given (non-integer) linear programming problem using the Simplex Algorithm (Simplex Method).
  226. What we mean here by "linear programming problem" (or LP problem, for short) can be formulated as:
  227. \f[\mbox{Maximize } c\cdot x\\
  228. \mbox{Subject to:}\\
  229. Ax\leq b\\
  230. x\geq 0\f]
  231. Where \f$c\f$ is fixed `1`-by-`n` row-vector, \f$A\f$ is fixed `m`-by-`n` matrix, \f$b\f$ is fixed `m`-by-`1`
  232. column vector and \f$x\f$ is an arbitrary `n`-by-`1` column vector, which satisfies the constraints.
  233. Simplex algorithm is one of many algorithms that are designed to handle this sort of problems
  234. efficiently. Although it is not optimal in theoretical sense (there exist algorithms that can solve
  235. any problem written as above in polynomial time, while simplex method degenerates to exponential
  236. time for some special cases), it is well-studied, easy to implement and is shown to work well for
  237. real-life purposes.
  238. The particular implementation is taken almost verbatim from **Introduction to Algorithms, third
  239. edition** by T. H. Cormen, C. E. Leiserson, R. L. Rivest and Clifford Stein. In particular, the
  240. Bland's rule <http://en.wikipedia.org/wiki/Bland%27s_rule> is used to prevent cycling.
  241. @param Func This row-vector corresponds to \f$c\f$ in the LP problem formulation (see above). It should
  242. contain 32- or 64-bit floating point numbers. As a convenience, column-vector may be also submitted,
  243. in the latter case it is understood to correspond to \f$c^T\f$.
  244. @param Constr `m`-by-`n+1` matrix, whose rightmost column corresponds to \f$b\f$ in formulation above
  245. and the remaining to \f$A\f$. It should contain 32- or 64-bit floating point numbers.
  246. @param z The solution will be returned here as a column-vector - it corresponds to \f$c\f$ in the
  247. formulation above. It will contain 64-bit floating point numbers.
  248. @return One of cv::SolveLPResult
  249. */
  250. CV_EXPORTS_W int solveLP(InputArray Func, InputArray Constr, OutputArray z);
  251. //! @}
  252. }// cv
  253. #endif