complex leja points on an unit circle

classical/newtonian.py

@@ -29,29 +29,48 @@ def divdiff_2(x, y, dtype=float, fpprec='fixed', prec=None):
     return np.diagonal(a)
 
 
-def leja_points_ga(number, fpprec='fixed', prec=None, ptype='real'):
+def trial_points(number, dtype=float, fpprec='fixed', prec=None):
+    """ construct the trial points for the algorithm from
+        G. Ashkenazi et al. JCP 103, 10005 (1995) """
+    zero = {'fixed': dtype(0), 'multi': sp.N(dtype(0), prec)}
+    one = {'fixed': dtype(1), 'multi': sp.N(dtype(1), prec)}
+    ci = {'fixed': 1j, 'multi': sp.N(1j, prec)}
+    pi = {'fixed': np.pi, 'multi': sp.N(sp.pi, prec)}
+    cosf = {'fixed': np.cos, 'multi': lambda x: sp.N(sp.cos(x), prec)}
+    sinf = {'fixed': np.sin, 'multi': lambda x: sp.N(sp.sin(x), prec)}
+    fspace = {
+        'fixed': lambda l, h, n: np.linspace(l, h, n),
+        'multi': lambda l, h, n: np.arange(l, h+(h-l)/(n-1), (h-l)/(n-1))
+    }
+    center = zero[fpprec]
+
+    if dtype == complex:
+        # equally spaced points on an unit circle on the complex plane
+        low, high = 0, 2*pi[fpprec]
+        phisf = np.arange(low, high, (high-low)/number)
+        circlef = lambda f: cosf[fpprec](f)+ci[fpprec]*sinf[fpprec](f)
+        zsfun = np.vectorize(circlef)
+        trials = zsfun(phisf)
+    elif dtype == float:
+        # equally spaced points in the interval [-1, 1]
+        low, high = -one[fpprec], one[fpprec]
+        trials = fspace[fpprec](low, high, number)
+    else:
+        raise ValueError('dtype must be complex or float')
+    return trials.tolist(), center
+
+
+def leja_points_ga(number, dtype=float, fpprec='fixed', prec=None):
     """ calculate Leja interpolation points based on an algorithm from
         G. Ashkenazi et al. JCP 103, 10005 (1995) """
 
-    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
+    isclose = {'fixed': np.isclose,
+               'multi': lambda x, y: abs(x-y) <= sp.N(10**(-prec/2), prec)}
+    ntrials = {float: 3*number, complex: 2*number}
+    trials, center = trial_points(ntrials[dtype], dtype, fpprec, prec)
 
     for index, trial in enumerate(trials):
-        if isclose(trial, center):
+        if isclose[fpprec](trial, center):
             del trials[index]
     chosen = [trials.pop(0)]
     while len(chosen) < number:
@@ -61,30 +80,23 @@ def leja_points_ga(number, fpprec='fixed', prec=None, ptype='real'):
     return np.array(chosen)/scaling, scaling
 
 
-def leja_points_jb(nflp, fpprec='fixed', prec=None, ptype='real'):
-    """ calculate fast Leja points based on an algorithm from
+def leja_points_jb(nflp, dtype=float, fpprec='fixed', prec=None):
+    """ calculate fast Leja points on an interval based on an algorithm from
         Baglama et al. Electronic Trans. Numer. Anal. 7, 124 (1998) """
 
-    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
+    assert dtype == float
+    ndtype = {'fixed': dtype, 'multi': np.object}
+    scaling = {'fixed': float(1), 'multi': sp.N(float(1), prec)}
+    left = {'fixed': dtype(-1), 'multi': sp.N(dtype(-1), prec)}
+    right = {'fixed': dtype(1), 'multi': sp.N(dtype(1), prec)}
 
-    zs = np.empty(nflp, dtype=dtype)
-    zt = np.empty(nflp, dtype=dtype)
-    zprod = np.empty(nflp, dtype=dtype)
+    zs = np.empty(nflp, dtype=ndtype[fpprec])
+    zt = np.empty(nflp, dtype=ndtype[fpprec])
+    zprod = np.empty(nflp, dtype=ndtype[fpprec])
     index = np.empty((nflp, 2), dtype=np.uint32)
 
-    zt[0] = left
-    zt[1] = right
+    zt[0] = left[fpprec]
+    zt[1] = right[fpprec]
     zt[2] = (zt[0] + zt[1])/2
     zs[0] = (zt[1] + zt[2])/2
     zs[1] = (zt[2] + zt[0])/2
@@ -105,7 +117,7 @@ def leja_points_jb(nflp, fpprec='fixed', prec=None, ptype='real'):
         zprod[maxi] = np.prod(zs[maxi]-zt[:i])
         zprod[i-1] = np.prod(zs[i-1]-zt[:i])
         zprod[:i] = zprod[:i]*(zs[:i]-zt[i])
-    return zt
+    return zt, scaling[fpprec]
 
 
 def chebyshev_points(number):
@@ -162,12 +174,12 @@ class Newtonian(Propagator):
         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'):
+    def ip_points(self, dtype=float, 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)
+        return algo_map[algo](self.nterms, dtype=dtype, fpprec=self.fpprec,
+                              prec=self.prec)
 
     def step(self, time, q, p):
         """ perform a time step with initial coordinate and momentum """

classical/tests.py

@@ -313,14 +313,14 @@ class PointsTest(unittest.TestCase):
         self.assertTrue(np.allclose(pol, func(vals), rtol=1e-3))
 
     def leja_points_exponent_test(self):
-        """ Leja points intepolation of exponential function """
+        """ Leja points interpolation of exponential function """
         vals = np.arange(-1, 1, 0.1)
         func = lambda x: 2*np.exp(3*x)
         lambdas, scaling = leja_points_ga(15)
         pol = newton_scalar(lambdas, func, vals, scale=scaling)
         self.assertTrue(np.allclose(pol, func(vals), rtol=1e-8))
-        lambdas = leja_points_jb(15)
-        pol = newton_scalar(lambdas, func, vals)
+        lambdas, scaling = leja_points_jb(15)
+        pol = newton_scalar(lambdas, func, vals, scale=scaling)
         self.assertTrue(np.allclose(pol, func(vals), rtol=1e-8))
 
     def leja_points_runge_test(self):
@@ -330,8 +330,8 @@ class PointsTest(unittest.TestCase):
         lambdas, scaling = leja_points_ga(60)
         pol = newton_scalar(lambdas, func, vals, scale=scaling)
         self.assertTrue(np.allclose(pol, func(vals), rtol=1e-4))
-        lambdas = leja_points_jb(60)
-        pol = newton_scalar(lambdas, func, vals)
+        lambdas, scaling = leja_points_jb(60)
+        pol = newton_scalar(lambdas, func, vals, scale=scaling)
         self.assertTrue(np.allclose(pol, func(vals), rtol=1e-4))
 
     def leja_points_runge_test_affine(self):
@@ -343,7 +343,7 @@ class PointsTest(unittest.TestCase):
         pol = newton_scalar(lambdas, func, vals, scale=scaling*(right-left)/2,
                             shift=(left+right)/2)
         self.assertTrue(np.allclose(pol, func(vals), atol=1e-3))
-        lambdas_jb = leja_points_jb(60)
+        lambdas_jb, scaling = leja_points_jb(60)
         pol = newton_scalar(lambdas_jb, func, vals, scale=(right-left)/2,
                             shift=(left+right)/2)
         self.assertTrue(np.allclose(pol, func(vals), atol=1e-3))