SP TanhSinh works

stocproc/__init__.py

@@ -59,5 +59,6 @@ if sys.version_info.major < 3:
 try:
     from .stocproc import StocProc_FFT
     from .stocproc import StocProc_KLE
+    from .stocproc import StocProc_TanhSinh
 except ImportError:
     pass

stocproc/method_fft.py

@@ -170,7 +170,7 @@ def fourier_integral_TanhSinh_with_feedback(integrand, w_max, tau, h, kmax):
     I = []
     feed_back = "ok"
     for ti in tau:
-        r = weights_k[0] * integrand(tmp1)*np.exp(-1j*tmp1*ti)
+        r = weights_k[0] * integrand(tmp1) * np.exp(-1j * tmp1 * ti)
         for i in range(1, kmax+1):
             if (y_k[i] * tmp1) == 0:
                 log.debug("y_k is 0")

stocproc/stocproc.py

@@ -462,7 +462,7 @@ class StocProc_TanhSinh(_absStocProc):
         h = 0.1
         k = 15
         kstep = 5
-        conv_fac = 2
+        conv_fac = 1
         old_d = None
         log.info("find h and kmax for TanhSinh integration ...")
         while True:
@@ -476,7 +476,7 @@ class StocProc_TanhSinh(_absStocProc):
             if fb == 'ok':
                 k += kstep
             else:
-                log.info("lowest diff with h {:.3e}: {:.3e} < tol ({:.3e}) -> new h {:.3e}".format(h, d, intgr_tol, h/2))
+                log.info("lowest diff with h {:.3e}: {:.3e} < tol ({:.3e}) -> new h {:.3e}".format(h, d[0], intgr_tol, h/2))
                 h /= 2
                 k = 15
                 kstep *= 2
@@ -511,41 +511,46 @@ class StocProc_TanhSinh(_absStocProc):
         wk = [method_fft.wk(h, ki) for ki in range(1, k+1)]
         wk = np.hstack((wk[::-1], [method_fft.wk(h, 0)], wk))
 
-        yk = [method_fft.yk(h, ki) for ki in range(1, k+1)]
+        yk = np.asarray([method_fft.yk(h, ki) for ki in range(1, k+1)])
+        tmp1 = wmax/2
+        self.omega_k = np.hstack( (yk[::-1] * tmp1, tmp1, (2 - yk) * tmp1))
+        self.fl = np.sqrt(tmp1*wk*spectral_density(self.omega_k)/np.pi)
 
-        tmp1 = wmax / 2
-        #sd = np.hstack((spectral_density(yk * tmp1, ))
+        super().__init__(t_max=t_max,
+                         num_grid_points=N,
+                         seed=seed,
+                         scale=scale)
 
 
+    @staticmethod
+    def get_key(t_max, bcf_ref, intgr_tol=1e-2, intpl_tol=1e-2):
+        return bcf_ref, t_max, intgr_tol, intpl_tol
 
 
-        self.yl = wk*spectral_density()
+    def __getstate__(self):
+        return self.fl, self.omega_k, self.num_grid_points, self.t_max, self._seed, self.scale, self.key
 
 
+    def __setstate__(self, state):
+        self.fl, self.omega_k, num_grid_points, t_max, seed, scale, self.key = state
+        super().__init__(t_max=t_max,
+                         num_grid_points=num_grid_points,
+                         seed=seed,
+                         scale=scale)
 
 
+    def calc_z(self, y):
+        r"""calculate
+    
+        .. math::
+            Z(t_l) = sum_k \sqrt{\frac{w_k J(\omega_k)}{\pi}} Y_k e^{-\i \omega_k t_l}
+        """
+        z = np.empty(shape=self.num_grid_points, dtype=np.complex128)
+        for i, ti in enumerate(self.t):
+            z[i] = np.sum(self.fl*y*np.exp(-1j*self.omega_k*ti))
+
+        return z
 
 
-
-        # assert abs(2 * np.pi - N * dx * dt) < 1e-12
-        #
-        # print("Fourier Integral Boundaries: [{:.3e}, {:.3e}]".format(a, b))
-        # print("Number of Nodes            : {}".format(N))
-        # print("yields dx                  : {:.3e}".format(dx))
-        # print("yields dt                  : {:.3e}".format(dt))
-        # print("yields t_max               : {:.3e}".format((N - 1) * dt))
-        #
-        # num_grid_points = int(np.ceil(t_max / dt)) + 1
-        #
-        # assert num_grid_points <= N
-        #
-        # t_max = (num_grid_points - 1) * dt
-        #
-        # super().__init__(t_max=t_max,
-        #                  num_grid_points=num_grid_points,
-        #                  seed=seed,
-        #                  scale=scale)
-        #
-        # self.yl = spectral_density(dx * np.arange(N) + a + dx / 2) * dx / np.pi
-        # self.yl = np.sqrt(self.yl)
-        # self.omega_min_correction = np.exp(-1j * (a + dx / 2) * self.t)  # self.t is from the parent class
+    def get_num_y(self):
+        return len(self.fl)

tests/test_method_fft.py

@@ -378,17 +378,41 @@ def test_calc_abN():
 
 
 
+def test_SP_TanhSinh():
+    s = -0.5
+    wc = 5
+    tmax = 25
+
+    sd = lambda w: w**s*np.exp(-w/wc)
+    bcf = lambda tau: (wc/(1+1j*wc*tau))**(s+1)*gamma_func(s+1)/np.pi
+
+    _sp = sp.StocProc_TanhSinh(spectral_density=sd, t_max=tmax, bcf_ref=bcf, intgr_tol=1e-2, intpl_tol=1e-2, seed=0)
+
+    d_tol = [0.09, 0.05]
+    for j, N in enumerate([1000, 5000]):
+        t = np.linspace(0, tmax, 500)
+        idx = 200
+        c = 0
+        for i in range(N):
+            _sp.new_process()
+            zt = _sp(t)
+            c += zt*np.conj(zt[idx])
+        c /= N
+
+        d = np.max(np.abs(c - bcf(t-t[idx])))
+        print(d)
+        assert d < d_tol[j]
 
-    
     
 if __name__ == "__main__":
     logging.basicConfig(level=logging.INFO)
     #logging.basicConfig(level=logging.DEBUG)
-    test_find_integral_boundary()
-    test_fourier_integral_finite_boundary()
-    test_fourier_integral_infinite_boundary(plot=False)
-    test_get_N_a_b_for_accurate_fourier_integral()
-    test_get_N_a_b_for_accurate_fourier_integral_b_only()
-    test_get_dt_for_accurate_interpolation()
-    test_sclicing()
-    test_calc_abN()
+    # test_find_integral_boundary()
+    # test_fourier_integral_finite_boundary()
+    # test_fourier_integral_infinite_boundary(plot=False)
+    # test_get_N_a_b_for_accurate_fourier_integral()
+    # test_get_N_a_b_for_accurate_fourier_integral_b_only()
+    # test_get_dt_for_accurate_interpolation()
+    # test_sclicing()
+    # test_calc_abN()
+    test_SP_TanhSinh()