Coverage for /usr/lib/python3/dist-packages/sympy/polys/matrices/linsolve.py: 12%

1# 

2# sympy.polys.matrices.linsolve module 

3# 

4# This module defines the _linsolve function which is the internal workhorse 

5# used by linsolve. This computes the solution of a system of linear equations 

6# using the SDM sparse matrix implementation in sympy.polys.matrices.sdm. This 

7# is a replacement for solve_lin_sys in sympy.polys.solvers which is 

8# inefficient for large sparse systems due to the use of a PolyRing with many 

9# generators: 

10# 

11# https://github.com/sympy/sympy/issues/20857 

12# 

13# The implementation of _linsolve here handles: 

14# 

15# - Extracting the coefficients from the Expr/Eq input equations. 

16# - Constructing a domain and converting the coefficients to 

17# that domain. 

18# - Using the SDM.rref, SDM.nullspace etc methods to generate the full 

19# solution working with arithmetic only in the domain of the coefficients. 

20# 

21# The routines here are particularly designed to be efficient for large sparse 

22# systems of linear equations although as well as dense systems. It is 

23# possible that for some small dense systems solve_lin_sys which uses the 

24# dense matrix implementation DDM will be more efficient. With smaller systems 

25# though the bulk of the time is spent just preprocessing the inputs and the 

26# relative time spent in rref is too small to be noticeable. 

27# 

28 

29from collections import defaultdict 

30 

31from sympy.core.add import Add 

32from sympy.core.mul import Mul 

33from sympy.core.singleton import S 

34 

35from sympy.polys.constructor import construct_domain 

36from sympy.polys.solvers import PolyNonlinearError 

37 

38from .sdm import ( 

39 SDM, 

40 sdm_irref, 

41 sdm_particular_from_rref, 

42 sdm_nullspace_from_rref 

43) 

44 

45from sympy.utilities.misc import filldedent 

46 

47 

48def _linsolve(eqs, syms): 

49 

50 """Solve a linear system of equations. 

51 

52 Examples 

53 ======== 

54 

55 Solve a linear system with a unique solution: 

56 

57 >>> from sympy import symbols, Eq 

58 >>> from sympy.polys.matrices.linsolve import _linsolve 

59 >>> x, y = symbols('x, y') 

60 >>> eqs = [Eq(x + y, 1), Eq(x - y, 2)] 

61 >>> _linsolve(eqs, [x, y]) 

62 {x: 3/2, y: -1/2} 

63 

64 In the case of underdetermined systems the solution will be expressed in 

65 terms of the unknown symbols that are unconstrained: 

66 

67 >>> _linsolve([Eq(x + y, 0)], [x, y]) 

68 {x: -y, y: y} 

69 

70 """ 

71 # Number of unknowns (columns in the non-augmented matrix) 

72 nsyms = len(syms) 

73 

74 # Convert to sparse augmented matrix (len(eqs) x (nsyms+1)) 

75 eqsdict, const = _linear_eq_to_dict(eqs, syms) 

76 Aaug = sympy_dict_to_dm(eqsdict, const, syms) 

77 K = Aaug.domain 

78 

79 # sdm_irref has issues with float matrices. This uses the ddm_rref() 

80 # function. When sdm_rref() can handle float matrices reasonably this 

81 # should be removed... 

82 if K.is_RealField or K.is_ComplexField: 

83 Aaug = Aaug.to_ddm().rref()[0].to_sdm() 

84 

85 # Compute reduced-row echelon form (RREF) 

86 Arref, pivots, nzcols = sdm_irref(Aaug) 

87 

88 # No solution: 

89 if pivots and pivots[-1] == nsyms: 

90 return None 

91 

92 # Particular solution for non-homogeneous system: 

93 P = sdm_particular_from_rref(Arref, nsyms+1, pivots) 

94 

95 # Nullspace - general solution to homogeneous system 

96 # Note: using nsyms not nsyms+1 to ignore last column 

97 V, nonpivots = sdm_nullspace_from_rref(Arref, K.one, nsyms, pivots, nzcols) 

98 

99 # Collect together terms from particular and nullspace: 

100 sol = defaultdict(list) 

101 for i, v in P.items(): 

102 sol[syms[i]].append(K.to_sympy(v)) 

103 for npi, Vi in zip(nonpivots, V): 

104 sym = syms[npi] 

105 for i, v in Vi.items(): 

106 sol[syms[i]].append(sym * K.to_sympy(v)) 

107 

108 # Use a single call to Add for each term: 

109 sol = {s: Add(*terms) for s, terms in sol.items()} 

110 

111 # Fill in the zeros: 

112 zero = S.Zero 

113 for s in set(syms) - set(sol): 

114 sol[s] = zero 

115 

116 # All done! 

117 return sol 

118 

119 

120def sympy_dict_to_dm(eqs_coeffs, eqs_rhs, syms): 

121 """Convert a system of dict equations to a sparse augmented matrix""" 

122 elems = set(eqs_rhs).union(*(e.values() for e in eqs_coeffs)) 

123 K, elems_K = construct_domain(elems, field=True, extension=True) 

124 elem_map = dict(zip(elems, elems_K)) 

125 neqs = len(eqs_coeffs) 

126 nsyms = len(syms) 

127 sym2index = dict(zip(syms, range(nsyms))) 

128 eqsdict = [] 

129 for eq, rhs in zip(eqs_coeffs, eqs_rhs): 

130 eqdict = {sym2index[s]: elem_map[c] for s, c in eq.items()} 

131 if rhs: 

132 eqdict[nsyms] = -elem_map[rhs] 

133 if eqdict: 

134 eqsdict.append(eqdict) 

135 sdm_aug = SDM(enumerate(eqsdict), (neqs, nsyms + 1), K) 

136 return sdm_aug 

137 

138 

139def _linear_eq_to_dict(eqs, syms): 

140 """Convert a system Expr/Eq equations into dict form, returning 

141 the coefficient dictionaries and a list of syms-independent terms 

142 from each expression in ``eqs```. 

143 

144 Examples 

145 ======== 

146 

147 >>> from sympy.polys.matrices.linsolve import _linear_eq_to_dict 

148 >>> from sympy.abc import x 

149 >>> _linear_eq_to_dict([2*x + 3], {x}) 

150 ([{x: 2}], [3]) 

151 """ 

152 coeffs = [] 

153 ind = [] 

154 symset = set(syms) 

155 for i, e in enumerate(eqs): 

156 if e.is_Equality: 

157 coeff, terms = _lin_eq2dict(e.lhs, symset) 

158 cR, tR = _lin_eq2dict(e.rhs, symset) 

159 # there were no nonlinear errors so now 

160 # cancellation is allowed 

161 coeff -= cR 

162 for k, v in tR.items(): 

163 if k in terms: 

164 terms[k] -= v 

165 else: 

166 terms[k] = -v 

167 # don't store coefficients of 0, however 

168 terms = {k: v for k, v in terms.items() if v} 

169 c, d = coeff, terms 

170 else: 

171 c, d = _lin_eq2dict(e, symset) 

172 coeffs.append(d) 

173 ind.append(c) 

174 return coeffs, ind 

175 

176 

177def _lin_eq2dict(a, symset): 

178 """return (c, d) where c is the sym-independent part of ``a`` and 

179 ``d`` is an efficiently calculated dictionary mapping symbols to 

180 their coefficients. A PolyNonlinearError is raised if non-linearity 

181 is detected. 

182 

183 The values in the dictionary will be non-zero. 

184 

185 Examples 

186 ======== 

187 

188 >>> from sympy.polys.matrices.linsolve import _lin_eq2dict 

189 >>> from sympy.abc import x, y 

190 >>> _lin_eq2dict(x + 2*y + 3, {x, y}) 

191 (3, {x: 1, y: 2}) 

192 """ 

193 if a in symset: 

194 return S.Zero, {a: S.One} 

195 elif a.is_Add: 

196 terms_list = defaultdict(list) 

197 coeff_list = [] 

198 for ai in a.args: 

199 ci, ti = _lin_eq2dict(ai, symset) 

200 coeff_list.append(ci) 

201 for mij, cij in ti.items(): 

202 terms_list[mij].append(cij) 

203 coeff = Add(*coeff_list) 

204 terms = {sym: Add(*coeffs) for sym, coeffs in terms_list.items()} 

205 return coeff, terms 

206 elif a.is_Mul: 

207 terms = terms_coeff = None 

208 coeff_list = [] 

209 for ai in a.args: 

210 ci, ti = _lin_eq2dict(ai, symset) 

211 if not ti: 

212 coeff_list.append(ci) 

213 elif terms is None: 

214 terms = ti 

215 terms_coeff = ci 

216 else: 

217 # since ti is not null and we already have 

218 # a term, this is a cross term 

219 raise PolyNonlinearError(filldedent(''' 

220 nonlinear cross-term: %s''' % a)) 

221 coeff = Mul._from_args(coeff_list) 

222 if terms is None: 

223 return coeff, {} 

224 else: 

225 terms = {sym: coeff * c for sym, c in terms.items()} 

226 return coeff * terms_coeff, terms 

227 elif not a.has_xfree(symset): 

228 return a, {} 

229 else: 

230 raise PolyNonlinearError('nonlinear term: %s' % a)