solveLP() infeasible results fix by rayonnant14 · Pull Request #19530 · opencv/opencv

rayonnant14 · 2021-02-15T08:59:25Z

fix #12343

Pull Request Readiness Checklist

See details at https://github.com/opencv/opencv/wiki/How_to_contribute#making-a-good-pull-request

I agree to contribute to the project under Apache 2 License.
To the best of my knowledge, the proposed patch is not based on a code under GPL or other license that is incompatible with OpenCV
The PR is proposed to proper branch
There is reference to original bug report and related work
There is accuracy test, performance test and test data in opencv_extra repository, if applicable
Patch to opencv_extra has the same branch name.
The feature is well documented and sample code can be built with the project CMake

alalek · 2021-02-17T19:45:22Z

Regression test case:

should fail without the patch
should pass with enabled patch.

Provided test case doesn't fail without the patch.

rayonnant14 · 2021-02-18T15:57:55Z

@alalek thanks for comment.
I made a change in the body of the test, it should be correct now

alalek

Thank you for update!

alalek · 2021-02-20T04:47:23Z

modules/core/src/lpsolver.cpp

        bool all_nonzero=true;
        for(pos_ptr=c.begin();pos_ptr!=c.end();pos_ptr++,pos_ctr++){
-            if(*pos_ptr==0){
+            if ((*pos_ptr >= -1e-12) && (*pos_ptr <= 1e-12))


avoid magic constants

use fabs() to compare near 0

const double SOLVELP_EPS = 1e-12; ... if (fabs(*pos_ptr) < SOLVELP_EPS)

alalek · 2021-02-20T04:49:41Z

modules/core/src/lpsolver.cpp

-            if(*pos_ptr==0){
+            if ((*pos_ptr >= -1e-12) && (*pos_ptr <= 1e-12))
+                *pos_ptr = 0.0;
+            if(*pos_ptr==0.0){


The second check can be eliminated:

if the first condition is true, then the second condition is true too.

if the first condition is false, then the second condition is false too.

For finite / infinite and NaN values.

alalek · 2021-02-20T04:54:23Z

modules/core/src/lpsolver.cpp

                double val=it[b.cols-1]/myite;
-                if(val<min || (val==min && B[row_it]<min_var)){
+                if((val<min) || ((val==min) && (B[row_it]<min_var))){


(not a part of this fix)
prefer multiplication over division.

(a/b < c) ==> (a < c*b) (works if b > 0)
(a/b == c) ==> (a == c*b)

alalek · 2021-02-20T05:01:12Z

modules/core/test/test_lpsolver.cpp

+                                       4., 0., 4., 1., 4.);
+
+    Mat z;
+    int result = cv::solveLP(A,B,z);


I believe, we still need some normalization or auto-tuning of SOLVELP_EPS value.

cv::solveLP(A,B,z); => OK

cv::solveLP(A,B*2,z); => OK

cv::solveLP(A,B*10,z); => OK

cv::solveLP(A,B*1e3,z); => OK

cv::solveLP(A,B*1e4,z); => FAIL

A can be multiplied by any positive scalar too.

BTW, use 2**k multiplier to minimize calculation errors in normalization (if normalization is going to be used).

alalek · 2021-02-20T05:03:00Z

modules/core/src/lpsolver.cpp

@@ -259,10 +268,12 @@ static int inner_simplex(Mat_<double>& c, Mat_<double>& b,double& v,vector<int>&
        int e=-1,pos_ctr=0,min_var=INT_MAX;


Please remove buggy "static" modifier from variable too

alalek · 2021-02-20T05:16:31Z

modules/core/src/lpsolver.cpp

    }
+    //check constraints feasibility
+    Mat prod = Constr(Rect(0, 0, Constr.cols - 1, Constr.rows)) * z;
+    Mat constr_check = (prod - 1e-12) > Constr.col(Constr.cols - 1);


Mat constr_check = Constr.col(Constr.cols - 1) - prod; double minValue = 0; minMaxIdx(constr_check, &minValue) // we need minValue only Analyze `minValue`

alalek · 2021-02-20T05:20:05Z

modules/core/include/opencv2/core/optim.hpp

 //! return codes for cv::solveLP() function
 enum SolveLPResult
 {
+    SOLVELP_LOST   = -3, //!< problem is feasible, but solver lost solution due to floating-point numbers representation error


"due to floating-point arithmetic errors"

alalek · 2021-06-07T18:43:29Z

modules/core/test/test_lpsolver.cpp

+    for (int i = 1; i < 1e5; i*=5)
+    {
+        Mat z;
+        int result = cv::solveLP(A,B*i,z);


In linear task both A/B should be multiplied to keep the same result.

No it's not true. 5x+10y < 25 is the same as x+2y < 5. So mathematically it's correct test, but the idea of test is not obvious. Is the provided matrix is bed defined for Simplex Method?

asmorkalov · 2021-08-03T06:55:38Z

@rayonnant14 Friendly reminder.

asmorkalov · 2021-08-23T07:41:11Z

@rayonnant14 Friendly reminder.

asmorkalov · 2021-10-05T10:33:58Z

@rayonnant14 Friednly reminder.

asmorkalov · 2022-05-23T09:11:48Z

@rayonnant14 Friednly reminder.

asmorkalov · 2023-05-24T09:57:00Z

modules/core/src/lpsolver.cpp

-            double myite=0;
-            //check constraints, select the tightest one, reinforcing Bland's rule
-            if((myite=it[e])>0){
+            if (fabs(it[b.cols-1]) <= SOLVELP_EPS_REL * fabs(it[b.cols-1]) + SOLVELP_EPS_ABS)


The check is equivalent to (1-SOLVELP_EPS_REL)*fabs(it[b.cols-1] <= SOLVELP_EPS_ABSorfabs(it[b.cols-1] < SOLVELP_EPS_ABS/0.2`. I do not think that it's expected.

fixed, added test

c762cfd

rayonnant14 added bug category: core labels Feb 15, 2021

asmorkalov requested a review from vpisarev February 15, 2021 10:37

fixed whitespaces

9d10454

test data change

065924d

rayonnant14 requested a review from alalek February 18, 2021 15:58

alalek reviewed Feb 20, 2021

View reviewed changes

change in error estimation

a4e2c6d

alalek reviewed Jun 7, 2021

View reviewed changes

vpisarev approved these changes May 19, 2023

View reviewed changes

asmorkalov reviewed May 24, 2023

View reviewed changes

asmorkalov mentioned this pull request May 24, 2023

Added final constrants check to solveLP to filter out flating-point numeric issues #23664

Merged

6 tasks

asmorkalov removed this from the 3.4.20 milestone May 24, 2023

asmorkalov closed this May 26, 2023

		@@ -259,10 +268,12 @@ static int inner_simplex(Mat_<double>& c, Mat_<double>& b,double& v,vector<int>&
		int e=-1,pos_ctr=0,min_var=INT_MAX;

Uh oh!

Conversation

rayonnant14 commented Feb 15, 2021

Pull Request Readiness Checklist

Uh oh!

alalek commented Feb 17, 2021

Uh oh!

rayonnant14 commented Feb 18, 2021

Uh oh!

alalek left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

asmorkalov commented Aug 3, 2021

Uh oh!

asmorkalov commented Aug 23, 2021

Uh oh!

asmorkalov commented Oct 5, 2021

Uh oh!

asmorkalov commented May 23, 2022

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants