Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

pretty printing gets stuck doing evalf #10800

Open
asmeurer opened this issue Mar 11, 2016 · 19 comments · May be fixed by #10814
Open

pretty printing gets stuck doing evalf #10800

asmeurer opened this issue Mar 11, 2016 · 19 comments · May be fixed by #10814

Comments

@asmeurer
Copy link
Member

In [147]: a = 2*integrate(sin(exp(x)), (x, 1, oo))+2*Si(E)

In [148]: a
^COut[148]: ---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-148-60b725f10c9c> in <module>()
----> 1 a

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/IPython/core/displayhook.py in __call__(self, result)
    244             self.start_displayhook()
    245             self.write_output_prompt()
--> 246             format_dict, md_dict = self.compute_format_data(result)
    247             self.update_user_ns(result)
    248             self.fill_exec_result(result)

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/IPython/core/displayhook.py in compute_format_data(self, result)
    150
    151         """
--> 152         return self.shell.display_formatter.format(result)
    153
    154     def write_format_data(self, format_dict, md_dict=None):

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/IPython/core/formatters.py in format(self, obj, include, exclude)
    175             md = None
    176             try:
--> 177                 data = formatter(obj)
    178             except:
    179                 # FIXME: log the exception

<decorator-gen-10> in __call__(self, obj)

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/IPython/core/formatters.py in catch_format_error(method, self, *args, **kwargs)
    220     """show traceback on failed format call"""
    221     try:
--> 222         r = method(self, *args, **kwargs)
    223     except NotImplementedError:
    224         # don't warn on NotImplementedErrors

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/IPython/core/formatters.py in __call__(self, obj)
    697                 type_pprinters=self.type_printers,
    698                 deferred_pprinters=self.deferred_printers)
--> 699             printer.pretty(obj)
    700             printer.flush()
    701             return stream.getvalue()

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/IPython/lib/pretty.py in pretty(self, obj)
    366                 if cls in self.type_pprinters:
    367                     # printer registered in self.type_pprinters
--> 368                     return self.type_pprinters[cls](obj, self, cycle)
    369                 else:
    370                     # deferred printer

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/interactive/printing.py in _print_plain(arg, p, cycle)
     66         """caller for pretty, for use in IPython 0.11"""
     67         if _can_print_latex(arg):
---> 68             p.text(stringify_func(arg))
     69         else:
     70             p.text(IPython.lib.pretty.pretty(arg))

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/printing/pretty/pretty.py in pretty(expr, **settings)
   2037
   2038     try:
-> 2039         return pp.doprint(expr)
   2040     finally:
   2041         pretty_use_unicode(uflag)

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/printing/pretty/pretty.py in doprint(self, expr)
     56
     57     def doprint(self, expr):
---> 58         return self._print(expr).render(**self._settings)
     59
     60     # empty op so _print(stringPict) returns the same

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/printing/printer.py in _print(self, expr, *args, **kwargs)
    255                 printmethod = '_print_' + cls.__name__
    256                 if hasattr(self, printmethod):
--> 257                     return getattr(self, printmethod)(expr, *args, **kwargs)
    258
    259             # Unknown object, fall back to the emptyPrinter.

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/printing/pretty/pretty.py in _print_Add(self, expr, order)
   1234             terms = list(expr.args)
   1235         else:
-> 1236             terms = self._as_ordered_terms(expr, order=order)
   1237         pforms, indices = [], []
   1238

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/printing/printer.py in _as_ordered_terms(self, expr, order)
    269             return sorted(Add.make_args(expr), key=cmp_to_key(Basic._compare_pretty))
    270         else:
--> 271             return expr.as_ordered_terms(order=order)

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/core/expr.py in as_ordered_terms(self, order, data)
    886         """
    887         key, reverse = self._parse_order(order)
--> 888         terms, gens = self.as_terms()
    889
    890         if not any(term.is_Order for term, _ in terms):

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/core/expr.py in as_terms(self)
    925                     if factor.is_number:
    926                         try:
--> 927                             coeff *= complex(factor)
    928                         except TypeError:
    929                             pass

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/core/expr.py in __complex__(self)
    227
    228     def __complex__(self):
--> 229         result = self.evalf()
    230         re, im = result.as_real_imag()
    231         return complex(float(re), float(im))

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/core/evalf.py in evalf(self, n, subs, maxn, chop, strict, quad, verbose)
   1384             options['quad'] = quad
   1385         try:
-> 1386             result = evalf(self, prec + 4, options)
   1387         except NotImplementedError:
   1388             # Fall back to the ordinary evalf

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/core/evalf.py in evalf(x, prec, options)
   1276     try:
   1277         rf = evalf_table[x.func]
-> 1278         r = rf(x, prec, options)
   1279     except KeyError:
   1280         try:

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/core/evalf.py in evalf_integral(expr, prec, options)
   1006     maxprec = options.get('maxprec', INF)
   1007     while 1:
-> 1008         result = do_integral(expr, workprec, options)
   1009         accuracy = complex_accuracy(result)
   1010         if accuracy >= prec:  # achieved desired precision

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/core/evalf.py in do_integral(expr, prec, options)
    965             quadrature_error = MINUS_INF
    966         else:
--> 967             result, quadrature_error = quadts(f, [xlow, xhigh], error=1)
    968             quadrature_error = fastlog(quadrature_error._mpf_)
    969

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/calculus/quadrature.py in quadts(ctx, *args, **kwargs)
    784         """
    785         kwargs['method'] = 'tanh-sinh'
--> 786         return ctx.quad(*args, **kwargs)
    787
    788     def quadgl(ctx, *args, **kwargs):

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/calculus/quadrature.py in quad(ctx, f, *points, **kwargs)
    741             ctx.prec += 20
    742             if dim == 1:
--> 743                 v, err = rule.summation(f, points[0], prec, epsilon, m, verbose)
    744             elif dim == 2:
    745                 v, err = rule.summation(lambda x: \

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/calculus/quadrature.py in summation(self, f, points, prec, epsilon, max_degree, verbose)
    230                     print("Integrating from %s to %s (degree %s of %s)" % \
    231                         (ctx.nstr(a), ctx.nstr(b), degree, max_degree))
--> 232                 results.append(self.sum_next(f, nodes, degree, prec, results, verbose))
    233                 if degree > 1:
    234                     err = self.estimate_error(results, prec, epsilon)

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/calculus/quadrature.py in sum_next(self, f, nodes, degree, prec, previous, verbose)
    302         else:
    303             S = self.ctx.zero
--> 304         S += self.ctx.fdot((w,f(x)) for (x,w) in nodes)
    305         return h*S
    306

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/ctx_mp_python.py in fdot(ctx, A, B, conjugate)
    923         hasattr_ = hasattr
    924         types = (ctx.mpf, ctx.mpc)
--> 925         for a, b in A:
    926             if type(a) not in types: a = ctx.convert(a)
    927             if type(b) not in types: b = ctx.convert(b)

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/calculus/quadrature.py in <genexpr>(.0)
    302         else:
    303             S = self.ctx.zero
--> 304         S += self.ctx.fdot((w,f(x)) for (x,w) in nodes)
    305         return h*S
    306

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/core/evalf.py in f(t)
    938
    939         def f(t):
--> 940             re, im, re_acc, im_acc = evalf(func, mp.prec, {'subs': {x: t}})
    941
    942             have_part[0] = re or have_part[0]

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/core/evalf.py in evalf(x, prec, options)
   1276     try:
   1277         rf = evalf_table[x.func]
-> 1278         r = rf(x, prec, options)
   1279     except KeyError:
   1280         try:

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/core/evalf.py in evalf_trig(v, prec, options)
    775     if xsize >= 10:
    776         xprec = prec + xsize
--> 777         re, im, re_acc, im_acc = evalf(arg, xprec, options)
    778     # Need to repeat in case the argument is very close to a
    779     # multiple of pi (or pi/2), hitting close to a root

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/core/evalf.py in evalf(x, prec, options)
   1276     try:
   1277         rf = evalf_table[x.func]
-> 1278         r = rf(x, prec, options)
   1279     except KeyError:
   1280         try:

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/core/evalf.py in <lambda>(x, prec, options)
   1245
   1246         exp: lambda x, prec, options: evalf_pow(
-> 1247             Pow(S.Exp1, x.args[0], evaluate=False), prec, options),
   1248
   1249         cos: evalf_trig,

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/core/evalf.py in evalf_pow(v, prec, options)
    704             re, im = libmp.mpc_exp((yre or fzero, yim), prec)
    705             return finalize_complex(re, im, target_prec)
--> 706         return mpf_exp(yre, target_prec), None, target_prec, None
    707
    708     xre, xim, _, _ = evalf(base, prec + 5, options)

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/libmp/libelefun.py in mpf_exp(x, prec, rnd)
   1174             else:
   1175                 t = man >> (-offset)
-> 1176             lg2 = ln2_fixed(wpmod)
   1177             n, t = divmod(t, lg2)
   1178             n = int(n)

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/libmp/libelefun.py in g(prec, **kwargs)
     97             return f.memo_val >> (memo_prec-prec)
     98         newprec = int(prec*1.05+10)
---> 99         f.memo_val = f(newprec, **kwargs)
    100         f.memo_prec = newprec
    101         return f.memo_val >> (newprec-prec)

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/libmp/libelefun.py in ln2_fixed(prec)
    166     with binary splitting at high precision.
    167     """
--> 168     return machin([(18, 26), (-2, 4801), (8, 8749)], prec, True)
    169
    170 @constant_memo

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/libmp/libelefun.py in machin(coefs, prec, hyperbolic)
    154     s = MPZ_ZERO
    155     for a, b in coefs:
--> 156         s += MPZ(a) * acot_fixed(MPZ(b), prec+extraprec, hyperbolic)
    157     return (s >> extraprec)
    158

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/libmp/libelefun.py in acot_fixed(a, prec, hyperbolic)
    141     """
    142     N = int(0.35 * prec/math.log(a) + 20)
--> 143     p, q, r = bsp_acot(a, 0,N, hyperbolic)
    144     return ((p+q)<<prec)//(q*a)
    145

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/libmp/libelefun.py in bsp_acot(q, a, b, hyperbolic)
    129             return -MPZ_ONE, a1 * q**2, a1
    130     m = (a+b)//2
--> 131     p1, q1, r1 = bsp_acot(q, a, m, hyperbolic)
    132     p2, q2, r2 = bsp_acot(q, m, b, hyperbolic)
    133     return q2*p1 + r1*p2, q1*q2, r1*r2

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/libmp/libelefun.py in bsp_acot(q, a, b, hyperbolic)
    129             return -MPZ_ONE, a1 * q**2, a1
    130     m = (a+b)//2
--> 131     p1, q1, r1 = bsp_acot(q, a, m, hyperbolic)
    132     p2, q2, r2 = bsp_acot(q, m, b, hyperbolic)
    133     return q2*p1 + r1*p2, q1*q2, r1*r2

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/libmp/libelefun.py in bsp_acot(q, a, b, hyperbolic)
    129             return -MPZ_ONE, a1 * q**2, a1
    130     m = (a+b)//2
--> 131     p1, q1, r1 = bsp_acot(q, a, m, hyperbolic)
    132     p2, q2, r2 = bsp_acot(q, m, b, hyperbolic)
    133     return q2*p1 + r1*p2, q1*q2, r1*r2

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/libmp/libelefun.py in bsp_acot(q, a, b, hyperbolic)
    129             return -MPZ_ONE, a1 * q**2, a1
    130     m = (a+b)//2
--> 131     p1, q1, r1 = bsp_acot(q, a, m, hyperbolic)
    132     p2, q2, r2 = bsp_acot(q, m, b, hyperbolic)
    133     return q2*p1 + r1*p2, q1*q2, r1*r2

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/libmp/libelefun.py in bsp_acot(q, a, b, hyperbolic)
    129             return -MPZ_ONE, a1 * q**2, a1
    130     m = (a+b)//2
--> 131     p1, q1, r1 = bsp_acot(q, a, m, hyperbolic)
    132     p2, q2, r2 = bsp_acot(q, m, b, hyperbolic)
    133     return q2*p1 + r1*p2, q1*q2, r1*r2

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/libmp/libelefun.py in bsp_acot(q, a, b, hyperbolic)
    129             return -MPZ_ONE, a1 * q**2, a1
    130     m = (a+b)//2
--> 131     p1, q1, r1 = bsp_acot(q, a, m, hyperbolic)
    132     p2, q2, r2 = bsp_acot(q, m, b, hyperbolic)
    133     return q2*p1 + r1*p2, q1*q2, r1*r2

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/libmp/libelefun.py in bsp_acot(q, a, b, hyperbolic)
    129             return -MPZ_ONE, a1 * q**2, a1
    130     m = (a+b)//2
--> 131     p1, q1, r1 = bsp_acot(q, a, m, hyperbolic)
    132     p2, q2, r2 = bsp_acot(q, m, b, hyperbolic)
    133     return q2*p1 + r1*p2, q1*q2, r1*r2

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/libmp/libelefun.py in bsp_acot(q, a, b, hyperbolic)
    130     m = (a+b)//2
    131     p1, q1, r1 = bsp_acot(q, a, m, hyperbolic)
--> 132     p2, q2, r2 = bsp_acot(q, m, b, hyperbolic)
    133     return q2*p1 + r1*p2, q1*q2, r1*r2
    134

/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/libmp/libelefun.py in bsp_acot(q, a, b, hyperbolic)
    131     p1, q1, r1 = bsp_acot(q, a, m, hyperbolic)
    132     p2, q2, r2 = bsp_acot(q, m, b, hyperbolic)
--> 133     return q2*p1 + r1*p2, q1*q2, r1*r2
    134
    135 # the acoth(x) series converges like the geometric series for x^2

KeyboardInterrupt:

The printer shouldn't be calling evalf at all for any expression. We really need to fix this.

@ThomasHickman
Copy link
Contributor

Intergral(sin(exp(x)), (x, 1, oo)) is being evaluated when considering the order of the terms printed in 2*integrate(sin(exp(x)), (x, 1, oo))+2*Si(E), as the values of coefficients are taken into account when ordering terms. This would be fixed if only complex terms were compared and therefore integrate(sin(exp(x)), (x, 1, oo)) would not be evaluated.

@smichr
Copy link
Member

smichr commented Mar 20, 2016

Why don't we just use Add(*list(ordered(Add.make_args(constant_expr)))) where constant_expr is the part of the expression containing constants?

@asmeurer
Copy link
Member Author

Does ordered also use evalf? I've noticed that it can also be slow.

@smichr
Copy link
Member

smichr commented Apr 23, 2016

Does ordered also use evalf? I've noticed that it can also be slow

By default it only uses _nodes and default_sort_key. I would be interested to know of any case where it is slow.

@smichr
Copy link
Member

smichr commented Apr 23, 2016

Here's another simple printing failure:

>>> Sum(x,(x,1,oo))+3
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "sympy\core\basic.py", line 393, in __repr__
    return sstr(self, order=None)
  File "sympy\printing\str.py", line 753, in sstr
    s = p.doprint(expr)
  File "sympy\printing\printer.py", line 233, in doprint
    return self._str(self._print(expr))
  File "sympy\printing\printer.py", line 257, in _print
    return getattr(self, printmethod)(expr, *args, **kwargs)
  File "sympy\printing\str.py", line 51, in _print_Add
    terms = self._as_ordered_terms(expr, order=order)
  File "sympy\printing\printer.py", line 271, in _as_ordered_terms
    return expr.as_ordered_terms(order=order)
  File "sympy\core\expr.py", line 890, in as_ordered_terms
    terms, gens = self.as_terms()
  File "sympy\core\expr.py", line 929, in as_terms
    coeff *= complex(factor)
  File "sympy\core\expr.py", line 231, in __complex__
    result = self.evalf()
  File "sympy\core\evalf.py", line 1386, in evalf
    result = evalf(self, prec + 4, options)
  File "sympy\core\evalf.py", line 1278, in evalf
    r = rf(x, prec, options)
  File "sympy\core\evalf.py", line 1165, in evalf_sum
    v = hypsum(func, n, int(a), prec2)
  File "sympy\core\evalf.py", line 1111, in hypsum
    raise ValueError("Sum diverges like n^%i" % (-p))
ValueError: Sum diverges like n^1

@asmeurer
Copy link
Member Author

asmeurer commented Apr 23, 2016

It does get called. An easy way to check is to put a stack print in _eval_evalf.

In [4]: class Test(Function):
   ...:     def _eval_evalf(self, prec):
   ...:         import traceback
   ...:         traceback.print_stack()
   ...:         return S(1)
   ...:

In [7]: Test(1)
Out[7]: Test(1)

In [8]: Test(1) + 1
  File "./bin/isympy", line 357, in <module>
    main()
  File "./bin/isympy", line 354, in main
    init_session(ipython, **args)
  File "/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/interactive/session.py", line 475, in init_session
    mainloop(message)
  File "/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/IPython/terminal/interactiveshell.py", line 548, in mainloop
    self.interact(display_banner=display_banner)
  File "/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/IPython/terminal/interactiveshell.py", line 672, in interact
    self.run_cell(source_raw, store_history=True)
  File "/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/IPython/core/interactiveshell.py", line 2723, in run_cell
    interactivity=interactivity, compiler=compiler, result=result)
  File "/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/IPython/core/interactiveshell.py", line 2831, in run_ast_nodes
    if self.run_code(code, result):
  File "/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/IPython/core/interactiveshell.py", line 2885, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-8-6471e92dcaf1>", line 1, in <module>
    Test(1) + 1
  File "/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/IPython/core/displayhook.py", line 246, in __call__
    format_dict, md_dict = self.compute_format_data(result)
  File "/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/IPython/core/displayhook.py", line 152, in compute_format_data
    return self.shell.display_formatter.format(result)
  File "/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/IPython/core/formatters.py", line 177, in format
    data = formatter(obj)
  File "<decorator-gen-10>", line 2, in __call__
  File "/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/IPython/core/formatters.py", line 222, in catch_format_error
    r = method(self, *args, **kwargs)
  File "/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/IPython/core/formatters.py", line 699, in __call__
    printer.pretty(obj)
  File "/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/IPython/lib/pretty.py", line 368, in pretty
    return self.type_pprinters[cls](obj, self, cycle)
  File "/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/interactive/printing.py", line 68, in _print_plain
    p.text(stringify_func(arg))
  File "/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/printing/pretty/pretty.py", line 2049, in pretty
    return pp.doprint(expr)
  File "/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/printing/pretty/pretty.py", line 59, in doprint
    return self._print(expr).render(**self._settings)
  File "/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/printing/printer.py", line 257, in _print
    return getattr(self, printmethod)(expr, *args, **kwargs)
  File "/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/printing/pretty/pretty.py", line 1247, in _print_Add
    terms = self._as_ordered_terms(expr, order=order)
  File "/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/printing/printer.py", line 271, in _as_ordered_terms
    return expr.as_ordered_terms(order=order)
  File "/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/core/expr.py", line 890, in as_ordered_terms
    terms, gens = self.as_terms()
  File "/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/core/expr.py", line 929, in as_terms
    coeff *= complex(factor)
  File "/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/core/expr.py", line 231, in __complex__
    result = self.evalf()
  File "/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/core/evalf.py", line 1386, in evalf
    result = evalf(self, prec + 4, options)
  File "/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/core/evalf.py", line 1284, in evalf
    xe = x._eval_evalf(prec)
  File "<ipython-input-4-a3c4975d39eb>", line 4, in _eval_evalf
    traceback.print_stack()
Out[8]: 1 + Test(1)

So it looks like as_terms calls complex(), which uses evalf (this is the same instance in all three tracebacks here). I think for testing this, we should create a dummy Function subclass whose _eval_evalf raises an exception and print some expressions containing it.

Also worth noting: _should_evalf is ignored (see https://stackoverflow.com/questions/35752150/how-to-avoid-automatic-evalf-on-new-sympy-functions).

@oscarbenjamin
Copy link
Contributor

The original examples in this thread now works fine:

In [28]: a = 2*integrate(sin(exp(x)), (x, 1, oo))+2*Si(E)                                                                                      

In [29]: a                                                                                                                                     
Out[29]: π

I guess that's because it now returns a different result. Whatever result it previously returned may still not be printable...

@asmeurer
Copy link
Member Author

Replace the integral with one that can't be computed. For example, a = 2*integrate(sin(exp(x**2)), (x, 1, oo))+2*Si(E). Or, more simply, replace integrate with Integral in the original example.

@asmeurer
Copy link
Member Author

asmeurer commented Mar 20, 2019

@smichr's example works fine now though:

>>> pprint(Sum(x, (x, 1, oo)) + 3)
  ∞
 ___
 ╲
  ╲   x
  ╱     + 3
 ╱
 ‾‾‾
x = 1

Although there's an obvious centering issue there #16346.

@czgdp1807
Copy link
Member

@asmeurer
It looks like the issue has been fixed in master as of this comment,

Python 3.6.9 (default, Nov  7 2019, 10:44:02) 
[GCC 8.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
(SymPyConsole)
>>> a = 2*integrate(sin(exp(x)), (x, 1, oo))+2*Si(E)
>>> a
π

Attaching a Could close label.

@oscarbenjamin
Copy link
Contributor

That's just because the integral now evaluates. If you try to print the unevaluated integral you will see the error:

a = 2*Integral(sin(exp(x)), (x, 1, oo))+2*Si(E)

@oscarbenjamin
Copy link
Contributor

The problem is that the printing code tries to sort some expressions (not sure why). The assumptions system attempts to numerically evaluate the expressions as part of a comparison like expra < exprb. Then apparently evalf fails for whatever expression is evaluated. So there are 3 possible things to fix:

  1. The printing code should not sort or should sort in some way that doesn't involve evaluating mathematical inequalities.
  2. The assumptions code for Expr should not attempt numerical evaluation when comparing objects that have is_number=True (this is heavily depended on though).
  3. There is some object for which evalf hangs and that should be fixed.

@asmeurer
Copy link
Member Author

Certainly evalf should be improved, but we shouldn't assume that it will be fast for every object.

Here is the issue

>>> pi + E # Sorted from smallest to largest
E + pi
>>> (pi - E).is_positive
True

So if we remove evalf from the assumptions, it would remove a lot of functionality and would probably break things. So I'm not sure how to fix it. Clearly it shouldn't try to evaluate very slow to evaluate things, but I don't know if it should try to remove evaluation entirely. Maybe for certain objects like Sum and Integral the evaluation shouldn't be enabled by default for such things.

I don't know if the printing sorting is that important, but the assumptions calling evalf could cause issues in other contexts as well.

I don't know how much performance is affected in general just from "normal" functions evaluating (throughout SymPy, not just in the printers). Maybe we should try removing the evalf calls from the assumptions and seeing 1) what fails, and 2) if it affects performance benchmarks outside of the printers.

@oscarbenjamin
Copy link
Contributor

See also #17609

@asmeurer
Copy link
Member Author

Oh that's a good point. Even "normal" expressions can cause issues with evalf. I think we should have a "fast" evalf that gives up on anything that could potentially be slow (like integration) or produce excessively large expressions. evalf() itself should just do the evaluation, but anything used in library code to get a numerical estimate should use the safer mode.

The question is to what degree can we disable only the "bad" things, without disabling entire classes of expressions wholesale. It would be nice if mpmath itself had some sort of mode to do a fail if the internal values get too large. I think the problem is that mpmath represents things as man * 2^exp (http://mpmath.org/doc/current/technical.html#representation-of-numbers), where man is limited by the precision, but unlike machine floats, exp can be any Python long integer. According to that doc page mpmath doesn't have such a mode.

It would also help if the slower mpmath functions like quad had a timeout ability. Maybe those flags do exist (it looks like we can use maxdegree http://mpmath.org/doc/current/calculus/integration.html#mpmath.quad). We have to investigate them.

@oscarbenjamin
Copy link
Contributor

The question is to what degree can we disable only the "bad" things, without disabling entire classes of expressions wholesale. It would be nice if mpmath itself had some sort of mode to do a fail if the internal values get too large.

In evalf there is already a limit on the mantissa so we can have one for the exponent as well. As I said in #17609 the obvious place to check this would be evalf_pow but there are probably others. Most floating point systems limit the exponent e.g.:

In [2]: import decimal                                                                                                                         

In [3]: decimal.getcontext()                                                                                                                   
Out[3]: Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[], traps=[InvalidOperation, DivisionByZero, Overflow])

Emax there is the maximum value that any exponent can have in the default decimal context. You can change the context but you are still limited by
https://docs.python.org/3/library/decimal.html#decimal.MAX_EMAX
which gives a maximum exponent of 999999999999999999. That's plenty big enough for normal calculations. I doubt that much "real" usage of sympy needs more than that.

@czgdp1807
Copy link
Member

Current status,

>>> a = 2*integrate(sin(exp(x)), (x, 1, oo))+2*Si(E)
>>> a
π
>>> a = 2*Integral(sin(exp(x)), (x, 1, oo))+2*Si(E)
>>> a
^CTraceback (most recent call last):
  File "/usr/lib/python3.6/code.py", line 91, in runcode
    exec(code, self.locals)
  File "<console>", line 1, in <module>
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/interactive/printing.py", line 29, in _displayhook
    print(stringify_func(arg, **settings))
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/printing/pretty/pretty.py", line 2647, in pretty
    return pp.doprint(expr)
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/printing/pretty/pretty.py", line 65, in doprint
    return self._print(expr).render(**self._settings)
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/printing/printer.py", line 289, in _print
    return getattr(self, printmethod)(expr, **kwargs)
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/printing/pretty/pretty.py", line 1707, in _print_Add
    terms = self._as_ordered_terms(expr, order=order)
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/printing/printer.py", line 308, in _as_ordered_terms
    return expr.as_ordered_terms(order=order)
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/core/expr.py", line 1149, in as_ordered_terms
    terms, gens = self.as_terms()
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/core/expr.py", line 1188, in as_terms
    coeff *= complex(factor)
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/core/expr.py", line 332, in __complex__
    result = self.evalf()
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/core/evalf.py", line 1462, in evalf
    result = evalf(self, prec + 4, options)
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/core/evalf.py", line 1314, in evalf
    r = rf(x, prec, options)
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/core/evalf.py", line 1044, in evalf_integral
    result = do_integral(expr, workprec, options)
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/core/evalf.py", line 1003, in do_integral
    result, quadrature_error = quadts(f, [xlow, xhigh], error=1)
  File "/usr/local/lib/python3.6/dist-packages/mpmath-1.1.0-py3.6.egg/mpmath/calculus/quadrature.py", line 785, in quadts
    return ctx.quad(*args, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/mpmath-1.1.0-py3.6.egg/mpmath/calculus/quadrature.py", line 742, in quad
    v, err = rule.summation(f, points[0], prec, epsilon, m, verbose)
  File "/usr/local/lib/python3.6/dist-packages/mpmath-1.1.0-py3.6.egg/mpmath/calculus/quadrature.py", line 232, in summation
    results.append(self.sum_next(f, nodes, degree, prec, results, verbose))
  File "/usr/local/lib/python3.6/dist-packages/mpmath-1.1.0-py3.6.egg/mpmath/calculus/quadrature.py", line 304, in sum_next
    S += self.ctx.fdot((w,f(x)) for (x,w) in nodes)
  File "/usr/local/lib/python3.6/dist-packages/mpmath-1.1.0-py3.6.egg/mpmath/ctx_mp_python.py", line 944, in fdot
    for a, b in A:
  File "/usr/local/lib/python3.6/dist-packages/mpmath-1.1.0-py3.6.egg/mpmath/calculus/quadrature.py", line 304, in <genexpr>
    S += self.ctx.fdot((w,f(x)) for (x,w) in nodes)
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/core/evalf.py", line 976, in f
    re, im, re_acc, im_acc = evalf(func, mp.prec, {'subs': {x: t}})
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/core/evalf.py", line 1314, in evalf
    r = rf(x, prec, options)
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/core/evalf.py", line 813, in evalf_trig
    re, im, re_acc, im_acc = evalf(arg, xprec, options)
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/core/evalf.py", line 1314, in evalf
    r = rf(x, prec, options)
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/core/evalf.py", line 1283, in <lambda>
    Pow(S.Exp1, x.args[0], evaluate=False), prec, options),
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/core/evalf.py", line 742, in evalf_pow
    return mpf_exp(yre, target_prec), None, target_prec, None
  File "/usr/local/lib/python3.6/dist-packages/mpmath-1.1.0-py3.6.egg/mpmath/libmp/libelefun.py", line 1176, in mpf_exp
    lg2 = ln2_fixed(wpmod)
  File "/usr/local/lib/python3.6/dist-packages/mpmath-1.1.0-py3.6.egg/mpmath/libmp/libelefun.py", line 99, in g
    f.memo_val = f(newprec, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/mpmath-1.1.0-py3.6.egg/mpmath/libmp/libelefun.py", line 168, in ln2_fixed
    return machin([(18, 26), (-2, 4801), (8, 8749)], prec, True)
  File "/usr/local/lib/python3.6/dist-packages/mpmath-1.1.0-py3.6.egg/mpmath/libmp/libelefun.py", line 156, in machin
    s += MPZ(a) * acot_fixed(MPZ(b), prec+extraprec, hyperbolic)
  File "/usr/local/lib/python3.6/dist-packages/mpmath-1.1.0-py3.6.egg/mpmath/libmp/libelefun.py", line 143, in acot_fixed
    p, q, r = bsp_acot(a, 0,N, hyperbolic)
  File "/usr/local/lib/python3.6/dist-packages/mpmath-1.1.0-py3.6.egg/mpmath/libmp/libelefun.py", line 131, in bsp_acot
    p1, q1, r1 = bsp_acot(q, a, m, hyperbolic)
  File "/usr/local/lib/python3.6/dist-packages/mpmath-1.1.0-py3.6.egg/mpmath/libmp/libelefun.py", line 131, in bsp_acot
    p1, q1, r1 = bsp_acot(q, a, m, hyperbolic)
  File "/usr/local/lib/python3.6/dist-packages/mpmath-1.1.0-py3.6.egg/mpmath/libmp/libelefun.py", line 131, in bsp_acot
    p1, q1, r1 = bsp_acot(q, a, m, hyperbolic)
  [Previous line repeated 2 more times]
  File "/usr/local/lib/python3.6/dist-packages/mpmath-1.1.0-py3.6.egg/mpmath/libmp/libelefun.py", line 132, in bsp_acot
    p2, q2, r2 = bsp_acot(q, m, b, hyperbolic)
  File "/usr/local/lib/python3.6/dist-packages/mpmath-1.1.0-py3.6.egg/mpmath/libmp/libelefun.py", line 133, in bsp_acot
    return q2*p1 + r1*p2, q1*q2, r1*r2
KeyboardInterrupt

This issue is still open for work.

@smichr
Copy link
Member

smichr commented Sep 22, 2023

  1. printing code should not sort or should sort in some way that doesn't involve evaluating mathematical inequalities

Use ordered to put things in canonical order.

@asmeurer
Copy link
Member Author

asmeurer commented Sep 22, 2023

It shouldn't sort at all. The args of an Add are already sorted, and even if they weren't, it doesn't matter if something prints as x + y or y + x. The only ordering that really matters is printing the monomials of a polynomial in order, so that x**2 + x + 1 doesn't print as x + 1 + x**2. But there's some idea in the printer now that numerical constants need to be sorted numerically:

>>> pi + E + 3
E + 3 + pi

which is something that I don't think anyone actually expects, or even realizes is there. But it incurs a huge cost in the printing, which should be an extremely cheap operation.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
5 participants