Cython

gistfile1.pyx

import numpy as np
cimport cython
cimport numpy as np
 
from libc.stdint cimport uint32_t, int32_t
from libc.math cimport sqrt
from libc.math cimport fabs
from libc.math cimport pow
 
 
ctypedef float     real_t
ctypedef uint32_t  uint_t
ctypedef int32_t   int_t
 
cdef:
    real_t f1 = .1
    real_t g1 = .01
    real_t g2 = .1
 
cdef real_t f(real_t X) nogil:
    return -f1*X
cdef real_t g(real_t X) nogil:
    return sqrt((g1 + g2*X*X))
cdef real_t scale(real_t X) nogil:
    return 1/fabs(g1 + g2*X*X)
 
@cython.boundscheck(False)
@cython.wraparound(False)
def simulation_c(real_t L = 0, uint_t N = 100000, real_t dt = .001, real_t init = .1):
    
    cdef:
        uint_t t
        np.ndarray[real_t, ndim=1] dW = np.asarray(np.random.randn(N)*np.sqrt(dt), dtype = np.float32)
        np.ndarray[real_t, ndim=1] X = np.zeros(N, dtype = np.float32)
        np.ndarray[real_t, ndim=1] T = np.zeros(N, dtype = np.float32)
    
    X[0] = init
    
    for t in range(N - 1):
        X[t+1] = X[t] + f(X[t])*dt*pow(scale(X[t]), L) +  g(X[t])*dW[t]*sqrt(pow(scale(X[t]), L))
        T[t+1] = T[t] + dt*pow(scale(X[t]), L)
    
    return X, T