Ask Your Question
1

find_root fails in for loop with type coercion error

asked 2021-01-15 17:54:16 +0100

akhhhh gravatar image

Hi there, I'm trying to reproduce a model in the literature. Here is the code that causes the error (description+stack trace follows):

electron_chg = 1.6*10**(-19) # Coulombs
eps0 = 8.85418782e-12 # m^-3 kg^-1 s^4 A^2 = F/m
angs_per_meter = 1*10**10
T1 = 108
T11 = 120
T111 = 120
a10 = 3.34e5
a110 = 4.69e6
a1110 = -5.52e7
a12 = 3.23e8
a112 = 4.47e9
a123 = 4.91e9
m110 = -202e6
m1110 = 276e7

var('Px,Py,Pz, T')
a1(T) = a10*(T - T1)
a11(T) = a110*(T - T11)+m110
a111(T) = a1110*(T - T111)+m1110
Glan(Px,Py,Pz, T) = (a1(T))*(Px**2+Py**2+Pz**2) + (a11(T))*(Px**4+Py**4+Pz**4) + (a111(T))*(Px**6+Py**6+Pz**6) + a12*((Px**2)*(Py**2)+(Pz**2)*(Py**2)+(Px**2)*(Pz**2)) + a112*((Px**2)*(Py**4+Pz**4)+(Py**2)*(Px**4+Pz**4)+(Pz**2)*(Px**4+Py**4)) + a123*(Px**2)*(Py**2)*(Pz**2)

Ex(Px,Py,Pz, T) = derivative(Glan,Px)
Ey(Px,Py,Pz, T) = derivative(Glan,Py)
Ez(Px,Py,Pz, T) = derivative(Glan,Pz)

# this works
var('Ptet')
find_root(Ez(0,0,Ptet,32),0.0001,1.0)

# this doesn't, even for i = 0, which is the case that works above
var('Ptet')
Ptetx = zero_vector(101)
Ptety = zero_vector(101)
Ptetz = zero_vector(101)
for i in range(101):
    print(i)
    Ptetx[i] = find_root(Ex(0,0,Ptet,32+i),0.0001,1.0)
    Ptety[i] = find_root(Ey(0,0,Ptet,32+i),0.0001,1.0)
    Ptetz[i] = find_root(Ez(0,0,Ptet,32+i),0.0001,1.0)

Ptetx

Here is the full trace of the error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-9-3d8c5561fbc3> in <module>
      5 for i in range(Integer(101)):
      6     print(i)
----> 7     Ptetx[i] = find_root(Ex(Integer(0),Integer(0),Ptet,Integer(32)+i),RealNumber('0.0001'),RealNumber('1.0'))
      8     Ptety[i] = find_root(Ey(Integer(0),Integer(0),Ptet,Integer(32)+i),RealNumber('0.0001'),RealNumber('1.0'))
      9     Ptetz[i] = find_root(Ez(Integer(0),Integer(0),Ptet,Integer(32)+i),RealNumber('0.0001'),RealNumber('1.0'))

/opt/sagemath-9.2/local/lib/python3.7/site-packages/sage/modules/free_module_element.pyx in sage.modules.free_module_element.FreeModuleElement.__setitem__ (build/cythonized/sage/modules/free_module_element.c:13856)()
   1813             if n < 0 or n >= d:
   1814                 raise IndexError("vector index out of range")
-> 1815             self.set_unsafe(n, R(value))
   1816 
   1817     cdef int set_unsafe(self, Py_ssize_t i, value) except -1:

/opt/sagemath-9.2/local/lib/python3.7/site-packages/sage/structure/parent.pyx in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9337)()
    898         if mor is not None:
    899             if no_extra_args:
--> 900                 return mor._call_(x)
    901             else:
    902                 return mor._call_with_args(x, args, kwds)

/opt/sagemath-9.2/local/lib/python3.7/site-packages/sage/structure/coerce_maps.pyx in sage.structure.coerce_maps.NamedConvertMap._call_ (build/cythonized/sage/structure/coerce_maps.c:6043)()
    285             raise TypeError("Cannot coerce {} to {}".format(x, C))
    286         cdef Map m
--> 287         cdef Element e = method(C)
    288         if e is None:
    289             raise RuntimeError("BUG in coercion model: {} method of {} returned None".format(self.method_name, type(x)))

/opt/sagemath-9.2/local/lib/python3.7/site-packages/sage/rings/real_mpfr.pyx in sage.rings.real_mpfr.RealNumber._integer_ (build/cythonized/sage/rings/real_mpfr.c:17319)()
   2227             return n
   2228 
-> 2229         raise TypeError("Attempt to coerce non-integral RealNumber to Integer")
   2230 
   2231     def integer_part(self):

TypeError: Attempt to coerce non-integral RealNumber to Integer

The error describes itself as a type coercion, but I can't see what the origin is, given that the 0 case works outside the loop. Could it be related to how range() works, is it different than the python implementation?

edit retag flag offensive close merge delete

1 Answer

Sort by ยป oldest newest most voted
1

answered 2021-01-15 19:37:58 +0100

Emmanuel Charpentier gravatar image

The error you reported is attributable to your initial allocation, which, by default, creates vectors of Integers, in which one cannot store floats. Try :

Ptetx = zero_vector(RR, 101)
Ptety = zero_vector(RR, 101)
Ptetz = zero_vector(RR 101)

(you can also try RDF : look up this in previous questions...).

However, this is not sufficient : it happens that your functions have no root for some i values, and this is not caught by your code, which stops at the first occurrence. You should bracket your root findings in try: except constructs (or, alternatively, use utility functions wrapping find_root in such constructs).

An interesting problem is to record this "no solution" occurrence in your results vectors. A possibility is RR("NaN"), but this depends in a large part on what you plan to do with these results...

HTH,

edit flag offensive delete link more

Comments

Thanks so much Emmanuel! And thank you for the suggestion about using try and except constructs for the instances with no root--as implemented the model is actually incomplete, so hopefully when I finish setting it up it will have roots everywhere as expected, but in the mean time the way I have it set up should allow me to get a better understanding of its behavior before building it out fully.

akhhhh gravatar imageakhhhh ( 2021-01-15 19:54:18 +0100 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2021-01-15 17:54:16 +0100

Seen: 502 times

Last updated: Jan 15 '21