From ddc0010e7da6daa4d8b6a543e853ef0b0ff4ff33 Mon Sep 17 00:00:00 2001 From: Alex Leontiev Date: Fri, 28 Jun 2013 15:28:57 +0300 Subject: [PATCH] The first draft of simplex algorithm, simple tests. What we have now corresponds to "formal simplex algorithm", described in Cormen's "Intro to Algorithms". It will work *only* if the initial problem has (0,0,0,...,0) as feasible solution (consequently, it will work unpredictably if problem was unfeasible or did not have zero-vector as feasible solution). Moreover, it might cycle. TODO (first priority) 1. Implement initialize_simplex() procedure, that shall check for feasibility and generate initial feasible solution. (in particular, code should pass all 4 tests implemented at the moment) 2. Implement Bland's rule to avoid cycling. 3. Make the code more clear. 4. Implement several non-trivial tests (??) and check algorithm against them. Debug if necessary. TODO (second priority) 1. Concentrate on stability and speed (make difficult tests) --- modules/optim/include/opencv2/optim.hpp | 15 +- modules/optim/src/lpsolver.cpp | 257 +++++++++++++++++++++++++++++++- modules/optim/test/test_lpsolver.cpp | 61 ++++++++ modules/optim/test/test_main.cpp | 3 + modules/optim/test/test_precomp.cpp | 1 + modules/optim/test/test_precomp.hpp | 15 ++ 6 files changed, 338 insertions(+), 14 deletions(-) create mode 100644 modules/optim/test/test_lpsolver.cpp create mode 100644 modules/optim/test/test_main.cpp create mode 100644 modules/optim/test/test_precomp.cpp create mode 100644 modules/optim/test/test_precomp.hpp diff --git a/modules/optim/include/opencv2/optim.hpp b/modules/optim/include/opencv2/optim.hpp index 8464cf8..f4a639f 100644 --- a/modules/optim/include/opencv2/optim.hpp +++ b/modules/optim/include/opencv2/optim.hpp @@ -82,12 +82,12 @@ class CV_EXPORTS LPSolver : public Solver public: class CV_EXPORTS LPFunction:public Solver::Function { - cv::Mat z; + 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(cv::Mat z_in):z(z_in){} + LPFunction(Mat z_in):z(z_in){} ~LPFunction(){}; - const cv::Mat& getz()const{return z;} + const Mat& getz()const{return z;} double calc(InputArray args)const; }; @@ -96,18 +96,19 @@ public: //!this form and **we shall create various constructors for this class that will perform these conversions**. class CV_EXPORTS LPConstraints:public Solver::Constraints { - cv::Mat A,b; + 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(cv::Mat A_in, cv::Mat b_in):A(A_in),b(b_in){} - const cv::Mat& getA()const{return A;} - const cv::Mat& getb()const{return b;} + 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; }; +CV_EXPORTS_W int solveLP(const Mat& Func, const Mat& Constr, Mat& z); }}// cv #endif diff --git a/modules/optim/src/lpsolver.cpp b/modules/optim/src/lpsolver.cpp index 56dec12..b9cabc1 100644 --- a/modules/optim/src/lpsolver.cpp +++ b/modules/optim/src/lpsolver.cpp @@ -1,15 +1,12 @@ +#include "opencv2/ts.hpp" #include "precomp.hpp" +#include +#include namespace cv{namespace optim{ +using std::vector; double LPSolver::solve(const Function& F,const Constraints& C, OutputArray result)const{ - printf("call to solve\n"); - - //TODO: sanity check and throw exception, if appropriate - - //TODO: copy A,b,z - - //TODO: run simplex algo return 0.0; } @@ -18,5 +15,251 @@ double LPSolver::LPFunction::calc(InputArray args)const{ printf("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); + for(int i=0;i(i,j)); + } + printf("]\n"); + } +} +namespace solveLP_aux{ + //return -1 if problem is unfeasible, 0 if feasible + //in latter case it returns feasible solution in z with homogenised b's and v + int initialize_simplex(const Mat& c, Mat& b, Mat& z,double& v); +} +int solveLP(const Mat& Func, const Mat& Constr, Mat& z){ + printf("call to solveLP\n");//-3(incorrect),-2 (no_sol - unbdd),-1(no_sol - unfsbl), 0(single_sol), 1(multiple_sol=>least_l2_norm) + + //sanity check (size, type, no. of channels) (and throw exception, if appropriate) + if(Func.type()!=CV_64FC1 || Constr.type()!=CV_64FC1){ + printf("both Func and Constr should be one-channel matrices of double's\n"); + return -3; + } + if(Func.rows!=1){ + printf("Func should be row-vector\n"); + return -3; + } + vector N(Func.cols); + N[0]=1; + for (std::vector::iterator it = N.begin()+1 ; it != N.end(); ++it){ + *it=it[-1]+1; + } + if((Constr.cols-1)!=Func.cols){ + printf("Constr should have one more column when compared to Func\n"); + return -3; + } + vector B(Constr.rows); + B[0]=Func.cols+1; + for (std::vector::iterator it = B.begin()+1 ; it != B.end(); ++it){ + *it=it[-1]+1; + } + + //copy arguments for we will shall modify them + Mat c=Func.clone(), + b=Constr.clone(); + double v=0; + + solveLP_aux::initialize_simplex(c,b,z,v); + + int count=0; + while(1){ + printf("iteration #%d\n",count++); + + MatIterator_ pos_ptr; + int e=0; + for(pos_ptr=c.begin();(*pos_ptr<=0) && pos_ptr!=c.end();pos_ptr++,e++); + if(pos_ptr==c.end()){ + break; + } + printf("offset of first nonneg coef is %d\n",e);//TODO: choose the var with the smallest index + + int l=-1; + double min=DBL_MAX; + int row_it=0; + double ite=0; + MatIterator_ min_row_ptr=b.begin(); + for(MatIterator_ it=b.begin();it!=b.end();it+=b.cols,row_it++){ + double myite=0; + //check constraints, select the tightest one, TODO: smallest index + if((myite=it[e])>0){ + double val=it[b.cols-1]/myite; + if(val it=min_row_ptr;col_count it=b.begin();row_count()); + if(row_count==l){ + it+=b.cols; + }else{ + //remaining constraints + double coef=it[e]; + MatIterator_ shadow_it=min_row_ptr; + for(int col_it=0;col_it<(b.cols);col_it++,it++,shadow_it++){ + if(col_it==e){ + *it=-coef*(*shadow_it); + }else{ + *it-=coef*(*shadow_it); + } + } + } + } + //objective function + double coef=*pos_ptr; + MatIterator_ shadow_it=min_row_ptr; + MatIterator_ it=c.begin(); + for(int col_it=0;col_it<(b.cols-1);col_it++,it++,shadow_it++){ + if(col_it==e){ + *it=-coef*(*shadow_it); + }else{ + *it-=coef*(*shadow_it); + } + } + v+=coef*(*shadow_it); + + //new basis and nonbasic sets + int tmp=N[e]; + N[e]=B[l]; + B[l]=tmp; + + printf("objective, v=%g\n",v); + print_matrix(c); + printf("constraints\n"); + print_matrix(b); + printf("non-basic: "); + for (std::vector::iterator it = N.begin() ; it != N.end(); ++it){ + printf("%d, ",*it); + } + printf("\nbasic: "); + for (std::vector::iterator it = B.begin() ; it != B.end(); ++it){ + printf("%d, ",*it); + } + printf("\n"); + } + + //return the optimal solution + //z=cv::Mat_(1,c.cols,0); + MatIterator_ it=z.begin(); + for(int i=1;i<=c.cols;i++,it++){ + std::vector::iterator pos=B.begin(); + if((pos=std::find(B.begin(),B.end(),i))==B.end()){ + *it+=0; + }else{ + *it+=b.at(pos-B.begin(),b.cols-1); + } + } + + return 0; +} +int solveLP_aux::initialize_simplex(const Mat& c, Mat& b, Mat& z,double& v){//TODO + z=Mat_(1,c.cols,0.0); + v=0; + return 0; + + cv::Mat mod_b=(cv::Mat_(1,b.rows)); + bool gen_new_sol_flag=false,hom_sol_given=false; + if(z.type()!=CV_64FC1 || z.rows!=1 || z.cols!=c.cols || (hom_sol_given=(countNonZero(z)==0))){ + printf("line %d\n",__LINE__); + if(hom_sol_given==false){ + printf("line %d, %d\n",__LINE__,hom_sol_given); + z=cv::Mat_(1,c.cols,(double)0); + } + //check homogeneous solution + printf("line %d\n",__LINE__); + for(MatIterator_ b_it=b.begin()+b.cols-1,mod_b_it=mod_b.begin();mod_b_it!=mod_b.end(); + b_it+=b.cols,mod_b_it++){ + if(*b_it<0){ + //if no - we need feasible solution + gen_new_sol_flag=true; + break; + } + } + printf("line %d, gen_new_sol_flag=%d - I've got here!!!\n",__LINE__,gen_new_sol_flag); + //if yes - we have feasible solution! + }else{ + //check for feasibility + MatIterator_ it=b.begin(); + for(MatIterator_ mod_b_it=mod_b.begin();it!=b.end();mod_b_it++){ + double sum=0; + for(MatIterator_ z_it=z.begin();z_it!=z.end();z_it++,it++){ + sum+=(*it)*(*z_it); + } + if((*mod_b_it=(*it-sum))<0){ + break; + } + it++; + } + if(it==b.end()){ + //z contains feasible solution - just homogenise b's TODO: and v + gen_new_sol_flag=false; + for(MatIterator_ b_it=b.begin()+b.cols-1,mod_b_it=mod_b.begin();mod_b_it!=mod_b.end(); + b_it+=b.cols,mod_b_it++){ + *b_it=*mod_b_it; + } + }else{ + //if no - we need feasible solution + gen_new_sol_flag=true; + } + } + if(gen_new_sol_flag==true){ + //we should generate new solution - TODO + printf("we should generate new solution\n"); + Mat new_c=Mat_(1,c.cols+1,0.0), + new_b=Mat_(b.rows,b.cols+1,-1.0), + new_z=Mat_(1,c.cols+1,0.0); + + new_c.end()[-1]=-1; + c.copyTo(new_c.colRange(0,new_c.cols-1)); + + b.col(b.cols-1).copyTo(new_b.col(new_b.cols-1)); + b.colRange(0,b.cols-1).copyTo(new_b.colRange(0,new_b.cols-2)); + + Mat b_slice=b.col(b.cols-1); + new_z.end()[-1]=-*(std::min_element(b_slice.begin(),b_slice.end())); + + /*printf("matrix new_c\n"); + print_matrix(new_c); + printf("matrix new_b\n"); + print_matrix(new_b); + printf("matrix new_z\n"); + print_matrix(new_z);*/ + + printf("run for the second time!\n"); + solveLP(new_c,new_b,new_z); + printf("original z was\n"); + print_matrix(z); + printf("that's what I've got\n"); + print_matrix(new_z); + printf("for the constraints\n"); + print_matrix(b); + return 0; + } + +} }} diff --git a/modules/optim/test/test_lpsolver.cpp b/modules/optim/test/test_lpsolver.cpp new file mode 100644 index 0000000..c95edf5 --- /dev/null +++ b/modules/optim/test/test_lpsolver.cpp @@ -0,0 +1,61 @@ +#include "test_precomp.hpp" +#include "opencv2/optim.hpp" + +TEST(Optim_LpSolver, regression) +{ + cv::Mat A,B,z,etalon_z; + + if(true){ + //cormen's example #1 + A=(cv::Mat_(1,3)<<3,1,2); + B=(cv::Mat_(3,4)<<1,1,3,30,2,2,5,24,4,1,2,36); + std::cout<<"here A goes\n"<(1,3)<<8,4,0); + ASSERT_EQ(cv::countNonZero(z!=etalon_z),0); + } + + if(true){ + //cormen's example #2 + A=(cv::Mat_(1,2)<<18,12.5); + B=(cv::Mat_(3,3)<<1,1,20,1,0,20,0,1,16); + std::cout<<"here A goes\n"<(1,2)<<20,0); + ASSERT_EQ(cv::countNonZero(z!=etalon_z),0); + } + + if(true){ + //cormen's example #3 + A=(cv::Mat_(1,2)<<5,-3); + B=(cv::Mat_(2,3)<<1,-1,1,2,1,2); + std::cout<<"here A goes\n"<(1,2)<<1,0); + ASSERT_EQ(cv::countNonZero(z!=etalon_z),0); + + } + if(false){ + //cormen's example #4 - unfeasible + A=(cv::Mat_(1,3)<<-1,-1,-1); + B=(cv::Mat_(2,4)<<-2,-7.5,-3,-10000,-20,-5,-10,-30000); + std::cout<<"here A goes\n"<(1,2)<<1,0); + ASSERT_EQ(cv::countNonZero(z!=etalon_z),0); + } +} + +//TODO +// get optimal solution from initial (0,0,...,0) - DONE +// milestone: pass first test (wo initial solution) - DONE + // learn how to get initial solution + // Blands_rule + // 1_more_test & make_more_clear + // -> **contact_Vadim**: min_l2_norm, init_optional_fsbl_check, error_codes, comment_style-too_many?, copyTo temp headers +// ??how to get smallest l2 norm +// FUTURE: compress&debug-> more_tests(Cormen) -> readNumRecipes-> fast&stable || hill_climbing diff --git a/modules/optim/test/test_main.cpp b/modules/optim/test/test_main.cpp new file mode 100644 index 0000000..6b24993 --- /dev/null +++ b/modules/optim/test/test_main.cpp @@ -0,0 +1,3 @@ +#include "test_precomp.hpp" + +CV_TEST_MAIN("cv") diff --git a/modules/optim/test/test_precomp.cpp b/modules/optim/test/test_precomp.cpp new file mode 100644 index 0000000..5956e13 --- /dev/null +++ b/modules/optim/test/test_precomp.cpp @@ -0,0 +1 @@ +#include "test_precomp.hpp" diff --git a/modules/optim/test/test_precomp.hpp b/modules/optim/test/test_precomp.hpp new file mode 100644 index 0000000..9a86cab --- /dev/null +++ b/modules/optim/test/test_precomp.hpp @@ -0,0 +1,15 @@ +#ifdef __GNUC__ +# pragma GCC diagnostic ignored "-Wmissing-declarations" +# if defined __clang__ || defined __APPLE__ +# pragma GCC diagnostic ignored "-Wmissing-prototypes" +# pragma GCC diagnostic ignored "-Wextra" +# endif +#endif + +#ifndef __OPENCV_TEST_PRECOMP_HPP__ +#define __OPENCV_TEST_PRECOMP_HPP__ + +#include "opencv2/ts.hpp" +#include "opencv2/optim.hpp" + +#endif -- 2.7.4