Repository navigation
Expand file tree
/
Copy pathsolvers.py
More file actions
435 lines (341 loc) · 13.2 KB
/
Copy pathsolvers.py
File metadata and controls
435 lines (341 loc) · 13.2 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
"""Low-level linear systems solver. """
from sympy.utilities.exceptions import sympy_deprecation_warning
from sympy.utilities.iterables import connected_components
from sympy.core.sympify import sympify
from sympy.core.numbers import Integer, Rational
from sympy.matrices.dense import MutableDenseMatrix
from sympy.polys.domains import ZZ, QQ
from sympy.polys.domains import EX
from sympy.polys.rings import sring
from sympy.polys.polyerrors import NotInvertible
from sympy.polys.domainmatrix import DomainMatrix
class PolyNonlinearError(Exception):
"""Raised by solve_lin_sys for nonlinear equations"""
pass
class RawMatrix(MutableDenseMatrix):
"""
.. deprecated:: 1.9
This class fundamentally is broken by design. Use ``DomainMatrix`` if
you want a matrix over the polys domains or ``Matrix`` for a matrix
with ``Expr`` elements. The ``RawMatrix`` class will be removed/broken
in future in order to reestablish the invariant that the elements of a
Matrix should be of type ``Expr``.
"""
_sympify = staticmethod(lambda x: x)
def __init__(self, *args, **kwargs):
sympy_deprecation_warning(
"""
The RawMatrix class is deprecated. Use either DomainMatrix or
Matrix instead.
""",
deprecated_since_version="1.9",
active_deprecations_target="deprecated-rawmatrix",
)
domain = ZZ
for i in range(self.rows):
for j in range(self.cols):
val = self[i,j]
if getattr(val, 'is_Poly', False):
K = val.domain[val.gens]
val_sympy = val.as_expr()
elif hasattr(val, 'parent'):
K = val.parent()
val_sympy = K.to_sympy(val)
elif isinstance(val, (int, Integer)):
K = ZZ
val_sympy = sympify(val)
elif isinstance(val, Rational):
K = QQ
val_sympy = val
else:
for K in ZZ, QQ:
if K.of_type(val):
val_sympy = K.to_sympy(val)
break
else:
raise TypeError
domain = domain.unify(K)
self[i,j] = val_sympy
self.ring = domain
def eqs_to_matrix(eqs_coeffs, eqs_rhs, gens, domain):
"""Get matrix from linear equations in dict format.
Explanation
===========
Get the matrix representation of a system of linear equations represented
as dicts with low-level DomainElement coefficients. This is an
*internal* function that is used by solve_lin_sys.
Parameters
==========
eqs_coeffs: list[dict[Symbol, DomainElement]]
The left hand sides of the equations as dicts mapping from symbols to
coefficients where the coefficients are instances of
DomainElement.
eqs_rhs: list[DomainElements]
The right hand sides of the equations as instances of
DomainElement.
gens: list[Symbol]
The unknowns in the system of equations.
domain: Domain
The domain for coefficients of both lhs and rhs.
Returns
=======
The augmented matrix representation of the system as a DomainMatrix.
Examples
========
>>> from sympy import symbols, ZZ
>>> from sympy.polys.solvers import eqs_to_matrix
>>> x, y = symbols('x, y')
>>> eqs_coeff = [{x:ZZ(1), y:ZZ(1)}, {x:ZZ(1), y:ZZ(-1)}]
>>> eqs_rhs = [ZZ(0), ZZ(-1)]
>>> eqs_to_matrix(eqs_coeff, eqs_rhs, [x, y], ZZ)
DomainMatrix([[1, 1, 0], [1, -1, 1]], (2, 3), ZZ)
See also
========
solve_lin_sys: Uses :func:`~eqs_to_matrix` internally
"""
sym2index = {x: n for n, x in enumerate(gens)}
nrows = len(eqs_coeffs)
ncols = len(gens) + 1
rows = [[domain.zero] * ncols for _ in range(nrows)]
for row, eq_coeff, eq_rhs in zip(rows, eqs_coeffs, eqs_rhs):
for sym, coeff in eq_coeff.items():
row[sym2index[sym]] = domain.convert(coeff)
row[-1] = -domain.convert(eq_rhs)
return DomainMatrix(rows, (nrows, ncols), domain)
def sympy_eqs_to_ring(eqs, symbols):
"""Convert a system of equations from Expr to a PolyRing
Explanation
===========
High-level functions like ``solve`` expect Expr as inputs but can use
``solve_lin_sys`` internally. This function converts equations from
``Expr`` to the low-level poly types used by the ``solve_lin_sys``
function.
Parameters
==========
eqs: List of Expr
A list of equations as Expr instances
symbols: List of Symbol
A list of the symbols that are the unknowns in the system of
equations.
Returns
=======
Tuple[List[PolyElement], Ring]: The equations as PolyElement instances
and the ring of polynomials within which each equation is represented.
Examples
========
>>> from sympy import symbols
>>> from sympy.polys.solvers import sympy_eqs_to_ring
>>> a, x, y = symbols('a, x, y')
>>> eqs = [x-y, x+a*y]
>>> eqs_ring, ring = sympy_eqs_to_ring(eqs, [x, y])
>>> eqs_ring
[x - y, x + a*y]
>>> type(eqs_ring[0])
<class 'sympy.polys.rings.PolyElement'>
>>> ring
ZZ(a)[x,y]
With the equations in this form they can be passed to ``solve_lin_sys``:
>>> from sympy.polys.solvers import solve_lin_sys
>>> solve_lin_sys(eqs_ring, ring)
{y: 0, x: 0}
"""
try:
K, eqs_K = sring(eqs, symbols, field=True, extension=True)
except NotInvertible:
# https://github.com/sympy/sympy/issues/18874
K, eqs_K = sring(eqs, symbols, domain=EX)
return eqs_K, K.to_domain()
def solve_lin_sys(eqs, ring, _raw=True):
"""Solve a system of linear equations from a PolynomialRing
Explanation
===========
Solves a system of linear equations given as PolyElement instances of a
PolynomialRing. The basic arithmetic is carried out using instance of
DomainElement which is more efficient than :class:`~sympy.core.expr.Expr`
for the most common inputs.
While this is a public function it is intended primarily for internal use
so its interface is not necessarily convenient. Users are suggested to use
the :func:`sympy.solvers.solveset.linsolve` function (which uses this
function internally) instead.
Parameters
==========
eqs: list[PolyElement]
The linear equations to be solved as elements of a
PolynomialRing (assumed equal to zero).
ring: PolynomialRing
The polynomial ring from which eqs are drawn. The generators of this
ring are the unknowns to be solved for and the domain of the ring is
the domain of the coefficients of the system of equations.
_raw: bool
If *_raw* is False, the keys and values in the returned dictionary
will be of type Expr (and the unit of the field will be removed from
the keys) otherwise the low-level polys types will be returned, e.g.
PolyElement: PythonRational.
Returns
=======
``None`` if the system has no solution.
dict[Symbol, Expr] if _raw=False
dict[Symbol, DomainElement] if _raw=True.
Examples
========
>>> from sympy import symbols
>>> from sympy.polys.solvers import solve_lin_sys, sympy_eqs_to_ring
>>> x, y = symbols('x, y')
>>> eqs = [x - y, x + y - 2]
>>> eqs_ring, ring = sympy_eqs_to_ring(eqs, [x, y])
>>> solve_lin_sys(eqs_ring, ring)
{y: 1, x: 1}
Passing ``_raw=False`` returns the same result except that the keys are
``Expr`` rather than low-level poly types.
>>> solve_lin_sys(eqs_ring, ring, _raw=False)
{x: 1, y: 1}
See also
========
sympy_eqs_to_ring: prepares the inputs to ``solve_lin_sys``.
linsolve: ``linsolve`` uses ``solve_lin_sys`` internally.
sympy.solvers.solvers.solve: ``solve`` uses ``solve_lin_sys`` internally.
"""
as_expr = not _raw
assert ring.domain.is_Field
eqs_dict = [dict(eq) for eq in eqs]
one_monom = ring.one.monoms()[0]
zero = ring.domain.zero
eqs_rhs = []
eqs_coeffs = []
for eq_dict in eqs_dict:
eq_rhs = eq_dict.pop(one_monom, zero)
eq_coeffs = {}
for monom, coeff in eq_dict.items():
if sum(monom) != 1:
msg = "Nonlinear term encountered in solve_lin_sys"
raise PolyNonlinearError(msg)
eq_coeffs[ring.gens[monom.index(1)]] = coeff
if not eq_coeffs:
if not eq_rhs:
continue
else:
return None
eqs_rhs.append(eq_rhs)
eqs_coeffs.append(eq_coeffs)
result = _solve_lin_sys(eqs_coeffs, eqs_rhs, ring)
if result is not None and as_expr:
def to_sympy(x):
as_expr = getattr(x, 'as_expr', None)
if as_expr:
return as_expr()
else:
return ring.domain.to_sympy(x)
tresult = {to_sympy(sym): to_sympy(val) for sym, val in result.items()}
# Remove 1.0x
result = {}
for k, v in tresult.items():
if k.is_Mul:
c, s = k.as_coeff_Mul()
result[s] = v/c
else:
result[k] = v
return result
def _solve_lin_sys(eqs_coeffs, eqs_rhs, ring):
"""Solve a linear system from dict of PolynomialRing coefficients
Explanation
===========
This is an **internal** function used by :func:`solve_lin_sys` after the
equations have been preprocessed. The role of this function is to split
the system into connected components and pass those to
:func:`_solve_lin_sys_component`.
Examples
========
Setup a system for $x-y=0$ and $x+y=2$ and solve:
>>> from sympy import symbols, sring
>>> from sympy.polys.solvers import _solve_lin_sys
>>> x, y = symbols('x, y')
>>> R, (xr, yr) = sring([x, y], [x, y])
>>> eqs = [{xr:R.one, yr:-R.one}, {xr:R.one, yr:R.one}]
>>> eqs_rhs = [R.zero, -2*R.one]
>>> _solve_lin_sys(eqs, eqs_rhs, R)
{y: 1, x: 1}
See also
========
solve_lin_sys: This function is used internally by :func:`solve_lin_sys`.
"""
V = ring.gens
E = []
for eq_coeffs in eqs_coeffs:
syms = list(eq_coeffs)
E.extend(zip(syms[:-1], syms[1:]))
G = V, E
components = connected_components(G)
sym2comp = {}
for n, component in enumerate(components):
for sym in component:
sym2comp[sym] = n
subsystems = [([], []) for _ in range(len(components))]
for eq_coeff, eq_rhs in zip(eqs_coeffs, eqs_rhs):
sym = next(iter(eq_coeff), None)
sub_coeff, sub_rhs = subsystems[sym2comp[sym]]
sub_coeff.append(eq_coeff)
sub_rhs.append(eq_rhs)
sol = {}
for subsystem in subsystems:
subsol = _solve_lin_sys_component(subsystem[0], subsystem[1], ring)
if subsol is None:
return None
sol.update(subsol)
return sol
def _solve_lin_sys_component(eqs_coeffs, eqs_rhs, ring):
"""Solve a linear system from dict of PolynomialRing coefficients
Explanation
===========
This is an **internal** function used by :func:`solve_lin_sys` after the
equations have been preprocessed. After :func:`_solve_lin_sys` splits the
system into connected components this function is called for each
component. The system of equations is solved using Gauss-Jordan
elimination with division followed by back-substitution.
Examples
========
Setup a system for $x-y=0$ and $x+y=2$ and solve:
>>> from sympy import symbols, sring
>>> from sympy.polys.solvers import _solve_lin_sys_component
>>> x, y = symbols('x, y')
>>> R, (xr, yr) = sring([x, y], [x, y])
>>> eqs = [{xr:R.one, yr:-R.one}, {xr:R.one, yr:R.one}]
>>> eqs_rhs = [R.zero, -2*R.one]
>>> _solve_lin_sys_component(eqs, eqs_rhs, R)
{y: 1, x: 1}
See also
========
solve_lin_sys: This function is used internally by :func:`solve_lin_sys`.
"""
# transform from equations to matrix form
matrix = eqs_to_matrix(eqs_coeffs, eqs_rhs, ring.gens, ring.domain)
# convert to a field for rref
if not matrix.domain.is_Field:
matrix = matrix.to_field()
# solve by row-reduction
echelon, pivots = matrix.rref()
# construct the returnable form of the solutions
keys = ring.gens
if pivots and pivots[-1] == len(keys):
return None
if len(pivots) == len(keys):
sol = []
for s in [row[-1] for row in echelon.rep.to_ddm()]:
a = s
sol.append(a)
sols = dict(zip(keys, sol))
else:
sols = {}
g = ring.gens
# Extract ground domain coefficients and convert to the ring:
if hasattr(ring, 'ring'):
convert = ring.ring.ground_new
else:
convert = ring.ground_new
echelon = echelon.rep.to_ddm()
vals_set = {v for row in echelon for v in row}
vals_map = {v: convert(v) for v in vals_set}
echelon = [[vals_map[eij] for eij in ei] for ei in echelon]
for i, p in enumerate(pivots):
v = echelon[i][-1] - sum(echelon[i][j]*g[j] for j in range(p+1, len(g)) if echelon[i][j])
sols[keys[p]] = v
return sols