uqfoundation / dill

serialize all of Python
http://dill.rtfd.io
Other
2.27k stars 181 forks source link

pickling a lambdified expression causes sympy function to throw ValueError #201

Closed MarcoMagl closed 7 years ago

MarcoMagl commented 7 years ago

I have described the issue here Here is a minimalist example: from pdb import set_trace import sympy as sp import dill as pickle

x, y = sp.symbols("x, y") 

f_example = 2*x**2 + 2*x*y
f_lbda= sp.lambdify((x, y),f_example ) 
pickle.settings['recurse'] = True 
# works 
fileW = open('file_where_I_dump.dill', 'wb')
pickle.dump([f_lbda, f_lbda], fileW)
fileR = open('file_where_I_dump.dill', 'rb')
f_lbda_loaded = pickle.load(fileR)

And a more complicated one where dill is working nicely to import sympy expressions, but not the corresponding lambda functions:

import sympy as sp
import dill
import pickle

# declare multiply symbolic variable at the same time 
x, y , Alpha, Ksi, x1x, x1y, x2x, x2y, x3x, x3y, N, nx, ny  = sp.symbols("x, y , Alpha, Ksi, x1x, x1y, x2x, x2y, x3x, x3y, N, nx, ny ") 
# gap function and penalty stiffness 
gN, kPenN = sp.symbols("gN, kPenN")

"""
Explanations for the symbolic variable declared
"""
S = sp.Matrix([x, y])
x1 = sp.Matrix([x1x, x1y])
x2 = sp.Matrix([x2x, x2y])
x3 = sp.Matrix([x3x, x3y])
# normal to the curve taken at Ksi = Ksi_bar
N  = sp.Matrix([nx, ny])
# control points
b0 = x2 + 0.5 * (x1 - x2)
#b3
b3 = x2 + 0.5 * (x3 - x2)
#b1
b1 = b0 + (x2- b0) * Alpha
#b2
b2 = x2 + (b3 - x2) * ( 1 - Alpha)

# Berstein polynomials
# B1 
B0 =  (1./8.) * ((1 - Ksi)**3)
# B1
B1 = (3./8.) * (1 - Ksi)**2 *  (1 + Ksi)
# B3
B2 =  (3./8.) * (1 - Ksi) * (1 + Ksi)**2
#B4
B3 = (1./8.) *(1 + Ksi)**3

# Berstein polynomials first order derivative
B0_Ksi = sp.diff(B0, Ksi)
B1_Ksi = sp.diff(B1, Ksi)
B2_Ksi = sp.diff(B2, Ksi)
B3_Ksi = sp.diff(B3, Ksi)

# expression of the position vector along the 
# Bezier polynomial depending on the convective 
# coordinate Ksi
x_interp = b0 * B0 +  b1 * B1 + b2 * B2 + b3 * B3
x_Ksi = sp.diff(x_interp, Ksi)
#second order derivative
x_KsiKsi = sp.diff(x_interp, Ksi, 2)

#check the shapes
for el in [x_Ksi, x_KsiKsi]:
    if el.shape != (2,1):
        raise Error('incorrect size')

"""
Important: all these matrixes will be evaluated numerically
for Ksi = Ksi_bar in the main file
"""
# Matrices used to compute different contribution to contact terms and so on 
Bn  = sp.flatten([N, -B0 * N, -B1 * N, -B2 * N, -B3 * N])
Bn_Ksi  = sp.flatten([sp.zeros(2,1), B0_Ksi * N, B1_Ksi * N, B2_Ksi * N, B3_Ksi * N])
Bt      = sp.flatten([x_Ksi, -B0_Ksi * x_Ksi, -B1_Ksi * x_Ksi, -B2_Ksi * x_Ksi, -B3_Ksi * x_Ksi]) 
Bn = sp.Matrix(Bn)
Bn_Ksi = sp.Matrix(Bn_Ksi)
Bt = sp.Matrix(Bt)
H_KsiKsi = x_Ksi.dot(x_Ksi) - gN * N.dot(x_KsiKsi)
B_Ksi  = (H_KsiKsi**(-1)) * ( Bt + gN*Bn_Ksi )
b_KsiKsi = x_KsiKsi.dot(N)
#check the shapes
for el in [Bn_Ksi, Bt, B_Ksi]:
    if el.shape != (10,1):
        raise Error('incorrect size')

# We now have all the elements to compute the contribution of ONE contact element to the residuum
# Note that I did not take into account the vector of variations because I do not use the 'test function'
# in my code
dWN = kPenN * gN * Bn_Ksi

"""
# just to check that the 2 expressions are equal
# except that norm() adds ABS that are not convenient to manipulate ..
expr1 = x_Ksi.norm() ** 2
expr2 =  x_Ksi.dot(x_Ksi)
set_trace()
"""

# computing the contribution to the stiffness locally. I used directly
# the form provided in 9.239
DdWN = kPenN * (   Bn.dot(Bn.T) - gN * ( Bn_Ksi.dot(B_Ksi.T) + B_Ksi.dot(Bn_Ksi.T) + b_KsiKsi * B_Ksi.dot(B_Ksi.T) - (gN / (x_Ksi.dot(x_Ksi)) * (Bn_Ksi + b_KsiKsi * B_Ksi).dot((Bn_Ksi + b_KsiKsi * B_Ksi).T)))) # only the terms corresponding to the variation of the dof of the slave node are
# interesting to me in the case of contact with a rigid seg 

# Saving the objects:
with open('normal_contact_Bezier.dill', 'wb') as f:  # Python 3: open(..., 'wb'):
    dill.dump([dWN , DdWN ] , f)

#to load the objects as they are
with open('/Users/marco.magliulo/repos/COMMON_ROOT_MYPROG/Python/OOP_2DLattice/Symbolic_computations/normal_contact_Bezier.dill', 'rb') as f:  # Python 3: open(..., 'rb')
    dWN , DdWN = pickle.load(f)

dWN_lbda= sp.lambdify((x, y , Alpha, Ksi, x1x, x1y, x2x, x2y, x3x, x3y, N, nx, ny, gN, kPenN), dWN) 
DdWN_lbda= sp.lambdify((x, y , Alpha, Ksi, x1x, x1y, x2x, x2y, x3x, x3y, N, nx, ny, gN, kPenN), DdWN) 

with open('normal_contact_Bezier.dill', 'wb') as f:  # Python 3: open(..., 'wb'):
    dill.dump([dWN_lbda, DdWN_lbda] , f)

with open('/Users/marco.magliulo/repos/COMMON_ROOT_MYPROG/Python/OOP_2DLattice/Symbolic_computations/normal_contact_Bezier.dill', 'rb') as f:  # Python 3: open(..., 'rb')
    dWN_lbda, DdWN_lbda= pickle.load(f)

please do not hesitate to contact me if you require more info. I use Python 3.5.2

mmckerns commented 7 years ago

The above is confusing as to what you are asking, and there's no traceback, and your link just points to this question… however, since I know of your original posting here, I'm going to assume I know what you are asking. Please correct me if I'm mistaken.

Here's the minimalist example of what happens:

>>> import sympy as sp
>>> x,y = sp.symbols('x,y')
>>> f_ex = 2*x**2 + 2*x*y
>>> f_la = sp.lambdify((x,y),f_ex)
>>> f_ex
2*x**2 + 2*x*y
>>> import dill
>>> dill.settings['recurse'] = True
>>> la = dill.loads(dill.dumps(f_la))
>>> la
<function <lambda> at 0x10b4ed668>
>>> f_ex
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sympy/core/basic.py", line 392, in __repr__
    return sstr(self, order=None)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sympy/printing/str.py", line 748, in sstr
    s = p.doprint(expr)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sympy/printing/printer.py", line 233, in doprint
    return self._str(self._print(expr))
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sympy/printing/printer.py", line 257, in _print
    return getattr(self, printmethod)(expr, *args, **kwargs)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sympy/printing/str.py", line 51, in _print_Add
    terms = self._as_ordered_terms(expr, order=order)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sympy/printing/printer.py", line 271, in _as_ordered_terms
    return expr.as_ordered_terms(order=order)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sympy/core/expr.py", line 888, in as_ordered_terms
    terms, gens = self.as_terms()
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sympy/core/expr.py", line 920, in as_terms
    coeff = complex(coeff)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sympy/core/expr.py", line 229, in __complex__
    result = self.evalf()
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sympy/core/evalf.py", line 1377, in evalf
    prec = dps_to_prec(n)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/mpmath/libmp/libmpf.py", line 67, in dps_to_prec
    return max(1, int(round((int(n)+1)*3.3219280948873626)))
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/core/fromnumeric.py", line 2293, in amax
    out=out, **kwargs)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/core/_methods.py", line 26, in _amax
    return umr_maximum(a, axis, None, out, keepdims)
ValueError: 'axis' entry is out of bounds

As I said on SO, I think this means that the serialization of one object is destroying a pointer or a weakref on another object. It's pretty weird, and I think most definitely a sympy issue. I'm willing to lend a hand, but it should be a PR to sympy and not dill.

MarcoMagl commented 7 years ago

You're right, sorry. I will close this thread and open a new one in sympy. I am going to take your title which is probably more accurate than the one I gave :)

mmckerns commented 7 years ago

I've just done it.