Python Monte Carlo Optimisierung

Keine Ahnung, wie gut Threading in Python so ist... aber neben der Optimierung der Zufallsfunktionen wäre mein nächster Schritt definitiv die Parallelisierung der Schleife. Sehe bei deinen ganzen Variablen nur so mittelmäßig durch, aber scheinbar schreibst du "global" nur in die Matrizen ab und an? Das kostet ja nicht viel Zeit, also kann man das fein lock {} -en.
 
Hallo,

ich bin mir nicht sicher, ob die Parallelisierung so einfach durchzuführen ist. Beim Ising-Model entscheiden ja die umliegenden Spins ob ein bestimmter Spin geflippt wird. Wenn man jetzt das komplette Gitter in 4 Teilgitter unterteilt und die parallel rechnet, braucht man zumindest für die Spins im Randbereich die Information von den anderen Teilgittern. Machbar ist das schon, aber da muss man vermutlich schon relativ viel Arbeit reinstecken.

Ich würde erstmal versuchen alles aus der Schleife rauszuholen, ohne an Parallelisierung zu denken..
 
@kuddlmuddl:
ich glaube nicht das sowas ausreichen würde. meine wahrscheinlichkeiten ändern sich ja mit jedem teilchen/spin der auf dem platz ist, also kann ich von vornerein keine ausrechnen.
was ich mir überlegt hatte wäre ein array für den zufallswurf von random.random zu machen den ich momentan 2 mal benötige, was ja genau die vektorisierung sein sollte die numpy bevorzugt. falls ich aber wirklich die simulation mit ca 100mrd. durchläufen machen muss, wäre das array größer als mein RAM... und ich bräuchte 2 davon :D

@Enurian:
das threading bzw multiprocessing habe ich mir auch angeschaut, bin mir aber nicht sicher ob ich das hier realisieren kann. wenn zb ein platz 1 von prozess 1 und platz 2 von prozess ausgewählt wird und ein teilchen von 1 nach rechts auf 2 springt, ändern sich auch die warhscheinlichkeiten mit der ein teilchen von platz 2 eine aktion macht.

@stwe:
fast :freaky:
ich habe hier ein aktives ising modell, da sind die spinflips nicht von den nachbarn sondern von der masse/magnetisierung am jeweiligen platz abhängig.


habe gestern das programm mal durchlaufen lassen und es hat ca. 6h gedauert (~600mio durchläufe). in einem beispiel das mir vorliegt wurden 3,6 billionen durchläufe gemacht; bei mir würde das ca 4 jahre dauern.


außerdem habe ich heute mal Cython installiert und eine setup.py datei erstellt mit

Code:
import pyximport; pyximport.install()
import Bachelorarbeit

wenn ich das laufen lasse erstellt mir cython in einem seperaten ordner die datei Bachelorarbeit.cpython-36 was den typ Compiled Python File hat. dieses läuft jetzt auch seit heute morgen und ich hoffe das es etwas schneller ist, bis jetzt scheint das aber nicht der fall zu sein.

außerdem habe ich noch gesehen das man den code verbessern kann indem man feste datentypen definiert.
 
Erstmal vorneweg: ich habe auch, in der guten alten Zeit, vor vielleicht 20 Jahren, mit Ising-Systemen herumgespielt. Allerdings "normale" Spin-Flip-Isingmodelle und Potts-Modelle.

Achso, noch zwei Vorbemerkungen, die Du nicht sooo ernst zu nehmen brauchst, wie sie gemeint sind:
  • die Überschrift müßte heißen: "Optimierung eines Monte-Carlo-Programms in Python". "Optimisierung"gibt es nicht und "Monte-Carlo" ist ein Eigenname für ein statistisches Verfahren.
  • Eine Art "Programmschema" in Python zu Lehrzwecken zu basteln halte ich für eine ausgezeichnete Idee, sich dann damit "zu beschäftigen", das "zu optimieren" um es produktiv zu nutzen ist (imho) reine Zeitverschwendung - es sei denn, Du möchtest etwas über Python lernen und die (nebensächliche) Simulation als Benchmark verwenden.

Battletoad schrieb:
ich glaube nicht das sowas ausreichen würde. meine wahrscheinlichkeiten ändern sich ja mit jedem teilchen/spin der auf dem platz ist, also kann ich von vornerein keine ausrechnen.

beim normalen Ising geht das gut, da Du die "Muster" leicht vorher identifizieren kannst, bei dem andere Spins (-1/+1) um den aktuellen Spin herum angeordnet sind. Dann brauchst Du eine vorher berechnete Liste mit den 2^4 Werten der Umgebung und sparst die exp-Funktion. "Aktives Ising" habe ich noch nicht gerechnet, entweder Du erklärst es mir hier kurz oder ich lese noch mal nach.

was ich mir überlegt hatte wäre ein array für den zufallswurf von random.random zu machen den ich momentan 2 mal benötige, was ja genau die vektorisierung sein sollte die numpy bevorzugt. falls ich aber wirklich die simulation mit ca 100mrd. durchläufen machen muss, wäre das array größer als mein RAM... und ich bräuchte 2 davon

Vorsicht! Jedes statistische Verfahren, und so auch die Ising-Modelle, sind immer sowohl von der Genauigkeit als auch von der Geschwindigkeit wesentlich von der verwendeten Methode der Pseudozufallszahlenerzeugung abhängig. Da kann man sooo viel falsch machen und erhält manchmal ganz lustige Ergebnisse. Been there, done that ;) Hier musst Du viel Gehirnschmalz reinstecken. Du kannst nicht einfach irgend ein random() von irgendeiner Bibliothek verwenden, es sei denn, Du kennst die statistischen Eigenschaften der Zahlenfolge ganz genau (Triplett-Korrelation usw.).

das threading bzw multiprocessing habe ich mir auch angeschaut, bin mir aber nicht sicher ob ich das hier realisieren kann. wenn zb ein platz 1 von prozess 1 und platz 2 von prozess ausgewählt wird und ein teilchen von 1 nach rechts auf 2 springt, ändern sich auch die warhscheinlichkeiten mit der ein teilchen von platz 2 eine aktion macht.
Ising-Systeme lassen sich parallelisieren, indem Du die Flächen in genug Unterflächen aufteilst und dann parallel immer in zwei Unterflächen, die sich nicht berühren, den Spin veränderst. Allerdings ist das in der Literatur seit 20 Jahren beschrieben und ein Literaturstudium sollte für einen Masterstudenten nichts unbekanntes mehr sein ;)

habe gestern das programm mal durchlaufen lassen und es hat ca. 6h gedauert (~600mio durchläufe). in einem beispiel das mir vorliegt wurden 3,6 billionen durchläufe gemacht; bei mir würde das ca 4 jahre dauern.

Unnötige, wenn nicht vergebliche Mühe. Das bringt meiner Meinung nach nichts. Selbst wenn Du es schaffst, mit deiner Schleife nahezu komplett in numpy-Arrays zu bleiben, wird der Code hinterher deutlich komplizierter aussehen, als wenn Du ein einfaches C-Programm schreibst. Es ist übrigens eine wunderbare Sache, wenn man C (oder einfache C++-Konzepte) anhand eines Algorithmus lernt, den man schon gut kennt. Du wirst als Physiker sowieso kaum um C/C++ herumkommen, also könntest Du das nutzen.

Aber: zuerst studierst Du die statistischen Eigenschaften und die relative Performance von Pseudozufallsalgorithmen - in C!

Hier mal was zum kauen für Dich. Was hältst Du von dem:
Code:
/*
 SEMI-RANDOM ALGORITHM "A5"
 Returns uniform random numbers in the
 range if 0. to 1. (excluding 1.)
 In C++ source, simply do:

  main() {
   ...
   ...
   double next_uniform_number = fastRNG::draw();
   ...
   ...
  }
*/
 #if defined (__cplusplus)
 #include <cmath>

 class fastRNG {
   static double a, b, c;
 public:
   static double draw(void) {   // static access
      double r = a + b + c;
      while(r >= 1.0) r -= 1.0; // r = fmod(a+b+c, 1.0);
      return a=b, b=c, c=r;
   }
 };
  double fastRNG::a=0.3022777;
 double fastRNG::b=0.2371781;
 double fastRNG::c=0.5263637;
 
#else
  #include <math.h>
 double fastRNG (void)
{
 static double a = 0.3022777, b = 0.2371781, c = 0.5263637;
 double r = a + b + c;
 while(r >= 1.0) r -= 1.0;
 return a=b, b=c, c=r;
}
 #endif
Kannst Du als header in einen C-Quelltext einbinden.
 
der name passt ja mal garnicht zum text :p

aktives ising modell unterscheidet sich dadrin, dass die teilchen aufeinandre sitzen und jetzt nach dichte/magnetisierung an diesem platz ein spinflip machen, oder nach oben/unten/rechts/links springen. damit ergibt sich nach einiger zeit das sich alle teilchen in der mitte des feldes ansammeln.

zu den zufallsgeneratoren: die random-fkt in numpy benutzt den mersenne twister. den habe ich früher mal genau angeschaut (und einen eigenen geschrieben) inklusive chi-quadrat text etc.

python habe ich mir deshalb ausgesucht, da ich später in die astrophysik will und das dort oft verwendet wird. außerdem macht es mir mehr spaß als C/C++ :D. mithilfe von jemanden von stackoverflow konnte die laufzeit um den factor 435 verringert werden. dazu wurde im numba kompilliert was sich nicht mehr groß von der laufzeit in C unterscheiden sollte (google numba vs C).
 
kannst du bitte die finale Version des Source posten, würde gerne mal einen Blick darauf werfen. Nehme mal an, dass einige Decorators verwendet worden.
 
Zurück
Oben