search

JuliaPy / SymPy.jl

Julia interface to SymPy via PyCall
http://juliapy.github.io/SymPy.jl/
MIT License
269 stars 61 forks source link

Incorrect substitution of multivariate derivatives #375

Open Gyoshi opened 4 years ago

Gyoshi commented 4 years ago 
using SymPy
@vars(u, w)
@symfuns(x)

dww = diff(x(u,w), (u, 1), (w, 2))(w=>0)
u0 = (
    u=>0,
    diff(x(u,w), (u, 2), (w, 4))(0,0)=>0,
)
diff(dww, (u, 3))(u0...)

for N in 1:6
    t = diff(x(u, w), (u, N), (w, 1))
    println(t, ':', t(0,0)(t(0,0)=>0))
end

gives the output

Derivative(x(u, w), u, w):0
Derivative(x(u, w), (u, 2), w):0
Derivative(x(u, w), (u, 3), w):0
Derivative(x(u, w), (u, 4), w):0
Derivative(x(u, w), (u, 5), w):Subs(Derivative(x(u, w), (u, 4), (w, 2)), (u, w), (0, 0))
Derivative(x(u, w), (u, 6), w):0

I expect them all to be 0, but for some reason the evaluation of a (u, 5), (w, 1) derivative and other multivariate derivatives of total order 6 (not shown here) do not self-substitute correctly after evaluating the line diff(dww, (u, 3))(u0...).

Gyoshi commented 4 years ago

Ok, so I realised the above is not an actual MWE; it only works after running some of my code. Will take some time to trace what is actually causing this

mzaffalon commented 4 years ago

Have you seen that Derivative(x(u, w), (u, 4), (w, 2)) shows it has taken the second derivative along w not the first?

The result look OK to me:

julia> @vars u w
(u, w)
julia> @symfuns(x, real=true)
(x,)
julia> diff(x(u, w), (u, 5), (w, 1))
   6
  d
------(x(u, w))
     5
dw du
Gyoshi commented 4 years ago

It's a strange one, but I updated the OP to include a working example now...as minimal as I could get it. It is annoying because the problematic code 'contaminates' the REPL session, so you have to restart it everytime you want to check if the result changed.

Gyoshi commented 4 years ago

@mzaffalon Yes, that does look OK, but the problem comes when doing substitution on the derivatives. Please see the OP again; I have updated it with a working example.

jverzani commented 4 years ago

I can't replicate, I get 0 for all the cases. Did you replicate this behaviour @mzaffalon?

mzaffalon commented 4 years ago

It's 0 in all cases for me as well, with SymPy 1.0.30.

Gyoshi commented 4 years ago

I tried updating SymPy 1.0.28->1.0.30 to no effect. \:( This sucks because it means I cannot really trust any of the results I get.

What is likely to be the problem then? PyCall? Julia? Python?

jverzani commented 4 years ago

My guess sympy. What is the output of sympy.__version__? (Mine is 1.5.1, but 1.6.2 is already a month old)

Gyoshi commented 4 years ago 
julia> sympy.__version__
"1.6.2"
jverzani commented 4 years ago

I'll have to update and check. I'll try and do so sometime soon.

jverzani commented 4 years ago

Well, this is odd. I have two different instances running each with SymPy 1.0.30, PyCall with 1.91.4, sympy.version at 1.6.2 and different answers for:

@vars u v
@symfuns x
phi = sympy.Derivative(x(u,v), (u,5), (v,1))
phi00 = phi.subs(u,0).subs(v,0)
phi00.subs(phi00, 0)

The different one is my local master, but that should be the same, and if not, certainly not what is seen here.

One thing I tried, we to create a new project, activate that, add SymPy@1.0.30 and run these commands. This worked. Does the above work as expected on your machine if you try that?

Gyoshi commented 4 years ago

create a new project, activate that, add SymPy@1.0.30 and run these commands

Neither your code nor mine in the OP gave the correct result after activating the fresh project. But after restarting the REPL both gave the correct results!? ... After some testing (restarting the REPL, then activating the project or not) it seems like the only thing that matters (regardless of project) is which block of code I run first: yours or mine. If I run mine first, both are wrong. If I run yours first, both are right.

no idea why

jverzani commented 4 years ago

Wierd, I'll have to play around with that to see what might trigger this. Definitely worthy of the "!?" punctuation!

jverzani commented 4 years ago

Yeah, wierd. I couldn't replicate this AM; then I upgraded to 1.5.2 and voila -- it was back. I then tried with a reduced environment and it went away; now I can't make the issue come back. I'm at a loss for now. If it is an issue for you, can you try this PyCall-only test:

using PyCall
sympy = PyCall.pyimport_conda("sympy", "sympy")

u,v = sympy.sympify("u,v")
x = sympy.Function("x")
ϕ = sympy.Derivative(x(u,v), (u,5), (v, 1))
ϕ₀₀ = sympy.Derivative(x(u,v), (u,5), (v, 1)).subs(u,0).subs(v,0)
ϕ₀₀.subs(ϕ₀₀, 0)
jverzani commented 4 years ago

I'm still curious if the issue appears with the PyCall only version of the problem I posted. If you have time, I'd love to hear.

On Fri, Sep 25, 2020 at 10:21 AM Johannes van den Berg < notifications@github.com> wrote:

I noticed this is sometimes an issue for other orders of derivatives too, so the workaround I'm using right now is to run this at the start of the session:

using SymPy@vars u w@symfuns x y zfor ui in 0:6 # greater or equal to the order of derivative which will be used for wi in 0:6 # idem phi = sympy.Derivative(x(u,w), (u,ui), (w,wi)) phi00 = phi.subs(u,0).subs(w,0) @assert phi00.subs(phi00, 0) == Sym(0)

  phi = sympy.Derivative(y(u,w), (u,ui), (w,wi))
  phi00 = phi.subs(u,0).subs(w,0)
  @assert phi00.subs(phi00, 0) == Sym(0)

  phi = sympy.Derivative(z(u,w), (u,ui), (w,wi))
  phi00 = phi.subs(u,0).subs(w,0)
  @assert phi00.subs(phi00, 0) == Sym(0)

endend

With this everything seems to work as expected. Very far from an elegant solution, but it will do for my purposes I think. Thanks for the help :)

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/JuliaPy/SymPy.jl/issues/375#issuecomment-698957992, or unsubscribe https://github.com/notifications/unsubscribe-auth/AADG6TD5GQ5TKP4TLFIP6ATSHSRPVANCNFSM4RXDNAWA .

-- John Verzani Department of Mathematics College of Staten Island, CUNY verzani@math.csi.cuny.edu

Gyoshi commented 4 years ago

By the way, the I had deleted the comment that you replied to. The "workaround" fixes the output of the MWE, but then when I run some of my other code it gives me wrong answers and the MWE returns the wrong value again. \:/

So...apparently PyCall was not installed. I figured it would be a dependency of SymPy.jl, but seems not. The PyCall test you wrote works fine, but for me the SymPy.jl equivalent never failed, so I'm not convinced it tells us much.

I rewrote my MWE in python, and it gave the same erroneous result. So I suppose it is not really a problem with SymPy.jl then.

jverzani commented 4 years ago

Thanks!

BTW, regardless of where the issue lies, this is a weird one with its intermittency.

jverzani commented 4 years ago

I'm trying to get an example to file a bug report, but this is working for me. Perhaps it is os specific. I'm on a Mac. Does this run for you?

julia> using PyCall

julia> py"""
import sympy
from sympy import *
u,v = sympify("u,v")
x = Function("x")
ϕ = Derivative(x(u,v), (u,5), (v, 1))
ϕ00 = Derivative(x(u,v), (u,5), (v, 1)).subs(u,0).subs(v,0)
print(sympy.__version__)
print(ϕ00.subs(ϕ00, 0))
"""
1.6.2
0
Gyoshi commented 4 years ago

Yes that runs for me no problem. (I am running on Windows.) Try this:

py"""
import sympy
from sympy import *
u,v = sympify("u,v")
x = Function("x")
f0 = diff(x(u,v), u, 4, v, 2).subs(v, 0)
g00 = diff(x(u,v), u, 2, v, 4).subs(u, 0).subs(v, 0)
f0.subs([(u, 0), (g00, 0)])
ϕ00 = diff(x(u,v), u, 5, v, 1).subs(u,0).subs(v,0)
print(sympy.__version__)
print(ϕ00.subs(ϕ00, 0))
"""
jverzani commented 4 years ago

So the straight python code runs, but the pycall interface doesn't always?

That is 

using PyCall
sympy = PyCall.pyimport_conda("sympy", "sympy")

u,v = sympy.sympify("u,v")
x = sympy.Function("x")
ϕ = sympy.Derivative(x(u,v), (u,5), (v, 1))
ϕ₀₀ = sympy.Derivative(x(u,v), (u,5), (v, 1)).subs(u,0).subs(v,0)
ϕ₀₀.subs(ϕ₀₀, 0)

can fail, but

py"""
import sympy
from sympy import *
u,v = sympify("u,v")
x = Function("x")
f0 = diff(x(u,v), u, 4, v, 2).subs(v, 0)
g00 = diff(x(u,v), u, 2, v, 4).subs(u, 0).subs(v, 0)
f0.subs([(u, 0), (g00, 0)])
ϕ00 = diff(x(u,v), u, 5, v, 1).subs(u,0).subs(v,0)
print(sympy.__version__)
print(ϕ00.subs(ϕ00, 0))

doesn't. If that is so, it would seem that perhaps PyCall is the bottleneck.

Gyoshi commented 4 years ago

Ah no, I added a few lines there. The code gives the same result with Python vs PyCall.

2020年9月28日(月) 19:39 john verzani notifications@github.com:

So the straight python code runs, but the pycall interface doesn't always?

That is

using PyCall

sympy = PyCall.pyimport_conda("sympy", "sympy")

u,v = sympy.sympify("u,v")

x = sympy.Function("x")

ϕ = sympy.Derivative(x(u,v), (u,5), (v, 1))

ϕ₀₀ = sympy.Derivative(x(u,v), (u,5), (v, 1)).subs(u,0).subs(v,0)

ϕ₀₀.subs(ϕ₀₀, 0)

can fail, but

py"""

import sympy

from sympy import *

u,v = sympify("u,v")

x = Function("x")

f0 = diff(x(u,v), u, 4, v, 2).subs(v, 0)

g00 = diff(x(u,v), u, 2, v, 4).subs(u, 0).subs(v, 0)

f0.subs([(u, 0), (g00, 0)])

ϕ00 = diff(x(u,v), u, 5, v, 1).subs(u,0).subs(v,0)

print(sympy.version)

print(ϕ00.subs(ϕ00, 0))

doesn't. If that is so, it would seem that perhaps PyCall is the bottleneck.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaPy/SymPy.jl/issues/375#issuecomment-700180783, or unsubscribe https://github.com/notifications/unsubscribe-auth/AJ6PMBIWRYGCDYVPVUOUDY3SIDC37ANCNFSM4RXDNAWA .

jverzani commented 4 years ago

Alright, now I'm confused. Both are working or none?

Gyoshi commented 4 years ago

Neither Python nor PyCall gives the correct result when I include the lines:


f0 = diff(x(u,v), u, 4, v, 2).subs(v, 0)

g00 = diff(x(u,v), u, 2, v, 4).subs(u, 0).subs(v, 0)

f0.subs([(u, 0), (g00, 0)])

2020年9月28日(月) 19:44 john verzani notifications@github.com:

Alright, now I'm confused. Both are working or none?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaPy/SymPy.jl/issues/375#issuecomment-700183382, or unsubscribe https://github.com/notifications/unsubscribe-auth/AJ6PMBLWFCVHHTD6JF6A32TSIDDOXANCNFSM4RXDNAWA .