Coverage for /usr/lib/python3/dist-packages/sympy/plotting/experimental_lambdify.py: 18%

237 statements  

« prev     ^ index     » next       coverage.py v7.9.1, created at 2025-06-14 15:55 +0200

1""" rewrite of lambdify - This stuff is not stable at all. 

2 

3It is for internal use in the new plotting module. 

4It may (will! see the Q'n'A in the source) be rewritten. 

5 

6It's completely self contained. Especially it does not use lambdarepr. 

7 

8It does not aim to replace the current lambdify. Most importantly it will never 

9ever support anything else than SymPy expressions (no Matrices, dictionaries 

10and so on). 

11""" 

12 

13 

14import re 

15from sympy.core.numbers import (I, NumberSymbol, oo, zoo) 

16from sympy.core.symbol import Symbol 

17from sympy.utilities.iterables import numbered_symbols 

18 

19# We parse the expression string into a tree that identifies functions. Then 

20# we translate the names of the functions and we translate also some strings 

21# that are not names of functions (all this according to translation 

22# dictionaries). 

23# If the translation goes to another module (like numpy) the 

24# module is imported and 'func' is translated to 'module.func'. 

25# If a function can not be translated, the inner nodes of that part of the 

26# tree are not translated. So if we have Integral(sqrt(x)), sqrt is not 

27# translated to np.sqrt and the Integral does not crash. 

28# A namespace for all this is generated by crawling the (func, args) tree of 

29# the expression. The creation of this namespace involves many ugly 

30# workarounds. 

31# The namespace consists of all the names needed for the SymPy expression and 

32# all the name of modules used for translation. Those modules are imported only 

33# as a name (import numpy as np) in order to keep the namespace small and 

34# manageable. 

35 

36# Please, if there is a bug, do not try to fix it here! Rewrite this by using 

37# the method proposed in the last Q'n'A below. That way the new function will 

38# work just as well, be just as simple, but it wont need any new workarounds. 

39# If you insist on fixing it here, look at the workarounds in the function 

40# sympy_expression_namespace and in lambdify. 

41 

42# Q: Why are you not using Python abstract syntax tree? 

43# A: Because it is more complicated and not much more powerful in this case. 

44 

45# Q: What if I have Symbol('sin') or g=Function('f')? 

46# A: You will break the algorithm. We should use srepr to defend against this? 

47# The problem with Symbol('sin') is that it will be printed as 'sin'. The 

48# parser will distinguish it from the function 'sin' because functions are 

49# detected thanks to the opening parenthesis, but the lambda expression won't 

50# understand the difference if we have also the sin function. 

51# The solution (complicated) is to use srepr and maybe ast. 

52# The problem with the g=Function('f') is that it will be printed as 'f' but in 

53# the global namespace we have only 'g'. But as the same printer is used in the 

54# constructor of the namespace there will be no problem. 

55 

56# Q: What if some of the printers are not printing as expected? 

57# A: The algorithm wont work. You must use srepr for those cases. But even 

58# srepr may not print well. All problems with printers should be considered 

59# bugs. 

60 

61# Q: What about _imp_ functions? 

62# A: Those are taken care for by evalf. A special case treatment will work 

63# faster but it's not worth the code complexity. 

64 

65# Q: Will ast fix all possible problems? 

66# A: No. You will always have to use some printer. Even srepr may not work in 

67# some cases. But if the printer does not work, that should be considered a 

68# bug. 

69 

70# Q: Is there same way to fix all possible problems? 

71# A: Probably by constructing our strings ourself by traversing the (func, 

72# args) tree and creating the namespace at the same time. That actually sounds 

73# good. 

74 

75from sympy.external import import_module 

76import warnings 

77 

78#TODO debugging output 

79 

80 

81class vectorized_lambdify: 

82 """ Return a sufficiently smart, vectorized and lambdified function. 

83 

84 Returns only reals. 

85 

86 Explanation 

87 =========== 

88 

89 This function uses experimental_lambdify to created a lambdified 

90 expression ready to be used with numpy. Many of the functions in SymPy 

91 are not implemented in numpy so in some cases we resort to Python cmath or 

92 even to evalf. 

93 

94 The following translations are tried: 

95 only numpy complex 

96 - on errors raised by SymPy trying to work with ndarray: 

97 only Python cmath and then vectorize complex128 

98 

99 When using Python cmath there is no need for evalf or float/complex 

100 because Python cmath calls those. 

101 

102 This function never tries to mix numpy directly with evalf because numpy 

103 does not understand SymPy Float. If this is needed one can use the 

104 float_wrap_evalf/complex_wrap_evalf options of experimental_lambdify or 

105 better one can be explicit about the dtypes that numpy works with. 

106 Check numpy bug http://projects.scipy.org/numpy/ticket/1013 to know what 

107 types of errors to expect. 

108 """ 

109 def __init__(self, args, expr): 

110 self.args = args 

111 self.expr = expr 

112 self.np = import_module('numpy') 

113 

114 self.lambda_func_1 = experimental_lambdify( 

115 args, expr, use_np=True) 

116 self.vector_func_1 = self.lambda_func_1 

117 

118 self.lambda_func_2 = experimental_lambdify( 

119 args, expr, use_python_cmath=True) 

120 self.vector_func_2 = self.np.vectorize( 

121 self.lambda_func_2, otypes=[complex]) 

122 

123 self.vector_func = self.vector_func_1 

124 self.failure = False 

125 

126 def __call__(self, *args): 

127 np = self.np 

128 

129 try: 

130 temp_args = (np.array(a, dtype=complex) for a in args) 

131 results = self.vector_func(*temp_args) 

132 results = np.ma.masked_where( 

133 np.abs(results.imag) > 1e-7 * np.abs(results), 

134 results.real, copy=False) 

135 return results 

136 except ValueError: 

137 if self.failure: 

138 raise 

139 

140 self.failure = True 

141 self.vector_func = self.vector_func_2 

142 warnings.warn( 

143 'The evaluation of the expression is problematic. ' 

144 'We are trying a failback method that may still work. ' 

145 'Please report this as a bug.') 

146 return self.__call__(*args) 

147 

148 

149class lambdify: 

150 """Returns the lambdified function. 

151 

152 Explanation 

153 =========== 

154 

155 This function uses experimental_lambdify to create a lambdified 

156 expression. It uses cmath to lambdify the expression. If the function 

157 is not implemented in Python cmath, Python cmath calls evalf on those 

158 functions. 

159 """ 

160 

161 def __init__(self, args, expr): 

162 self.args = args 

163 self.expr = expr 

164 self.lambda_func_1 = experimental_lambdify( 

165 args, expr, use_python_cmath=True, use_evalf=True) 

166 self.lambda_func_2 = experimental_lambdify( 

167 args, expr, use_python_math=True, use_evalf=True) 

168 self.lambda_func_3 = experimental_lambdify( 

169 args, expr, use_evalf=True, complex_wrap_evalf=True) 

170 self.lambda_func = self.lambda_func_1 

171 self.failure = False 

172 

173 def __call__(self, args): 

174 try: 

175 #The result can be sympy.Float. Hence wrap it with complex type. 

176 result = complex(self.lambda_func(args)) 

177 if abs(result.imag) > 1e-7 * abs(result): 

178 return None 

179 return result.real 

180 except (ZeroDivisionError, OverflowError): 

181 return None 

182 except TypeError as e: 

183 if self.failure: 

184 raise e 

185 

186 if self.lambda_func == self.lambda_func_1: 

187 self.lambda_func = self.lambda_func_2 

188 return self.__call__(args) 

189 

190 self.failure = True 

191 self.lambda_func = self.lambda_func_3 

192 warnings.warn( 

193 'The evaluation of the expression is problematic. ' 

194 'We are trying a failback method that may still work. ' 

195 'Please report this as a bug.', stacklevel=2) 

196 return self.__call__(args) 

197 

198 

199def experimental_lambdify(*args, **kwargs): 

200 l = Lambdifier(*args, **kwargs) 

201 return l 

202 

203 

204class Lambdifier: 

205 def __init__(self, args, expr, print_lambda=False, use_evalf=False, 

206 float_wrap_evalf=False, complex_wrap_evalf=False, 

207 use_np=False, use_python_math=False, use_python_cmath=False, 

208 use_interval=False): 

209 

210 self.print_lambda = print_lambda 

211 self.use_evalf = use_evalf 

212 self.float_wrap_evalf = float_wrap_evalf 

213 self.complex_wrap_evalf = complex_wrap_evalf 

214 self.use_np = use_np 

215 self.use_python_math = use_python_math 

216 self.use_python_cmath = use_python_cmath 

217 self.use_interval = use_interval 

218 

219 # Constructing the argument string 

220 # - check 

221 if not all(isinstance(a, Symbol) for a in args): 

222 raise ValueError('The arguments must be Symbols.') 

223 # - use numbered symbols 

224 syms = numbered_symbols(exclude=expr.free_symbols) 

225 newargs = [next(syms) for _ in args] 

226 expr = expr.xreplace(dict(zip(args, newargs))) 

227 argstr = ', '.join([str(a) for a in newargs]) 

228 del syms, newargs, args 

229 

230 # Constructing the translation dictionaries and making the translation 

231 self.dict_str = self.get_dict_str() 

232 self.dict_fun = self.get_dict_fun() 

233 exprstr = str(expr) 

234 newexpr = self.tree2str_translate(self.str2tree(exprstr)) 

235 

236 # Constructing the namespaces 

237 namespace = {} 

238 namespace.update(self.sympy_atoms_namespace(expr)) 

239 namespace.update(self.sympy_expression_namespace(expr)) 

240 # XXX Workaround 

241 # Ugly workaround because Pow(a,Half) prints as sqrt(a) 

242 # and sympy_expression_namespace can not catch it. 

243 from sympy.functions.elementary.miscellaneous import sqrt 

244 namespace.update({'sqrt': sqrt}) 

245 namespace.update({'Eq': lambda x, y: x == y}) 

246 namespace.update({'Ne': lambda x, y: x != y}) 

247 # End workaround. 

248 if use_python_math: 

249 namespace.update({'math': __import__('math')}) 

250 if use_python_cmath: 

251 namespace.update({'cmath': __import__('cmath')}) 

252 if use_np: 

253 try: 

254 namespace.update({'np': __import__('numpy')}) 

255 except ImportError: 

256 raise ImportError( 

257 'experimental_lambdify failed to import numpy.') 

258 if use_interval: 

259 namespace.update({'imath': __import__( 

260 'sympy.plotting.intervalmath', fromlist=['intervalmath'])}) 

261 namespace.update({'math': __import__('math')}) 

262 

263 # Construct the lambda 

264 if self.print_lambda: 

265 print(newexpr) 

266 eval_str = 'lambda %s : ( %s )' % (argstr, newexpr) 

267 self.eval_str = eval_str 

268 exec("MYNEWLAMBDA = %s" % eval_str, namespace) 

269 self.lambda_func = namespace['MYNEWLAMBDA'] 

270 

271 def __call__(self, *args, **kwargs): 

272 return self.lambda_func(*args, **kwargs) 

273 

274 

275 ############################################################################## 

276 # Dicts for translating from SymPy to other modules 

277 ############################################################################## 

278 ### 

279 # builtins 

280 ### 

281 # Functions with different names in builtins 

282 builtin_functions_different = { 

283 'Min': 'min', 

284 'Max': 'max', 

285 'Abs': 'abs', 

286 } 

287 

288 # Strings that should be translated 

289 builtin_not_functions = { 

290 'I': '1j', 

291# 'oo': '1e400', 

292 } 

293 

294 ### 

295 # numpy 

296 ### 

297 

298 # Functions that are the same in numpy 

299 numpy_functions_same = [ 

300 'sin', 'cos', 'tan', 'sinh', 'cosh', 'tanh', 'exp', 'log', 

301 'sqrt', 'floor', 'conjugate', 

302 ] 

303 

304 # Functions with different names in numpy 

305 numpy_functions_different = { 

306 "acos": "arccos", 

307 "acosh": "arccosh", 

308 "arg": "angle", 

309 "asin": "arcsin", 

310 "asinh": "arcsinh", 

311 "atan": "arctan", 

312 "atan2": "arctan2", 

313 "atanh": "arctanh", 

314 "ceiling": "ceil", 

315 "im": "imag", 

316 "ln": "log", 

317 "Max": "amax", 

318 "Min": "amin", 

319 "re": "real", 

320 "Abs": "abs", 

321 } 

322 

323 # Strings that should be translated 

324 numpy_not_functions = { 

325 'pi': 'np.pi', 

326 'oo': 'np.inf', 

327 'E': 'np.e', 

328 } 

329 

330 ### 

331 # Python math 

332 ### 

333 

334 # Functions that are the same in math 

335 math_functions_same = [ 

336 'sin', 'cos', 'tan', 'asin', 'acos', 'atan', 'atan2', 

337 'sinh', 'cosh', 'tanh', 'asinh', 'acosh', 'atanh', 

338 'exp', 'log', 'erf', 'sqrt', 'floor', 'factorial', 'gamma', 

339 ] 

340 

341 # Functions with different names in math 

342 math_functions_different = { 

343 'ceiling': 'ceil', 

344 'ln': 'log', 

345 'loggamma': 'lgamma' 

346 } 

347 

348 # Strings that should be translated 

349 math_not_functions = { 

350 'pi': 'math.pi', 

351 'E': 'math.e', 

352 } 

353 

354 ### 

355 # Python cmath 

356 ### 

357 

358 # Functions that are the same in cmath 

359 cmath_functions_same = [ 

360 'sin', 'cos', 'tan', 'asin', 'acos', 'atan', 

361 'sinh', 'cosh', 'tanh', 'asinh', 'acosh', 'atanh', 

362 'exp', 'log', 'sqrt', 

363 ] 

364 

365 # Functions with different names in cmath 

366 cmath_functions_different = { 

367 'ln': 'log', 

368 'arg': 'phase', 

369 } 

370 

371 # Strings that should be translated 

372 cmath_not_functions = { 

373 'pi': 'cmath.pi', 

374 'E': 'cmath.e', 

375 } 

376 

377 ### 

378 # intervalmath 

379 ### 

380 

381 interval_not_functions = { 

382 'pi': 'math.pi', 

383 'E': 'math.e' 

384 } 

385 

386 interval_functions_same = [ 

387 'sin', 'cos', 'exp', 'tan', 'atan', 'log', 

388 'sqrt', 'cosh', 'sinh', 'tanh', 'floor', 

389 'acos', 'asin', 'acosh', 'asinh', 'atanh', 

390 'Abs', 'And', 'Or' 

391 ] 

392 

393 interval_functions_different = { 

394 'Min': 'imin', 

395 'Max': 'imax', 

396 'ceiling': 'ceil', 

397 

398 } 

399 

400 ### 

401 # mpmath, etc 

402 ### 

403 #TODO 

404 

405 ### 

406 # Create the final ordered tuples of dictionaries 

407 ### 

408 

409 # For strings 

410 def get_dict_str(self): 

411 dict_str = dict(self.builtin_not_functions) 

412 if self.use_np: 

413 dict_str.update(self.numpy_not_functions) 

414 if self.use_python_math: 

415 dict_str.update(self.math_not_functions) 

416 if self.use_python_cmath: 

417 dict_str.update(self.cmath_not_functions) 

418 if self.use_interval: 

419 dict_str.update(self.interval_not_functions) 

420 return dict_str 

421 

422 # For functions 

423 def get_dict_fun(self): 

424 dict_fun = dict(self.builtin_functions_different) 

425 if self.use_np: 

426 for s in self.numpy_functions_same: 

427 dict_fun[s] = 'np.' + s 

428 for k, v in self.numpy_functions_different.items(): 

429 dict_fun[k] = 'np.' + v 

430 if self.use_python_math: 

431 for s in self.math_functions_same: 

432 dict_fun[s] = 'math.' + s 

433 for k, v in self.math_functions_different.items(): 

434 dict_fun[k] = 'math.' + v 

435 if self.use_python_cmath: 

436 for s in self.cmath_functions_same: 

437 dict_fun[s] = 'cmath.' + s 

438 for k, v in self.cmath_functions_different.items(): 

439 dict_fun[k] = 'cmath.' + v 

440 if self.use_interval: 

441 for s in self.interval_functions_same: 

442 dict_fun[s] = 'imath.' + s 

443 for k, v in self.interval_functions_different.items(): 

444 dict_fun[k] = 'imath.' + v 

445 return dict_fun 

446 

447 ############################################################################## 

448 # The translator functions, tree parsers, etc. 

449 ############################################################################## 

450 

451 def str2tree(self, exprstr): 

452 """Converts an expression string to a tree. 

453 

454 Explanation 

455 =========== 

456 

457 Functions are represented by ('func_name(', tree_of_arguments). 

458 Other expressions are (head_string, mid_tree, tail_str). 

459 Expressions that do not contain functions are directly returned. 

460 

461 Examples 

462 ======== 

463 

464 >>> from sympy.abc import x, y, z 

465 >>> from sympy import Integral, sin 

466 >>> from sympy.plotting.experimental_lambdify import Lambdifier 

467 >>> str2tree = Lambdifier([x], x).str2tree 

468 

469 >>> str2tree(str(Integral(x, (x, 1, y)))) 

470 ('', ('Integral(', 'x, (x, 1, y)'), ')') 

471 >>> str2tree(str(x+y)) 

472 'x + y' 

473 >>> str2tree(str(x+y*sin(z)+1)) 

474 ('x + y*', ('sin(', 'z'), ') + 1') 

475 >>> str2tree('sin(y*(y + 1.1) + (sin(y)))') 

476 ('', ('sin(', ('y*(y + 1.1) + (', ('sin(', 'y'), '))')), ')') 

477 """ 

478 #matches the first 'function_name(' 

479 first_par = re.search(r'(\w+\()', exprstr) 

480 if first_par is None: 

481 return exprstr 

482 else: 

483 start = first_par.start() 

484 end = first_par.end() 

485 head = exprstr[:start] 

486 func = exprstr[start:end] 

487 tail = exprstr[end:] 

488 count = 0 

489 for i, c in enumerate(tail): 

490 if c == '(': 

491 count += 1 

492 elif c == ')': 

493 count -= 1 

494 if count == -1: 

495 break 

496 func_tail = self.str2tree(tail[:i]) 

497 tail = self.str2tree(tail[i:]) 

498 return (head, (func, func_tail), tail) 

499 

500 @classmethod 

501 def tree2str(cls, tree): 

502 """Converts a tree to string without translations. 

503 

504 Examples 

505 ======== 

506 

507 >>> from sympy.abc import x, y, z 

508 >>> from sympy import sin 

509 >>> from sympy.plotting.experimental_lambdify import Lambdifier 

510 >>> str2tree = Lambdifier([x], x).str2tree 

511 >>> tree2str = Lambdifier([x], x).tree2str 

512 

513 >>> tree2str(str2tree(str(x+y*sin(z)+1))) 

514 'x + y*sin(z) + 1' 

515 """ 

516 if isinstance(tree, str): 

517 return tree 

518 else: 

519 return ''.join(map(cls.tree2str, tree)) 

520 

521 def tree2str_translate(self, tree): 

522 """Converts a tree to string with translations. 

523 

524 Explanation 

525 =========== 

526 

527 Function names are translated by translate_func. 

528 Other strings are translated by translate_str. 

529 """ 

530 if isinstance(tree, str): 

531 return self.translate_str(tree) 

532 elif isinstance(tree, tuple) and len(tree) == 2: 

533 return self.translate_func(tree[0][:-1], tree[1]) 

534 else: 

535 return ''.join([self.tree2str_translate(t) for t in tree]) 

536 

537 def translate_str(self, estr): 

538 """Translate substrings of estr using in order the dictionaries in 

539 dict_tuple_str.""" 

540 for pattern, repl in self.dict_str.items(): 

541 estr = re.sub(pattern, repl, estr) 

542 return estr 

543 

544 def translate_func(self, func_name, argtree): 

545 """Translate function names and the tree of arguments. 

546 

547 Explanation 

548 =========== 

549 

550 If the function name is not in the dictionaries of dict_tuple_fun then the 

551 function is surrounded by a float((...).evalf()). 

552 

553 The use of float is necessary as np.<function>(sympy.Float(..)) raises an 

554 error.""" 

555 if func_name in self.dict_fun: 

556 new_name = self.dict_fun[func_name] 

557 argstr = self.tree2str_translate(argtree) 

558 return new_name + '(' + argstr 

559 elif func_name in ['Eq', 'Ne']: 

560 op = {'Eq': '==', 'Ne': '!='} 

561 return "(lambda x, y: x {} y)({}".format(op[func_name], self.tree2str_translate(argtree)) 

562 else: 

563 template = '(%s(%s)).evalf(' if self.use_evalf else '%s(%s' 

564 if self.float_wrap_evalf: 

565 template = 'float(%s)' % template 

566 elif self.complex_wrap_evalf: 

567 template = 'complex(%s)' % template 

568 

569 # Wrapping should only happen on the outermost expression, which 

570 # is the only thing we know will be a number. 

571 float_wrap_evalf = self.float_wrap_evalf 

572 complex_wrap_evalf = self.complex_wrap_evalf 

573 self.float_wrap_evalf = False 

574 self.complex_wrap_evalf = False 

575 ret = template % (func_name, self.tree2str_translate(argtree)) 

576 self.float_wrap_evalf = float_wrap_evalf 

577 self.complex_wrap_evalf = complex_wrap_evalf 

578 return ret 

579 

580 ############################################################################## 

581 # The namespace constructors 

582 ############################################################################## 

583 

584 @classmethod 

585 def sympy_expression_namespace(cls, expr): 

586 """Traverses the (func, args) tree of an expression and creates a SymPy 

587 namespace. All other modules are imported only as a module name. That way 

588 the namespace is not polluted and rests quite small. It probably causes much 

589 more variable lookups and so it takes more time, but there are no tests on 

590 that for the moment.""" 

591 if expr is None: 

592 return {} 

593 else: 

594 funcname = str(expr.func) 

595 # XXX Workaround 

596 # Here we add an ugly workaround because str(func(x)) 

597 # is not always the same as str(func). Eg 

598 # >>> str(Integral(x)) 

599 # "Integral(x)" 

600 # >>> str(Integral) 

601 # "<class 'sympy.integrals.integrals.Integral'>" 

602 # >>> str(sqrt(x)) 

603 # "sqrt(x)" 

604 # >>> str(sqrt) 

605 # "<function sqrt at 0x3d92de8>" 

606 # >>> str(sin(x)) 

607 # "sin(x)" 

608 # >>> str(sin) 

609 # "sin" 

610 # Either one of those can be used but not all at the same time. 

611 # The code considers the sin example as the right one. 

612 regexlist = [ 

613 r'<class \'sympy[\w.]*?.([\w]*)\'>$', 

614 # the example Integral 

615 r'<function ([\w]*) at 0x[\w]*>$', # the example sqrt 

616 ] 

617 for r in regexlist: 

618 m = re.match(r, funcname) 

619 if m is not None: 

620 funcname = m.groups()[0] 

621 # End of the workaround 

622 # XXX debug: print funcname 

623 args_dict = {} 

624 for a in expr.args: 

625 if (isinstance(a, Symbol) or 

626 isinstance(a, NumberSymbol) or 

627 a in [I, zoo, oo]): 

628 continue 

629 else: 

630 args_dict.update(cls.sympy_expression_namespace(a)) 

631 args_dict.update({funcname: expr.func}) 

632 return args_dict 

633 

634 @staticmethod 

635 def sympy_atoms_namespace(expr): 

636 """For no real reason this function is separated from 

637 sympy_expression_namespace. It can be moved to it.""" 

638 atoms = expr.atoms(Symbol, NumberSymbol, I, zoo, oo) 

639 d = {} 

640 for a in atoms: 

641 # XXX debug: print 'atom:' + str(a) 

642 d[str(a)] = a 

643 return d