Numpy hohe Präzision

Ich bin mit numpy und pyfits zu manipulieren Spektren und I erfordern eine hohe Präzision (sowas wie 8-10 Dezimalstellen auf einen Wert, der gehen könnte so hoch wie 10^12). Für die, die den Datentyp "decimal" wäre perfekt (float64 ist nicht gut genug), aber unfortunalely numpy.interp nicht mag:

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

Eine vereinfachte version des Codes, die ich benutze:

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

wo die "Welle" und "new_wave" sind ein ein-dimensionales numpy-array repräsentiert ein 1D-Spektrum. Dieser code wird benötigt, um meine shift-Spektren entlang der x-Achse (die Wellenlänge)

Mein größtes Problem ist, dass weiter unten der code, den ich Teile meine Spektren durch eine template-Spektrum konstruiert sich aus der Summe aller meiner Spektren werden und die eine analyse der Unterschiede und da ich nicht genügend Nachkommastellen ich bin immer Rundungsfehler. Irgendwelche Ideen?

Dank!

UPDATE:

Test Beispiel:

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

Der Fehler ist:

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

dies ist der nächste, den ich bieten kann, ohne mit den realen Spektren im fits-format.

UPDATE 2:
Einige typische Werte von meinem Spektren, gedruckt mit Dezimal:

  18786960689.118938446044921875
  18473926205.282184600830078125
  18325454516.792461395263671875
  18400241010.149127960205078125
2577901751996.03857421875
2571812230557.63330078125
2567431795280.80712890625

dem problem bin ich immer, wenn ich make-Vorgänge zwischen Ihnen, bekomme ich die Rundung up-Fehler. Zum Beispiel habe ich eine Vorlage erstellen, die für alle Spektren, die durch addieren aller von Ihnen. Dann habe ich diese Vorlage benutzen, zu normalisieren alle Spektren. Ein Beispiel:

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

NormSpectra = Spectra/Template

Dieser zurückkehren sollte mir nur das Rauschen in den Spektren (vorausgesetzt, die Vorlage ist eine gute Darstellung der Sterne). Ich habe versucht, die Normalisierung jedes Spektrum seiner gesamten flux 

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

ebenso wie die Vorlage, aber würde noch schlimmer kommen, aufrunden Fehler.

Verwenden von Decimal funktionieren würde, gut für mich, aber ich muss die "shift" - mein Spektren, so dass alle spektralen Merkmale/Linien ausgerichtet sind.

Hoffe, das macht Sinn?

  • Haben Sie versucht, mit numpy.longdouble als "dtype"? mail.scipy.org/pipermail/scipy-dev/2008-March/008562.html
  • habe es gerade ausprobiert, und ich bin stil-gettin "TypeError: array kann nicht sicher gewirkt erforderlich Typ" 🙁
  • können Sie uns einen kleinen Selbstversorger-Beispiel demonstriert das problem?
  • "so etwas wie 8-10 Dezimalstellen auf einen Wert, der gehen könnte so hoch wie 10^12" könnte problematisch sein für float32, aber es ist weit von den Grenzen der float64. Kannst du ein Beispiel zeigen, wie float64 ist nicht genug? Sie könnten versuchen, die Skalierung Ihres Problems zu.
  • Ihre aktualisierte code, der sich beschwert, dass die Konvertierung von float in eine Dezimalzahl ist nicht erlaubt. Jedenfalls Stimme ich mit @jordeca, dass ein hilfreicher Beispiel wäre zu zeigen, wie float64 ist nicht genug. Außerdem würde ich prüfen, (i) ob longdouble ist nicht genug, und (ii) ob ein problem besteht, wenn Sie scipy.interp1d statt numpy.interp.
  • Dezimal ist kein numpy dtype Objekt oder convertable zu eins, so dass Sie nie in der Lage sein, es zu benutzen als ein dtype. Was Numpy macht, ist zu sehen, 'oh, Sie eine beliebige Klasse, besser, erstellen Sie ein Python-Objekt array. 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: Cannot cast-array-Daten von dtype('O') zu dtype('float64') nach der Regel "sicher". Es unterstützt diese Theorie. Ich wäre sehr überrascht, wenn scipy.interp1d arbeitete für dieses Objekt-array, das Sie erstellt, da ein array..
  • .. seit ein array von dtype 'O' ist nicht als sich wie ein numerischer Typ, und NumPy werde nicht versuchen, führen Sie numerische Berechnungen mit.
  • Ich aktualisiert die post mit ein paar typische Werte für meine Spektren. Das Problem kommt, wenn ich versuche, Sie zu teilen jeder-Spektren von seiner Vorlage, zu bekommen, was sieht aus wie ein Rundungsfehler.
  • Ich dachte darüber nach, dem bus, da scipy verwendet numpy, ich dachte, es könnte die gleichen Probleme haben? jedenfalls habe ich beschlossen zu gehen und meine eigenen Interpolations-Funktion, ich brauche nichts besonderes, eine lineare interpolation ausreichend.
  • hum, ich verstehe das. Also besser meine eigenen Interpolations-routine in reinem python und so kann ich die decimal.
  • Tut mir Leid, aber ich weiß nicht, wie zur Reproduktion Ihrer Rundungsfehler aus den Daten, die Sie zur Verfügung gestellt. Übrigens, haben Sie versucht, ob die zusätzlichen bits dtype=np.longdouble sind genug? Auch, mpmath für Sie interessant sein könnten (es ist eine Bibliothek, die für multiprecision floating-point-Arithmetik in pure Python, als auch schnell zusammengestellt backends).
  • Hey danke, ich denke mal, dass mpmath helfen könnte :). Ich war auch auf der Suche auf das problem der falsche Weg: die meiste Zeit werde ich nicht benötigen diese Art der Präzision, die auf meine numpy Operationen, nur, wenn ich Normalisiere mein Spektren (dh, die Aufteilung von Sachen wie 18400241010.149127960205078125 von 18786960689.118938446044921875), so kann ich nur umgehen, der interp problem insgesamt. Danke Euch allen für die Hilfe!!!!
  • Wie viel signifikante Ziffern hat Ihre Daten? Ist das problem in der Rundungsfehler nur? Wenn Sie weniger als 14 dezimal signifikanten Ziffern, und Sie müssen die genaue Summe, die Sie verwenden können, Karatsuba-Algorithmus.

InformationsquelleAutor jorgehumberto | 2013-05-15
  1. 5

    Wie können Sie sicher sein np.float64? In typischer usecases, kann man erwarten, ~15 wichtige Fakten aus einer doppelten.

    Wenn Sie sicher sind, dass dies nicht genug ist, können Sie versuchen np.float128 (aka np.longdouble).

    Aber dein problem scheint tiefer als das: es scheint eine schlecht gestellte problem (division von großen zahlen durch kleine zahlen in der Regel sind). Und das ist nicht etwas, was Sie wollen. Die Erhöhung der Genauigkeit sollte das problem lösen bis zu einem gewissen Ausmaß, aber Sie begegnen werden, einige Daten, die erforderlich wäre float256/float512/etc. um zu vermeiden, pathologische Rundungsfehler.

    Ich würde Ihnen raten, zu erklären Ihr problem, nicht Ihre Lösung, so dass wir hoffentlich einen anderen Weg finden, es zu lösen in allen Fällen (XY-Problem).

    InformationsquelleAutor
Schreibe einen Kommentar