Merge pull request #5 from prmiles/extend_precision

Extend precision

Author: Paul Miles

.gitignore

@@ -1,3 +1,3 @@
 __pycache__
 *.spyproject
-
+*.coverage

pyfod/GaussLaguerre.py

@@ -8,14 +8,26 @@
 
 import numpy as np
 import sys
+from sympy.integrals.quadrature import gauss_gen_laguerre as sp_gauss_laguerre
+import sympy as sp
 
 
 class GaussLaguerre:
 
-    def __init__(self, N=5, start=0.0, finish=1.0, alpha=0.0, f=None):
+    def __init__(self, N=5, start=0.0, finish=1.0, alpha=0.0,
+                 f=None, extend_precision=True, n_digits=30):
         self.description = 'Gaussian-Laguerre Quadrature'
-        points, weights = np.polynomial.laguerre.laggauss(deg=N)
-        self.points = 1 - np.exp(-points)
+        if extend_precision is False:
+            points, weights = np.polynomial.laguerre.laggauss(deg=N)
+            self.points = 1 - np.exp(-points)
+        else:
+            points, weights = sp_gauss_laguerre(
+                    n=N, n_digits=n_digits, alpha=0)
+            points = [-p for p in points]
+            points = sp.Array(points)
+            self.points = sp.Array(
+                    np.ones(shape=len(points))) - points.applyfunc(sp.exp)
+            weights = sp.Array(weights)
         self.weights = weights
         self.start = start
         self.finish = finish
@@ -31,23 +43,59 @@ def integrate(self, f=None, alpha=None):
             sys.exit('No function defined... provide function f')
         # transform kernel
         span = self.finish - self.start
-        return (self.weights*(span**(1-alpha)*f(span*self.points
-                              + self.start)*(1-self.points)**(-alpha))).sum()
+        # check if sympy
+        if isinstance(self.points, sp.Array):
+            evalpoints = self.points.applyfunc(
+                    lambda x: span*x + self.start)
+            feval = evalpoints.applyfunc(f)
+            coef = self.points.applyfunc(
+                    lambda x: span**(1-alpha)*(1-x)**(-alpha))
+            s = 0
+            for ii, (w, f, t) in enumerate(zip(self.weights, feval, coef)):
+                s += w*f*t
+            return s
+        else:
+            coef = span**(1-alpha)*(1-self.points)**(-alpha)
+            s = (self.weights*(coef*f(span*self.points + self.start))).sum()
+            return s
 
 
 if __name__ == '__main__':  # pragma: no cover
 
-    def f(t):
-        return np.cos(2*t)
+    n_digits = 50
+    alpha = 0.9
+    start = 0.0
+    finish = 1.0
+    dt = 1e-6
+    N = 24
 
-    GLag = GaussLaguerre(N=100, start=1.0, finish=12.0)
+    def fcosnp(x):
+        return np.cos(2*x) + 3
 
-    a = GLag.integrate(f)
-    print('a = {}'.format(a))
+    GLag = GaussLaguerre(N=N, start=start, finish=finish,
+                         f=fcosnp, alpha=0, extend_precision=False)
+    F1 = GLag.integrate()
+    print('Int(cos(2t)+3) = {} ({})'.format(F1, 3.4546))
+
+    def fexp(x):
+        return sp.exp(2*x)
+    GLag = GaussLaguerre(N=N, start=start, finish=finish,
+                         f=fexp, alpha=alpha, n_digits=n_digits)
+    F1 = GLag.integrate()
+    GLag = GaussLaguerre(N=N, start=start, finish=finish-dt,
+                         f=fexp, alpha=alpha, n_digits=n_digits)
+    F2 = GLag.integrate()
+    print('D[f(t)=exp(2t)] = {} ({})'.format(
+            (F1-F2)/(dt*sp.gamma(1-0.9)), 13.815))
+
+    def fcos(x):
+        return sp.cos(2*x)
 
-    GLag = GaussLaguerre(N=8, start=0.0, finish=1.0, f=f, alpha=0.9)
+    GLag = GaussLaguerre(N=N, start=start, finish=finish,
+                         f=fcos, alpha=alpha, n_digits=n_digits)
     F1 = GLag.integrate()
-    dt = 1e-4
-    GLag = GaussLaguerre(N=8, start=0.0, finish=1.0-dt, f=f, alpha=0.9)
+    GLag = GaussLaguerre(N=N, start=start, finish=finish-dt,
+                         f=fcos, alpha=alpha, n_digits=n_digits)
     F2 = GLag.integrate()
-    print('D[f(t)] = {}'.format((F1-F2)/dt))
+    print('D[f(t)=cos(2t)] = {} ({})'.format(
+            (F1-F2)/(dt*sp.gamma(1-0.9)), -1.779))

requirements.txt

@@ -1,3 +1,4 @@
 numpy>=1.14
 scipy>=1.0
+sympy>=1.3
 coveralls

test/test_GaussLaguerre.py

@@ -5,6 +5,7 @@
 @author: prmiles
 """
 
+import sympy as sp
 import numpy as np
 import unittest
 from pyfod.GaussLaguerre import GaussLaguerre
@@ -13,12 +14,15 @@
 def f(t):
     return np.cos(t)
 
+def fsp(t):
+    return sp.cos(t)
+
 
 # --------------------------
-class Initialization(unittest.TestCase):
+class Initialization_non_extended(unittest.TestCase):
 
     def test_init(self):
-        GQ = GaussLaguerre(N=10, start=1.0, finish=12.0)
+        GQ = GaussLaguerre(N=10, start=1.0, finish=12.0, extend_precision=False)
         self.assertTrue(hasattr(GQ, 'f'), msg='Expect attribute f to exist')
         self.assertEqual(GQ.f, None, msg='Expect value of None')
         self.assertEqual(GQ.points.size, 10, msg='Expect 10 nodes')
@@ -28,7 +32,7 @@ def test_init(self):
 
     def test_init_with_f(self):
 
-        GQ = GaussLaguerre(N=10, start=1.0, finish=12.0, alpha=0.5, f=f)
+        GQ = GaussLaguerre(N=10, start=1.0, finish=12.0, alpha=0.5, f=f, extend_precision=False)
         self.assertTrue(hasattr(GQ, 'f'), msg='Expect attribute f to exist')
         self.assertEqual(GQ.f, f, msg='Expect function f')
         self.assertEqual(GQ.points.size, 10, msg='Expect 10 nodes')
@@ -37,19 +41,60 @@ def test_init_with_f(self):
 
 
 # --------------------------
-class Integrate(unittest.TestCase):
+class Initialization_with_extended(unittest.TestCase):
+
+    def test_init(self):
+        GQ = GaussLaguerre(N=10, start=1.0, finish=12.0)
+        self.assertTrue(hasattr(GQ, 'f'), msg='Expect attribute f to exist')
+        self.assertEqual(GQ.f, None, msg='Expect value of None')
+        self.assertEqual(len(GQ.points), 10, msg='Expect 10 nodes')
+        self.assertEqual(len(GQ.weights), 10, msg='Expect 10 weights')
+        self.assertEqual(GQ.alpha, 0, msg='Expect alpha eq 0')
+
+
+    def test_init_with_f(self):
+        
+        GQ = GaussLaguerre(N=10, start=1.0, finish=12.0, alpha=0.5, f=fsp)
+        self.assertTrue(hasattr(GQ, 'f'), msg='Expect attribute f to exist')
+        self.assertEqual(GQ.f, fsp, msg='Expect function fsp')
+        self.assertEqual(len(GQ.points), 10, msg='Expect 10 nodes')
+        self.assertEqual(len(GQ.weights), 10, msg='Expect 10 weights')
+        self.assertEqual(GQ.alpha, 0.5, msg='Expect alpha eq 0.5')
+
+
+# --------------------------
+class Integrate_non_extended(unittest.TestCase):
 
     def test_no_f(self):
-        GQ = GaussLaguerre()
+        GQ = GaussLaguerre(extend_precision=False)
         with self.assertRaises(SystemExit):
             GQ.integrate()
 
     def test_with_f(self):
-        GQ = GaussLaguerre()
+        GQ = GaussLaguerre(extend_precision=False)
         a = GQ.integrate(f=f)
         self.assertEqual(a.size, 1, msg='Expect float return')
 
     def test_with_alpha(self):
-        GQ = GaussLaguerre()
+        GQ = GaussLaguerre(extend_precision=False)
         a = GQ.integrate(f=f, alpha=0.5)
-        self.assertEqual(a.size, 1, msg='Expect float return')
+        self.assertEqual(a.size, 1, msg='Expect float return')
+
+
+# --------------------------
+class Integrate_with_extended(unittest.TestCase):
+    
+    def test_no_f(self):
+        GQ = GaussLaguerre()
+        with self.assertRaises(SystemExit):
+            GQ.integrate()
+
+    def test_with_f(self):
+        GQ = GaussLaguerre()
+        a = GQ.integrate(f=fsp)
+        self.assertTrue(a.is_Float, msg='Expect float return')
+
+    def test_with_alpha(self):
+        GQ = GaussLaguerre()
+        a = GQ.integrate(f=fsp, alpha=0.5)
+        self.assertTrue(a.is_Float, msg='Expect float return')