generic-library/opencv
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Cleaning the code of simplex method

Browse Source 
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.
This commit is contained in:
Alex Leontiev 2013-07-10 20:11:52 +03:00
parent a4a5e98cc0
commit a95650111f
5 changed files with 193 additions and 118 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

48
modules/optim/doc/linear_programming.rst Normal file
View File

@ -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
};

10
modules/optim/doc/optim.rst Normal file
View File

@ -0,0 +1,10 @@
**************************************
optim. Generic numerical optimization
**************************************
.. highlight:: cpp
.. toctree::
:maxdepth: 2
linear_programming

94
modules/optim/include/opencv2/optim.hpp
View File

@ -47,67 +47,77 @@
#include "opencv2/core.hpp" #include "opencv2/core.hpp"
#include "opencv2/core/mat.hpp" #include "opencv2/core/mat.hpp"
/*! \namespace cv //uncomment the next line to print the debug info
Namespace where all the C++ OpenCV functionality resides //#define ALEX_DEBUG
*/
namespace cv{namespace optim namespace cv{namespace optim
{ {
//! generic class for optimization algorithms */ //! generic class for optimization algorithms */
class CV_EXPORTS Solver : public Algorithm /* Algorithm is the base OpenCV class */ class CV_EXPORTS Solver : public Algorithm /* Algorithm is the base OpenCV class */
{ {
public: public:
class CV_EXPORTS Function class CV_EXPORTS Function
{ {
public: public:
virtual ~Function(){} virtual ~Function(){}
virtual double calc(InputArray args) const = 0; virtual double calc(InputArray args) const = 0;
}; };
class CV_EXPORTS Constraints class CV_EXPORTS Constraints
{ {
public: public:
virtual ~Constraints(){} virtual ~Constraints(){}
}; };
//! could be reused for all the generic algorithms like downhill simplex. Return value is the maximum value of a function*/ //! 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 double solve(const Function& F,const Constraints& C, OutputArray result) const = 0;
/*virtual void setTermCriteria(const TermCriteria& criteria) = 0; /*virtual void setTermCriteria(const TermCriteria& criteria) = 0;
virtual TermCriteria getTermCriteria() = 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 class CV_EXPORTS LPSolver : public Solver
{ {
public: public:
class CV_EXPORTS LPFunction:public Solver::Function class CV_EXPORTS LPFunction:public Solver::Function
{ {
Mat z; Mat z;
public: public:
//! Note, that this class is supposed to be immutable, so it's ok to make only a shallow copy of z_in.*/ //! 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(Mat z_in):z(z_in){}
~LPFunction(){}; ~LPFunction(){};
const Mat& getz()const{return z;} const Mat& getz()const{return z;}
double calc(InputArray args)const; 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. //!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 //!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**. //!this form and **we shall create various constructors for this class that will perform these conversions**.
class CV_EXPORTS LPConstraints:public Solver::Constraints class CV_EXPORTS LPConstraints:public Solver::Constraints
{ {
Mat A,b; Mat A,b;
public: public:
~LPConstraints(){}; ~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.*/ //! 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){} LPConstraints(Mat A_in, Mat b_in):A(A_in),b(b_in){}
const Mat& getA()const{return A;} const Mat& getA()const{return A;}
const Mat& getb()const{return b;} const Mat& getb()const{return b;}
}; };
LPSolver(){} LPSolver(){}
double solve(const Function& F,const Constraints& C, OutputArray result)const; 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_EXPORTS_W int solveLP(const Mat& Func, const Mat& Constr, Mat& z);
}}// cv }}// cv

151
modules/optim/src/lpsolver.cpp
View File

@ -2,73 +2,87 @@
#include "precomp.hpp" #include "precomp.hpp"
#include <climits> #include <climits>
#include <algorithm> #include <algorithm>
#include <cstdarg>
namespace cv{namespace optim{ namespace cv{namespace optim{
using std::vector; 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; return 0.0;
} }
double LPSolver::LPFunction::calc(InputArray args)const{ double LPSolver::LPFunction::calc(InputArray args)const{
printf("call to LPFunction::calc()\n"); dprintf("call to LPFunction::calc()\n");
return 0.0; 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++){ for(int i=0;i<X.rows;i++){
printf("\t["); dprintf("\t[");
for(int j=0;j<X.cols;j++){ 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); print_matrix(c);
printf("non-basic: "); dprintf("non-basic: ");
for (std::vector<int>::const_iterator it = N.begin() ; it != N.end(); ++it){ 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); print_matrix(b);
printf("basic: "); dprintf("basic: ");
for (std::vector<int>::const_iterator it = B.begin() ; it != B.end(); ++it){ 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:
/**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
*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
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
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.
* @return -1 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);
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);
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.
/**@return -2 means the problem is unbdd, 1 means multiple solutions, 0 means successful. */
*/ const int inner_simplex(Mat_<double>& c, Mat_<double>& b,double& v,vector<int>& N,vector<int>& B);
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);
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) //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){ 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(Func.type()==CV_64FC1);
CV_Assert(Constr.type()==CV_64FC1); CV_Assert(Constr.type()==CV_64FC1);
CV_Assert(Func.rows==1); CV_Assert(Func.rows==1);
@ -82,15 +96,15 @@ int solveLP(const Mat& Func, const Mat& Constr, Mat& z){
double v=0; double v=0;
vector<int> N,B; vector<int> N,B;
if(solveLP_aux::initialize_simplex(bigC,bigB,v,N,B)==-1){ if(initialize_simplex(bigC,bigB,v,N,B)==SOLVELP_UNFEASIBLE){
return -1; return SOLVELP_UNFEASIBLE;
} }
Mat_<double> c=bigC.colRange(1,bigC.cols), Mat_<double> c=bigC.colRange(1,bigC.cols),
b=bigB.colRange(1,bigB.cols); b=bigB.colRange(1,bigB.cols);
int res=0; int res=0;
if((res=solveLP_aux::inner_simplex(c,b,v,N,B))==-2){ if((res=inner_simplex(c,b,v,N,B))==SOLVELP_UNBOUNDED){
return -2; return SOLVELP_UNBOUNDED;
} }
//return the optimal solution //return the optimal solution
@ -109,7 +123,7 @@ int solveLP(const Mat& Func, const Mat& Constr, Mat& z){
return res; 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.resize(c.cols);
N[0]=0; N[0]=0;
for (std::vector<int>::iterator it = N.begin()+1 ; it != N.end(); ++it){ 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); 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); pivot(c,b,v,N,B,k,0);
print_simplex_state(c,b,v,N,B); print_simplex_state(c,b,v,N,B);
inner_simplex(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); print_simplex_state(c,b,v,N,B);
vector<int>::iterator it=std::find(B.begin(),B.end(),0); vector<int>::iterator it=std::find(B.begin(),B.end(),0);
if(it!=B.end()){ if(it!=B.end()){
int it_offset=it-B.begin(); int it_offset=it-B.begin();
if(b(it_offset,b.cols-1)>0){ if(b(it_offset,b.cols-1)>0){
return -1; return SOLVELP_UNFEASIBLE;
} }
pivot(c,b,v,N,B,it_offset,0); 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(c,it_offset,0);
swap_columns(b,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); print_simplex_state(c,b,v,N,B);
//start from 1, because we ignore x_0 //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; v=0;
for(int i=1;i<old_c.cols;i++){ for(int i=1;i<old_c.cols;i++){
if((it=std::find(N.begin(),N.end(),i))!=N.end()){ 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); fflush(stdout);
int it_offset=it-N.begin(); int it_offset=it-N.begin();
c(0,it_offset)+=old_c(0,i); c(0,it_offset)+=old_c(0,i);
print_matrix(c); print_matrix(c);
}else{ }else{
//cv::Mat_ //cv::Mat_
printf("i=%d from basic\n",i); dprintf("i=%d from basic\n",i);
fflush(stdout); fflush(stdout);
int it_offset=std::find(B.begin(),B.end(),i)-B.begin(); 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); 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); print_simplex_state(c,b,v,N,B);
N.erase(N.begin()); N.erase(N.begin());
return 0; 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; int count=0;
while(1){ 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; int e=-1,pos_ctr=0,min_var=INT_MAX;
bool all_nonzero=true; bool all_nonzero=true;
for(pos_ptr=c.begin();pos_ptr!=c.end();pos_ptr++,pos_ctr++){ 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){ if(e==-1){
printf("hello from e==-1\n"); dprintf("hello from e==-1\n");
print_matrix(c); print_matrix(c);
if(all_nonzero==true){ if(all_nonzero==true){
return 0; return SOLVELP_SINGLE;
}else{ }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; int l=-1;
min_var=INT_MAX; min_var=INT_MAX;
double min=DBL_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){ if(l==-1){
//unbounded return SOLVELP_UNBOUNDED;
return -2;
} }
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); print_matrix(c);
printf("constraints\n"); dprintf("constraints\n");
print_matrix(b); print_matrix(b);
printf("non-basic: "); dprintf("non-basic: ");
for (std::vector<int>::iterator it = N.begin() ; it != N.end(); ++it){ 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){ 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); double coef=b(leaving_index,entering_index);
for(int i=0;i<b.cols;i++){ for(int i=0;i<b.cols;i++){
if(i==entering_index){ 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); 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); v+=coef*b(leaving_index,b.cols-1);
int tmp=N[entering_index]; 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; 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++){ for(int i=0;i<A.rows;i++){
double tmp=A(i,col1); double tmp=A(i,col1);
A(i,col1)=A(i,col2); A(i,col1)=A(i,col2);

8
modules/optim/test/test_lpsolver.cpp
View File

@ -1,8 +1,7 @@
#include "test_precomp.hpp" #include "test_precomp.hpp"
#include "opencv2/optim.hpp" #include "opencv2/optim.hpp"
TEST(Optim_LpSolver, regression_basic) TEST(Optim_LpSolver, regression_basic){
{
cv::Mat A,B,z,etalon_z; cv::Mat A,B,z,etalon_z;
if(true){ if(true){
@ -120,9 +119,10 @@ TEST(Optim_LpSolver, regression_cycling){
// //
// ??how_check_multiple_solutions & pass_test - DONE // ??how_check_multiple_solutions & pass_test - DONE
// Blands_rule - DONE // Blands_rule - DONE
// (&1tests on cycling) // (assert, assign) - DONE
// //
// make_more_clear (assert, assign) // (&1tests on cycling)
// make_more_clear
// wrap in OOP // wrap in OOP
// //
// non-trivial tests // non-trivial tests