Brownian

quantlib/_time_grid.pxd

@@ -1,7 +1,7 @@
-include 'types.pxi'
+from quantlib.types cimport Size, Time
 from libcpp.vector cimport vector
 
-cdef extern from 'ql/timegrid.hpp' namespace 'QuantLib':
+cdef extern from 'ql/timegrid.hpp' namespace 'QuantLib' nogil:
     cdef cppclass TimeGrid:
         TimeGrid()
         TimeGrid(Time end, Size steps)

quantlib/math/_array.pxd

@@ -7,7 +7,7 @@
  FOR A PARTICULAR PURPOSE.  See the license for more details.
 """
 
-include '../types.pxi'
+from quantlib.types cimport Real, Size
 
 cdef extern from 'ql/math/array.hpp' namespace 'QuantLib':
     cdef cppclass Array:
@@ -17,3 +17,4 @@ cdef extern from 'ql/math/array.hpp' namespace 'QuantLib':
         Real& at(Size i) except +IndexError
         Size size()
         Real& operator[](Size)
+        const Real* begin()

quantlib/math/array.pyx

@@ -7,15 +7,19 @@
  FOR A PARTICULAR PURPOSE.  See the license for more details.
 """
 
-include '../types.pxi'
+from quantlib.types cimport Real, Size
+from cpython.ref cimport Py_INCREF
 from libcpp.utility cimport move
+cimport numpy as np
+
+np.import_array()
 
 cdef class Array:
     """
     1D array for linear algebra
     """
 
-    def __init__(self, Size size, value=None):
+    def __init__(self, Size size=0, value=None):
         if value is None:
             self._thisptr = move[_arr.Array](_arr.Array(size))
         else:
@@ -31,9 +35,16 @@ cdef class Array:
             raise IndexError("index {} is larger than the size of the array {}".
                                format(key, self._thisptr.size()))
 
-    property size:
-        def __get__(self):
-            return self._thisptr.size()
+    def __len__(self):
+        return self._thisptr.size()
+
+    def to_ndarray(self):
+        cdef np.npy_intp[1] dims
+        dims[0] = self._thisptr.size()
+        cdef arr = np.PyArray_SimpleNewFromData(1, &dims[0], np.NPY_DOUBLE, <void*>(self._thisptr.begin()))
+        Py_INCREF(self)
+        np.PyArray_SetBaseObject(arr, self)
+        return arr
 
 cpdef qlarray_from_pyarray(p):
     cdef Array x = Array(len(p))

quantlib/math/randomnumbers/_inverse_cumulative_rsg.pxd

@@ -1,4 +1,4 @@
-include '../../types.pxi'
+from quantlib.types cimport Real, Size
 from quantlib.methods.montecarlo._sample cimport Sample
 from libcpp.vector cimport vector
 

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()

quantlib/methods/montecarlo/_multipathgenerator.pxd

@@ -0,0 +1,18 @@
+from libcpp cimport bool
+from quantlib.handle cimport shared_ptr
+
+from quantlib._time_grid cimport TimeGrid
+from libcpp.vector cimport vector
+from ._multipath cimport MultiPath
+from ._sample cimport Sample
+from quantlib._stochastic_process cimport StochasticProcess
+
+cdef extern from 'ql/methods/montecarlo/multipathgenerator.hpp' namespace 'QuantLib' nogil:
+    cdef cppclass MultiPathGenerator[GSG]:
+        ctypedef Sample[MultiPath] sample_type
+        MultiPathGenerator(const shared_ptr[StochasticProcess]&,
+                           const TimeGrid&,
+                           GSG generator,
+                           bool brownianBridge) # = false)
+        const sample_type& next()
+        const sample_type& antithetic()

quantlib/methods/montecarlo/_sample.pxd

@@ -1,4 +1,4 @@
-include '../../types.pxi'
+from quantlib.types cimport Real
 cdef extern from 'ql/methods/montecarlo/sample.hpp' namespace 'QuantLib' nogil:
     cdef cppclass Sample[T]:
         T value

quantlib/methods/montecarlo/multipathgenerator.pxd

@@ -0,0 +1,13 @@
+from . cimport _multipathgenerator as _mpg
+from quantlib.math.randomnumbers._randomsequencegenerator cimport RandomSequenceGenerator
+from quantlib.math.randomnumbers._rngtraits cimport PseudoRandom as QlPseudoRandom, LowDiscrepancy as QlLowDiscrepancy, ZigguratPseudoRandom
+
+
+cdef class PseudoRandomMultiPathGenerator:
+    cdef _mpg.MultiPathGenerator[QlPseudoRandom.rsg_type]* _thisptr
+
+cdef class LowDiscrepancyMultiPathGenerator:
+    cdef _mpg.MultiPathGenerator[QlLowDiscrepancy.rsg_type]* _thisptr
+
+cdef class ZigguratMultiPathGenerator:
+    cdef _mpg.MultiPathGenerator[ZigguratPseudoRandom]* _thisptr

quantlib/methods/montecarlo/multipathgenerator.pyx

@@ -0,0 +1,50 @@
+from cython.operator cimport dereference as deref
+from libcpp cimport bool
+from quantlib.stochastic_process cimport StochasticProcess
+from quantlib.time_grid cimport TimeGrid
+from quantlib.math.randomnumbers.rngtraits cimport PseudoRandom, LowDiscrepancy, FastPseudoRandom
+from ._sample cimport Sample
+from .multipath cimport MultiPath
+from . cimport _multipath as _mp
+
+cdef class PseudoRandomMultiPathGenerator:
+    def __init__(self, StochasticProcess process, TimeGrid time_grid, PseudoRandom gen, bool brownian_bridge):
+        self._thisptr = new _mpg.MultiPathGenerator[QlPseudoRandom.rsg_type](process._thisptr, time_grid._thisptr, deref(gen._thisptr), brownian_bridge)
+
+    def __next__(self):
+        cdef Sample[_mp.MultiPath]* r = <Sample[_mp.MultiPath]*>&self._thisptr.next()
+        cdef MultiPath mp = MultiPath.__new__(MultiPath)
+        mp._thisptr = new _mp.MultiPath(r.value)
+        return (r.weight, mp)
+
+    def __dealloc__(self):
+        del self._thisptr
+
+
+
+cdef class LowDiscrepancyMultiPathGenerator:
+    def __init__(self, StochasticProcess process, TimeGrid time_grid, LowDiscrepancy gen, bool brownian_bridge):
+        self._thisptr = new _mpg.MultiPathGenerator[QlLowDiscrepancy.rsg_type](process._thisptr, time_grid._thisptr, deref(gen._thisptr), brownian_bridge)
+
+    def __next__(self):
+        cdef Sample[_mp.MultiPath]* r = <Sample[_mp.MultiPath]*>&self._thisptr.next()
+        cdef MultiPath mp = MultiPath.__new__(MultiPath)
+        mp._thisptr = new _mp.MultiPath(r.value)
+        return (r.weight, mp)
+
+    def __dealloc__(self):
+        del self._thisptr
+
+
+cdef class ZigguratMultiPathGenerator:
+    def __init__(self, StochasticProcess process, TimeGrid time_grid, FastPseudoRandom gen, bool brownian_bridge):
+        self._thisptr = new _mpg.MultiPathGenerator[ZigguratPseudoRandom](process._thisptr, time_grid._thisptr, deref(gen._thisptr), brownian_bridge)
+
+    def __next__(self):
+        cdef Sample[_mp.MultiPath]* r = <Sample[_mp.MultiPath]*>&self._thisptr.next()
+        cdef MultiPath mp = MultiPath.__new__(MultiPath)
+        mp._thisptr = new _mp.MultiPath(r.value)
+        return (r.weight, mp)
+
+    def __dealloc__(self):
+        del self._thisptr

quantlib/models/marketmodels/_brownian_generator.pxd

@@ -0,0 +1,14 @@
+from libcpp.vector cimport vector
+from quantlib.types cimport Real, Size
+from quantlib.handle cimport shared_ptr
+
+cdef extern from 'ql/models/marketmodels/browniangenerator.hpp' namespace 'QuantLib' nogil:
+    cdef cppclass BrownianGenerator:
+        Real nextStep(vector[Real]&) except +
+        Real nextPath()
+
+        Size numberOfFactors() const
+        Size numberOfSteps() const
+
+    cdef cppclass BrownianGeneratorFactory:
+        shared_ptr[BrownianGenerator] create(Size factor, Size steps) const

quantlib/models/marketmodels/brownian_generator.pxd

@@ -0,0 +1,8 @@
+from quantlib.handle cimport shared_ptr
+from . cimport _brownian_generator as _bg
+
+cdef class BrownianGenerator:
+    cdef shared_ptr[_bg.BrownianGenerator] gen
+
+cdef class BrownianGeneratorFactory:
+    cdef _bg.BrownianGeneratorFactory* factory