diff --git a/.gitignore b/.gitignore index 0a19f3ceeb..caaebed0e3 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,4 @@ OpenCV4Tegra/ .sw[a-z] .*.swp tags +build/ diff --git a/modules/optim/CMakeLists.txt b/modules/optim/CMakeLists.txt new file mode 100644 index 0000000000..c36c24d9de --- /dev/null +++ b/modules/optim/CMakeLists.txt @@ -0,0 +1,2 @@ +set(the_description "Generic optimization") +ocv_define_module(optim opencv_core) diff --git a/modules/optim/doc/linear_programming.rst b/modules/optim/doc/linear_programming.rst new file mode 100644 index 0000000000..946df9e955 --- /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). It should contain 32- or 64-bit floating point numbers. As a convenience, column-vector may be also submitted, in the latter case it is understood to correspond to :math:`c^T`. + + :param Constr: *m*-by-*n\+1* matrix, whose rightmost column corresponds to :math:`b` in formulation above and the remaining to :math:`A`. It should containt 32- or 64-bit floating point numbers. + + :param z: The solution will be returned here as a column-vector - it corresponds to :math:`c` in the formulation above. It will contain 64-bit floating point numbers. + + :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 0000000000..bead2122a4 --- /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 new file mode 100644 index 0000000000..0bbedad38f --- /dev/null +++ b/modules/optim/include/opencv2/optim.hpp @@ -0,0 +1,62 @@ +/*M/////////////////////////////////////////////////////////////////////////////////////// +// +// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. +// +// By downloading, copying, installing or using the software you agree to this license. +// If you do not agree to this license, do not download, install, +// copy or use the software. +// +// +// License Agreement +// For Open Source Computer Vision Library +// +// Copyright (C) 2000-2008, Intel Corporation, all rights reserved. +// Copyright (C) 2008-2012, Willow Garage Inc., all rights reserved. +// Third party copyrights are property of their respective owners. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// * Redistribution's of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// * Redistribution's in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// * The name of the copyright holders may not be used to endorse or promote products +// derived from this software without specific prior written permission. +// +// This software is provided by the copyright holders and contributors "as is" and +// any express or implied warranties, including, but not limited to, the implied +// warranties of merchantability and fitness for a particular purpose are disclaimed. +// In no event shall the Intel Corporation or contributors be liable for any direct, +// indirect, incidental, special, exemplary, or consequential damages +// (including, but not limited to, procurement of substitute goods or services; +// loss of use, data, or profits; or business interruption) however caused +// and on any theory of liability, whether in contract, strict liability, +// or tort (including negligence or otherwise) arising in any way out of +// the use of this software, even if advised of the possibility of such damage. +// +//M*/ + +#ifndef __OPENCV_OPTIM_HPP__ +#define __OPENCV_OPTIM_HPP__ + +#include "opencv2/core.hpp" + +namespace cv{namespace optim +{ +//!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 + +#endif diff --git a/modules/optim/include/opencv2/optim/optim.hpp b/modules/optim/include/opencv2/optim/optim.hpp new file mode 100644 index 0000000000..b5a9ebf794 --- /dev/null +++ b/modules/optim/include/opencv2/optim/optim.hpp @@ -0,0 +1,48 @@ +/*M/////////////////////////////////////////////////////////////////////////////////////// +// +// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. +// +// By downloading, copying, installing or using the software you agree to this license. +// If you do not agree to this license, do not download, install, +// copy or use the software. +// +// +// License Agreement +// For Open Source Computer Vision Library +// +// Copyright (C) 2000-2008, Intel Corporation, all rights reserved. +// Copyright (C) 2009, Willow Garage Inc., all rights reserved. +// Copyright (C) 2013, OpenCV Foundation, all rights reserved. +// Third party copyrights are property of their respective owners. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// * Redistribution's of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// * Redistribution's in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// * The name of the copyright holders may not be used to endorse or promote products +// derived from this software without specific prior written permission. +// +// This software is provided by the copyright holders and contributors "as is" and +// any express or implied warranties, including, but not limited to, the implied +// warranties of merchantability and fitness for a particular purpose are disclaimed. +// In no event shall the Intel Corporation or contributors be liable for any direct, +// indirect, incidental, special, exemplary, or consequential damages +// (including, but not limited to, procurement of substitute goods or services; +// loss of use, data, or profits; or business interruption) however caused +// and on any theory of liability, whether in contract, strict liability, +// or tort (including negligence or otherwise) arising in any way out of +// the use of this software, even if advised of the possibility of such damage. +// +//M*/ + +#ifdef __OPENCV_BUILD +#error this is a compatibility header which should not be used inside the OpenCV library +#endif + +#include "opencv2/optim.hpp" diff --git a/modules/optim/src/lpsolver.cpp b/modules/optim/src/lpsolver.cpp new file mode 100644 index 0000000000..72db13a0c4 --- /dev/null +++ b/modules/optim/src/lpsolver.cpp @@ -0,0 +1,322 @@ +#include "precomp.hpp" +#include +#include +#include + +namespace cv{namespace optim{ +using std::vector; + +#ifdef ALEX_DEBUG +#define dprintf(x) printf x +static void print_matrix(const Mat& x){ + print(x); + printf("\n"); +} +static void print_simplex_state(const Mat& c,const Mat& b,double v,const std::vector N,const std::vector B){ + printf("\tprint simplex state\n"); + + printf("v=%g\n",v); + + printf("here c goes\n"); + print_matrix(c); + + printf("non-basic: "); + print(Mat(N)); + printf("\n"); + + printf("here b goes\n"); + print_matrix(b); + printf("basic: "); + + print(Mat(B)); + printf("\n"); +} +#else +#define dprintf(x) +#define print_matrix(x) +#define print_simplex_state(c,b,v,N,B) +#endif + +/**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. +*/ +static int initialize_simplex(Mat_& c, Mat_& b,double& v,vector& N,vector& B,vector& indexToRow); +static inline void pivot(Mat_& c,Mat_& b,double& v,vector& N,vector& B,int leaving_index, + int entering_index,vector& indexToRow); +/**@return SOLVELP_UNBOUNDED means the problem is unbdd, SOLVELP_MULTI means multiple solutions, SOLVELP_SINGLE means one solution. + */ +static int inner_simplex(Mat_& c, Mat_& b,double& v,vector& N,vector& B,vector& indexToRow); +static void swap_columns(Mat_& A,int col1,int col2); +#define SWAP(type,a,b) {type tmp=(a);(a)=(b);(b)=tmp;} + +//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){ + dprintf(("call to solveLP\n")); + + //sanity check (size, type, no. of channels) + CV_Assert(Func.type()==CV_64FC1 || Func.type()==CV_32FC1); + CV_Assert(Constr.type()==CV_64FC1 || Constr.type()==CV_32FC1); + CV_Assert((Func.rows==1 && (Constr.cols-Func.cols==1))|| + (Func.cols==1 && (Constr.cols-Func.rows==1))); + + //copy arguments for we will shall modify them + Mat_ bigC=Mat_(1,(Func.rows==1?Func.cols:Func.rows)+1), + bigB=Mat_(Constr.rows,Constr.cols+1); + if(Func.rows==1){ + Func.convertTo(bigC.colRange(1,bigC.cols),CV_64FC1); + }else{ + Mat FuncT=Func.t(); + FuncT.convertTo(bigC.colRange(1,bigC.cols),CV_64FC1); + } + Constr.convertTo(bigB.colRange(1,bigB.cols),CV_64FC1); + double v=0; + vector N,B; + vector indexToRow; + + if(initialize_simplex(bigC,bigB,v,N,B,indexToRow)==SOLVELP_UNFEASIBLE){ + return SOLVELP_UNFEASIBLE; + } + Mat_ c=bigC.colRange(1,bigC.cols), + b=bigB.colRange(1,bigB.cols); + + int res=0; + if((res=inner_simplex(c,b,v,N,B,indexToRow))==SOLVELP_UNBOUNDED){ + return SOLVELP_UNBOUNDED; + } + + //return the optimal solution + z.create(c.cols,1,CV_64FC1); + MatIterator_ it=z.begin(); + for(int i=1;i<=c.cols;i++,it++){ + if(indexToRow[i](indexToRow[i]-N.size(),b.cols-1); + } + } + + return res; +} + +static int initialize_simplex(Mat_& c, Mat_& b,double& v,vector& N,vector& B,vector& indexToRow){ + N.resize(c.cols); + N[0]=0; + for (std::vector::iterator it = N.begin()+1 ; it != N.end(); ++it){ + *it=it[-1]+1; + } + B.resize(b.rows); + B[0]=N.size(); + for (std::vector::iterator it = B.begin()+1 ; it != B.end(); ++it){ + *it=it[-1]+1; + } + indexToRow.resize(c.cols+b.rows); + indexToRow[0]=0; + for (std::vector::iterator it = indexToRow.begin()+1 ; it != indexToRow.end(); ++it){ + *it=it[-1]+1; + } + v=0; + + int k=0; + { + double min=DBL_MAX; + for(int i=0;i=0){ + N.erase(N.begin()); + for (std::vector::iterator it = indexToRow.begin()+1 ; it != indexToRow.end(); ++it){ + --(*it); + } + return 0; + } + + Mat_ old_c=c.clone(); + c=0; + c(0,0)=-1; + for(int i=0;i=N.size()){ + int iterator_offset=indexToRow[0]-N.size(); + if(b(iterator_offset,b.cols-1)>0){ + return SOLVELP_UNFEASIBLE; + } + pivot(c,b,v,N,B,iterator_offset,0,indexToRow); + } + + vector::iterator iterator; + { + int iterator_offset=indexToRow[0]; + iterator=N.begin()+iterator_offset; + std::iter_swap(iterator,N.begin()); + SWAP(int,indexToRow[*iterator],indexToRow[0]); + swap_columns(c,iterator_offset,0); + swap_columns(b,iterator_offset,0); + } + + dprintf(("after swaps\n")); + print_simplex_state(c,b,v,N,B); + + //start from 1, because we ignore x_0 + c=0; + v=0; + for(int I=1;I::iterator it = indexToRow.begin()+1 ; it != indexToRow.end(); ++it){ + --(*it); + } + return 0; +} + +static int inner_simplex(Mat_& c, Mat_& b,double& v,vector& N,vector& B,vector& indexToRow){ + int count=0; + for(;;){ + dprintf(("iteration #%d\n",count)); + count++; + + 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++){ + if(*pos_ptr==0){ + all_nonzero=false; + } + if(*pos_ptr>0){ + if(N[pos_ctr] 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, reinforcing Bland's rule + if((myite=it[e])>0){ + double val=it[b.cols-1]/myite; + if(val& c,Mat_& b,double& v,vector& N,vector& B, + int leaving_index,int entering_index,vector& indexToRow){ + double Coef=b(leaving_index,entering_index); + for(int i=0;i& A,int col1,int col2){ + for(int i=0;i + +TEST(Optim_LpSolver, regression_basic){ + cv::Mat A,B,z,etalon_z; + +#if 1 + //cormen's example #1 + A=(cv::Mat_(3,1)<<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"<(3,1)<<8,4,0); + ASSERT_EQ(cv::countNonZero(z!=etalon_z),0); +#endif + +#if 1 + //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"<(2,1)<<20,0); + ASSERT_EQ(cv::countNonZero(z!=etalon_z),0); +#endif + +#if 1 + //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"<(2,1)<<1,0); + ASSERT_EQ(cv::countNonZero(z!=etalon_z),0); +#endif +} + +TEST(Optim_LpSolver, regression_init_unfeasible){ + cv::Mat A,B,z,etalon_z; + +#if 1 + //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"<(3,1)<<1250,1000,0); + ASSERT_EQ(cv::countNonZero(z!=etalon_z),0); +#endif +} + +TEST(Optim_LpSolver, regression_absolutely_unfeasible){ + cv::Mat A,B,z,etalon_z; + +#if 1 + //trivial absolutely unfeasible example + A=(cv::Mat_(1,1)<<1); + B=(cv::Mat_(2,2)<<1,-1); + std::cout<<"here A goes\n"<(2,1)<<1,1); + B=(cv::Mat_(1,3)<<1,1,1); + std::cout<<"here A goes\n"<(4,1)<<10,-57,-9,-24); + B=(cv::Mat_(3,5)<<0.5,-5.5,-2.5,9,0,0.5,-1.5,-0.5,1,0,1,0,0,0,1); + std::cout<<"here A goes\n"<