add Ziggurat Gaussian random sequence generator

Author: thrasibule

quantlib/math/randomnumbers/_mt19937uniformrng.pxd

@@ -0,0 +1,3 @@
+cdef extern from 'ql/math/randomnumbers/mt19937uniformrng.hpp' namespace 'QuantLib' nogil:
+    cdef cppclass MersenneTwisterUniformRng:
+        pass

quantlib/math/randomnumbers/_randomsequencegenerator.pxd

@@ -0,0 +1,13 @@
+from libcpp.vector cimport vector
+from quantlib.types cimport Real, Size
+from quantlib._stochastic_process cimport StochasticProcess
+from quantlib.methods.montecarlo._sample cimport Sample
+
+cdef extern from 'ql/math/randomnumbers/randomsequencegenerator.hpp' namespace 'QuantLib' nogil:
+    cdef cppclass RandomSequenceGenerator[RNG]:
+        ctypedef Sample[vector[Real]] sample_type
+        RandomSequenceGenerator(Size dimensionality,
+                                const RNG& rgn)
+        const sample_type& nextSequence() const
+        const sample_type& lastSequence() const
+        Size dimension() const

quantlib/math/randomnumbers/_rngtraits.pxd

@@ -1,12 +1,25 @@
 from quantlib.types cimport BigNatural, Size
 from ._inverse_cumulative_rsg cimport InverseCumulativeRsg
 from ._sobol_rsg cimport SobolRsg
+from ._randomsequencegenerator cimport RandomSequenceGenerator
+from ._mt19937uniformrng cimport MersenneTwisterUniformRng
+from ._zigguratgaussianrng cimport ZigguratGaussianRng
+from ._xoshiro256starstaruniformrng cimport Xoshiro256StarStarUniformRng
 
 cdef extern from 'ql/math/distributions/normaldistribution.hpp' namespace 'QuantLib' nogil:
     cdef cppclass InverseCumulativeNormal:
         pass
 
 cdef extern from 'ql/math/randomnumbers/rngtraits.hpp' namespace 'QuantLib' nogil:
+    cdef cppclass GenericPseudoRandom[URNG, IC]:
+        ctypedef RandomSequenceGenerator[URNG] ursg_type
+        ctypedef InverseCumulativeRsg[ursg_type, IC] rsg_type
+
+        @staticmethod
+        rsg_type make_sequence_generator(Size dimension,
+                                         BigNatural seed)
+    ctypedef GenericPseudoRandom[MersenneTwisterUniformRng, InverseCumulativeNormal] PseudoRandom
+
     cdef cppclass GenericLowDiscrepancy[URSG, IC]:
         ctypedef InverseCumulativeRsg[URSG, IC] rsg_type
 
@@ -15,3 +28,5 @@ cdef extern from 'ql/math/randomnumbers/rngtraits.hpp' namespace 'QuantLib' nogi
                                          BigNatural seed)
 
     ctypedef GenericLowDiscrepancy[SobolRsg, InverseCumulativeNormal] LowDiscrepancy
+
+ctypedef RandomSequenceGenerator[ZigguratGaussianRng[Xoshiro256StarStarUniformRng]] ZigguratPseudoRandom

quantlib/math/randomnumbers/_xoshiro256starstaruniformrng.pxd

@@ -0,0 +1,11 @@
+from libc.stdint cimport uint64_t
+from quantlib.types cimport Real
+from quantlib.methods.montecarlo._sample cimport Sample
+
+cdef extern from 'ql/math/randomnumbers/xoshiro256starstaruniformrng.hpp' namespace 'QuantLib' nogil:
+    cdef cppclass Xoshiro256StarStarUniformRng:
+        ctypedef Sample[Real] sample_type
+        Xoshiro256StarStarUniformRng(uint64_t seed) #=0
+        sample_type next() const
+        Real nextReal() const
+        uint64_t nextInt64()

quantlib/math/randomnumbers/_zigguratgaussianrng.pxd

@@ -0,0 +1,9 @@
+from quantlib.types cimport Real
+from quantlib.methods.montecarlo._sample cimport Sample
+
+cdef extern from 'ql/math/randomnumbers/zigguratgaussianrng.hpp' namespace 'QuantLib' nogil:
+    cdef cppclass ZigguratGaussianRng[RNG]:
+        ctypedef Sample[Real] sample_type
+        ZigguratGaussianRng(const RNG& uint64Generator)
+        sample_type next() const
+        Real nextReal() const

quantlib/math/randomnumbers/rngtraits.pxd

@@ -1,7 +1,16 @@
-from . cimport _rngtraits
 from ._inverse_cumulative_rsg cimport InverseCumulativeRsg
+from ._randomsequencegenerator cimport RandomSequenceGenerator
 from ._sobol_rsg cimport SobolRsg
 from ._rngtraits cimport InverseCumulativeNormal
+from ._mt19937uniformrng cimport MersenneTwisterUniformRng
+from ._zigguratgaussianrng cimport ZigguratGaussianRng
+from ._xoshiro256starstaruniformrng cimport Xoshiro256StarStarUniformRng
 
 cdef class LowDiscrepancy:
     cdef InverseCumulativeRsg[SobolRsg, InverseCumulativeNormal]* _thisptr
+
+cdef class PseudoRandom:
+    cdef InverseCumulativeRsg[RandomSequenceGenerator[MersenneTwisterUniformRng], InverseCumulativeNormal]* _thisptr
+
+cdef class FastPseudoRandom:
+    cdef RandomSequenceGenerator[ZigguratGaussianRng[Xoshiro256StarStarUniformRng]]* _thisptr

quantlib/math/randomnumbers/rngtraits.pyx

@@ -1,15 +1,51 @@
-from quantlib.types cimport BigNatural, Size
+from libc.string cimport memcpy
+from quantlib.types cimport BigNatural, Real, Size
+from cython.operator cimport dereference as deref
+from . cimport _rngtraits
 cimport cython
 cimport numpy as np
 from cython.operator cimport dereference as deref
 from libcpp.vector cimport vector
+from quantlib.methods.montecarlo._sample cimport Sample
 np.import_array()
 
 cdef class LowDiscrepancy:
 
     def __init__(self, Size dimension, BigNatural seed=0):
         self._thisptr = new InverseCumulativeRsg[SobolRsg, InverseCumulativeNormal](
-                _rngtraits.LowDiscrepancy.make_sequence_generator(dimension, seed))
+                _rngtraits.LowDiscrepancy.make_sequence_generator(dimension, seed)
+        )
+
+    def __dealloc__(self):
+        if self._thisptr is not NULL:
+            del self._thisptr
+
+    def __iter__(self):
+        return self
+
+    @cython.boundscheck(False)
+    def __next__(self):
+        cdef np.npy_intp[1] dims
+        dims[0] = self._thisptr.dimension()
+        cdef arr = np.PyArray_SimpleNew(1, &dims[0], np.NPY_DOUBLE)
+        cdef Sample[vector[Real]]* s = <Sample[vector[Real]]*>&self._thisptr.nextSequence()
+        memcpy(np.PyArray_DATA(arr), s.value.data(), self._thisptr.dimension() * sizeof(Real))
+        return (s.weight, arr)
+
+    @property
+    def dimension(self):
+        return self._thisptr.dimension()
+
+    def last_sequence(self):
+        return self._thisptr.lastSequence().value
+
+
+cdef class PseudoRandom:
+
+    def __init__(self, Size dimension, BigNatural seed=0):
+        self._thisptr = new _rngtraits.PseudoRandom.rsg_type(
+            _rngtraits.PseudoRandom.make_sequence_generator(dimension, seed)
+        )
 
     def __dealloc__(self):
         if self._thisptr is not NULL:
@@ -23,16 +59,43 @@ cdef class LowDiscrepancy:
         cdef np.npy_intp[1] dims
         dims[0] = self._thisptr.dimension()
         cdef arr = np.PyArray_SimpleNew(1, &dims[0], np.NPY_DOUBLE)
-        cdef double[:] r = arr
-        cdef size_t i
-        cdef vector[double] s = self._thisptr.nextSequence().value
-        for i in range(dims[0]):
-            r[i] = s[i]
-        return arr
+        cdef Sample[vector[Real]]* s = <Sample[vector[Real]]*>&self._thisptr.nextSequence()
+        memcpy(np.PyArray_DATA(arr), s.value.data(), self._thisptr.dimension() * sizeof(Real))
+        return (s.weight, arr)
 
     @property
     def dimension(self):
         return self._thisptr.dimension()
 
     def last_sequence(self):
         return self._thisptr.lastSequence().value
+
+
+cdef class FastPseudoRandom:
+
+    def __init__(self, Size dimension, BigNatural seed=0):
+        cdef Xoshiro256StarStarUniformRng* uniform_random = new Xoshiro256StarStarUniformRng(seed)
+        cdef ZigguratGaussianRng[Xoshiro256StarStarUniformRng]* rng = new ZigguratGaussianRng[Xoshiro256StarStarUniformRng](deref(uniform_random))
+        self._thisptr = new RandomSequenceGenerator[ZigguratGaussianRng[Xoshiro256StarStarUniformRng]](dimension, deref(rng))
+        del uniform_random
+        del rng
+
+    def __dealloc__(self):
+        if self._thisptr is not NULL:
+            del self._thisptr
+
+    def __iter__(self):
+        return self
+
+    @cython.boundscheck(False)
+    def __next__(self):
+        cdef np.npy_intp[1] dims
+        dims[0] = self._thisptr.dimension()
+        cdef arr = np.PyArray_SimpleNew(1, &dims[0], np.NPY_DOUBLE)
+        cdef Sample[vector[Real]]* s = <Sample[vector[Real]]*>&self._thisptr.nextSequence()
+        memcpy(np.PyArray_DATA(arr), s.value.data(), self._thisptr.dimension() * sizeof(Real))
+        return (s.weight, arr)
+
+    @property
+    def dimension(self):
+        return self._thisptr.dimension()

test/test_cms_spread.py

@@ -168,7 +168,7 @@ def mc_reference_value(cpn1, cpn2, vol, correlation):
         it = 0
         z = np.empty((samples, 2))
         for i in range(samples):
-            z[i] = next(g)
+            _, z[i] = next(g)
         z = z.dot(C.T) + avg
         if vol.volatility_type == VolatilityType.ShiftedLognormal:
             z = (atm_rate + vol_shift) * np.exp(z) - vol_shift

test/test_randomnumbers.py

@@ -39,7 +39,7 @@ def test_mean_variance(self):
         samples = 100000
         X = np.empty((samples, self.g.dimension))
         for i in range(samples):
-            X[i] = next(self.g)
+            _, X[i] = next(self.g)
         C = np.cov(X.T)
         self.assertAlmostEqual(np.linalg.norm(X.mean(axis=0)), 0., 3)
         self.assertAlmostEqual(np.linalg.norm(C - np.eye(10)), 0., 2)