From da812bae05edb045a85fe77ef47b9ad916ce1fd7 Mon Sep 17 00:00:00 2001
From: Ivan Kondov <ivan.kondov@kit.edu>
Date: Sat, 19 Jan 2019 21:32:20 +0100
Subject: [PATCH] multi precision for scalar test, refactoring of besself()
---
classical/chebyshev.py | 20 ++++++++++----------
classical/tests.py | 11 +++++++----
2 files changed, 17 insertions(+), 14 deletions(-)
diff --git a/classical/chebyshev.py b/classical/chebyshev.py
index d98bf1d..c0f89fd 100644
--- a/classical/chebyshev.py
+++ b/classical/chebyshev.py
@@ -6,15 +6,18 @@ from scipy.special import jn, iv
from propagators import PolynomialPropagator
from liouville import LiouvillePowersQP
+
def cheb_coef(bf, nterms, alpha):
""" chebyshev expansion coefficients for the function exp(z*t) """
coef = bf(range(nterms), np.full(nterms, alpha))
coef[1:] = coef[1:]*2.
return coef
-def besf(fpprec, signed, prec=None):
+
+def besself(fpprec, signed, prec=None):
""" select the proper bessel function according to requested precision
- and selected axis: if signed if True [-1, +1] otherwise [-i, +i] """
+ and selected axis: if signed is True [-1, +1] otherwise [-i, +i] """
+ assert fpprec in ['fixed', 'multi']
if fpprec == 'fixed':
bessel = iv if signed else jn
elif fpprec == 'multi':
@@ -23,16 +26,14 @@ def besf(fpprec, signed, prec=None):
def bessel(orders, args):
bfun = [bess(o, a).evalf(*fargs) for o, a in zip(orders, args)]
return np.array(bfun)
- else:
- bessel = None
return bessel
+
def chebyshev_scalar(z, t, delta=None, lambda_min=None, nterms=None,
signed=False, fpprec='fixed', prec=None):
""" Chebyshev polynomial expansion of the scalar function exp(z*t) """
-
assert isinstance(z, np.ndarray)
- delta = max(abs(z)) if delta is None else delta
+ delta = max(abs(z)) if delta is None else delta
alpha = delta*t/2.
lambda_min = -delta/2. if lambda_min is None else lambda_min
nterms = int(np.ceil(abs(alpha)))+1 if nterms is None else nterms
@@ -41,10 +42,9 @@ def chebyshev_scalar(z, t, delta=None, lambda_min=None, nterms=None,
shift = -(2.*lambda_min/delta+1.)
sign = -1 if signed else 1
z = z*scale+shift
- exp = np.exp if fpprec == 'fixed' else lambda x: sp.exp(x).evalf(prec)
+ exp = np.exp if fpprec == 'fixed' else sp.exp
scaling = exp((lambda_min+delta/2.)*t)
- bessel_func = besf(fpprec, signed, prec=prec)
- ch_coef = cheb_coef(bessel_func, nterms, alpha)*scaling
+ ch_coef = cheb_coef(besself(fpprec, signed, prec), nterms, alpha)*scaling
ch_poly = np.empty((nterms, len(z)), dtype=z.dtype)
ch_poly[0, :] = 1.
ch_poly[1, :] = z
@@ -68,7 +68,7 @@ class Chebyshev(PolynomialPropagator):
self.lambda_min = -self.delta/2. if lambda_min is None else lambda_min
self.scale = 2./self.delta
self.shift = -(2.*self.lambda_min/self.delta+1.)
- bf = besf(fpprec, signed, prec=prec)
+ bf = besself(fpprec, signed, prec=prec)
sign = -1 if signed else 1
exp = np.exp if fpprec == 'fixed' else lambda x: sp.exp(x).evalf(prec)
diff --git a/classical/tests.py b/classical/tests.py
index 4f3d7c2..f6c79d7 100644
--- a/classical/tests.py
+++ b/classical/tests.py
@@ -470,10 +470,13 @@ class ChebyshevAssessmentTest(unittest.TestCase):
sample_length = 10
lambdas = (np.random.rand((sample_length))-0.5)*6*np.random.rand()
- for sign, fact in product((False, True), (1., 1.j)):
- apr = chebyshev_scalar(lambdas*fact, 1.1, nterms=16, signed=sign)
- ref = np.exp(lambdas*1.1*fact)
- self.assertTrue(np.allclose(apr, ref))
+ pars = [(False, True), (1., 1.j), ('fixed', 'multi')]
+ for sign, fact, fpprec in product(*pars):
+ values = lambdas*fact
+ apr = chebyshev_scalar(values, 1.1, nterms=24, signed=sign,
+ fpprec=fpprec)
+ ref = np.exp(values*1.1)
+ self.assertTrue(np.allclose(np.array(apr, dtype=ref.dtype), ref))
def chebyshev_symplectic_1p_test(self):
""" with one-particle system: anharmonic oscillator """
--
GitLab