From a95650111f813d99f82f672b3caaee7f9fa3df04 Mon Sep 17 00:00:00 2001 From: Alex Leontiev Date: Wed, 10 Jul 2013 20:11:52 +0300 Subject: [PATCH] Cleaning the code of simplex method In particular, the following things are done: *) Consistent tabulation of 4 spaces is ensured *) New function dprintf() is introduced, so now printing of the debug information can be turned on/off via the ALEX_DEBUG macro *) Removed solveLP_aux namespace *) All auxiliary functions are declared as static *) The return codes of solveLP() are encapsulated in enum. --- modules/optim/doc/linear_programming.rst | 48 ++++++++++ modules/optim/doc/optim.rst | 10 ++ modules/optim/include/opencv2/optim.hpp | 94 ++++++++++--------- modules/optim/src/lpsolver.cpp | 151 ++++++++++++++++--------------- modules/optim/test/test_lpsolver.cpp | 8 +- 5 files changed, 193 insertions(+), 118 deletions(-) create mode 100644 modules/optim/doc/linear_programming.rst create mode 100644 modules/optim/doc/optim.rst diff --git a/modules/optim/doc/linear_programming.rst b/modules/optim/doc/linear_programming.rst new file mode 100644 index 0000000..fbc1d04 --- /dev/null +++ b/modules/optim/doc/linear_programming.rst @@ -0,0 +1,48 @@ +Linear Programming +================== + +.. highlight:: cpp + +optim::solveLP +-------------------- +Solve given (non-integer) linear programming problem using the Simplex Algorithm (Simplex Method). +What we mean here by "linear programming problem" (or LP problem, for short) can be +formulated as: + +.. math:: + \mbox{Maximize } c\cdot x\\ + \mbox{Subject to:}\\ + Ax\leq b\\ + x\geq 0 + +Where :math:`c` is fixed *1*-by-*n* row-vector, :math:`A` is fixed *m*-by-*n* matrix, :math:`b` is fixed *m*-by-*1* column vector and +:math:`x` is an arbitrary *n*-by-*1* column vector, which satisfies the constraints. + +Simplex algorithm is one of many algorithms that are designed to handle this sort of problems efficiently. Although it is not optimal in theoretical +sense (there exist algorithms that can solve any problem written as above in polynomial type, while simplex method degenerates to exponential time +for some special cases), it is well-studied, easy to implement and is shown to work well for real-life purposes. + +The particular implementation is taken almost verbatim from **Introduction to Algorithms, third edition** +by T. H. Cormen, C. E. Leiserson, R. L. Rivest and Clifford Stein. In particular, the Bland's rule +(`http://en.wikipedia.org/wiki/Bland%27s\_rule `_) is used to prevent cycling. + +.. ocv:function:: int optim::solveLP(const Mat& Func, const Mat& Constr, Mat& z) + + :param Func: This row-vector corresponds to :math:`c` in the LP problem formulation (see above). + + :param Constr: *m*-by-*n\+1* matrix, whose rightmost column corresponds to :math:`b` in formulation above and the remaining to :math:`A`. + + :param z: The solution will be returned here as a row-vector - it corresponds to (transposed) :math:`c` in the formulation above. + + :return: One of the return codes: + +:: + + //!the return codes for solveLP() function + enum + { + SOLVELP_UNBOUNDED = -2, //problem is unbounded (target function can achieve arbitrary high values) + SOLVELP_UNFEASIBLE = -1, //problem is unfeasible (there are no points that satisfy all the constraints imposed) + SOLVELP_SINGLE = 0, //there is only one maximum for target function + SOLVELP_MULTI = 1 //there are multiple maxima for target function - the arbitrary one is returned + }; diff --git a/modules/optim/doc/optim.rst b/modules/optim/doc/optim.rst new file mode 100644 index 0000000..bead212 --- /dev/null +++ b/modules/optim/doc/optim.rst @@ -0,0 +1,10 @@ +************************************** +optim. Generic numerical optimization +************************************** + +.. highlight:: cpp + +.. toctree:: + :maxdepth: 2 + + linear_programming diff --git a/modules/optim/include/opencv2/optim.hpp b/modules/optim/include/opencv2/optim.hpp index f4a639f..f40456c 100644 --- a/modules/optim/include/opencv2/optim.hpp +++ b/modules/optim/include/opencv2/optim.hpp @@ -47,67 +47,77 @@ #include "opencv2/core.hpp" #include "opencv2/core/mat.hpp" -/*! \namespace cv - Namespace where all the C++ OpenCV functionality resides - */ +//uncomment the next line to print the debug info +//#define ALEX_DEBUG + namespace cv{namespace optim { //! generic class for optimization algorithms */ class CV_EXPORTS Solver : public Algorithm /* Algorithm is the base OpenCV class */ { - public: - class CV_EXPORTS Function - { - public: + public: + class CV_EXPORTS Function + { + public: virtual ~Function(){} virtual double calc(InputArray args) const = 0; - }; - class CV_EXPORTS Constraints - { - public: + }; + class CV_EXPORTS Constraints + { + public: virtual ~Constraints(){} - }; + }; - //! could be reused for all the generic algorithms like downhill simplex. Return value is the maximum value of a function*/ - virtual double solve(const Function& F,const Constraints& C, OutputArray result) const = 0; + //! could be reused for all the generic algorithms like downhill simplex. Return value is the maximum value of a function*/ + virtual double solve(const Function& F,const Constraints& C, OutputArray result) const = 0; - /*virtual void setTermCriteria(const TermCriteria& criteria) = 0; - virtual TermCriteria getTermCriteria() = 0;*/ + /*virtual void setTermCriteria(const TermCriteria& criteria) = 0; + virtual TermCriteria getTermCriteria() = 0;*/ - // more detailed API to be defined later ... + // more detailed API to be defined later ... }; class CV_EXPORTS LPSolver : public Solver { public: - class CV_EXPORTS LPFunction:public Solver::Function - { - Mat z; - public: - //! Note, that this class is supposed to be immutable, so it's ok to make only a shallow copy of z_in.*/ - LPFunction(Mat z_in):z(z_in){} - ~LPFunction(){}; - const Mat& getz()const{return z;} - double calc(InputArray args)const; - }; + class CV_EXPORTS LPFunction:public Solver::Function + { + Mat z; + public: + //! Note, that this class is supposed to be immutable, so it's ok to make only a shallow copy of z_in.*/ + LPFunction(Mat z_in):z(z_in){} + ~LPFunction(){}; + const Mat& getz()const{return z;} + double calc(InputArray args)const; + }; - //!This class represents constraints for linear problem. There are two matrix stored: m-by-n matrix A and n-by-1 column-vector b. - //!What this represents is the set of constraints Ax\leq b and x\geq 0. It can be shown that any set of linear constraints can be converted - //!this form and **we shall create various constructors for this class that will perform these conversions**. - class CV_EXPORTS LPConstraints:public Solver::Constraints - { - Mat A,b; - public: - ~LPConstraints(){}; - //! Note, that this class is supposed to be immutable, so it's ok to make only a shallow copy of A_in and b_in.*/ - LPConstraints(Mat A_in, Mat b_in):A(A_in),b(b_in){} - const Mat& getA()const{return A;} - const Mat& getb()const{return b;} - }; + //!This class represents constraints for linear problem. There are two matrix stored: m-by-n matrix A and n-by-1 column-vector b. + //!What this represents is the set of constraints Ax\leq b and x\geq 0. It can be shown that any set of linear constraints can be converted + //!this form and **we shall create various constructors for this class that will perform these conversions**. + class CV_EXPORTS LPConstraints:public Solver::Constraints + { + Mat A,b; + public: + ~LPConstraints(){}; + //! Note, that this class is supposed to be immutable, so it's ok to make only a shallow copy of A_in and b_in.*/ + LPConstraints(Mat A_in, Mat b_in):A(A_in),b(b_in){} + const Mat& getA()const{return A;} + const Mat& getb()const{return b;} + }; - LPSolver(){} - double solve(const Function& F,const Constraints& C, OutputArray result)const; + LPSolver(){} + double solve(const Function& F,const Constraints& C, OutputArray result)const; }; + +//!the return codes for solveLP() function +enum +{ + SOLVELP_UNBOUNDED = -2, //problem is unbounded (target function can achieve arbitrary high values) + SOLVELP_UNFEASIBLE = -1, //problem is unfeasible (there are no points that satisfy all the constraints imposed) + SOLVELP_SINGLE = 0, //there is only one maximum for target function + SOLVELP_MULTI = 1 //there are multiple maxima for target function - the arbitrary one is returned +}; + CV_EXPORTS_W int solveLP(const Mat& Func, const Mat& Constr, Mat& z); }}// cv diff --git a/modules/optim/src/lpsolver.cpp b/modules/optim/src/lpsolver.cpp index b56dc3e..1e80dfa 100644 --- a/modules/optim/src/lpsolver.cpp +++ b/modules/optim/src/lpsolver.cpp @@ -2,73 +2,87 @@ #include "precomp.hpp" #include #include +#include namespace cv{namespace optim{ using std::vector; -double LPSolver::solve(const Function& F,const Constraints& C, OutputArray result)const{ +const void dprintf(const char* format,...){ +#ifdef ALEX_DEBUG + va_list args; + va_start (args,format); + vprintf(format,args); + va_end(args); +#endif +} +double LPSolver::solve(const Function& F,const Constraints& C, OutputArray result)const{ return 0.0; } double LPSolver::LPFunction::calc(InputArray args)const{ - printf("call to LPFunction::calc()\n"); + dprintf("call to LPFunction::calc()\n"); return 0.0; } -void print_matrix(const Mat& X){ - printf("\ttype:%d vs %d,\tsize: %d-on-%d\n",X.type(),CV_64FC1,X.rows,X.cols); + + +void const print_matrix(const Mat& X){ +#ifdef ALEX_DEBUG + dprintf("\ttype:%d vs %d,\tsize: %d-on-%d\n",X.type(),CV_64FC1,X.rows,X.cols); for(int i=0;i(i,j)); + dprintf("%g, ",X.at(i,j)); } - printf("]\n"); + dprintf("]\n"); } +#endif } -void print_simplex_state(const Mat& c,const Mat&b,double v,const vector& N,const vector& B){ - printf("\tprint simplex state\n"); - printf("v=%g\n",v); +void const print_simplex_state(const Mat& c,const Mat&b,double v,const vector& N,const vector& B){ +#ifdef ALEX_DEBUG + dprintf("\tprint simplex state\n"); - printf("here c goes\n"); + dprintf("v=%g\n",v); + + dprintf("here c goes\n"); print_matrix(c); - printf("non-basic: "); + dprintf("non-basic: "); for (std::vector::const_iterator it = N.begin() ; it != N.end(); ++it){ - printf("%d, ",*it); + dprintf("%d, ",*it); } - printf("\n"); + dprintf("\n"); - printf("here b goes\n"); + dprintf("here b goes\n"); print_matrix(b); - printf("basic: "); + dprintf("basic: "); for (std::vector::const_iterator it = B.begin() ; it != B.end(); ++it){ - printf("%d, ",*it); + dprintf("%d, ",*it); } - printf("\n"); + dprintf("\n"); +#endif } -namespace solveLP_aux{ - /**Due to technical considerations, the format of input b and c is somewhat special: - *both b and c should be one column bigger than corresponding b and c of linear problem and the leftmost column will be used internally - by this procedure - it should not be cleaned before the call to procedure and may contain mess after - it also initializes N and B and does not make any assumptions about their init values - * @return -1 if problem is unfeasible, 0 if feasible. - */ - int initialize_simplex(Mat_& c, Mat_& b,double& v,vector& N,vector& B); - inline void pivot(Mat_& c,Mat_& b,double& v,vector& N,vector& B, int leaving_index,int entering_index); - /**@return -2 means the problem is unbdd, 1 means multiple solutions, 0 means successful. - */ - int inner_simplex(Mat_& c, Mat_& b,double& v,vector& N,vector& B); - void swap_columns(Mat_& A,int col1,int col2); -} +/**Due to technical considerations, the format of input b and c is somewhat special: + *both b and c should be one column bigger than corresponding b and c of linear problem and the leftmost column will be used internally + by this procedure - it should not be cleaned before the call to procedure and may contain mess after + it also initializes N and B and does not make any assumptions about their init values + * @return SOLVELP_UNFEASIBLE if problem is unfeasible, 0 if feasible. +*/ +const int initialize_simplex(Mat_& c, Mat_& b,double& v,vector& N,vector& B); +const inline void pivot(Mat_& c,Mat_& b,double& v,vector& N,vector& B, int leaving_index,int entering_index); +/**@return SOLVELP_UNBOUNDED means the problem is unbdd, SOLVELP_MULTI means multiple solutions, SOLVELP_SINGLE means one solution. + */ +const int inner_simplex(Mat_& c, Mat_& b,double& v,vector& N,vector& B); +const void swap_columns(Mat_& A,int col1,int col2); //return codes:-2 (no_sol - unbdd),-1(no_sol - unfsbl), 0(single_sol), 1(multiple_sol=>least_l2_norm) int solveLP(const Mat& Func, const Mat& Constr, Mat& z){ - printf("call to solveLP\n"); + dprintf("call to solveLP\n"); - //sanity check (size, type, no. of channels) (and throw exception, if appropriate) + //sanity check (size, type, no. of channels) CV_Assert(Func.type()==CV_64FC1); CV_Assert(Constr.type()==CV_64FC1); CV_Assert(Func.rows==1); @@ -82,15 +96,15 @@ int solveLP(const Mat& Func, const Mat& Constr, Mat& z){ double v=0; vector N,B; - if(solveLP_aux::initialize_simplex(bigC,bigB,v,N,B)==-1){ - return -1; + if(initialize_simplex(bigC,bigB,v,N,B)==SOLVELP_UNFEASIBLE){ + return SOLVELP_UNFEASIBLE; } Mat_ c=bigC.colRange(1,bigC.cols), b=bigB.colRange(1,bigB.cols); int res=0; - if((res=solveLP_aux::inner_simplex(c,b,v,N,B))==-2){ - return -2; + if((res=inner_simplex(c,b,v,N,B))==SOLVELP_UNBOUNDED){ + return SOLVELP_UNBOUNDED; } //return the optimal solution @@ -109,7 +123,7 @@ int solveLP(const Mat& Func, const Mat& Constr, Mat& z){ return res; } -int solveLP_aux::initialize_simplex(Mat_& c, Mat_& b,double& v,vector& N,vector& B){ +const int initialize_simplex(Mat_& c, Mat_& b,double& v,vector& N,vector& B){ N.resize(c.cols); N[0]=0; for (std::vector::iterator it = N.begin()+1 ; it != N.end(); ++it){ @@ -147,21 +161,21 @@ int solveLP_aux::initialize_simplex(Mat_& c, Mat_& b,double& v,v print_simplex_state(c,b,v,N,B); - printf("\tWE MAKE PIVOT\n"); + dprintf("\tWE MAKE PIVOT\n"); pivot(c,b,v,N,B,k,0); print_simplex_state(c,b,v,N,B); inner_simplex(c,b,v,N,B); - printf("\tAFTER INNER_SIMPLEX\n"); + dprintf("\tAFTER INNER_SIMPLEX\n"); print_simplex_state(c,b,v,N,B); vector::iterator it=std::find(B.begin(),B.end(),0); if(it!=B.end()){ int it_offset=it-B.begin(); if(b(it_offset,b.cols-1)>0){ - return -1; + return SOLVELP_UNFEASIBLE; } pivot(c,b,v,N,B,it_offset,0); } @@ -172,7 +186,7 @@ int solveLP_aux::initialize_simplex(Mat_& c, Mat_& b,double& v,v swap_columns(c,it_offset,0); swap_columns(b,it_offset,0); - printf("after swaps\n"); + dprintf("after swaps\n"); print_simplex_state(c,b,v,N,B); //start from 1, because we ignore x_0 @@ -180,14 +194,14 @@ int solveLP_aux::initialize_simplex(Mat_& c, Mat_& b,double& v,v v=0; for(int i=1;i& c, Mat_& b,double& v,v } } - printf("after restore\n"); + dprintf("after restore\n"); print_simplex_state(c,b,v,N,B); N.erase(N.begin()); return 0; } -int solveLP_aux::inner_simplex(Mat_& c, Mat_& b,double& v,vector& N,vector& B){ +const int inner_simplex(Mat_& c, Mat_& b,double& v,vector& N,vector& B){ int count=0; while(1){ - printf("iteration #%d\n",count++); + dprintf("iteration #%d\n",count++); - MatIterator_ pos_ptr; + static MatIterator_ pos_ptr; int e=-1,pos_ctr=0,min_var=INT_MAX; bool all_nonzero=true; for(pos_ptr=c.begin();pos_ptr!=c.end();pos_ptr++,pos_ctr++){ @@ -223,21 +237,15 @@ int solveLP_aux::inner_simplex(Mat_& c, Mat_& b,double& v,vector } } if(e==-1){ - printf("hello from e==-1\n"); + dprintf("hello from e==-1\n"); print_matrix(c); if(all_nonzero==true){ - return 0; + return SOLVELP_SINGLE; }else{ - return 1; + return SOLVELP_MULTI; } } - - /*for(pos_ptr=c.begin();(*pos_ptr<=0) && pos_ptr!=c.end();pos_ptr++,e++);//TODO: select the smallest index var w/ pos coef - if(pos_ptr==c.end()){ - return 0; - }*/ - int l=-1; min_var=INT_MAX; double min=DBL_MAX; @@ -259,30 +267,29 @@ int solveLP_aux::inner_simplex(Mat_& c, Mat_& b,double& v,vector } } if(l==-1){ - //unbounded - return -2; + return SOLVELP_UNBOUNDED; } - printf("the tightest constraint is in row %d with %g\n",l,min); + dprintf("the tightest constraint is in row %d with %g\n",l,min); - solveLP_aux::pivot(c,b,v,N,B,l,e); + pivot(c,b,v,N,B,l,e); - printf("objective, v=%g\n",v); + dprintf("objective, v=%g\n",v); print_matrix(c); - printf("constraints\n"); + dprintf("constraints\n"); print_matrix(b); - printf("non-basic: "); + dprintf("non-basic: "); for (std::vector::iterator it = N.begin() ; it != N.end(); ++it){ - printf("%d, ",*it); + dprintf("%d, ",*it); } - printf("\nbasic: "); + dprintf("\nbasic: "); for (std::vector::iterator it = B.begin() ; it != B.end(); ++it){ - printf("%d, ",*it); + dprintf("%d, ",*it); } - printf("\n"); + dprintf("\n"); } } -inline void solveLP_aux::pivot(Mat_& c,Mat_& b,double& v,vector& N,vector& B, int leaving_index,int entering_index){ +const inline void pivot(Mat_& c,Mat_& b,double& v,vector& N,vector& B, int leaving_index,int entering_index){ double coef=b(leaving_index,entering_index); for(int i=0;i& c,Mat_& b,double& v,vector< c(0,i)-=coef*b(leaving_index,i); } } - printf("v was %g\n",v); + dprintf("v was %g\n",v); v+=coef*b(leaving_index,b.cols-1); int tmp=N[entering_index]; @@ -322,7 +329,7 @@ inline void solveLP_aux::pivot(Mat_& c,Mat_& b,double& v,vector< B[leaving_index]=tmp; } -void solveLP_aux::swap_columns(Mat_& A,int col1,int col2){ +const inline void swap_columns(Mat_& A,int col1,int col2){ for(int i=0;i