Typeerror: Cannot Cast Array Data From Dtype('o') To Dtype('float64') According To The Rule 'safe'

I need to make an integral of the type g(u)jn(u) where g(u) is a smooth function without zeros and jn(u) in the Bessel function with infinity zeros, but I got the following error:

Solution 1:

The full traceback (the part you sniped) shows that the error is in the __call__ method of the univariatespline object. So indeed the problem is that the mpmath integration routine feeds in its mpf decimals, and scipy has no way of dealing with them.

A simplest fix is then to manually cast the offending part of the argument of the integrand to a float:

integrand = lambda U: dxdu_u(float(U)) * mp.besselj(n,U)

In general this is prone to numerical errors (mpmath uses its high-precision variables on purpose!) so proceed with caution. In this specific case it might be OK, because the interpolation is actually done in double precision. Still, best check the results.

A possible alternative might be to avoid mpmath and use the weights argument to scipy.integrate.quad, see the docs (scroll down to weights="sin" part)

Another alternative is to stick with mpmath all the way and implement the interpolation yourselves in pure python (this way, mpf objects are probably fine since they should support usual arithmetics). It's likely a simple linear interpolation is enough. If it's not, it's not too big of a deal to code up your own cubic spline interpolator.

Solution 2:

The full traceback:

In [443]: f(0)                                                                  
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent calllast)
<ipython-input-443-6bfbdbfff9c4>in<module>----> 1 f(0)<ipython-input-440-7ebeff3611f6>in f(n)
      2     integrand = lambda U: dxdu_u(U) * mp.besselj(n,U)
      3     bjz = lambda nth: mp.besseljzero(n, nth)
----> 4     return mp.quadosc(integrand, [0,mp.inf], zeros=bjz)5/usr/local/lib/python3.6/dist-packages/mpmath/calculus/quadrature.py in quadosc(ctx, f, interval, omega, period, zeros)
    998         #    raise ValueError("zeros do not appear to be correctly indexed")
    999         n =1->1000         s = ctx.quadgl(f, [a, zeros(n)])
   1001         def term(k):
   1002return ctx.quadgl(f, [zeros(k), zeros(k+1)])

/usr/local/lib/python3.6/dist-packages/mpmath/calculus/quadrature.py in quadgl(ctx, *args, **kwargs)
    807         """
    808         kwargs['method'] = 'gauss-legendre'
--> 809         return ctx.quad(*args, **kwargs)
    810 
    811     def quadosc(ctx, f, interval, omega=None, period=None, zeros=None):

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

/usr/local/lib/python3.6/dist-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)

/usr/local/lib/python3.6/dist-packages/mpmath/calculus/quadrature.py in sum_next(self, f, nodes, degree, prec, previous, verbose)
    252         case the quadrature rule is able to reuse them.
    253         """
--> 254         return self.ctx.fdot((w, f(x)) for (x,w) in nodes)255256/usr/local/lib/python3.6/dist-packages/mpmath/ctx_mp_python.py in fdot(ctx, A, B, conjugate)
    942         hasattr_ = hasattr
    943         types = (ctx.mpf, ctx.mpc)
--> 944         for a, b in A:945             if type(a) notin types: a = ctx.convert(a)
    946             if type(b) notin types: b = ctx.convert(b)

/usr/local/lib/python3.6/dist-packages/mpmath/calculus/quadrature.py in<genexpr>(.0)
    252case the quadrature rule is able to reuse them.
    253         """
--> 254         return self.ctx.fdot((w, f(x)) for (x,w) in nodes)
    255 
    256 

<ipython-input-440-7ebeff3611f6> in <lambda>(U)
      1 def f(n):
----> 2     integrand = lambda U: dxdu_u(U) * mp.besselj(n,U)
      3     bjz = lambda nth: mp.besseljzero(n, nth)
      4     return mp.quadosc(integrand, [0,mp.inf], zeros=bjz)
      5 

at this point it starts using the scipy interpolation code

/usr/local/lib/python3.6/dist-packages/scipy/interpolate/fitpack2.py in __call__(self, x, nu, ext)
    310             except KeyError:
    311                 raise ValueError("Unknown extrapolation mode %s." % ext)
--> 312         return fitpack.splev(x, self._eval_args, der=nu, ext=ext)
    313314     def get_knots(self):

/usr/local/lib/python3.6/dist-packages/scipy/interpolate/fitpack.py in splev(x, tck, der, ext)
    366         return tck(x, der, extrapolate=extrapolate)
    367     else:
--> 368         return _impl.splev(x, tck, der, ext)
    369370 

/usr/local/lib/python3.6/dist-packages/scipy/interpolate/_fitpack_impl.py in splev(x, tck, der, ext)
    596         shape = x.shape
    597         x = atleast_1d(x).ravel()
--> 598         y, ier = _fitpack._spl_(x, der, t, c, k, ext)
    599600         if ier == 10:

TypeError: Cannot cast array data from dtype('O') to dtype('float64') according to the rule 'safe'

_fitpack._spl_ probably is compiled code (for speed). It can't take the mpmath objects directly; it has to pass their values as C compatible doubles.

To illustrate the problem, make a numpy array of mpmath objects:

In [444]: one,two = mp.mpmathify(1), mp.mpmathify(2)                            
In [445]: arr = np.array([one,two])                                             
In [446]: arr                                                                   
Out[446]: array([mpf('1.0'), mpf('2.0')], dtype=object)

In [447]: arr.astype(float)    # default'unsafe' casting                                                     
Out[447]: array([1., 2.])
In [448]: arr.astype(float, casting='safe')                                     
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent calllast)
<ipython-input-448-4860036bcca8>in<module>----> 1 arr.astype(float, casting='safe')

TypeError: Cannot cast arrayfrom dtype('O') to dtype('float64') according to the rule 'safe'

With integrand = lambda U: dxdu_u(float(U)) * mp.besselj(n,U),

In [453]: f(0)      # a minute or so later                                                                
Out[453]: mpf('0.61060303588231069')