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

N + Subs + Zeta + Derivative = maximum recursion depth exceeded #11802

Open
latot opened this issue Nov 2, 2016 · 57 comments
Open

N + Subs + Zeta + Derivative = maximum recursion depth exceeded #11802

latot opened this issue Nov 2, 2016 · 57 comments

Comments

@latot
Copy link
Contributor

latot commented Nov 2, 2016

Hi, here the problem:

>>> x=symbols('x')
>>> p=zeta(x)
>>> p=p.diff(x, 3)
>>> p=p.subs(x, 4)
>>> N(p)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib64/python3.4/site-packages/sympy/core/evalf.py", line 1481, in N
    return sympify(x).evalf(n, **options)
#The next 2 lines x490 times
  File "/usr/lib64/python3.4/site-packages/sympy/core/function.py", line 561, in _eval_derivative
    l.append(df * da)
  File "/usr/lib64/python3.4/site-packages/sympy/core/decorators.py", line 77, in __sympifyit_wrapper
    return func(a, b)
  File "/usr/lib64/python3.4/site-packages/sympy/core/decorators.py", line 118, in binary_op_wrapper
    return func(self, other)
  File "/usr/lib64/python3.4/site-packages/sympy/core/expr.py", line 141, in __mul__
    return Mul(self, other)
  File "/usr/lib64/python3.4/site-packages/sympy/core/cache.py", line 93, in wrapper
    retval = cfunc(*args, **kwargs)
  File "/usr/lib64/python3.4/functools.py", line 458, in wrapper
    key = make_key(args, kwds, typed)
  File "/usr/lib64/python3.4/functools.py", line 382, in _make_key
    return _HashedSeq(key)
  File "/usr/lib64/python3.4/functools.py", line 351, in __init__
    self.hashvalue = hash(tup)
  File "/usr/lib64/python3.4/site-packages/sympy/core/basic.py", line 110, in __hash__
    h = hash((type(self).__name__,) + self._hashable_content())
RuntimeError: maximum recursion depth exceeded

Thx. Cya.

@cbm755
Copy link
Contributor

cbm755 commented Apr 29, 2019

For some reason that I haven't tracked down, this can actually crash Python on Fedora 30, Intel hardware, SymPy 1.4

Python 3.7.3 (default, Mar 27 2019, 13:36:35) 
[GCC 9.0.1 20190227 (Red Hat 9.0.1-0.8)] on linux

>>> A = zeta(x).diff(x,2).subs(x, 2)
⎛  2      ⎞│   
⎜ d       ⎟│   
⎜───(ζ(x))⎟│   
⎜  2      ⎟│   
⎝dx       ⎠│x=2

>>> N(A)
Fatal Python error: Cannot recover from stack overflow.

Current thread 0x00007f5c7a9c3680 (most recent call first):
  File "/usr/lib/python3.7/site-packages/sympy/core/sympify.py", line 293 in sympify
  File "/usr/lib/python3.7/site-packages/sympy/core/containers.py", line 51 in <genexpr>
  File "/usr/lib/python3.7/site-packages/sympy/core/containers.py", line 52 in __new__
  File "/usr/lib/python3.7/site-packages/sympy/core/containers.py", line 146 in <lambda>
  File "/usr/lib/python3.7/site-packages/sympy/core/sympify.py", line 293 in sympify
  File "/usr/lib/python3.7/site-packages/sympy/core/containers.py", line 51 in <genexpr>
  File "/usr/lib/python3.7/site-packages/sympy/core/containers.py", line 52 in __new__
  File "/usr/lib/python3.7/site-packages/sympy/core/containers.py", line 146 in <lambda>
  File "/usr/lib/python3.7/site-packages/sympy/core/sympify.py", line 293 in sympify
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 1218 in __new__
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 1556 in _eval_derivative
  File "/usr/lib/python3.7/site-packages/sympy/core/basic.py", line 1663 in _visit_eval_derivative_scalar
  File "/usr/lib/python3.7/site-packages/sympy/core/basic.py", line 1658 in _accept_eval_derivative
  File "/usr/lib/python3.7/site-packages/sympy/core/basic.py", line 1682 in _eval_derivative_n_times
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 1392 in __new__
  File "/usr/lib/python3.7/site-packages/sympy/core/expr.py", line 3095 in diff
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2115 in doit
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  File "/usr/lib/python3.7/site-packages/sympy/core/function.py", line 2136 in evalf
  ...
Aborted (core dumped)

That is, not just recursion in Python but a proper stack overflow... Ouch.

The first derivative still causes a "maximum recusion depth exceeded" exception as noted above.

@latot
Copy link
Contributor Author

latot commented Apr 30, 2019

maybe now pthon3.7 have less limit than the previos versions, maybe test sys.setrecursionlimit(99999)?

@oscarbenjamin
Copy link
Collaborator

A here is an object that doits to itself:

In [1]: A = zeta(x).diff(x,2).subs(x, 2)                                                                                          

In [2]: A                                                                                                                         
Out[2]:2      ⎞│   
⎜ d       ⎟│   
⎜───(ζ(x))⎟│   
⎜  2      ⎟│   
⎝dx       ⎠│x=2

In [3]: A.doit()                                                                                                                  
Out[3]:2      ⎞│   
⎜ d       ⎟│   
⎜───(ζ(x))⎟│   
⎜  2      ⎟│   
⎝dx       ⎠│x=2

So this will give infinite recursion:

sympy/sympy/core/function.py

Lines 2135 to 2136 in c162866

def evalf(self, prec=None, **options):
return self.doit().evalf(prec, **options)

@latot
Copy link
Contributor Author

latot commented Apr 30, 2019

mm, that sounds more like a sympy bug.

@sachin-4099
Copy link
Contributor

I would like to take up this issue. Please guide me as to what needs to be done.
@oscarbenjamin @latot @cbm755

@latot
Copy link
Contributor Author

latot commented Nov 30, 2019

Ouuu, this problem seems to be..., this really need a test for sympy.

Checking deeper, this seems to not belong to Zeta function, other example stieltjes, i don't really know where the bug is, but, that line from @oscarbenjamin , the doit(), seems to always be returning the same object, this will happen with any function that doesn't have a solution with elementary functions, or maybe implemented functions?

My big question to know why this happen, Why evalf works with Zeta, and fails with the derivative?

It almost seems that Zeta, have a different evalf with their derivative.

Maybe track in both cases how works evalf, and check in what section do differents things.

Hope this can help, and if other one can say other idea how work on this would be great, i don't know very well how works sympy.

Thx.

@oscarbenjamin
Copy link
Collaborator

The idea is that a Subs.evalf() can only really work by doing the substitution so doit forces the substitution to happen. Normally Subs.doit will return a non-Subs object. Certain cases like derivatives that can't be evaluated can only be represented as Subs so doit can't do anything.

Perhaps the fallback for evalf here could be to numerically evaluate the derivative:

In [81]: h = S(10)**-10                                                                                                           

In [82]: h                                                                                                                        
Out[82]: 1/10000000000

In [83]: differentiate_finite(zeta(x).diff(x), x, points=[2-h, 2, 2+h])                                                           
Out[83]: 
                         2                                                                              
  100000000000000000000π                           ⎛20000000001⎞                          ⎛19999999999- ──────────────────────── + 100000000000000000000ζ⎜───────────⎟ + 100000000000000000000ζ⎜───────────⎟
             310000000000⎠                          ⎝10000000000⎠

In [84]: _.evalf()                                                                                                                
Out[84]: 1.98928023429890

@latot
Copy link
Contributor Author

latot commented Nov 30, 2019

Yes, that wold works.

I have a question.
I don't get the why Zeta can be evaluated and their derivative not.

@oscarbenjamin
Copy link
Collaborator

I don't get the why Zeta can be evaluated and their derivative not.

I'm not sure exactly what you mean by "evaluated" there. SymPy does not have defined symbolic functions that can represent the derivatives of the zeta function. Those could be added and if they were then they could use explicit algorithms for numerical approximate evaluation:
https://en.wikipedia.org/wiki/Riemann_zeta_function#Numerical_algorithms

Otherwise if you mean the fact that zeta(x).diff(x).subs(x, 1).evalf() doesn't work then that's the whole point of this issue.

I think that a reasonable solution here is to add something to Subs.evalf that checks whether doit has returned a Subs and whether the expression in that Subs is a Derivative and then uses numerical differentiation (around the point to be substituted) to approximate the result.

@latot
Copy link
Contributor Author

latot commented Dec 1, 2019

Is there a way to know this, before?, i mean, intead check if doit() return a Subs, when the function is constructed, check this and set it as a propperty in a variable, so we will not need check this every time.

@oscarbenjamin
Copy link
Collaborator

There isn't a reasonable way to know if doit will return a Subs without calling doit.

@latot
Copy link
Contributor Author

latot commented Dec 2, 2019

In that case, when the expression is constructed, run doit and check there isn't vars (if there is a function f(x, y) and we only set x, doit will return the same), and if return the same save it in a var in the object.

Is to just avoid call doit every time that we want get some type of aproximation.

@oscarbenjamin
Copy link
Collaborator

Calling doit can be expensive so we don't want to call it if it can be avoided.

@latot
Copy link
Contributor Author

latot commented Dec 2, 2019

In that case, when we run doit save the result in the object, if not just keep it empty until we need it, that is other way

@oscarbenjamin
Copy link
Collaborator

That is probably worth doing but it should be handled as part of a broader overhaul of evaluation (#17280)

@latot
Copy link
Contributor Author

latot commented Dec 2, 2019

Okis, i think can works check if Subs return a derivative, i'm little worried if this happen with other object more than derivative, i think check if the expression have a symbolic function representation can work for all cases.

Now, it feels is similar to #17280, that would be better to store in a var to know this and don't check it every time.

In my case i don't know very well how check that in Sympy.

I don't know if that can be considerated part of this issue.

Maybe a simpler solution for now, is make a function that request this data, like exist_sympy_func_repr()? and for now only implement to recognize derivative, and in the new issue follow the complete implementation.

@latot
Copy link
Contributor Author

latot commented Dec 3, 2019

Hi, @oscarbenjamin if your are right with what we have until now, i think @sachin-4099 can work on this.

@sachin-4099
Copy link
Contributor

I am ready to start the work, but I want to know what all exactly needs to be changed. @oscarbenjamin @latot

@latot
Copy link
Contributor Author

latot commented Dec 3, 2019

Hi, writing for now, @oscarbenjamin can confirm this please.

First make a new function to recognize if a function have a sympy representation, store this in a var to don't run all the check every time for the same object, if you can find a nice way to this would be great, if not you can check if the return of doit is a Subs and a derivative, in that case use a numerical aproximation, if not, run evalf, in this case would be necesary open a new issue for a full implementation.

Thx.

@oscarbenjamin
Copy link
Collaborator

We don't need to store anything.

The Subs.evalf method should be changed to check if doit returns the object unchanged. Then if it is a Subs of a Derivative numerical differentiation can be attempted. Otherwise the object should just be returned unevaluated.

@latot
Copy link
Contributor Author

latot commented Dec 3, 2019

Hi, i think that can cause some problems, #11802 (comment), if you can comment that i would appreciate it.

@oscarbenjamin
Copy link
Collaborator

I don't understand what problems you are referring to

@latot
Copy link
Contributor Author

latot commented Dec 3, 2019

mm, the problem is not, the Subs is returning a derivative and that causing the recursion.

The problem is that the code doesn't detect if the expression have a symbolic representation, the derivative in Subs is a particular case of this, and when this happen there is a recursion.

To can handle this point i think we need check that, but how the check can be expensive, i think would be better sotore the value, similar to what we talk of doit.

And with that check (symbolic representation), define if use numeric aproximation or not.

Hope i could explain better this.

@oscarbenjamin
Copy link
Collaborator

Calling doit is supposed to find the symbolic representation if it exists. In this case the symbolic representation is not there so doit returns the result with unevaluated Derivative and Subs.

@latot
Copy link
Contributor Author

latot commented Dec 3, 2019

:O okis.

To follow both things, i'm still worried there to be other cases than derivative.

Ex:

>>> zeta(zeta(y)).subs(y, 2).doit()
ζ(π26)

I think we need a good way to detect this.

Edited:

For now, i think checking if there is a derivative, in any part of the result can works, but i really really would like a method that can says if there is a representative function or not in sympy, not start testing in this way, is too expensive while we would be able to know this type in a easly way.

@oscarbenjamin
Copy link
Collaborator

We can just check if doit returned the same result. If it did we don't want to call evalf because that will lead to infinite recursion:

newself = self.doit()
if newself != self:
    return self.evalf()

Then if doit did return self we can check if the Subs is for a Derivative:

if isinstance(self.args[0], Derivative):
    # attempt numerical differentiation

@latot
Copy link
Contributor Author

latot commented Dec 3, 2019

Here, i think a complete example would be better:

k=zeta(zeta(y)).diff().subs(y, 2)
>>> y.doit() == y
True
>>> isinstance(k, Derivative)
False
>>> k.evalf()
ddξ1ζ(ξ1)∣∣∣ξ1=π26ddyζ(y)∣∣∣y=2

@sachin-4099
Copy link
Contributor

Here itself @oscarbenjamin??

sympy/sympy/core/function.py

Lines 2135 to 2136 in c162866

def evalf(self, prec=None, **options):
return self.doit().evalf(prec, **options)

@oscarbenjamin
Copy link
Collaborator

Yes

@jmig5776
Copy link
Member

@sachin-4099 Please start a PR if you had understand what needs to be done and CC @oscarbenjamin and me to it for further review

@asmeurer
Copy link
Member

asmeurer commented Jan 8, 2020

mpmath has numerical differentiation routines http://mpmath.org/doc/current/calculus/differentiation.html?highlight=diff#mpmath.diff

@sachin-4099
Copy link
Contributor

sachin-4099 commented Jan 16, 2020

Then if doit did return self we can check if the Subs is for a Derivative:

if isinstance(self.args[0], Derivative):
    # attempt numerical differentiation

@oscarbenjamin can you please be specific as to what to code in the numerical differentiation part.

@latot
Copy link
Contributor Author

latot commented Jan 26, 2020

Hi, the idea is replace that #attemp with the differentiation from mpmath, the link of asmeurer have the function, the idea is send the actual derivarive with the right arguments (order of the derivateve, params, etc) to get the final value, then sympy will finish joining all the values.

@sachin-4099
Copy link
Contributor

def evalf(self, prec=None, **options):
        newself = self.doit()
        if newself != self:
            return self.evalf()
        if isinstance(self.args[0],Derivative):
            return diff(self.args[0].args[0],self.args[0].args[1].args[0],self.args[0].args[1].args[1])

Are you talking about this?

@latot
Copy link
Contributor Author

latot commented Jan 26, 2020

Yes, something like that, but be careful, you need use diff from mpmath for the aproximations, and need to work with mixed derivatives (not just one variable).
If you can, change the code and test this in your fork, i think in any case we need add tests to sympy to check this.
If you need values you can use wolframalpha to get some aproximations from test and compare with what you get.

@sachin-4099
Copy link
Contributor

sachin-4099 commented Jan 26, 2020

My current code looks like:

def evalf(self, prec=None, **options):
        newself = self.doit()
        if newself != self:
            return self.evalf()
        if isinstance(self.args[0],Derivative):
            symbol = self.args[0].args[1].args[0]
            return mpmath.diff(lambda symbol:self.args[0].args[0],self.args[2].args[0],self.args[0].args[1].args[1])

where self.args[0].args[0]-->zeta(x)
self.args[0].args[1].args[0]-->symbol-->x (the variable in the expression)
self.args[2].args[0]--> 4 (we have to substitute x=4 finally)
self.args[0].args[1].args[1]--> 3 (we have to to find the third derivative)

But N(p) is finally returning 0 instead of returning −0.0726408498913214
@oscarbenjamin @latot need your suggestions

@sachin-4099
Copy link
Contributor

sachin-4099 commented Jan 26, 2020

Moreover this 1 testcase in test_function.py is giving an exception/error:

assert Subs(f(x)*cos(y) + z, (x, y), (0, pi/3)).n(2) == \
       Subs(f(x)*cos(y) + z, (x, y), (0, pi/3)).evalf(2) == \
       z + Rational('1/2').n(2)*f(0)

RecursionError: maximum recursion depth exceeded while calling a Python object

@latot
Copy link
Contributor Author

latot commented Jan 26, 2020

For the first case, than can happend for the precision, it can be changed, in the docs explains, with addprec, relative and h.

To can do that test in a easy way, you can assign every value tu one variable, and next run the assert, so, if fails with the recursion you can know the line with the error (is failing in n, evalf or both?, maybe use try to can test both but we will not know what line inside sympy is the error).
For now, try this example testcase, at least while debugging, is just an idea:

try:
 a = Subs(f(x)*cos(y) + z, (x, y), (0, pi/3)).n(2)
except:
 print("failing a")
try:
 b = Subs(f(x)*cos(y) + z, (x, y), (0, pi/3)).evalf(2)
except:
 print("failing b")
try:
 c = z + Rational('1/2').n(2)*f(0)
except:
 print("failing c")
#yes, this now only to get the code line inside sympy
a = Subs(f(x)*cos(y) + z, (x, y), (0, pi/3)).n(2)
b = Subs(f(x)*cos(y) + z, (x, y), (0, pi/3)).evalf(2)
c = z + Rational('1/2').n(2)*f(0)

assert a == b == c

I know this isn't the most pretty test case :D but i think helps while debugging.
When at least the vales are returned, even of they are wrong, can be simplified.

Hope this can help.

@sachin-4099
Copy link
Contributor

sachin-4099 commented Jan 26, 2020

For the first case, than can happend for the precision, it can be changed, in the docs explains, with addprec, relative and h.

No success, it still returns 0. I have already tried these.

@latot
Copy link
Contributor Author

latot commented Jan 26, 2020

mm, can you test with other eq that the result is not 0? (and not close to it)

@sachin-4099
Copy link
Contributor

suggest some equations which will make use of the function we are trying to improve.

@latot
Copy link
Contributor Author

latot commented Jan 26, 2020

Hi, here you can try several examples like this:

https://www.wolframalpha.com/input/?i=diff%28zeta%28zeta%28x%29%29%29%2C+x%3D4

@oscarbenjamin
Copy link
Collaborator

I see that Derivative already has a doit_numerically which uses mpmath.diff although it only works for first derivatives:

In [1]: Derivative(zeta(x)).doit_numerically(2)                                                                                                
Out[1]: -0.937548254315844

We should adapt that so that it works for higher-order derivatives and then use it in Subs.evalf.

@sachin-4099
Copy link
Contributor

I see that Derivative already has a doit_numerically which uses mpmath.diff although it only works for first derivatives:

In [1]: Derivative(zeta(x)).doit_numerically(2)                                                                                                
Out[1]: -0.937548254315844

We should adapt that so that it works for higher-order derivatives and then use it in Subs.evalf.

Should we change Derivative - doit_numerically to work for higher order derivatives? Is that what you mean?

@oscarbenjamin
Copy link
Collaborator

Yes doit_numerically should be made to work for higher-order derivatives.

@sachin-4099
Copy link
Contributor

Ok so what are the required changes then?

@asmeurer
Copy link
Member

There are a few things that need to be done here:

  • Supporting for doing evalf on derivatives should be added, using http://mpmath.org/doc/current/calculus/differentiation.html?highlight=diff#mpmath.diff
  • For zeta specifically, the mpmath implementation supports evaluating the derivative via a third argument. We should support this directly, as it is likely mpmath uses a specialized derivative implementation that is more accurate and/or faster than diff(zeta(z)) (ditto for any other functions in mpmath that have a similar argument). http://mpmath.org/doc/current/functions/zeta.html
  • We should investigate the recursion error. It likely indicates a bug. Even though the above things aren't implemented, evalf should produce some sort of NotImplementedError, not a RecursionError.

I expect the second point will be the easiest.

@asmeurer
Copy link
Member

Another thing is we should make sure lambdify(modules='mpmath') also works correctly for derivatives. Depending on how we implement them in evalf, this may happen automatically, or not.

@asmeurer
Copy link
Member

Derivative support should be added for numpy as well using https://docs.scipy.org/doc/scipy/reference/generated/scipy.misc.derivative.html

@oscarbenjamin
Copy link
Collaborator

  • We should investigate the recursion error.

The recursion error is explained above: #11802 (comment)

@asmeurer
Copy link
Member

For lambdify, I think Derivative support has to be done in the printer. We can't do it via a replacement, because the function inside the derivative will be replaced with a numeric function which will evaluate to a numeric value. This also makes this very difficult to work around. You can only pass a custom Derivative to the modules if you know ahead of time what the Derivative will be, e.g., {"Derivative": lambda expr, x: mpmath.zeta(x, derivative=1)} if you know the Derivatives in your expression are Derivative(zeta(x), x).

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

No branches or pull requests

6 participants