tulip-control / polytope

Geometric operations on polytopes of any dimension
https://pypi.org/project/polytope
Other
74 stars 19 forks source link

Zero Volume for 14D polytope #69

Closed alhaas closed 3 years ago

alhaas commented 3 years ago

Dear polytope devs,

I hope this is the right place for this. I have created an example similar to the one in the tutorial for a convex polytope in H-representation, however it gives me a volume of 0.0 and also, there is no warning in case the dimensions do not match. I think I am missing something here, as the same problem gives a nonzero volume with the Volesti package in R.

My Python version is 3.7.6.

import numpy as np
import polytope as pt

A = np.array([[0.01281, -0.00435, -0.00303,  0.00307, -0.00161, -0.00153,
         0.00728, -0.00249,  0.02289,  0.00303,  0.00186, -0.01623,
        -0.00289, -0.00246],
       [ 0.01246, -0.00457, -0.0033 ,  0.00266, -0.00236, -0.00184,
         0.0068 , -0.00278,  0.0233 ,  0.00258,  0.00129,  0.11378,
        -0.00316, -0.00274],
       [ 0.03503, -0.0095 , -0.013  , -0.01412,  0.15911, -0.01669,
        -0.00333, -0.01422,  0.06059, -0.02439, -0.07807,  0.1101 ,
        -0.01338, -0.01452],
       [-0.03503,  0.0095 ,  0.013  ,  0.01412, -0.15911,  0.01669,
         0.00333,  0.01422, -0.06059,  0.02439,  0.07807, -0.1101 ,
         0.01338,  0.01452],
       [-0.05537,  0.01528,  0.01336, -0.00266,  0.19761,  0.01101,
        -0.0184 ,  0.01247, -0.09187,  0.00186,  0.02847, -0.12141,
         0.01314,  0.01252]])

b = np.array([62.2,   62.2, 1416. ,  978. , 1059.])

p = pt.Polytope(A, b)       
print(p.volume)
johnyf commented 3 years ago

Thank you for reporting this error. I confirm that zero volume is computed using commit f981f434a270ee005076b9c534383aeaba314f40. The reason appears to be a missing else branch at (within the computation of the bounding box, which is in turn used to compute the polytope's volume):

https://github.com/tulip-control/polytope/blob/f981f434a270ee005076b9c534383aeaba314f40/polytope/polytope.py#L1304-L1306

The consequence of this missing branch is that if the optimization does not succeed, then the value of 0 is used for all bounds, because this is how the bounds are initialized within that function:

https://github.com/tulip-control/polytope/blob/f981f434a270ee005076b9c534383aeaba314f40/polytope/polytope.py#L1296

Running the lines of code from that area in ipython, when using scipy I get:

>>> sol
{'status': 3,
 'x': array([-6.37083e+07,  4.27238e+07,  4.38813e+06, -1.10456e+08,
        -1.77619e+08, -3.76229e+07, -1.46560e+08, -1.02404e+07,
        -1.70243e+08, -1.55426e+08, -3.99459e+08, -6.16730e+07,
         3.36853e+05, -7.07112e+08]),
 'fun': -707111604.4008977}

and when using cvxopt, I get:

>>> sol
{'status': 3, 'x': None, 'fun': None}

Quoting from the documentation of the function scipy.optimize.linprog:

3 : Problem appears to be unbounded.

and the documentation of the function cvxopt.solvers.lp says:

With the GLPK option, the solver does not return certificates of primal or dual infeasibility: if the status is 'primal infeasible' or 'dual infeasible', all entries of the output dictionary are None.

which is relevant to the lines:

https://github.com/tulip-control/polytope/blob/f981f434a270ee005076b9c534383aeaba314f40/polytope/solvers.py#L107

https://github.com/tulip-control/polytope/blob/f981f434a270ee005076b9c534383aeaba314f40/polytope/solvers.py#L115-L116

So:

I had been aware that unbounded polytopes are currently unsupported by the package polytope, and intended to consider that topic in a separate issue at some point in the future.

In conclusion, the simplest approach for now seems to be to raise an exception from within the function polytope.polytope.bounding_box if no solution is found (until unbounded polytopes become supported).

Regarding the polytope given above (https://github.com/tulip-control/polytope/issues/69#issue-902614042), if indeed it is unbounded, then its volume is infinite.

alhaas commented 3 years ago

Thank you very much for the exhaustive answer. Indeed, the problem was unboundend, sorry for the confusion, I tried to prepare a minimal working example which ended up being too minimal as the bounds stayed in the pulp variable declaration.

I now created a file including bounds which results in a nonzero volume, however a different value is computed each time the script is executed. Is this behavior problem specific and/or is there a way to change the solving algorithm or it's parameters in order to better "suit" this problem? example.txt

necozay commented 3 years ago

Volume algorithm implemented in polytope package is a randomized algorithm that tries to give an estimate of the volume. It might be unreliable especially in higher dimensions. You can try to increase the number of samples on the following line (this number should increase exponentially with the dimension to get similar accuracy): https://github.com/tulip-control/polytope/blob/f981f434a270ee005076b9c534383aeaba314f40/polytope/polytope.py#L1467 to improve accuracy. This should reduce the variation but there will still be some variation due to randomness. Increasing the samples will slow down the computation. We currently don't have an alternative deterministic volume computation algorithm implemented.

johnyf commented 3 years ago

Thank you @alhaas for the example. Thank you @necozay for the comment.

I fixed the bug due to the missing else branch mentioned above (https://github.com/tulip-control/polytope/issues/69#issuecomment-849040906) in commit 7f5964597d1bdc3bf59b8f18bdcce9087c928b21, which is now on branch master.

I added a parameter nsamples to the function polytope.polytope.volume in commit 2bd0a5bb14ed2a526aebf532f81f5cce0c74f750, which is now on branch dev. This can be used to set the number N mentioned above (https://github.com/tulip-control/polytope/issues/69#issuecomment-849696661), which controls the behavior of the randomized algorithm. Also, I changed the function polytope.polytope.volume in commit bfc6b0bb87ca6ff190fdf85c546e348a3e5ced18 to not use the precomputed value of volume (this is already done within the property volume of the classes Polytope and Region).

So now the function volume can be called multiple times, and also after modifying a polytope (not recommended to modify a Polytope or Region, though possible because these classes are mutable), to recompute the volume (without the need to set _volume to None to trigger recomputation of the volume). There are a number of other changes too on branch dev.

I confirm that the example code given above (https://github.com/tulip-control/polytope/issues/69#issuecomment-849667752) returns varying volume values in each run. Nonetheless, some of the values computed by different runs are: 8.071350273957078e+59, 7.371284178766923e+59, 9.059678878931413e+59. This is quite large an exponent. Even though the value varies, the variation is still a relatively small percentage of the computed value.

Using commit 5fb8a70e7c4f0ca8c16a46d0719b48e8f2ebb099 and changing the example's (https://github.com/tulip-control/polytope/issues/69#issuecomment-849667752) last lines to:

p = pt.Polytope(A, b)
pt.volume(p, nsamples=10**6)  # the default value is 10**4
print(p.volume)

returns the value 8.524334217903648e+59 with variability < 2 %. So the changes on branch dev appear to address the issue for user code.

alhaas commented 3 years ago

Thanks again for the feedback, from my side the issue is solved. If possible, maybe it would be an enhancement for the future to also give a seed to the volume function, this way results should become reproducible.

johnyf commented 3 years ago

Thank you for the suggestion @alhaas. I implemented this approach by adding in commit 99e209039ac1c3499aba7eb081be2f0316caae57 a parameter named seed to the function polytope.polytope.volume.

In more detail, the randomness within the function polytope.polytope.volume originates from the function numpy.random.rand:

https://github.com/tulip-control/polytope/blob/7f5964597d1bdc3bf59b8f18bdcce9087c928b21/polytope/polytope.py#L1501

The documentation of the function numpy.random.rand reads:

This is a convenience function for users porting code from Matlab, and wraps random_sample.

The function numpy.random.rand does not include a parameter for defining a seed. The documentation of the function numpy.random.random_sample reads:

New code should use the random method of a default_rng() instance instead;

The function numpy.random.random_sample does not include a parameter for defining a seed. The documentation of the function numpy.random.default_rng does include a parameter seed that can be an integer, among other types. Also relevant is the documentation of the class numpy.random.SeedSequence.

So the line of code linked-to above has been replaced by the line:

https://github.com/tulip-control/polytope/blob/99e209039ac1c3499aba7eb081be2f0316caae57/polytope/polytope.py#L1562

where seed is provided as an argument to the function polytope.polytope.volume.

Aside

The CI tests for Python 2.7 for commit 99e209039ac1c3499aba7eb081be2f0316caae57 failed, because the function numpy.random.default_rng does not exist, as it was added to numpy in a recent version.

The reason is that the last numpy version compatible with Python 2.7 was released about 1.5 years ago, on Dec 29, 2019 numpy == 1.16.6. Larger versions of numpy do not support Python 2.7 since Jul 26, 2019 (numpy == 1.17.0). I do not think there is any reason for polytope to remain compatible with Python 2.7. For that reason I opened #70.

johnyf commented 3 years ago

As of commit 9ee3f9634089c61353dfa4aec2d3d9b4cfb6905f, support for Python 2.7 has been removed, and also support for Python 3.5 and 3.6 has been removed. For motivation, please read the commit messages on branch dev.

The latest setuptools and numpy are now required on branch dev, the supported Python versions are 3.7, 3.8, 3.9 (, and soon 3.10), and the CI tests pass.

johnyf commented 3 years ago

As of commit a23ddcfd0deb42f6caf2c34bfec23ad39e89289a, the changes discussed above have been merged to the mainline branch, which addresses this issue.

Support for unbounded polytopes is an orthogonal issue, and will be considered in a separate issue that I plan to open in time.