Numpy高精度

Question

我使用numpy和pyfits来操纵光谱，并且我需要高精度（如可能高达10^12的值的8-10小数位）。对于数据类型“十进制”将是完美的（float64不够好），但unfortunalely 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轴（这是波长）移动我的光谱

我最大的问题是，进一步下来的代码我的光谱由我的所有光谱的总和构建的模板谱分析差异，并且由于我没有足够的小数位，所以我正在舍入错误。有任何想法吗？

谢谢！

UPDATE：

试验例：

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), ...) 

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

使用小数对我来说很好，但我需要“移动”我的光谱，以便所有光谱特征/行都对齐。

希望这有意义吗？

Answer 1

你怎么确定np.float64？在典型的使用情况下，人们可以期望从双倍数中得到〜15个有效数字。

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

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

我会建议你解释你的问题，而不是你的解决方案，这样我们就可以在每种情况下找到另一种解决方法（XY Problem）。