some cython stuff to boost memory nice version of KLE, some reorganization

setup.py

@@ -0,0 +1,8 @@
+from distutils.core import setup
+from Cython.Build import cythonize
+import numpy as np  
+
+setup(
+  ext_modules = cythonize(["stocproc/stocproc_c.pyx"]),
+  include_dirs = [np.get_include()],   
+)

stocproc/__init__.py

@@ -4,7 +4,4 @@ from .class_stocproc_kle import *
 from .class_stocproc import StocProc_FFT
 from .class_stocproc import StocProc_KLE
 
-
-
-
 import gquad

stocproc/class_stocproc.py

@@ -1,6 +1,6 @@
 import numpy as np
 from scipy.interpolate import InterpolatedUnivariateSpline
-from . import StocProc
+from .class_stocproc_kle import StocProc
 
 class ComplexInterpolatedUnivariateSpline(object):
     r"""
@@ -37,12 +37,15 @@ class _absStocProc(object):
             :param k: polynomial degree used for spline interpolation  
         """
         self._verbose = verbose
+        self.t_max = t_max
+        self.num_grid_points = num_grid_points
         self.t = np.linspace(0, t_max, num_grid_points)
         self._z = None
         self._interpolator = None
         self._k = k
         if seed is not None:
             np.random.seed(seed)
+        self._one_over_sqrt_2 = 1/np.sqrt(2)
 
     def __call__(self, t):
         r"""
@@ -104,7 +107,7 @@ class _absStocProc(object):
             #random complex normal samples
             if self._verbose > 1:
                 print("generate samples ...")
-            y = np.random.normal(size = 2*self.get_num_y()).view(np.complex)
+            y = np.random.normal(scale=self._one_over_sqrt_2, size = 2*self.get_num_y()).view(np.complex)
             if self._verbose > 1:
                 print("done")        
 
@@ -119,26 +122,18 @@ class StocProc_KLE(_absStocProc):
         - Calculate discrete stochastic process (using interpolation solution of fredholm equation) with num_grid_points nodes
         - invoke spline interpolator when calling 
     """
-    def __init__(self, r_tau, t_max, num_grid_points, ng_fredholm=None, seed=None, sig_min=1e-5, verbose=1, k=3):
+    def __init__(self, r_tau, t_max, ng_fredholm, ng_fac=4, seed=None, sig_min=1e-5, verbose=1, k=3):
         r"""
             :param r_tau: auto correlation function of the stochastic process
             :param t_max: specify time axis as [0, t_max]
-            :param num_grid_points: number of equidistant times on that axis
-            :param ng_fredholm: number of discrete time used to solve the underlying fredholm equation, ``None`` set ``ng_fredholm = num_grid_points``
             :param seed: if not ``None`` set seed to ``seed``
             :param sig_min: eigenvalue threshold (see KLE method to generate stochastic processes)
             :param verbose: 0: no output, 1: informational output, 2: (eventually some) debug info
             :param k: polynomial degree used for spline interpolation  
              
         """
-        if ng_fredholm is None:
+        self.ng_fac = ng_fac
-            ng_fredholm = num_grid_points
+        if ng_fac == 1:
-            
-        if num_grid_points < ng_fredholm:
-            print("WARNING: found 'num_grid_points < ng_fredholm' -> set 'num_grid_points == ng_fredholm'")
-
-        # if ng_fredholm == num_grid_points -> no kle_fredholm interpolation is needed
-        if ng_fredholm == num_grid_points:
             self.kle_interp = False
         else:
             self.kle_interp = True
@@ -150,7 +145,9 @@ class StocProc_KLE(_absStocProc):
                                                                        seed    = seed,
                                                                        verbose = verbose)
         
-        super().__init__(t_max=t_max, num_grid_points=num_grid_points, seed=seed, verbose=verbose, k=k)
+        ng = ng_fac * (ng_fredholm - 1) + 1  
+        
+        super().__init__(t_max=t_max, num_grid_points=ng, seed=seed, verbose=verbose, k=k)
                 
         # this is only needed access the underlaying stocproc KLE class
         # in a convenient fashion 
@@ -161,14 +158,15 @@ class StocProc_KLE(_absStocProc):
         
         self.verbose = verbose
 
-    def _cal_z(self, y):
+    def _calc_z(self, y):
         r"""
         uses the underlaying stocproc class to generate the process (see :py:class:`StocProc` for details) 
         """
         self.stocproc.new_process(y)
         
         if self.kle_interp:
-            return self.stocproc.x_t_array(self.t)
+            #return self.stocproc.x_t_array(np.linspace(0, self.t_max, self.num_grid_points))
+            return self.stocproc.x_t_mem_save(delta_t_fac = self.ng_fac)
         else:
             return self.stocproc.x_for_initial_time_grid()
         
@@ -176,7 +174,7 @@ class StocProc_KLE(_absStocProc):
         return self.num_y
         
 
-class StocProc_FFT(object):
+class StocProc_FFT(_absStocProc):
     r"""
         Simulate Stochastic Process using FFT method 
     """
@@ -187,16 +185,16 @@ class StocProc_FFT(object):
                          verbose         = verbose,
                          k               = k)
         
-        self.n_dft           = self.num_grid_points * 2 - 1
+        self.n_dft           = num_grid_points * 2 - 1
-        delta_t              = self.t_max / (self.num_grid_points-1)
+        delta_t              = t_max / (num_grid_points-1)
         self.delta_omega     = 2 * np.pi / (delta_t * self.n_dft)
           
         #omega axis
         omega = self.delta_omega*np.arange(self.n_dft)
         #reshape for multiplication with matrix xi
-        self.sqrt_spectral_density_times_sqrt_delta_omega_over_sqrt_2 = np.sqrt(spectral_density(omega)) * np.sqrt(self.delta_omega) / np.sqrt(2) 
+        self.sqrt_spectral_density_times_sqrt_delta_omega = np.sqrt(spectral_density(omega)) * np.sqrt(self.delta_omega) 
         
-        if self.verbose > 0:
+        if self._verbose > 0:
             print("stoc proc fft, spectral density sampling information:")
             print("  t_max      :", (t_max))
             print("  ng         :", (num_grid_points))
@@ -204,13 +202,13 @@ class StocProc_FFT(object):
             print("  omega_max  :", (self.delta_omega * self.n_dft))
             print("  delta_omega:", (self.delta_omega))
             
-    def _cal_z(self, y):
+    def _calc_z(self, y):
-        weighted_integrand = self.sqrt_spectral_density_times_sqrt_delta_omega_over_sqrt_2 * y 
+        weighted_integrand = self.sqrt_spectral_density_times_sqrt_delta_omega * y 
         #compute integral using fft routine
-        if self.verbose > 1:
+        if self._verbose > 1:
             print("calc process via fft ...")
         z = np.fft.fft(weighted_integrand)[0:self.num_grid_points]
-        if self.verbose > 1:
+        if self._verbose > 1:
             print("done")
         return z
 

stocproc/class_stocproc_kle.py

@@ -2,10 +2,12 @@ import functools
 import numpy as np
 import pickle
 
-from . import solve_hom_fredholm
+from .stocproc import solve_hom_fredholm
-from . import get_mid_point_weights
+from .stocproc import get_mid_point_weights
-from . import get_trapezoidal_weights_times
+from .stocproc import get_trapezoidal_weights_times
-from . import gquad
+import gquad
+
+from . import stocproc_c
 
 
 class StocProc(object):
@@ -91,6 +93,7 @@ class StocProc(object):
                  verbose    = 1):
         
         self.verbose = verbose
+        self._one_over_sqrt_2 = 1/np.sqrt(2)
         if fname is None:
             
             assert r_tau is not None
@@ -205,7 +208,7 @@ class StocProc(object):
         if yi is None:
             if self.verbose > 1:
                 print("generate samples ...")
-            self._Y = np.random.normal(size=2*self._num_ev).view(np.complex).reshape(self._num_ev,1)
+            self._Y = np.random.normal(scale = self._one_over_sqrt_2, size=2*self._num_ev).view(np.complex).reshape(self._num_ev,1)
             if self.verbose > 1:
                 print("done!")
         else:
@@ -220,7 +223,7 @@ class StocProc(object):
         equation. This is equivalent to calling :py:func:`stochastic_process_kle` with the
         same weights :math:`w_i` and time grid points :math:`s_i`.
         """
-        tmp = self._Y / np.sqrt(2) * self._sqrt_eig_val.reshape(self._num_ev,1) 
+        tmp = self._Y * self._sqrt_eig_val.reshape(self._num_ev,1) 
         if self.verbose > 1:
             print("calc process via matrix prod ...")
         res = np.tensordot(tmp, self._eig_vec, axes=([0],[1])).flatten()
@@ -265,27 +268,40 @@ class StocProc(object):
         if self.verbose > 1:
             print("done!")
         return res
+    
+    def x_t_mem_save(self, delta_t_fac):
+        a = delta_t_fac        
+        N1 = len(self._s)
+        N2 = a * (N1 - 1) + 1        
+        T = self._s[-1]
+        alpha_k = self._r_tau(np.linspace(-T, T, 2*N2 - 1))        
+        return stocproc_c.z_t(delta_t_fac = delta_t_fac,
+                              time_axis   = self._s,
+                              alpha_k     = alpha_k,
+                              Y           = self._Y[:,0],
+                              A           = self._A)
 
     def u_i_mem_save(self, delta_t_fac, i):
-        r"""
+        a = delta_t_fac        
-        """
-        a = delta_t_fac
-        assert isinstance(a, int)
-        T = self._s[-1]
         N1 = len(self._s)
-        N2 = a * (N1 - 1) + 1
+        N2 = a * (N1 - 1) + 1        
+        T = self._s[-1]
         alpha_k = self._r_tau(np.linspace(-T, T, 2*N2 - 1))
-        
+        return stocproc_c.eig_func_interp(delta_t_fac,
-        print(alpha_k[N2-1])
+                                          self._s,
-        
+                                          alpha_k,
-        print("WARNING! this needs to be cythonized")
+                                          self._w,
-        u_res = np.zeros(shape=N2, dtype=np.complex)
+                                          self._eig_val[i],
-        for j in range(N2):
+                                          self._eig_vec[:, i])
-            for l in range(N1):
+         
-                k = j - a*l + N2-1
+#         print("WARNING! this needs to be cythonized")
-                u_res[j] += self._w[l] * alpha_k[k] * self._eig_vec[l, i]      
+#         u_res = np.zeros(shape=N2, dtype=np.complex)
-
+#         for j in range(N2):
-        return u_res / self._eig_val[i]
+#             for l in range(N1):
+#                 k = j - a*l + N2-1
+#                 u_res[j] += self._w[l] * alpha_k[k] * self._eig_vec[l, i]      
+# 
+#         return u_res / self._eig_val[i]
         
     
     def u_i(self, t_array, i):
@@ -330,6 +346,32 @@ class StocProc(object):
         # A_j,i = w_j / sqrt(lambda_i) u_i(s_j)
         return np.tensordot(tmp, 1/self._sqrt_eig_val.reshape(1,self._num_ev) * self._A, axes=([0],[0]))
 
+    def u_i_all_mem_save(self, delta_t_fac):
+        a = delta_t_fac        
+        N1 = len(self._s)
+        N2 = a * (N1 - 1) + 1        
+        T = self._s[-1]
+        alpha_k = self._r_tau(np.linspace(-T, T, 2*N2 - 1))
+        
+        return stocproc_c.eig_func_all_interp(delta_t_fac = delta_t_fac,
+                                              time_axis   = self._s,
+                                              alpha_k     = alpha_k, 
+                                              weights     = self._w,
+                                              eigen_val   = self._eig_val,
+                                              eigen_vec   = self._eig_vec)
+        
+#         print("WARNING! this needs to be cythonized")
+#         u_res = np.zeros(shape=(N2, self.num_ev()), dtype=np.complex)
+#         for i in range(self.num_ev()):
+#             for j in range(N2):
+#                 for l in range(N1):
+#                     k = j - a*l + N2-1
+#                     u_res[j, i] += self._w[l] * alpha_k[k] * self._eig_vec[l, i]      
+#  
+#             u_res[:, i] /= self._eig_val[i]
+#             
+#         return u_res
+
     def eigen_vector_i(self, i):
         r"""Returns the i-th eigenvector (solution of the discrete Fredhom equation)"""
         return self._eig_vec[:,i]

stocproc/stocproc.py

@@ -194,9 +194,9 @@ def stochastic_process_kle(r_tau, t, w, num_samples, seed = None, sig_min = 1e-4
 
     if verbose > 0:
         print("generate samples ...")
-    x = np.random.normal(size=(2*num_samples*num_of_functions)).view(np.complex).reshape(num_of_functions, num_samples).T \
+    x = np.random.normal(scale=1/np.sqrt(2), size=(2*num_samples*num_of_functions)).view(np.complex).reshape(num_of_functions, num_samples).T
-         / np.sqrt(2) * sig                                           # random quantities all aligned for num_samples samples
+                                                                      # random quantities all aligned for num_samples samples
-    x_t_array = np.tensordot(x, eig_vec, axes=([1],[1]))              # multiplication with the eigenfunctions == base of Karhunen-Loève expansion
+    x_t_array = np.tensordot(x*sig, eig_vec, axes=([1],[1]))              # multiplication with the eigenfunctions == base of Karhunen-Loève expansion
     
     if verbose > 0:
         print("done!")
@@ -386,9 +386,9 @@ def stochastic_process_fft(spectral_density, t_max, num_grid_points, num_samples
         print("  delta_omega: {:.2}".format(delta_omega))
         print("generate samples ...")
     #random complex normal samples
-    xi = (np.random.normal(size = (2*num_samples*n_dft)).view(np.complex)).reshape(num_samples, n_dft)
+    xi = (np.random.normal(scale=1/np.sqrt(2), size = (2*num_samples*n_dft)).view(np.complex)).reshape(num_samples, n_dft)
     #each row contain a different integrand
-    weighted_integrand = sqrt_spectral_density * np.sqrt(delta_omega) / np.sqrt(2) * xi 
+    weighted_integrand = sqrt_spectral_density * np.sqrt(delta_omega) * xi 
     #compute integral using fft routine
     z_ast = np.fft.fft(weighted_integrand, axis = 1)[:, 0:num_grid_points]
     #corresponding time axis

stocproc/stocproc_c.pyx

@@ -0,0 +1,101 @@
+import numpy as np
+cimport numpy as np
+
+DTYPE_CPLX = np.complex128
+ctypedef np.complex128_t DTYPE_CPLX_t
+
+DTYPE_DBL = np.float64
+ctypedef np.float64_t DTYPE_DBL_t
+
+cpdef np.ndarray[DTYPE_CPLX_t, ndim=1] eig_func_interp(unsigned int                     delta_t_fac,
+                                                       np.ndarray[DTYPE_DBL_t,  ndim=1] time_axis,
+                                                       np.ndarray[DTYPE_CPLX_t, ndim=1] alpha_k, 
+                                                       np.ndarray[DTYPE_DBL_t,  ndim=1] weights,
+                                                       double                           eigen_val,
+                                                       np.ndarray[DTYPE_CPLX_t, ndim=1] eigen_vec):
+
+    cdef unsigned int N1
+    N1 = len(time_axis)
+    
+    cdef unsigned int N2
+    N2 = delta_t_fac * (N1 - 1) + 1
+    
+    cdef np.ndarray[DTYPE_CPLX_t, ndim=1] u_res 
+    u_res = np.zeros(shape=N2, dtype=DTYPE_CPLX)
+    
+    cdef unsigned int j
+    cdef unsigned int l
+    cdef unsigned int k
+    for j in range(N2):
+        for l in range(N1):
+            k = j - delta_t_fac*l + N2-1
+            u_res[j] = u_res[j] + weights[l] * alpha_k[k] * eigen_vec[l]      
+
+    return u_res / eigen_val      
+
+
+
+cpdef np.ndarray[DTYPE_CPLX_t, ndim=2] eig_func_all_interp(unsigned int                     delta_t_fac,
+                                                           np.ndarray[DTYPE_DBL_t,  ndim=1] time_axis,
+                                                           np.ndarray[DTYPE_CPLX_t, ndim=1] alpha_k, 
+                                                           np.ndarray[DTYPE_DBL_t,  ndim=1] weights,
+                                                           np.ndarray[DTYPE_DBL_t,  ndim=1] eigen_val,
+                                                           np.ndarray[DTYPE_CPLX_t, ndim=2] eigen_vec):
+
+    cdef unsigned int N1
+    N1 = len(time_axis)
+    
+    cdef unsigned int N2
+    N2 = delta_t_fac * (N1 - 1) + 1
+    
+    cdef unsigned int num_ev
+    num_ev = len(eigen_val)
+    
+    cdef np.ndarray[DTYPE_CPLX_t, ndim=2] u_res 
+    u_res = np.zeros(shape=(N2,num_ev), dtype=DTYPE_CPLX)
+    
+    cdef unsigned int i
+    cdef unsigned int j
+    cdef unsigned int l
+    cdef unsigned int k
+    for i in range(num_ev):
+        for j in range(N2):
+            for l in range(N1):
+                k = j - delta_t_fac*l + N2-1
+                u_res[j, i] = u_res[j, i] + weights[l] * alpha_k[k] * eigen_vec[l, i]
+                      
+            u_res[j, i] = u_res[j, i]/eigen_val[i]
+            
+    return u_res   
+
+
+
+cpdef np.ndarray[DTYPE_CPLX_t, ndim=1] z_t(unsigned int                     delta_t_fac,
+                                           np.ndarray[DTYPE_DBL_t,  ndim=1] time_axis,
+                                           np.ndarray[DTYPE_CPLX_t, ndim=1] alpha_k,
+                                           np.ndarray[DTYPE_CPLX_t, ndim=1] Y,
+                                           np.ndarray[DTYPE_CPLX_t, ndim=2] A):
+    
+    cdef unsigned int N1
+    N1 = len(time_axis)
+    
+    cdef unsigned int N2
+    N2 = delta_t_fac * (N1 - 1) + 1
+    
+    cdef unsigned int num_ev
+    num_ev = len(Y)
+    
+    cdef np.ndarray[DTYPE_CPLX_t, ndim=1] z_t_res 
+    z_t_res = np.zeros(shape=N2, dtype=DTYPE_CPLX)
+ 
+    cdef unsigned int i
+    cdef unsigned int j
+    cdef unsigned int a
+    
+    for j in range(N2):
+        for i in range(N1):
+            k = j - delta_t_fac*i + N2-1
+            for a in range(num_ev):
+                z_t_res[j] = z_t_res[j] + Y[a]*alpha_k[k]*A[i,a]
+    
+    return z_t_res

tests/test_stocproc.py

@@ -197,7 +197,8 @@ def test_func_vs_class_KLE_FFT():
                                 num_grid_points  = ng,
                                 seed             = seed)
     
-    x_t_array_class = stoc_proc.get_z()    
+    stoc_proc.new_process()
+    x_t_array_class = stoc_proc.get_z()
     
     print(np.max(np.abs(x_t_array_func - x_t_array_class)))
     assert np.all(x_t_array_func == x_t_array_class), "stochastic_process_fft vs. StocProc Class not identical"
@@ -305,11 +306,16 @@ def test_stocproc_KLE_splineinterpolation(plot=False):
     # leads to N+1 grid points
     ng_fredholm   = 60
     ng_kle_interp = ng_fredholm*3 
-    ng_fine       = ng_fredholm*30  
+    ng_fine       = ng_fredholm*15
     
     seed = 0
     sig_min = 1e-5
-    stoc_proc = sp.StocProc_KLE(r_tau, t_max, ng_kle_interp, ng_fredholm, seed, sig_min)
+    stoc_proc = sp.StocProc_KLE(r_tau       = r_tau,
+                                t_max       = t_max, 
+                                ng_fredholm = ng_fredholm,
+                                ng_fac      = 3,
+                                seed        = seed, 
+                                sig_min     = sig_min)
   
     finer_t = np.linspace(0, t_max, ng_fine)
     
@@ -317,10 +323,12 @@ def test_stocproc_KLE_splineinterpolation(plot=False):
     
     ac_conj     = np.zeros(shape=(ng_fine, ng_fine), dtype=np.complex)
     ac_not_conj = np.zeros(shape=(ng_fine, ng_fine), dtype=np.complex)
+
     print("generate samples ...")    
     for n in range(ns):
         stoc_proc.new_process()
         x_t = stoc_proc(finer_t)
+
         ac_conj     += x_t.reshape(ng_fine, 1) * np.conj(x_t.reshape(1, ng_fine))
         ac_not_conj += x_t.reshape(ng_fine, 1) * x_t.reshape(1, ng_fine)
     print("done!")
@@ -652,17 +660,17 @@ def show_auto_grid_points_result():
     diff = sp.mean_error(r_t_s, r_t_s_exact)
     diff_max = sp.max_error(r_t_s, r_t_s_exact)
     
-    plt.plot(t_large, diff)
+#     plt.plot(t_large, diff)
-    plt.plot(t_large, diff_max)
+#     plt.plot(t_large, diff_max)
-    plt.yscale('log')
+#     plt.yscale('log')
-    plt.grid()
+#     plt.grid()
-    plt.show()
+#     plt.show()
     
 def test_ui_mem_save():
     s_param = 1
     gamma_s_plus_1 = gamma(s_param+1)
     r_tau = lambda tau : corr(tau, s_param, gamma_s_plus_1)
-    t_max = 1
+    t_max = 5
     
     N1 = 100
     a  = 5
@@ -674,6 +682,8 @@ def test_ui_mem_save():
 
     stoc_proc = sp.StocProc.new_instance_with_trapezoidal_weights(r_tau, t_max, ng=N1, sig_min = 1e-4)
     
+    ui_all_ms = stoc_proc.u_i_all_mem_save(delta_t_fac=a)
+    
     for i in range(stoc_proc.num_ev()):
 
         ui_ms = stoc_proc.u_i_mem_save(delta_t_fac=a, i=i)
@@ -692,16 +702,56 @@ def test_ui_mem_save():
 #         plt.show()
         
         assert np.allclose(ui_ms, ui), "{}".format(max(np.abs(ui_ms - ui)))
+        assert np.allclose(ui_all_ms[:, i], ui), "{}".format(max(np.abs(ui_all_ms[:, i] - ui)))
         
-            
+
+def test_z_t_mem_save():            
+    s_param = 1
+    gamma_s_plus_1 = gamma(s_param+1)
+    r_tau = lambda tau : corr(tau, s_param, gamma_s_plus_1)
+    t_max = 5
     
+    N1 = 100
+    a  = 5
+    N2 = a*(N1 - 1) + 1
+    sig_min = 0
+    
+    t_fine = np.linspace(0, t_max, N2)
+    
+    assert abs( (t_max/(N1-1)) - a*(t_fine[1]-t_fine[0]) ) < 1e-14, "{}".format(abs( (t_max/(N1-1)) - (t_fine[1]-t_fine[0]) ))
+
+    stoc_proc = sp.StocProc.new_instance_with_trapezoidal_weights(r_tau, t_max, ng=N1, sig_min=sig_min)
+    
+    z_t_mem_save = stoc_proc.x_t_mem_save(delta_t_fac = a)
+    z_t = stoc_proc.x_t_array(t_fine)
+    
+    z_t_rough = stoc_proc.x_for_initial_time_grid()
+        
+#     plt.plot(t_fine, np.real(z_t_mem_save), color='k')
+#     plt.plot(t_fine, np.imag(z_t_mem_save), color='k')
+#      
+#     plt.plot(t_fine, np.real(z_t), color='r')
+#     plt.plot(t_fine, np.imag(z_t), color='r')
+#      
+#     plt.plot(stoc_proc._s, np.real(z_t_rough), marker = 'o', ls='', color='b')
+#     plt.plot(stoc_proc._s, np.imag(z_t_rough), marker = 'o', ls='', color='b')
+#      
+#     plt.grid()
+#      
+#     plt.show()    
+    
+    assert np.allclose(z_t_mem_save, z_t), "{}".format(max(np.abs(z_t_mem_save - z_t)))
+    
+    
+    
+ 
         
 if __name__ == "__main__":
 #     test_stochastic_process_KLE_correlation_function(plot=False)
 #     test_stochastic_process_FFT_correlation_function(plot=False)
 #     test_func_vs_class_KLE_FFT()
 #     test_stochastic_process_KLE_interpolation(plot=False)
-#     test_stocproc_KLE_splineinterpolation(plot=False)
+    test_stocproc_KLE_splineinterpolation(plot=False)
 #     test_stochastic_process_FFT_interpolation(plot=False)
 #     test_stocProc_eigenfunction_extraction()
 #     test_orthonomality()
@@ -709,5 +759,7 @@ if __name__ == "__main__":
 #     show_auto_grid_points_result()
 #     test_chache()
 #     test_dump_load()
-    test_ui_mem_save()
+#     test_ui_mem_save()
+#     test_z_t_mem_save()
+
     pass