Röntgenoptik-Projekt

[Gigaz]

Hallo erstmal an alle

Ich bin Doktorand für experimentelle Physik. Bei der Diskussion mit einem Kollegen hatte ich vor kurzem die Idee, dass wir möglichweise ein Python-Paket schreiben und veröffentlichen sollten.

In der ersten Version soll man damit die Röntgenreflektion und Transmission von aufgedampften Filmen berechnen können. Dabei ist die Performance das absolut entscheidende Kriterium, damit die Leute ihre Daten auch fitten können. Daher ist der Plan, die aufwändigen Berechnungen nach C auszulagern. "Aufwändig" bedeutet in diesem Kontext dass man einige 10.000 Datenpunke hat und für jeden Datenpunkt (Energie, Einfallswinkel) ein paar hundert 2x2-Matrizen multipliziert werden müssen. Die C-Algorithmen haben wir weitestgehend.

So, jetzt kommen die schwierigen Fragen bei denen ihr mir vielleicht weiterhelfen könnt.

1. Ist es klüger, für jeden Datenpunkt das C-Programm aufzurufen, oder sollte man dem C-Programm alle Datenpunkte geben und es nur einmal aufrufen? Wäre die zweite Variante eine deutliche Zeitersparnis? :K

2. So oder so muss das C-Programm jede Menge Zahlen übernehmen. energieabhängige Brechungsindizees, Schichtdicken, Rauigkeiten, usw.
Die Rückgabe ist je nach Antwort auf Frage 1 ein großer Datensatz oder nur wenige Werte.
Wie genau mache ichs am besten. Über eine .dll? Rufe ich eine .exe auf? Oder was ganz Anderes? :K

3. Gibt es irgendwas zu beachten bei der Parallelisierung? Der Python-Benutzer sollte mehrere Instanzen der Röntgenoptik-Funktion parallel laufen lassen können ohne dass die sich verheddern. Daten in Files auslagern wäre wohl tabu, aber sonst?

4. Irgendwelche sonstigen guten Tips?

[BlackJack]

@Gigaz: Numpy reicht nicht? Und was wäre mit Numpy + Cython?

[Gigaz]

Wenn ichs probiert hätte wüsste ich es. Man kommt aber um eine Iteration durch alle Schichten nicht herum, und darin ist Python bekanntlich langsam. Numpy ist wahrscheinlich keine Option. Cython vielleicht schon eher, ich muss es mir mal ansehen.

[Sirius3]

@Gigaz: sobald Du einige 100000 Datenpunkte hast, ist es eigentlich egal, wie die äußere Schleife 100 mal aufgerufen wird. Zumindest solltest Du, statt etwas neues in C zu programmieren eine Erweiterung für numpy in Betracht ziehen. Zum Fitten bringt numpy/scipy eigentlich fast alles mit, was man sich denken kann.

[Gigaz]

@Sirius man muss für jeden einzelnen Datenpunkt durch die Schichten-Schleife, die ist leider nicht außen rum.

Die Algorithmen existieren in C alle, wir haben sie selbst geschrieben und optimiert. Die Problemstellung gibt es leider nicht her dass man mit mathematisch optimierten Algorithmen was rausholt. Man braucht eine Programmiersprache die für möglichst schnelle Berechnungen gedacht ist. Also Fortran oder C. Aber wir hätten eben auch gerne die Python-Funktionalität drumrum. Außerdem kann fast jeder Physiker Python, aber nur wenige können C...

Ich werde wohl zuerst mal ausprobieren wie die Performance eines Cython-Moduls im Vergleich zu den reinen C-Versionen ist.

[Hyperion]

Man braucht eine Programmiersprache die für möglichst schnelle Berechnungen gedacht ist. Also Fortran oder C.

Es sei denn, man kann die Berechnung parallelisieren - das geht zwar auch in genannten Sprachen, dürfte aber wesentlich einfacher mit Rahmenwerken sein, die primär Bindings für die JVM bieten (Hazelcast z.B. o.ä.)

[Gigaz]

@Hyperion Danke. Über die Parallelisierung habe ich mir weiter oben Gedanken gemacht. Da letztlich der Benutzer auch seine Daten fitten will ist es nicht sinnvoll auf dieser niedrigen Ebene den Code zu parallelisieren. Die Parallelisierung läuft dann so dass zB. jeder neue Punkt im Fitparameter-Raum seinen eigenen Thread bekommt und die Threads möglichst lange getrennt laufen.

[BlackJack]

@Gigaz: Das „global interpreter lock“ (GIL) bei CPython ist Dir ein Begriff?

[Gigaz]

@Gigaz: Das „global interpreter lock“ (GIL) bei CPython ist Dir ein Begriff?

Bisher nicht, scheint aber wichtig zu sein Ich habe bisher nie irgendwas in Python parallelisiert, aber ich habe einfach mal angenommen es ist nicht viel schwieriger als openmp für C.

[Gigaz]

So, erste Frage zur Benutzung von Cython

Ich habe eine Klasse namens "Layer" definiert, die alle Eigenschaften einer optischen Lage kennt.

Wenn ich versuche, Cython eine Funktion zu geben, die als Eingabewert ein Layerwertiges Array nimmt, dann kriege ich eine Fehlermeldung.

Code: Alles auswählen

def __Reflectivity(np.ndarray[Layer] HS):
                             ^
------------------------------------------------------------

Pythonreflectivity.pyx:118:30: dtype must be "object", numeric type or a struct

Prinzipiell kann man numpy-arrays aber mit Klassenartigen Objekten füllen. Hat jemand eine Idee wie ich kriege was ich will?

[BlackJack]

@Gigaz: Die Fehlermeldung sagt's ja eigentlich schon: Der Typ darf nur `object` (das wäre dann jedes beliebige Python-Objekt) ein numerischer Typ oder eine Struktur sein (für Record-Arrays). Wenn `Layer` aus Python-Sicht eine Klasse ist, dann wäre `object` richtig, allerdings machen Numpy-Arrays für Python-Objekte nicht so wirklich viel Sinn, weil es dafür nicht wirklich viele Numpy-Operationen gibt, die Vorteile gegenüber Listen bringen.

„Klassenartige Objekte“? Klassen sind in Python auch Objekte. Ich denke mal das meintest Du jetzt nicht‽

Warum hat der Name der Funktion zwei Unterstriche am Anfang? Das macht keinen Sinn.

[Gigaz]

@Gigaz: Die Fehlermeldung sagt's ja eigentlich schon: Der Typ darf nur `object` (das wäre dann jedes beliebige Python-Objekt) ein numerischer Typ oder eine Struktur sein (für Record-Arrays). Wenn `Layer` aus Python-Sicht eine Klasse ist, dann wäre `object` richtig, allerdings machen Numpy-Arrays für Python-Objekte nicht so wirklich viel Sinn, weil es dafür nicht wirklich viele Numpy-Operationen gibt, die Vorteile gegenüber Listen bringen.

„Klassenartige Objekte“? Klassen sind in Python auch Objekte. Ich denke mal das meintest Du jetzt nicht‽

Warum hat der Name der Funktion zwei Unterstriche am Anfang? Das macht keinen Sinn.

Danke fürs helfen.

Wenn "Layer" eine Klasse ist und damit auch ein Objekt, wieso kriege ich dann diese Fehlermeldung? :K
Die zwei Unterstriche haben in dem Kontext nicht viel zu bedeuten, ich teste bloß erstmal ein paar Versionen um die richtige Syntax rauszufinden die ich brauche. Wahrscheinlich teste ich noch eine Weile.

[BlackJack]

@Gigaz: Du bekommst die Fehlermeldung weil `Layer` nicht `object`, ein numerischer Numpy-Typ, oder eine Struktur für Record-Arrays ist. Du kannst kein Numpy-Array erstellen bei dem die Elemente einen festen Python-Datentyp haben, sondern nur ein Array wo die Elemente den Python-Datentyp `object` haben. Und von dem Typ ist letztendlich *jedes* Objekt in Python. Auch Objekte vom Typ `Layer`.

[Gigaz]

Ich bin nicht sicher, ob ich dich völlig verstehe.
Wenn ich statt Layer einfach "object" schreibe kriege ich jedenfalls eine andere Fehlermeldung.

Code: Alles auswählen

def __Reflectivity(np.ndarray[object] HS):
	PyRefl(&HS[0])
	Pythonreflectivity.pyx:121:11: Cannot take address of Python variable

Stünde da statt object "double" dann funktioniert diese Syntax aber.

[BlackJack]

@Gigaz: `double` ist ja aber auch ein Typ bei dem es Sinn macht einen Zeiger zu haben. Python-Objekte sind ja selbst schon Zeiger. Was versprichst Du Dir denn davon Python-Objekte in ein Numpy-Array zu speichern? Statt beispielsweise in einer Liste?

[Gigaz]

@Gigaz: `double` ist ja aber auch ein Typ bei dem es Sinn macht einen Zeiger zu haben. Python-Objekte sind ja selbst schon Zeiger. Was versprichst Du Dir denn davon Python-Objekte in ein Numpy-Array zu speichern? Statt beispielsweise in einer Liste?

Eigentlich verspreche ich mir davon gar nichts, ich habe bloß keinen Beispielcode dafür gefunden bisher
Ich werde es mal versuchen das Klassen-Array stattdessen als Liste zu übergeben. Danke.

[BlackJack]

@Gigaz: „Klassen-Array“ ist so — unpräzise. Was soll das sein? Da Klassen in Python wie schon gesagt Objekte/Werte sind, wäre ein Klassen-Array ein Array in dem Klassen als Elemente stecken, wenn man ein double-Array als ein Array in dem double-Werte als Elemente stecken, begreift. Was in Deinem Array steckt sind Python-Objekte. Genauer braucht man da nicht werden weil es ja keine statische Typisierung gibt.

Wofür hast Du bisher keinen Beispiel-Code gefunden? Ich verstehe ja bisher noch nicht welches Problem Du gerade versuchst zu lösen.

[Gigaz]

Ich habe fast nur autodidaktisch Programmieren gelernt, daher kann ich das mit dem Unpräzise leider nicht so ohne weiteres ändern.

Es gibt diverse Codeschnipsel für Cython wenn man googelt, aber je komplexer das Problem desto weniger findet man.
Ich hab jetzt eine detailliertere Anfrage auf python-forum.org gestellt mit einem Beispielcode:

http://python-forum.org/viewtopic.php?f=6&t=19689

Vielleicht finde ich auch morgen selbst die Lösung.

[BlackJack]

@Gigaz: Wenn Du ein Array oder eine Liste mit Elementen vom Typ `Myclass` hast, und die ein Attribut `ptr` haben, welches ein Zeiger auf einen Wert vom C++-Typ `_MyClass` ist, und Du dann eine C++-Funktion mit der Signatur ``int Take_A_Class_Pointer(_Myclass *A);`` aufrufen möchtest, dann musst Du wohl eher so etwas schreiben:

Code: Alles auswählen

def Take_An_Array_Of_Myclass(np.ndarray[object] input):
    cdef int L
    L = Take_A_Class_Pointer(input[0].ptr)
    return L

Das was Du da hattest (&input[0]) würde intuitiv betrachtet versuchen einen Zeiger auf einen Zeiger auf eine `PyObject`-C-Struktur zu bekommen, was Cython offenbar verweigert. Dein C++-Code könnte damit ja sowieso nichts anfangen.

[Gigaz]

@ Blackjack. Danke nochmal.
Der Vorschlag mit input[0].ptr hat nicht kompiliert, deine Idee mit einem zusätzlichen Pointer war aber richtig. Die funktionierende Syntax ist

Code: Alles auswählen

    cdef void *ptr
    ptr = <void *>input

    L=Take_A_Class_Pointer(ptr)

Und die C++-Funktion packt dann mit

Code: Alles auswählen

int Take_A_Class_Pointer(void *arg){
    _Myclass *A=(_Myclass *)arg;

    return 0;

}

wieder aus.