linear programing function added #21687
Conversation
|
✅ Hi, I am the SymPy bot. I'm here to help you write a release notes entry. Please read the guide on how to write release notes. Your release notes are in good order. Here is what the release notes will look like:
This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.13. Click here to see the pull request description that was parsed.
Update The release notes on the wiki have been updated. |
|
Hi Oscar, |
|
Okay the pull request has been opened correctly this time! This looks nice but changes will be needed. The main things needed are:
|
|
Also you will need to edit the OP at the top of this PR to add a release note. |
|
I updated the post in #20391 to include some usage information as much as I recall. |
|
Hi all, In addition, I would like to know what realese notes should I add? (that is my first time contributing..) |
We can deal with the release notes at the end but we need to make sure something is entered before this is merged. For now though there seem to be some test failures relating to unused imports: |
|
Hi, |
|
Benchmark results from GitHub Actions Lower numbers are good, higher numbers are bad. A ratio less than 1 Significantly changed benchmark results (PR vs master) Significantly changed benchmark results (master vs previous release) before after ratio
[907895ac] [68022011]
- 4.83±0.02s 295±6ms 0.06 polygon.PolygonArbitraryPoint.time_bench01
Full benchmark results can be found as artifacts in GitHub Actions |
|
The only thing I need to add now are the realese notes. |
|
The docstring doesn't seem to follow the style guide: Can you show timings for this method? How long does it take for larger systems? What is the expected big-O complexity and does that match up with the timings? Does the method here only work for matrices with rational elements? What if there are expressions like @sylee957 do you want to review the algorithm here? I don't know linear programming well enough to comment in detail. |
|
Hi, About the complexity, I believe it is polynomial in the number of variables and constraints in the general case. As in terms of big-O complexity I cant realy tell. In my opinion it would be better to check it with large sysems |
|
If you use isympy or ipython then you can use In [36]: M = randMatrix(4, 4)
In [37]: b = randMatrix(4, 1)
In [38]: %time linsolve((M, b))
CPU times: user 12.1 ms, sys: 87 µs, total: 12.2 ms
Wall time: 12.3 ms
Out[38]:
⎧⎛19756036 49730 -21978115 -6380890 ⎞⎫
⎨⎜────────, ───────, ──────────, ─────────⎟⎬
⎩⎝3856167 3856167 3856167 3856167 ⎠⎭
In [39]: for n in [1, 2, 4, 8, 16, 32, 64, 128]:
...: M = randMatrix(n, n)
...: b = randMatrix(n, 1)
...: print('\n ------------- n =', n)
...: %time ok = linsolve((M, b))
...:
------------- n = 1
CPU times: user 3.8 ms, sys: 46 µs, total: 3.85 ms
Wall time: 3.9 ms
------------- n = 2
CPU times: user 5.53 ms, sys: 362 µs, total: 5.9 ms
Wall time: 5.71 ms
------------- n = 4
CPU times: user 14.8 ms, sys: 77 µs, total: 14.9 ms
Wall time: 15 ms
------------- n = 8
CPU times: user 42.1 ms, sys: 411 µs, total: 42.5 ms
Wall time: 43 ms
------------- n = 16
CPU times: user 161 ms, sys: 912 µs, total: 162 ms
Wall time: 163 ms
------------- n = 32
CPU times: user 711 ms, sys: 4.63 ms, total: 716 ms
Wall time: 722 ms
------------- n = 64
CPU times: user 4.3 s, sys: 33.8 ms, total: 4.34 s
Wall time: 4.39 s
------------- n = 128
CPU times: user 37.7 s, sys: 222 ms, total: 37.9 s
Wall time: 38.2 sHere |
Can you add tests for these? Are you sure it works for symbols? Can you show an example? There must be some cases where it fails for symbols... |
sympy/solvers/inequalities.py
Outdated
| for i in range(len(argmax)): | ||
| argmax[i] = nsimplify(argmax[i]) | ||
| for i in range(len(argmin_dual)): | ||
| argmin_dual[i] = nsimplify(argmin_dual[i]) |
There was a problem hiding this comment.
I don't think that this should use nsimplify
There was a problem hiding this comment.
but if I won't use nsimplify the results will be really long and messy (although they are still correct).
still want me to remove it?
There was a problem hiding this comment.
Can you give an example of the difference?
There was a problem hiding this comment.
this one, for example-
Expected:
(33/4, [33/4, 0, 0], [0, 0, 3/2])
Got:
(33/4, [8.25000000000000, 0, 0], [0, 0, 1.50000000000000])
(for the second example on the docstring)
There was a problem hiding this comment.
I think that's fine. If floats are used then the output can be floats. I don't know whether the algorithm needs to consider floating point error in the pivoting. If you want exact output you should use exact input.
There was a problem hiding this comment.
ok. I'll remove it.
|
I tried this test to check the comlexity: for n in [1, 2, 4, 8, 16, 32, 64, 128]: and this is few of the results I got: ------------- n = 1 ------------- n = 1 ------------- n = 1 it is not constant complexity but I would say that in most of the cases it outperformed |
I added a test in the docstring showing that the user can use all this types. |
Doctests are for testing the documentation. Code needs to be tested in the main tests |
ok. I will add some more tests. |
sympy/solvers/inequalities.py
Outdated
| piv_rows.append((ratio, i)) | ||
|
|
||
| if not piv_rows: | ||
| assert False, 'Unbounded' |
There was a problem hiding this comment.
You should redesign these exceptions assert False, 'Unbounded', assert False, 'Infeasible' because they are too prototypal.
I just left the design issues to the community, but these should be resolved before this is merged in the end product because it will have backward compatibility issues that makes you cannot change them easier later.
For example, scipy wraps the result in an object called OptimizeResult
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linprog.html
we could do the similar work, but I would favor more minimal design than that
How do I squash the commits? never did it before.. |
| >>> linear_programming(A, B, C, D) | ||
| (3, [2, 1, 0], [0, 1, 1]) | ||
| """ | ||
| M = Matrix([[A, B], [-C, D]]) |
There was a problem hiding this comment.
Sorry for yet another comment, but wouldn't it make sense to check the dimensions of the matrices before joining them?
I assume that there will be an error that the matrices cannot be joined or something, but that will be quite confusing for the average user. .rows and .cols can be used to get the corresponding sizes.
There was a problem hiding this comment.
Also the docstring should explain the shape constraints.
|
|
||
|
|
||
| class UnboundedLinearProgrammingError(Exception): | ||
| pass |
There was a problem hiding this comment.
There should be a docstring that explains what this exception class is for, perhaps with an example doctest showing when it would be raised.
| >>> linear_programming(A, B, C, D) | ||
| (11/3, [0, 1/3, 2/3], [0, 2/3, 1]) | ||
|
|
||
| The method also works with inegers, floats and symbolic expressions (like sqrt(2)). |
There was a problem hiding this comment.
When you say "symbolic expressions" that makes me think of expressions involving symbols. This should make clear that the method does not work for expressions containing symbols.
There was a problem hiding this comment.
integers is misspelt
| *) Minimizing y^{T}B constrained to y^{T}A >= C^{T} and y >= 0. | ||
|
|
||
| The method returns a triplet of solutions optimum, argmax and argmax_dual where | ||
| optimum is the optimum solution, argmax is x and argmax_dual is y. |
There was a problem hiding this comment.
There should be explicit parameters and returns sections in the docstring:
https://docs.sympy.org/latest/documentation-style-guide.html#docstring-sections
| >>> C = Matrix([[1, 1, 5]]) | ||
| >>> D = Matrix([0]) | ||
| >>> linear_programming(A, B, C, D) | ||
| (11/3, [0, 1/3, 2/3], [0, 2/3, 1]) |
There was a problem hiding this comment.
This demonstration would be slightly clearer if you extract the return values from the tuple like:
>>> minimum, argmax, argmax_dual = linear_programming(A, B, C, D)
>>> minimum
11/3
>>> argmax
...
>>> argmax_dual
...
Also I think it would be better to list a simpler example with a 2x2 system of inequalities and to clearly state in words the problem that is being solved e.g.:
Minimize $2x + 3y$ subject to $x - y <= 2$ and $x + y = 1$ ...
| r = r_orig.copy() | ||
| s = s_orig.copy() | ||
|
|
||
| M, r, s = _simplex(M, r, s) |
There was a problem hiding this comment.
I don't think it's necessary to pass in r and s here. The _simplex function can just make a list of integers and permute it then it would return something like:
r = [2, 3, 1, 0]
Then below you can just do something like:
argmax = [M[ri, -1] for ri in r]
|
|
||
| Examples | ||
| ======== | ||
| >>> from sympy.matrices.dense import Matrix |
There was a problem hiding this comment.
There should be a clank line after =====
| >>> C = Matrix([[1, 1, sqrt(5)]]) | ||
| >>> linear_programming(A, B, C, D) | ||
| (3, [2, 1, 0], [0, 1, 1]) | ||
| """ |
There was a problem hiding this comment.
There should be a "see also" section with references to related functions like I guess solve_univariate_inequality or anything else that someone might want to use instead of this for their problem.
| D = Matrix([0]) | ||
| optimum, argmax, argmax_dual = linear_programming(A, B, C, D) | ||
| assert optimum == -2 | ||
| A = Matrix([[1, -1, 2], [-1, 2, -3], [2, 1, -7]]) |
There was a problem hiding this comment.
The tests should check what is returned for argmax and argmax_dual as well.
Add a blank line between the test cases because this function is quite dense.
| MM = Matrix.zeros(M.rows, M.cols) | ||
| for ii in range(M.rows): | ||
| for jj in range(M.cols): | ||
| if ii == i and jj == j: | ||
| MM[ii, jj] = 1 / M[i, j] | ||
| elif ii == i and jj != j: | ||
| MM[ii, jj] = M[ii, jj] / M[i, j] | ||
| elif ii != i and jj == j: | ||
| MM[ii, jj] = -M[ii, jj] / M[i, j] | ||
| else: | ||
| MM[ii, jj] = M[ii, jj] - M[ii, j] * M[i, jj] / M[i, j] |
There was a problem hiding this comment.
This can be more vectorised. Something like:
def _pivot2(M, i, j):
Mi, Mj, Mij = M[i,:], M[:,j], M[i,j]
MM = M - Mj * (Mi / Mij)
MM[i,:] = Mi / Mij
MM[:,j] = -Mj / Mij
MM[i,j] = 1 / Mij
return MMFor large matrices this is faster e.g.
In [25]: M = randMatrix(100)
In [26]: %time ok = _pivot(M, 2, 2)
CPU times: user 696 ms, sys: 5.19 ms, total: 701 ms
Wall time: 704 ms
In [27]: %time ok = _pivot2(M, 2, 2)
CPU times: user 93.6 ms, sys: 3.81 ms, total: 97.4 ms
Wall time: 98.8 ms| piv_rows.append((ratio, i)) | ||
|
|
||
| if not piv_rows: | ||
| raise UnboundedLinearProgrammingError('Objective function can assume arbitrarily large positive (resp. negative) values at feasible vectors!') |
There was a problem hiding this comment.
This line is too long
|
|
||
| if not piv_rows: | ||
| raise UnboundedLinearProgrammingError('Objective function can assume arbitrarily large positive (resp. negative) values at feasible vectors!') | ||
| piv_rows = sorted(piv_rows, key=lambda x: (x[0], x[1])) |
There was a problem hiding this comment.
Is the key function needed? Isn't that just how tuples automatically compare?
| if not piv_rows: | ||
| raise UnboundedLinearProgrammingError('Objective function can assume arbitrarily large positive (resp. negative) values at feasible vectors!') | ||
| piv_rows = sorted(piv_rows, key=lambda x: (x[0], x[1])) | ||
| piv_rows = [(ratio, i) for ratio, i in piv_rows if ratio == piv_rows[0][0]] |
There was a problem hiding this comment.
If piv_rows is sorted then it should be necessary to loop through the whole list. Only the first O(1) items are needed.
| piv_cols = [] | ||
| for j in range(A.cols): | ||
| if A[k, j] < 0: | ||
| piv_cols.append(j) |
There was a problem hiding this comment.
piv_cols = [j for j in range(A.cols) if A[k, j] < 0]| """ | ||
| Simplex method with randomized pivoting | ||
| http://web.tecnico.ulisboa.pt/mcasquilho/acad/or/ftp/FergusonUCLA_LP.pdf | ||
| """ |
There was a problem hiding this comment.
Include a proper reference: title, author etc rather than just a URL. A URL is nice but if the URL stops working it should still be possible to identify the document that you are referring to.
|
If you install pytest and pytest-cov then you can measure code coverage with: Opening |
| >>> linear_programming(A, B, C, D) | ||
| (11/3, [0, 1/3, 2/3], [0, 2/3, 1]) | ||
|
|
||
| The method also works with inegers, floats and symbolic expressions (like sqrt(2)). |
There was a problem hiding this comment.
| The method also works with inegers, floats and symbolic expressions (like sqrt(2)). | |
| The method also works with non-integer numeric coefficients. |
|
I would like to work on this. I think this could be a big part of my GSoC project. As @oscarbenjamin mentions in #20391 (comment), this can be used to handle relational assumptions which is what my GSoC project is about. |
This commit includes the following improvements: - Streamlined _simplex for improved conciseness and efficiency. - Add an option to skip the second phase of the simplex method, useful for the new find_feasible function that only requires the first phase. - Move column and row labeling management entirely within _simplex to avoid duplicate code between find_feasible and linear_programming. - Enhance docstrings for several functions. - Add new test cases including test cases covering non-square matrices and constraints with ">=" relations. - Extend _to_standard_form to support AppliedBinaryRelation from the new assumptions.
References to other Issues or PRs
cf #22389
Brief description of what is fixed or changed
Other comments
Release Notes
linear_programmingfunction added.