Python Monte Carlo Optimisierung

Battletoad · 27. Juli 2017

hallo,

ich habe für meine bachelor arbeit ein aktives-2D-ising-modell programmiert und möchte gerne den code etwas optimieren. die ist mein erstes größeres python projekt und google wurde auch benutzt, nur ist vieles davon nicht so verständlich für mich als anfänger.

zu anfang erstelle ich zwei matrizen: in einer stehen die teilchenzahlen pro platz und in der anderen die magnetisierung (summe von spin 1 und (-1) teilchen).
danach wird ein platz und ein teilchen ausgewählt und dann eine aktion ausgewürfelt (spinflip 1 <-> (-1), sprung nach oben, unten, links oder rechts, nichts tun). am ende werden die teilchendichte/magnetisierung der plätze angepasst.

ich habe die platzanpassungen jetzt direkt in die funktionen geschrieben damit das programm, sobald es eine aktion ausgeführt hat, mit dem nächsten schritt anfangen kann; damit sollte es schneller durchlaufen mMn.
vorher wurde die am ende der schleife ausgeführt (jetzt auskommentiert).

da dies als MC-simulation läuft komme ich hier auf mehrere mrd. durchläufe und würde gerne die laufzeit minimieren (momentan 10h+).

wäre toll wenn jemand ein paar vorschläge hätte.

Grimba · 27. Juli 2017

Mal http://cython.org/ probiert?

Battletoad · 27. Juli 2017

hatte davon gelesen.

mit pycharm habe ich cython installiert und oben im code

Code:

import pyximport; pyximport.install()

eingefügt.

langt das jetzt schon um cython zu benutzen?

Grimba · 27. Juli 2017

Battletoad schrieb:
hatte davon gelesen.

mit pycharm habe ich cython installiert und oben im code

Code:

import pyximport; pyximport.install()

eingefügt.

langt das jetzt schon um cython zu benutzen?

Erwischt...

ich hab selber nur ein bißchen gegooglet, mit dem Wissen, dass es doch sowas gab, mit dem man aus Python C-Code machen konnte, um ggf. Laufzeit zu sparen und/oder echte C-libs einzubinden, welche ggf. auch schneller sind.

Richtig auskennen tue ich mich damit nicht. Aber ich dachte bevor du keine Antwort bekommst, und das geht ja in etwa in die Richtung...

Ich kann dir nur empfehlen, dein Anliegen auch auf stackoverflow.com zu posten. Sehr häufig DIE Anlaufstelle

Battletoad · 27. Juli 2017

hatte ich auch vor dort zu posten, aber ob das dort jemand sieht ist was anderes

trotzdem danke

Mr.Seymour Buds · 27. Juli 2017

Die Optimierung solltest eigentlich Du machen. Das ist ja der Sinn der Bachelorarbeit.

Battletoad · 27. Juli 2017

Mr.Seymour Buds schrieb:
Die Optimierung solltest eigentlich Du machen. Das ist ja der Sinn der Bachelorarbeit.

nein, dass ist nicht sinn der bachelorarbeit...

kuddlmuddl · 27. Juli 2017

nein, dass ist nicht sinn der bachelorarbeit...

Sehe ich genau so. Deine Frage hört sich so an, als ob du nicht unbedingt Informatik/Mathematik studierst und da finde ich es absolut legitim eine FERTIGE und laufende Lösung mit Hilfe anderer Leute in der Lautzeit zu verbessern. Es ist natürlich nur nicht Sinn der Sache andere den Kram von 0 auf implementieren zu lassen.

Ich Programmiere auch gerne mal Python aber mache eigtl nur C++. Da ich immer wieder genervt bin, wenn lange Berechnungen ZU lange dauern in Python baue ich es halt in C++ nach und in meinen wenigen Fällen war es dadurch rund 100-10000 mal schneller. Das liegt aber natürlich nicht pauschal an der Sprache sondern auch daran, dass ich C++ besser kann und es dann auch gleich im Hinblick auf Laufzeit 'besser' neu gebaut habe. Falls das (mit Hilfe eines erfahrenen C++ Programmierers?) keine Option ist:

Was in deiner Python-Lösung sicherlich viel Zeit kostet ist der Inhalt von a, denn da sehe ich Strings. In C(++) würde man hier enums verwenden die als integer vom PC um ein vielfaches schneller verarbeitet werden können.
Ich meine insbesondere diese ganzen Zeilen: "if now == 'nope':", denn Stringvergleiche sind sehr teuer.
Entweder guckst du, wie man enum equivalent in python verwendet link1, link2 oder du machst es einfach mit einer Kodierung up=1, down=2, ... und 6 Variablen die am Anfang des scripts stehen.

Der allgemein sinnvolle Ansatz ist aber den PC messen zu lassen, WAS eigentlich genau lange dauert. Der Fachbegriff lautet 'profilen' und hier link1, link2 sind Sachen zum lesen. Beachte aber, dass Code während des profilings NOCH länger läuft. Du solltest also die Variablen (j?) so wählen, dass die Gesamtlaufzeit eher im Bereich von Minuten als 10h liegt bevor du das erste mal profilest.

Die Idee es mit soetwas wie CPython zu beschleunigen ist sicherlich auch sinnvoll zu verfolgen, weil es wohl wenig Aufwand bedeutet und viel bringen kann.

Battletoad · 27. Juli 2017

danke erstmal

ich studiere physik mit nebenfach mathe und die aufgabenstellung ist ein pyhsik. problem mithilfe vom pc zu berechnen. die optimierung wird empfohlen ist aber keineswegs pflicht.

das python an sich langsamer ist als C++ ist klar, durch numpy/scipy und cython lässt sich aber eine ähnliche geschwindigkeit erreichen wie in C++ (wurde mir so gesagt). in python hat man das zeug schneller programmiert, dauert dann aber länger zum rechnen im vgl zu C/C++. aber diese diskussion hat es schon zu genüge gegeben

das mit den string dachte ich mir auch, wollte es noch auf zahlen ändern. über den profiler habe ich auch gelesen, mich aber nocht nicht so richtig damit befasst.
das programm macht auf jeden fall was verlangt wurde und ich habe es auch schon mehrmals laufen lassen.

Mr.Seymour Buds · 27. Juli 2017

Lasse den Code gerade bei mir laufen und mein Taskmanager zeigt mir eine Auslastung von ca. 25% an (mit Ausreißern noch oben bis auf 40%).

Würde sagen, hier kann man schon einiges rausholen an Speedup. Ich bin kein Spezialist für Monte-Carlo Simulationen, aber es sieht so aus, als gibt es Ansätze, sie zu parallelisieren. Da würde ich ansetzen (neben den genannten Vorschlägen). *hust* pyCUDA.

convexus · 27. Juli 2017

es macht Sinn die Code-Stellen zu ermitteln in denen die meiste Zeit verbraucht wird. Dazu kannst du einen Python-Profiler verwenden. Z.B.

Code:

python -m cProfile -s cumtime <dein-skript>

liefert für (j=1):

         70720717 function calls (70706290 primitive calls) in 306.741 seconds

   Ordered by: cumulative time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1   64.947   64.947  306.742  306.742 mc_sim.py:2(<module>)
  3566185  137.431    0.000  192.637    0.000 {method 'choice' of 'mtrand.RandomState' objects}
  3566185   23.001    0.000   47.536    0.000 fromnumeric.py:2388(prod)
  3566185    2.194    0.000   24.535    0.000 _methods.py:34(_prod)
  3572657   22.370    0.000   22.370    0.000 {method 'reduce' of 'numpy.ufunc' objects}
  1783166    5.747    0.000   12.002    0.000 {map}
  1815250   10.350    0.000   10.350    0.000 {method 'randint' of 'mtrand.RandomState' objects}
  1783114    9.173    0.000    9.173    0.000 {zip}
  3566184    6.465    0.000    7.670    0.000 getlimits.py:96(__new__)
  3566255    6.266    0.000    6.266    0.000 {numpy.core.multiarray.arange}
  1783092    1.817    0.000    6.207    0.000 numeric.py:1927(array_str)
 26939057    5.682    0.000    5.682    0.000 {method 'item' of 'numpy.ndarray' objects}
  1783092    4.078    0.000    4.389    0.000 arrayprint.py:340(array2string)
  3570446    1.190    0.000    1.193    0.000 {method 'get' of 'dict' objects}

Die meiste Zeit wird demnach in random.choice verbraucht! Evtl. die p-Arrays mal außerhalb der Schleife definieren und nur die geänderte Elemente aktualisieren.

Battletoad · 28. Juli 2017

Mr.Seymour Buds schrieb:
Lasse den Code gerade bei mir laufen und mein Taskmanager zeigt mir eine Auslastung von ca. 25% an (mit Ausreißern noch oben bis auf 40%).

auf meinem heim pc schwankt es auch zw. 25% und 100% (in versch. threads). auf dem uni pc wird jedoch ein kern vom i7 2600 voll ausgelastet. ka woran das liegt.

convexus schrieb:
Die meiste Zeit wird demnach in random.choice verbraucht! Evtl. die p-Arrays mal außerhalb der Schleife definieren und nur die geänderte Elemente aktualisieren.

ich glaube das ist auch was mein prof gesagt hat. wie genau ich das machen soll ist mir aber noch nicht klar. werde es mal versuchen und später rückmeldung geben.

Mr.Seymour Buds · 28. Juli 2017

Why is random.choice so slow?

convexus · 28. Juli 2017

ersetze mal folgenden Source

Code:

spin = np.random.choice(s, 1, p=[1 - num / rgrid.item(x), num / rgrid.item(x)])  # wahrscheinlichkeit teilchen mit spin 1 auszuwählen P(+) = n(+)/rgrid

durch

Code:

spin = -1.0 if np.random.random() < 1.0-num / rgrid.item(x) else 1.0

Anstelle von int(spin) muss du dann spin schreiben. Das reduziert die Laufzeit um ein 1/3!

Battletoad · 28. Juli 2017

die laufzeit reduzierung ist sehr gut.

warum genau das jetzt schneller ist verstehe ich nicht so ganz. lese mich momentan durch die vektorisierung von python und versuche zu verstehen wie ich das bei mir anwenden kann.

edit:
wenn ich das also richtig verstanden habe ist random.choice gut wenn viele ereignisse auf einmal ausgewählt werden sollen (siehe link von Mr.Seymour Buds und hier der 1. link).
das geht bei mir aber nicht da sich die warscheinlichkeiten in jedem schritt ändern und neu berechnet werden müssen. für nur einen wurf ist die methode von convexus besser.

convexus · 28. Juli 2017

Probier mal diese Ersatzfunktion anstelle von np.random.choice

Code:

def mr_choice(a, _, p):
    r = np.random.random()
    i=0
    while r>p[i] and i<len(a)-1:
        r -= p[i]
        i += 1
    return a[i]

Battletoad · 28. Juli 2017

ich habe jetzt

Code:

W = np.exp(- spin * b * mgrid.item(x) / rgrid.item(x))
now = np.random.choice(a, 1, p=[W * dt, dt, dt, (1 - spin * e) * dt, (1 + spin * e) * dt, 1 - (4 + W) * dt])

durch

Code:

p = np.random.random()
p1 = np.exp(- spin * b * mgrid.item(x) / rgrid.item(x)) * dt
p2 = (1 - spin * e) * dt
p3 = (1 + spin * e) * dt

if p < p1:
    now = 0
elif p < p1 + dt:
    now = 1
elif p < p1 + 2 * dt:
    now = 2
elif p < p1 + 2 * dt + p2:
    now = 3
elif p < p1 + 2 * dt + p2 + p3:
     now = 4
else:
     now = 5

ersetzt und konnte dadurch die laufzeit um ~50-60% reduzieren.

@convexus:
habe deine funktion außerhalb der schleife definiert und den code abgeändert zu

Code:

W = np.exp(- spin * b * mgrid.item(x) / rgrid.item(x))
p = [W * dt, dt, dt, (1 - spin * e) * dt, (1 + spin * e) * dt, 1 - (4 + W) * dt]
mr_choice(a, _, p)

aber er gibt mir ne fehlermeldung aus das der underscore nicht definiert ist, was doch eigentlich garkeinen sinn ergibt, da der doch eine wegwerf-variable darstellen soll.

convexus · 28. Juli 2017

anstelle von mr_choice(a, _, p) schreibe mr_choice(a, 1, p). Der Underscore _ ist ein unbenannter Parameter.

Battletoad · 28. Juli 2017

ich war zu doof dafür, klappt jetzt

cProfile sagt das mr_choice minimal langsamer als die if-elif version ist.

habe außerdem deinen code von oben in

Code:

spin = 1.0 if np.random.random() < num / rgrid.item(x) else -1.0

abgeändert. macht es ca 0,1s schneller und sollte nichts am ergebnis ändern.

kuddlmuddl · 28. Juli 2017

Random würde ich im gesamten Quellcode zur Laufzeit des Algos niemals verwenden..
Mach doch einfach ne statische Liste von Zufallszahlen beim Start des Programms (zB 2000 Stück) und nimm reihum immer eine davon. Solange es genug sind und diese deine gewünschte Verteilung haben reicht das auch.

Python Monte Carlo Optimisierung

Ensign

Captain

Ensign

Captain

Ensign

Commodore

Ensign

Commander

Ensign

Commodore

Ensign

Ensign

Commodore

Ensign

Ensign

Ensign

Ensign

Ensign

Ensign

Commander

Ähnliche Themen