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
« 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.
3It is for internal use in the new plotting module.
4It may (will! see the Q'n'A in the source) be rewritten.
6It's completely self contained. Especially it does not use lambdarepr.
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"""
14import re
15from sympy.core.numbers import (I, NumberSymbol, oo, zoo)
16from sympy.core.symbol import Symbol
17from sympy.utilities.iterables import numbered_symbols
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.
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.
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.
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.
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.
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.
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.
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.
75from sympy.external import import_module
76import warnings
78#TODO debugging output
81class vectorized_lambdify:
82 """ Return a sufficiently smart, vectorized and lambdified function.
84 Returns only reals.
86 Explanation
87 ===========
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.
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
99 When using Python cmath there is no need for evalf or float/complex
100 because Python cmath calls those.
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')
114 self.lambda_func_1 = experimental_lambdify(
115 args, expr, use_np=True)
116 self.vector_func_1 = self.lambda_func_1
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])
123 self.vector_func = self.vector_func_1
124 self.failure = False
126 def __call__(self, *args):
127 np = self.np
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
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)
149class lambdify:
150 """Returns the lambdified function.
152 Explanation
153 ===========
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 """
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
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
186 if self.lambda_func == self.lambda_func_1:
187 self.lambda_func = self.lambda_func_2
188 return self.__call__(args)
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)
199def experimental_lambdify(*args, **kwargs):
200 l = Lambdifier(*args, **kwargs)
201 return l
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):
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
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
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))
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')})
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']
271 def __call__(self, *args, **kwargs):
272 return self.lambda_func(*args, **kwargs)
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 }
288 # Strings that should be translated
289 builtin_not_functions = {
290 'I': '1j',
291# 'oo': '1e400',
292 }
294 ###
295 # numpy
296 ###
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 ]
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 }
323 # Strings that should be translated
324 numpy_not_functions = {
325 'pi': 'np.pi',
326 'oo': 'np.inf',
327 'E': 'np.e',
328 }
330 ###
331 # Python math
332 ###
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 ]
341 # Functions with different names in math
342 math_functions_different = {
343 'ceiling': 'ceil',
344 'ln': 'log',
345 'loggamma': 'lgamma'
346 }
348 # Strings that should be translated
349 math_not_functions = {
350 'pi': 'math.pi',
351 'E': 'math.e',
352 }
354 ###
355 # Python cmath
356 ###
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 ]
365 # Functions with different names in cmath
366 cmath_functions_different = {
367 'ln': 'log',
368 'arg': 'phase',
369 }
371 # Strings that should be translated
372 cmath_not_functions = {
373 'pi': 'cmath.pi',
374 'E': 'cmath.e',
375 }
377 ###
378 # intervalmath
379 ###
381 interval_not_functions = {
382 'pi': 'math.pi',
383 'E': 'math.e'
384 }
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 ]
393 interval_functions_different = {
394 'Min': 'imin',
395 'Max': 'imax',
396 'ceiling': 'ceil',
398 }
400 ###
401 # mpmath, etc
402 ###
403 #TODO
405 ###
406 # Create the final ordered tuples of dictionaries
407 ###
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
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
447 ##############################################################################
448 # The translator functions, tree parsers, etc.
449 ##############################################################################
451 def str2tree(self, exprstr):
452 """Converts an expression string to a tree.
454 Explanation
455 ===========
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.
461 Examples
462 ========
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
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)
500 @classmethod
501 def tree2str(cls, tree):
502 """Converts a tree to string without translations.
504 Examples
505 ========
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
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))
521 def tree2str_translate(self, tree):
522 """Converts a tree to string with translations.
524 Explanation
525 ===========
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])
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
544 def translate_func(self, func_name, argtree):
545 """Translate function names and the tree of arguments.
547 Explanation
548 ===========
550 If the function name is not in the dictionaries of dict_tuple_fun then the
551 function is surrounded by a float((...).evalf()).
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
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
580 ##############################################################################
581 # The namespace constructors
582 ##############################################################################
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
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