diff --git a/wafo/integrate_oscillating.py b/wafo/integrate_oscillating.py index 7a62eca..4860c9d 100644 --- a/wafo/integrate_oscillating.py +++ b/wafo/integrate_oscillating.py @@ -116,11 +116,11 @@ def osc_weights(omega, g, dg, x, basis, ab, *args, **kwds): return np.asarray(w).ravel() -class QuadOsc(object): +class _Integrator(object): info = namedtuple('info', ['error_estimate', 'n']) - def __init__(self, f, g, dg=None, a=-1, b=1, basis=chebyshev_basis, s=15, - precision=10, endpoints=False, maxiter=17, full_output=False): + def __init__(self, f, g, dg=None, a=-1, b=1, basis=chebyshev_basis, s=1, + precision=10, endpoints=True, full_output=False): self.f = f self.g = g self.dg = nd.Derivative(g) if dg is None else dg @@ -129,10 +129,18 @@ class QuadOsc(object): self.b = b self.s = s self.endpoints = endpoints - self.maxiter = maxiter self.precision = precision self.full_output = full_output + +class QuadOsc(_Integrator): + def __init__(self, f, g, dg=None, a=-1, b=1, basis=chebyshev_basis, s=15, + precision=10, endpoints=False, full_output=False, maxiter=17): + self.maxiter = maxiter + super(QuadOsc, self).__init__(f, g, dg=dg, a=a, b=b, basis=basis, s=s, + precision=precision, endpoints=endpoints, + full_output=full_output) + @staticmethod def _change_interval_to_0_1(f, g, dg, a, b): def f1(t, *args, **kwds): @@ -344,22 +352,8 @@ def chebyshev_roots(M, delta=None): return x -class AdaptiveLevin(object): +class AdaptiveLevin(_Integrator): '''Return integral for the Levin-type and adaptive Levin-type methods''' - info = namedtuple('info', ['error_estimate', 'n']) - - def __init__(self, f, g, dg=None, a=-1, b=1, basis=chebyshev_basis, s=1, - precision=10, endpoints=True, full_output=False): - self.f = f - self.g = g - self.dg = nd.Derivative(g) if dg is None else dg - self.basis = basis - self.a = a - self.b = b - self.s = s - self.precision = precision - self.endpoints = endpoints - self.full_output = full_output @staticmethod def aLevinTQ(omega, ff, gg, dgg, x, s, basis, *args, **kwds): @@ -510,16 +504,10 @@ class EvansWebster(AdaptiveLevin): def __init__(self, f, g, dg=None, a=-1, b=1, basis=chebyshev_basis, s=8, precision=10, endpoints=False, full_output=False): - self.f = f - self.g = g - self.dg = nd.Derivative(g) if dg is None else dg - self.basis = basis - self.a = a - self.b = b - self.s = s - self.endpoints = endpoints - self.precision = precision - self.full_output = full_output + super(EvansWebster, + self).__init__(f, g, dg=dg, a=a, b=b, basis=basis, s=s, + precision=precision, endpoints=endpoints, + full_output=full_output) def aLevinTQ(self, omega, ff, gg, dgg, x, s, basis, *args, **kwds): w = evans_webster_weights(omega, gg, dgg, x, basis, *args, **kwds) @@ -531,55 +519,6 @@ class EvansWebster(AdaptiveLevin): def _get_num_points(self, s, prec, betam): return 8 if s > 1 else int(prec / max(np.log10(betam + 1), 1) + 1) -# def _QaL(self, s, a, b, omega, *args, **kwds): -# '''if s>1,the integral is computed by Q_s^L''' -# scale = (b - a) / 2 -# offset = (a + b) / 2 -# prec = self.precision # desired precision -# -# def ff(t, *args, **kwds): -# return scale * self.f(scale * t + offset, *args, **kwds) -# -# def gg(t, *args, **kwds): -# return self.g(scale * t + offset, *args, **kwds) -# -# def dgg(t, *args, **kwds): -# return scale * self.dg(scale * t + offset, *args, **kwds) -# dg_a = abs(omega * dgg(-1, *args, **kwds)) -# dg_b = abs(omega * dgg(1, *args, **kwds)) -# g_a = abs(omega * gg(-1, *args, **kwds)) -# g_b = abs(omega * gg(1, *args, **kwds)) -# delta, alpha = min(dg_a, dg_b), min(g_a, g_b) -# -# betam = delta # * scale -# if self.endpoints: -# if (delta < 10 or alpha <= 10 or s > 1): -# points = chebyshev_extrema -# else: -# points = adaptive_levin_points -# elif (delta < 10 or alpha <= 10 or s > 1): -# points = chebyshev_roots -# else: -# points = tanh_sinh_open_nodes -# -# m = 8 if s > 1 else int(prec / max(np.log10(betam + 1), 1) + 1) -# abseps = 10*10.0**-prec -# num_collocation_point_list = 2*m*np.arange(1, 6, 2) -# # range(num_points, num_points+3, 2) -# basis = self.basis -# Q = 1e+300 -# n = 0 -# for num_collocation_points in num_collocation_point_list: -# Q_old = Q -# x = points(num_collocation_points-1, betam) -# Q = self.aLevinTQ(omega, ff, gg, dgg, x, s, basis, *args, **kwds) -# n += num_collocation_points -# err = np.abs(Q-Q_old) -# if err <= abseps: -# break -# info = self.info(err, n) -# return Q, info - if __name__ == '__main__': tanh_sinh_nodes(16)