samedi 20 juillet 2019

3D plotting in python - random (prng - LCG)

I am trying to create a 3d chart in python for the following linear congruential generator:

import numpy as np

class LCG(object):

    UZERO: np.uint32 = np.uint32(0)
    UONE : np.uint32 = np.uint32(1)

    def __init__(self, seed: np.uint32, a: np.uint32, c: np.uint32) -> None:
        self._seed: np.uint32 = np.uint32(seed)
        self._a   : np.uint32 = np.uint32(a)
        self._c   : np.uint32 = np.uint32(c)

    def next(self) -> np.uint32:
        self._seed = self._a * self._seed + self._c
        return self._seed

    def seed(self) -> np.uint32:
        return self._seed

    def set_seed(self, seed: np.uint32) -> np.uint32:
        self._seed = seed

    def skip(self, ns: np.int32) -> None:
        """
        Signed argument - skip forward as well as backward

        The algorithm here to determine the parameters used to skip ahead is
        described in the paper F. Brown, "Random Number Generation with Arbitrary Stride,"
        Trans. Am. Nucl. Soc. (Nov. 1994). This algorithm is able to skip ahead in
        O(log2(N)) operations instead of O(N). It computes parameters
        A and C which can then be used to find x_N = A*x_0 + C mod 2^M.
        """

        nskip: np.uint32 = np.uint32(ns)

        a: np.uint32 = self._a
        c: np.uint32 = self._c

        a_next: np.uint32 = LCG.UONE
        c_next: np.uint32 = LCG.UZERO

        while nskip > LCG.UZERO:
            if (nskip & LCG.UONE) != LCG.UZERO:
                a_next = a_next * a
                c_next = c_next * a + c

            c = (a + LCG.UONE) * c
            a = a * a

            nskip = nskip >> LCG.UONE

        self._seed = a_next * self._seed + c_next


#%%
np.seterr(over='ignore')

a = np.uint32(1664525)
c = np.uint32(1013904223)
seed = np.uint32(2**32)

rng = LCG(seed, a, c)
q = [rng.next() for _ in range(0, 10000)]
# print(q)

I would like the chart to change every x seconds for a given n - like this (George-Marsaglia).

It's not necessary that it changes over time.

I tried to adapt to the following code butI don't know how to do it with my generator:

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
import numpy as np
import os

def lcg (X, a, c, m):
    return (a * X + c) % m;

X = []
Y = []
Z = []

n = int(input("N : "))

prev = 0
for i in range(n):
    prev = lcg(prev,43,5,4096)
    if i % 3 == 0:
        X.append(prev)
    elif i % 3 == 1:
        Y.append(prev)
    else:
        Z.append(prev)

X = np.array(X)
Y = np.array(Y)
Z = np.array(Z)

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.set_aspect('equal')
ax.set_title('N = ' + str(n))

scat = ax.scatter(X, Y, Z)

max_range = np.array([X.max()-X.min(), Y.max()-Y.min(), Z.max()-Z.min()]).max() / 2.0

mid_x = (X.max()+X.min()) * 0.5
mid_y = (Y.max()+Y.min()) * 0.5
mid_z = (Z.max()+Z.min()) * 0.5
ax.set_xlim(mid_x - max_range, mid_x + max_range)
ax.set_ylim(mid_y - max_range, mid_y + max_range)
ax.set_zlim(mid_z - max_range, mid_z + max_range)

i = 0
while os.path.exists('fig%s.png' % i):
    i += 1

plt.savefig('fig' + str(i) + '.png')
plt.show()

It is possible to adapt the following code to my generator?




Aucun commentaire:

Enregistrer un commentaire