【问题标题】:TypeError: Cannot cast array data from dtype('O') to dtype('float64') according to the rule 'safe'TypeError:无法根据规则“安全”将数组数据从 dtype('O') 转换为 dtype('float64')
【发布时间】:2020-01-08 19:39:39
【问题描述】:

我需要对 g(u)jn(u) 类型进行积分,其中 g(u) 是一个没有零的平滑函数,而贝塞尔函数中的 jn(u) 是具有无穷大零的,但出现以下错误:

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

首先我需要将变量 x 更改为变量 u 并在新变量 u 中进行积分,但是函数 u(x) 在解析上是不可逆的,所以我需要使用插值来进行数值求逆。

import numpy as np
from scipy.interpolate import InterpolatedUnivariateSpline

x = np.linspace(0.1, 100, 1000)
u = lambda x: x*np.exp(x)
dxdu_x = lambda x: 1/((1+x) * np.exp(x))               ## dxdu as function of x: not invertible
dxdu_u = InterpolatedUnivariateSpline(u(x), dxdu_x(x)) ## dxdu as function of u: change of variable

在此之后,积分为:

from mpmath import mp

def f(n):
    integrand = lambda U: dxdu_u(U) * mp.besselj(n,U)
    bjz = lambda nth: mp.besseljzero(n, nth)
    return mp.quadosc(integrand, [0,mp.inf], zeros=bjz)

我使用来自mpmathquadosc 而不是来自scipyquad,因为quadosc 更适合于对快速振荡函数(如贝塞尔函数)进行积分。但是,另一方面,这迫使我使用两个不同的包,scipy 通过插值计算dxdu_umpmath 计算贝塞尔函数mp.besselj(n,U) 和乘积的积分dxdu_u(U) * mp.bessel(n,U) 所以我怀疑两种不同软件包的混合可能会产生一些问题/冲突。所以当我制作时:

print(f(0))

我得到了错误:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-38-ac2976a6b736> in <module>
     12     return mp.quadosc(integrand, [0,mp.inf], zeros=bjz)
     13 
---> 14 f(0)

<ipython-input-38-ac2976a6b736> in f(n)
     10     integrand = lambda U: dxdu_u(U) * mp.besselj(n,U)
     11     bjz = lambda nth: mp.besseljzero(n, nth)
---> 12     return mp.quadosc(integrand, [0,mp.inf], zeros=bjz)
     13 
     14 f(0)

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

有谁知道我该如何解决这个问题? 谢谢

【问题讨论】:

  • dtype('O') 表示您有一个对象数组。我认为 mp.quadosc 返回某种形式的多精度类型,这是 numpy 无法识别的
  • 是的,我想这就是重点。我不知道是否存在某种方式使 mpmath 和 numyp 相互通信。
  • numpy/scipy 函数使用编译代码时,数字输入必须与C 兼容doublesmpmath 对象不是双精度对象。 numpy 不“知道”任何关于 mpmath 的信息,所以你真的无法让它们交流,除非将 mpmath 对象转换为常规数字。
  • @hpulj 确实,很好的发现!不过这里的问题更简单,因为 mpf/double 问题发生在插值中

标签: python scipy typeerror mpmath


【解决方案1】:

完整的回溯(您狙击的部分)显示错误出在 univariatespline 对象的 __call__ 方法中。所以确实问题在于 mpmath 集成例程输入其 mpf 小数,而 scipy 无法处理它们。

最简单的解决方法是手动将 integrand 的参数的违规部分转换为浮点数:

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

一般来说,这很容易出现数值错误(mpmath 故意使用其高精度变量!)因此请谨慎操作。在这种特定情况下,它可能没问题,因为插值实际上是以双精度完成的。不过,最好检查一下结果。

一种可能的替代方法可能是避免使用 mpmath 并使用 scipy.integrate.quadweights 参数,请参阅 docs(向下滚动到 weights="sin" 部分)

另一种选择是一直坚持使用 mpmath 并在纯 python 中自己实现插值(这样,mpf 对象可能很好,因为它们应该支持通常的算术)。一个简单的线性插值可能就足够了。如果不是,编写您自己的三次样条插值器也没什么大不了的。

【讨论】:

    【解决方案2】:

    完整的追溯:

    In [443]: f(0)                                                                  
    ---------------------------------------------------------------------------
    TypeError                                 Traceback (most recent call last)
    <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):
       1002             return 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)
        255 
        256 
    
    /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) not in types: a = ctx.convert(a)
        946             if type(b) not in types: b = ctx.convert(b)
    
    /usr/local/lib/python3.6/dist-packages/mpmath/calculus/quadrature.py in <genexpr>(.0)
        252         case 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 
    

    此时它开始使用scipy 插值代码

    /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)
        313 
        314     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)
        369 
        370 
    
    /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)
        599 
        600         if ier == 10:
    
    TypeError: Cannot cast array data from dtype('O') to dtype('float64') according to the rule 'safe'
    

    _fitpack._spl_ 可能是编译代码(为了速度)。不能直接取mpmath对象;它必须将它们的值作为 C 兼容的双精度值传递。

    为了说明问题,创建一个由 mpmath 对象组成的 numpy 数组:

    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 call last)
    <ipython-input-448-4860036bcca8> in <module>
    ----> 1 arr.astype(float, casting='safe')
    
    TypeError: Cannot cast array from dtype('O') to dtype('float64') according to the rule 'safe'
    

    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')
    

    【讨论】:

      猜你喜欢
      • 2019-02-26
      • 2021-01-26
      • 2019-03-11
      • 2020-09-06
      • 1970-01-01
      • 2018-08-15
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      相关资源
      最近更新 更多