Zufallsgenerator Mersenne-Twister und jumpahead()

[wollmich]

Habe ein Problem mit dem Zufallsgenerator von numpy.random.

Der von Python eingebaute Zufallsgenerator random kennt die Funktion jumpahead(n). Mit Hilfer dieser Funktion ist es sehr einfach mit Multi-Processing Zufallszahlen zu erzeugen.

Beispiel:

Code: Alles auswählen

# Prozess 0
import random
random.seed()
random_state = random.getstate()
# Nun random_state an Prozesse verteilen und je 100000 Zahlen generieren

# Prozess 1
random.setstate(random_state)
random.jumpahead(0)
return [random.random() for i in xrange(100000)]

# Prozess 2
random.setstate(random_state)
random.jumpahead(100000)
return [random.random() for i in xrange(100000)]

# Prozess 3
random.setstate(random_state)
random.jumpahead(200000)
return [random.random() for i in xrange(100000)]

# und so weiter

Dies funktioniert, braucht aber lange um 100000 oder mehr Zahlen zu würfeln, da jeder Prozess einen For Loop hat.

Mit NumPy kann ich Bereiche würfeln (numpy.random.rand(1000)), viel schneller als der eingebaute Python Zufallsgenerator. Ich vermisse aber die jumpahead Funktion bei numpy.random.

Ein möglicher Workaround wäre:

Code: Alles auswählen

def jumpahead(n):
    x = 100000
    # For Loop über Ganzahldivision von n/x wegen Speicherverbrauch
    for i in xrange(n/x):
        numpy.random.rand(x)
    # Nun noch der Rest
    numpy.random.rand(n%x)

Dieser ist aber langsam sobald n >> 100000

Gibt es einen besseren Workaround, irgendwelche Alternativen?

Bin für jeden Hinweis dankbar.

[CM]

Hm, das ist zwar nicht ganz dasselbe, aber: Was spricht gegen einen wiederholten Aufruf von numpy.random.seed(100000) (oder anderen, ggf. "zufälligen" Zahlen, gewonnen womöglich mit dem random-Modul)? Oder reicht es womöglich auch bei den versch. Prozessen den state mit numpy.random.get/set_state zu übergeben?

Ansonsten ist dies die jumpahead()-Funktion aus dem Python-random-Modul:

Code: Alles auswählen

def jumpahead(self, n):
        """Act as if n calls to random() were made, but quickly.

        n is an int, greater than or equal to 0.

        Example use:  If you have 2 threads and know that each will
        consume no more than a million random numbers, create two Random
        objects r1 and r2, then do
            r2.setstate(r1.getstate())
            r2.jumpahead(1000000)
        Then r1 and r2 will use guaranteed-disjoint segments of the full
        period.
        """

        if not n >= 0:
            raise ValueError("n must be >= 0")
        x, y, z = self._seed
        x = int(x * pow(171, n, 30269)) % 30269
        y = int(y * pow(172, n, 30307)) % 30307
        z = int(z * pow(170, n, 30323)) % 30323
        self._seed = x, y, z

Nicht schwer so etwas in den eigenen Code einzubauen.

Du hast aber auch Probleme ...

Gruß,
Christian

[mitsuhiko]

Ich würde jedem Prozess einen eigenen Seed verpassen.

[wollmich]

Ist leider beides nicht die Lösung.

Details zu meinem Problem:
Ich möchte zum Beispiel 1'000'000'000 Zahlen würfeln. Dafür merke ich mir den InitRandomState mit getstate(). Ich würfle 10'000 Blöcke à 100'000 Zahlen wobei die Blöcke auf verschiedene Prozesse verteilt werden sollen. Im Speicher habe ich nie mehr als 8 Blöcke gleichzeitig.
Ich muss in der Lage sein immer wieder die gleichen Zahlen zu würfeln.
Beispiel Objekt 'a' hat den InitRandomState 'IRS_a' und Objekt 'b' hat den InitRandomState 'IRS_b'. Nun muss ich in der Lage sein immer wieder die gleichen Zahlen von zum Beispiel Block 135 von Objekt 'a' zu würfeln. Ich kann mir unmöglich alle gewürfelten Zahlen merken (Speicherplatz (mehrere 100GB RAM)). Möchte mir aber auch nicht für jeden Block eines Objektes den RandomState merken, sondern nur den von Block 0. Würde ja alles so funktionieren, wenn es jumpahead für numpy.random gäbe.

Die von dir abgedruckte jumpahead()-Funktion aus dem Python-random-Modul gehört zur Subklasse WichmannHill, der eine für mich zu kleine Periode hat. Der normalle Zufallsgernator von random und numpy.random hat kein self._seed. Sein RadomState sieht komplexer aus.

Werd mich behühen auch mal ein einfache Frage zu stellen

[BlackJack]

Ich denke er möchte das sich die Zahlenfolgen nicht überlappen, das kann man mit zufälligen "seeds" nicht garantieren.

Das `jumpahead()` aus `random` muss nicht mit `numpy` funktionieren. Die konstanten Zahlen sind ja nicht zufällig, sondern passen genau zu der Implementierung im `random`-Modul in der Standardbibliothek.

[CM]

Das `jumpahead()` aus `random` muss nicht mit `numpy` funktionieren. Die konstanten Zahlen sind ja nicht zufällig, sondern passen genau zu der Implementierung im `random`-Modul in der Standardbibliothek.

Ja, aber es ist eine gutes Beispiel. numpy.random verwendet einen pyrex-wrapper um eine C++-Klasse. Rel. einfach zu lesen. Also, wenn die Performance so wichtig ist, würde ich mtrand.pyx editieren und jumpahead implementieren. Wenn ich das richtig verstehe, kann man hierzu den relevanten Teil aus numpys seed()-Methode und aus Python jumpahead() abschreiben.

Alternativ: Robert Kern anschreiben und fragen, ob er bereit ist eine solche Methode zu schreiben. Er schreibt regelmäßig auf der numpy-Liste und scheint sehr hilfsbereit.

Mich juckt es in den Fingern, aber ich muß bis Montag noch was wegschaffen. Ansonsten will ich mich gerne in der nächsten Woche dran versuchen.

Gruß,
Christian