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.

diff --git a/modules/optim/doc/linear_programming.rst b/modules/optim/doc/linear_programming.rst
new file mode 100644 (file)
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 <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 (file)
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 (file)
--- 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 (file)
--- a/modules/optim/src/lpsolver.cpp
+++ b/modules/optim/src/lpsolver.cpp
@@ -2,73 +2,87 @@
 #include "precomp.hpp"
 #include <climits>
 #include <algorithm>
+#include <cstdarg>
 
 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<X.rows;i++){
-      printf("\t[");
+      dprintf("\t[");
       for(int j=0;j<X.cols;j++){
-          printf("%g, ",X.at<double>(i,j));
+          dprintf("%g, ",X.at<double>(i,j));
       }
-      printf("]\n");
+      dprintf("]\n");
     } 
+#endif
 }
-void print_simplex_state(const Mat& c,const Mat&b,double v,const vector<int>& N,const vector<int>& 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<int>& N,const vector<int>& 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<int>::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<int>::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_<double>& c, Mat_<double>& b,double& v,vector<int>& N,vector<int>& B);
-    inline void pivot(Mat_<double>& c,Mat_<double>& b,double& v,vector<int>& N,vector<int>& 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_<double>& c, Mat_<double>& b,double& v,vector<int>& N,vector<int>& B);
-    void swap_columns(Mat_<double>& 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_<double>& c, Mat_<double>& b,double& v,vector<int>& N,vector<int>& B);
+const inline void pivot(Mat_<double>& c,Mat_<double>& b,double& v,vector<int>& N,vector<int>& 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_<double>& c, Mat_<double>& b,double& v,vector<int>& N,vector<int>& B);
+const void swap_columns(Mat_<double>& 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<int> 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_<double> 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_<double>& c, Mat_<double>& b,double& v,vector<int>& N,vector<int>& B){
+const int initialize_simplex(Mat_<double>& c, Mat_<double>& b,double& v,vector<int>& N,vector<int>& B){
     N.resize(c.cols);
     N[0]=0;
     for (std::vector<int>::iterator it = N.begin()+1 ; it != N.end(); ++it){
@@ -147,21 +161,21 @@ int solveLP_aux::initialize_simplex(Mat_<double>& c, Mat_<double>& 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<int>::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_<double>& c, Mat_<double>& 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_<double>& c, Mat_<double>& b,double& v,v
     v=0;
     for(int i=1;i<old_c.cols;i++){
         if((it=std::find(N.begin(),N.end(),i))!=N.end()){
-            printf("i=%d from nonbasic\n",i);
+            dprintf("i=%d from nonbasic\n",i);
             fflush(stdout);
             int it_offset=it-N.begin();
             c(0,it_offset)+=old_c(0,i);     
             print_matrix(c);
         }else{
             //cv::Mat_
-            printf("i=%d from basic\n",i);
+            dprintf("i=%d from basic\n",i);
             fflush(stdout);
             int it_offset=std::find(B.begin(),B.end(),i)-B.begin();
             c-=old_c(0,i)*b.row(it_offset).colRange(0,b.cols-1);
@@ -196,19 +210,19 @@ int solveLP_aux::initialize_simplex(Mat_<double>& c, Mat_<double>& 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_<double>& c, Mat_<double>& b,double& v,vector<int>& N,vector<int>& B){
+const int inner_simplex(Mat_<double>& c, Mat_<double>& b,double& v,vector<int>& N,vector<int>& B){
     int count=0;
     while(1){
-        printf("iteration #%d\n",count++);
+        dprintf("iteration #%d\n",count++);
 
-        MatIterator_<double> pos_ptr;
+        static MatIterator_<double> 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_<double>& c, Mat_<double>& 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_<double>& c, Mat_<double>& 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<int>::iterator it = N.begin() ; it != N.end(); ++it){
-            printf("%d, ",*it);
+            dprintf("%d, ",*it);
         }
-        printf("\nbasic: ");
+        dprintf("\nbasic: ");
         for (std::vector<int>::iterator it = B.begin() ; it != B.end(); ++it){
-            printf("%d, ",*it);
+            dprintf("%d, ",*it);
         }
-        printf("\n");
+        dprintf("\n");
     }
 }
 
-inline void solveLP_aux::pivot(Mat_<double>& c,Mat_<double>& b,double& v,vector<int>& N,vector<int>& B, int leaving_index,int entering_index){
+const inline void pivot(Mat_<double>& c,Mat_<double>& b,double& v,vector<int>& N,vector<int>& B, int leaving_index,int entering_index){
     double coef=b(leaving_index,entering_index);
     for(int i=0;i<b.cols;i++){
         if(i==entering_index){
@@ -314,7 +321,7 @@ inline void solveLP_aux::pivot(Mat_<double>& c,Mat_<double>& 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_<double>& c,Mat_<double>& b,double& v,vector<
     B[leaving_index]=tmp;
 }
 
-void solveLP_aux::swap_columns(Mat_<double>& A,int col1,int col2){
+const inline void swap_columns(Mat_<double>& A,int col1,int col2){
     for(int i=0;i<A.rows;i++){
         double tmp=A(i,col1);
         A(i,col1)=A(i,col2);

diff --git a/modules/optim/test/test_lpsolver.cpp b/modules/optim/test/test_lpsolver.cpp
index 9219f0a..8d14bf9 100644 (file)
--- a/modules/optim/test/test_lpsolver.cpp
+++ b/modules/optim/test/test_lpsolver.cpp
@@ -1,8 +1,7 @@
 #include "test_precomp.hpp"
 #include "opencv2/optim.hpp"
 
-TEST(Optim_LpSolver, regression_basic)
-{
+TEST(Optim_LpSolver, regression_basic){
     cv::Mat A,B,z,etalon_z;
 
     if(true){
@@ -120,9 +119,10 @@ TEST(Optim_LpSolver, regression_cycling){
     //
     // ??how_check_multiple_solutions & pass_test - DONE
     // Blands_rule - DONE
-    // (&1tests on cycling)
+    // (assert, assign) - DONE
     //
-    // make_more_clear (assert, assign)
+    // (&1tests on cycling)
+    // make_more_clear
     // wrap in OOP
     //
     // non-trivial tests