2013-05-15 15 views
9

Sto usando numpy e pyfits per manipolare gli spettri e ho bisogno di alta precisione (qualcosa come 8-10 posizioni decimali su un valore che potrebbe arrivare fino a 10^12). Per questo il tipo di dati "decimali" sarebbe perfetto (float64 non è abbastanza buono), ma numpy.interp unfortunalely non piace:Numpy alta precisione

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 

Una versione semplificata del codice che uso:

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 

dove 'wave' e 'new_wave' sono una matrice numpy a una dimensione che rappresenta uno spettro 1D. Questo codice è necessario per spostare i miei spettri lungo l'asse x (che è la lunghezza d'onda)

Il mio problema più grande è che più in basso nel codice divido i miei spettri con uno spettro modello costruito dalla somma di tutti i miei spettri per analizzare le differenze e poiché non ho abbastanza cifre decimali ho degli errori di arrotondamento. Qualche idea?

Grazie!

UPDATE:

test Esempio:

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 

L'errore è:

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 

questo è il più vicino posso fornire senza utilizzare il vero spettri in formato FITS.

UPDATE 2: Alcuni valori tipici della mia spettri, stampata utilizzando decimale:

18786960689.118938446044921875 
    18473926205.282184600830078125 
    18325454516.792461395263671875 
    18400241010.149127960205078125 
2577901751996.03857421875 
2571812230557.63330078125 
2567431795280.80712890625 

il problema che sto ottenendo è quando faccio le operazioni tra di loro, vengo arrotondamenti errori. Ad esempio, creo un modello per tutti gli spettri sommandoli tutti. Quindi uso questo modello per normalizzare ogni spettro. Un esempio:

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

NormSpectra = Spectra/Template 

Questo mi deve restituire solo il rumore del spettri (supponendo che il modello è una buona rappresentazione della stella). Ho provato normalizzare ogni spettro al suo flusso totale

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

così come template, ma sarebbe ottenere ancora peggiori errori di arrotondamento in su.

L'utilizzo di Decimal funzionerebbe correttamente per me, ma ho bisogno di "spostare" i miei spettri in modo che tutte le caratteristiche/linee spettrali siano allineate.

Spero che questo abbia senso?

+0

Hai provato a usare 'numpy.longdouble' come dtype? http://mail.scipy.org/pipermail/scipy-dev/2008-March/008562.html –

+0

@Zhenya: ho appena provato, e sto facendo lo stilin "TypeError: array non può essere lanciato in modo sicuro nel tipo richiesto" :( – jorgehumberto

+0

puoi fornire un piccolo esempio autosufficiente che dimostri il problema? –

risposta

5

Come si può essere sicuri np.float64? In casi tipici, ci si può aspettare ~ 15 cifre significative da un doppio.

Se si è certi che questo non è sufficiente, è possibile provare (alias np.longdouble).

Ma il il tuo problema sembra più profondo di quello: sembra essere un problema mal posto (la divisione di grandi numeri con numeri piccoli di solito lo sono). E questo non è qualcosa che vuoi.Aumentare la precisione dovrebbe risolvere il problema in una certa misura, ma si incontreranno alcuni dati che richiederebbero float256/float512/etc. per evitare errori di arrotondamento patologici.

Vorrei consigliarvi di spiegare il vostro problema, non la vostra soluzione in modo da poter sperare di trovare un altro modo per risolverlo in tutti i casi (XY Problem).

Problemi correlati