Numpy高精度答案

【问题标题】：Numpy high precisionNumpy高精度
【发布时间】：2013-05-15 15:18:24
【问题描述】：

我正在使用 numpy 和 pyfits 来操作光谱，并且我需要高精度（可能高达 10^12 的值的小数点后 8-10 位）。为此，“十进制”数据类型将是完美的（float64 不够好），但不幸的是 numpy.interp 不喜欢它：

File ".../modules/manip_fits.py", line 47, in get_shift
pix_shift = np.interp(x, xp, fp)-fp
File "/usr/lib/python2.7/dist-packages/numpy/lib/function_base.py", line 1053, in interp
return compiled_interp(x, xp, fp, left, right)
TypeError: array cannot be safely cast to required type

我使用的代码的简化版本：

fp = np.array(range(new_wave.shape[-1]),dtype=Decimal)
pix_shift = np.empty_like(wave,dtype=Decimal)
      x = wave
  xp = new_wave
 pix_shift = np.interp(x, xp, fp)-fp

其中 'wave' 和 'new_wave' 是表示一维频谱的一维 numpy 数组。需要此代码来沿 x 轴（即波长）移动我的光谱

我最大的问题是，在代码的后面，我将光谱除以由所有光谱的总和构建的模板光谱，以分析差异，并且由于我没有足够的小数位，我得到了舍入错误。有什么想法吗？

谢谢！

更新：

测试示例：

import numpy as np
from decimal import *
getcontext().prec = 12

wave = np.array([Decimal(xx*np.pi) for xx in range(0,10)],dtype=np.dtype(Decimal))
new_wave = np.array([Decimal(xx*np.pi+0.5) for xx in range(0,10)],dtype=np.dtype(Decimal))

fp = np.array(range(new_wave.shape[-1]),dtype=Decimal)
pix_shift = np.empty_like(wave,dtype=Decimal)

x = wave
xp = new_wave
pix_shift = np.interp(x, xp, fp)-fp

错误是：

Traceback (most recent call last):
  File "untitled.py", line 16, in <module>
    pix_shift = np.interp(x, xp, fp)-fp
  File "/usr/lib/python2.7/dist-packages/numpy/lib/function_base.py", line 1053, in interp
    return compiled_interp(x, xp, fp, left, right)
TypeError: array cannot be safely cast to required type

这是我在不使用拟合格式的真实光谱的情况下可以提供的最接近的值。

更新 2：我的光谱的一些典型值，使用十进制打印：

  18786960689.118938446044921875
  18473926205.282184600830078125
  18325454516.792461395263671875
  18400241010.149127960205078125
2577901751996.03857421875
2571812230557.63330078125
2567431795280.80712890625

我遇到的问题是，当我在它们之间进行操作时，会出现四舍五入的错误。例如，我通过对所有光谱求和来为所有光谱创建一个模板。然后我使用这个模板来标准化每个光谱。一个例子：

Spectra = np.array([Spectrum1, Spectrum2, ...])
Template = np.nansum(Spectra, axis= 0)

NormSpectra = Spectra/Template

这应该只返回光谱上的噪声（假设模板是恒星的良好表示）。我尝试将每个光谱归一化为其总通量

(Spectrum1 = Spectrum1/np.nansum(Spectrum1), ...)

以及模板，但会变得更糟糕的四舍五入错误。

使用 Decimal 对我来说效果很好，但我需要“移动”我的光谱，以便对齐所有光谱特征/线。

希望这有意义吗？

【问题讨论】：

您是否尝试过使用numpy.longdouble 作为数据类型？ mail.scipy.org/pipermail/scipy-dev/2008-March/008562.html
“可能高达 10^12 的值的小数点后 8-10 位”对于 float32 可能是有问题的，但它远未达到 float64 的限制。你能发布一个例子来说明 float64 是如何不够的吗？您也可以尝试扩展您的问题。
Decimal 不是 numpy dtype 对象或可转换为 1，因此您将永远无法将其用作 dtype。 Numpy 正在做的是看到'哦，你提供了一个任意类，最好创建一个 Python-Object 数组. Which works okay, until you have to cast back to a NumPy datatype to do the calculation.I get a more useful error message from your code, which is TypeError：无法根据规则将数组数据从 dtype('O') 转换为 dtype('float64')'安全'`。它支持这个理论。如果 scipy.interp1d 也适用于您创建的这个对象数组，我会感到非常惊讶，因为数组......
.. 因为 dtype 'O' 的数组不被视为数值类型，NumPy 不会尝试用它执行数值计算。
很抱歉，我不知道如何从您提供的数据中重现您的舍入错误。顺便问一下，dtype=np.longdouble 的额外位是否足够？此外，mpmath 对您来说可能很有趣（它是纯 Python 中的多精度浮点运算库，也可以使用快速编译的后端）。

标签： python arrays numpy decimal

【解决方案1】：

你怎么能确定np.float64？在典型的用例中，可以从 double 中获得约 15 个有效数字。

如果你确定这还不够，你可以试试np.float128（又名np.longdouble）。

但是您的问题似乎比这更深：这似乎是一个不适定问题（通常是大数除以小数）。这不是你想要的。提高精度应该可以在一定程度上解决问题，但是你会遇到一些需要float256/float512/等的数据。以避免病理性舍入误差。

我建议您解释您的问题，而不是您的解决方案，以便我们希望在每种情况下都能找到另一种解决方法 (XY Problem)。

【讨论】：