divided difference tests, arbitrary precision for newtonian related funcs

classical/newtonian.py

 """ Newtonian propagator for classical particle dynamics """
 import numpy as np
+import sympy as sp
 from propagators import Propagator
 from liouville import LiouvillePowersQP
 
@@ -27,14 +28,29 @@ def divdiff_2(x, y, dtype=float):
     return np.diagonal(a)
 
 
-def leja_points_ga(number):
-    """ calculate Leja interpolation points based on algorithm from 
+def leja_points_ga(number, fpprec='fixed', prec=None, ptype='real'):
+    """ calculate Leja interpolation points based on an algorithm from
         G. Ashkenazi et al. JCP 103, 10005 (1995) """
 
-    center = 0.0
-    trials = np.linspace(-1., 1., 3*number).tolist()
+    if ptype == 'real':
+        if fpprec == 'fixed':
+            dtype = np.float64
+            center = dtype(0)
+            trials = np.linspace(dtype(-1), dtype(1), 3*number).tolist()
+            isclose = np.isclose
+        else:
+            flt = sp.Float
+            tol = flt(10**(-prec/2), prec)
+            center, low, high = flt(0, prec), flt(-1, prec), flt(1, prec)
+            step = (high-low)/(3*number-1)
+            trials = np.arange(low, high+step, step).tolist()
+            isclose = lambda x, y: abs(x-y) <= tol
+    elif ptype == 'complex':
+        # fixme: complex trial points in fixed and arbitrary precision
+        pass
+
     for index, trial in enumerate(trials):
-        if np.isclose(trial, center):
+        if isclose(trial, center):
             del trials[index]
     chosen = [trials.pop(0)]
     while len(chosen) < number:
@@ -43,17 +59,31 @@ def leja_points_ga(number):
     scaling = np.prod(np.absolute(center-np.array(chosen)))**(1./number)
     return np.array(chosen)/scaling, scaling
 
-def leja_points_jb(nflp):
-    """ calculate fast Leja points based on algorithm from 
+
+def leja_points_jb(nflp, fpprec='fixed', prec=None, ptype='real'):
+    """ calculate fast Leja points based on an algorithm from
         Baglama et al. Electronic Trans. Numer. Anal. 7, 124 (1998) """
 
-    zs = np.empty(nflp)
-    zt = np.empty(nflp)
-    zprod = np.empty(nflp)
-    index = np.empty((nflp, 2), dtype=np.int)
+    if fpprec == 'fixed':
+        dtype = np.float64 if ptype == 'real' else np.complex128
+        left = dtype(-1)
+        right = dtype(1)
+    else:
+        dtype = np.object
+        if ptype == 'real':
+            left = sp.Float(-1, prec)
+            right = sp.Float(1, prec)
+        else:
+            # fixme: complex initial points in arbitrary precision
+            pass
+
+    zs = np.empty(nflp, dtype=dtype)
+    zt = np.empty(nflp, dtype=dtype)
+    zprod = np.empty(nflp, dtype=dtype)
+    index = np.empty((nflp, 2), dtype=np.uint32)
 
-    zt[0] = -1.
-    zt[1] = 1.
+    zt[0] = left
+    zt[1] = right
     zt[2] = (zt[0] + zt[1])/2
     zs[0] = (zt[1] + zt[2])/2
     zs[1] = (zt[2] + zt[0])/2
@@ -97,30 +127,47 @@ def newton_scalar(lambdas, func, vals, scale=1.0, shift=0.0):
 class Newtonian(Propagator):
     """ Newtonian propagator for classical particle dynamics """
     name = 'newtonian'
-    def __init__(self, delta=1.0, lambda_min=None, nterms=10,
-                 liouville=LiouvillePowersQP, **kwargs):
+    def __init__(self, delta=1.0, lambda_min=None, nterms=None,
+                 liouville=LiouvillePowersQP, liou_params={}, fpprec='fixed',
+                 prec=None, **kwargs):
         Propagator.__init__(self, **kwargs)
 
         assert delta > 0
         assert nterms > 0
+        self.fpprec = fpprec
+        self.prec = prec
         self.delta = delta
         self.nterms = nterms
         self.lambda_min = -self.delta/2. if lambda_min is None else lambda_min
-        lambdas, scaling = leja_points_ga(self.nterms)
+        lambdas, scaling = self.ip_points(algo='ga')
         self.scale = self.delta*scaling
         self.shift = (self.lambda_min+self.delta/2.)/self.scale
         lamb_sc = (lambdas*self.scale+self.shift)
         if self.substep is None:
-            self.ddiffs = divdiff_1(lambdas, np.exp(lamb_sc*self.tstep))
+            self.ddiffs = self.dd_coeffs(lambdas, lamb_sc, self.tstep)
         else:
-            self.ddiffs = np.array([
-                divdiff_1(lambdas, np.exp(lamb_sc*st)) for st in self.subtime
-            ])
+            self.ddiffs = np.array([self.dd_coeffs(lambdas, lamb_sc, st)
+                                    for st in self.subtime])
 
         self.iln = liouville(self.nterms, self.syst, scale=1./self.scale,
-                             shift=-self.shift, nodes=lambdas,
+                             shift=-self.shift, nodes=lambdas, fpprec=fpprec,
                              recursion='newton')
 
+    def dd_coeffs(self, lamb, lamb_sc, t, algo=1):
+        """ select the algorithm for divided differences """
+        algo_map = {1: divdiff_1, 2: divdiff_2}
+        if self.fpprec == 'fixed':
+            return algo_map[algo](lamb, np.exp(lamb_sc*t))
+        elif self.fpprec == 'multi':
+            return algo_map[algo](lamb, [sp.exp(i) for i in lamb_sc*t])
+
+    def ip_points(self, ptype='real', algo='ga'):
+        """ select the algorithm for interpolation points """
+        algo_map = {'ga': leja_points_ga, 'jb': leja_points_jb,
+                    'ch': chebyshev_points}
+        return algo_map[algo](self.nterms, fpprec=self.fpprec, prec=self.prec,
+                              ptype=ptype)
+
     def step(self, time, q, p):
         """ perform a time step with initial coordinate and momentum """
         qt = np.dot(self.ddiffs, self.iln.qobj.get_val_float(q, p))

classical/tests.py

@@ -349,6 +349,19 @@ class PointsTest(unittest.TestCase):
         self.assertTrue(np.allclose(pol, func(vals), atol=1e-3))
 
 
+class DividedDifferenceTest(unittest.TestCase):
+    """ Tests for finite difference functions """
+
+    def divided_difference_test(self):
+        """ compare two algorithms for divided differences """
+        from newtonian import divdiff_1, divdiff_2
+        nodes = np.arange(-1, 1, 0.1)
+        funcvals = lambda x: 2*np.exp(3*x)
+        diffs1 = divdiff_1(nodes, funcvals(nodes))
+        diffs2 = divdiff_2(nodes, funcvals(nodes))
+        self.assertTrue(np.allclose(diffs1, diffs2))
+
+
 class NewtonianPropagatorTest(unittest.TestCase):
     """ Tests for the Newtonian polynomial propagator """