Refactored franke function + added test error_estimate from padua_fit

wafo/padua.py

@@ -132,12 +132,14 @@ class _ExampleFunctions(object):
 
         Maple: 0.40696958949155611906
         '''
-        exp = np.exp
+        def _exp(x, y, loc, scale, p2=2):
+            return np.exp(- (x-loc[0])**2/scale[0] - (y-loc[1])**p2/scale[1])
+        # exp = np.exp
         x9, y9 = 9. * x, 9. * y
-        return (3. / 4 * exp(-((x9 - 2)**2 + (y9 - 2)**2) / 4) +
-                3. / 4 * exp(-(x9 + 1)**2 / 49 - (y9 + 1) / 10) +
-                1. / 2 * exp(-((x9 - 7)**2 + (y9 - 3)**2) / 4) -
-                1. / 5 * exp(-(x9 - 4)**2 - (y9 - 7)**2))
+        return (3. / 4 * _exp(x9, y9, [2, 2], [4, 4]) +
+                3. / 4 * _exp(x9, y9, [-1, -1], [49, 10], p2=1) +
+                1. / 2 * _exp(x9, y9, [7, 3], [4, 4]) -
+                1. / 5 * _exp(x9, y9, [4, 7], [1, 1]))
 
     @staticmethod
     def half_sphere(x, y):
@@ -449,7 +451,7 @@ def paduavals2coefs(f):
 
 # TODO: padua_fit2 does not work correctly yet.
 def padua_fit2(Pad, fun, *args):
-    N = np.shape(Pad)[1]
+    # N = np.shape(Pad)[1]
     # recover the degree n from N = (n+1)(n+2)/2
     # _n = int(round(-3 + np.sqrt(1 + 8 * N)) / 2)
     C0f = fun(Pad[0], Pad[1], *args)

wafo/tests/test_padua.py

@@ -6,7 +6,7 @@ from numpy import cos, pi
 import numpy.testing as npt
 from numpy.testing import assert_array_almost_equal
 from wafo.padua import (padua_points, example_functions, padua_fit,
-                        padua_fit2,
+                        # padua_fit2,
                         padua_cubature, padua_val)
 
 
@@ -40,17 +40,19 @@ class PaduaTestCase(unittest.TestCase):
 
     def test_padua_fit_even_degree(self):
         points = padua_points(10)
-        C0f, _abs_error = padua_fit(points, example_functions, 6)
+        C0f, abs_error = padua_fit(points, example_functions, 6)
         expected = np.zeros((11, 11))
         expected[0, 0] = 1
         assert_array_almost_equal(C0f, expected, 15)
+        assert_array_almost_equal(abs_error, 1.2168216554799264e-15)
 
     def test_padua_fit_odd_degree(self):
         points = padua_points(9)
-        C0f, _abs_error = padua_fit(points, example_functions, 6)
+        C0f, abs_error = padua_fit(points, example_functions, 6)
         expected = np.zeros((10, 10))
         expected[0, 0] = 1
         assert_array_almost_equal(C0f, expected, 15)
+        assert_array_almost_equal(abs_error, 4.509537093983535e-17)
 
 # TODO: padua_fit2 does not work correctly
 #    def test_padua_fit_odd_degree2(self):
@@ -63,33 +65,36 @@ class PaduaTestCase(unittest.TestCase):
     def test_padua_cubature(self):
         domain = [0, 1, 0, 1]
         points = padua_points(500, domain)
-        C0f, _abs_error = padua_fit(points, example_functions, 0)
+        C0f, abs_error = padua_fit(points, example_functions, 0)
         val = padua_cubature(C0f, domain)
         expected = 4.06969589491556e-01
         assert_array_almost_equal(val, expected, 15)
+        assert_array_almost_equal(abs_error,  3.66470417665e-16)
 
     def test_padua_val_unordered(self):
         domain = [0, 1, 0, 1]
         points = padua_points(20, domain)
-        C0f, _abs_error = padua_fit(points, example_functions, 0)
-        X = [0, 0.5, 1]
-
+        C0f, abs_error = padua_fit(points, example_functions, 0)
+        X = np.array([0, 0.5, 1])
+        # true_val = example_functions.franke(X, X)
         val = padua_val(X, X, C0f, domain)
-        expected = [7.664205912849228e-01, 3.2621734202884815e-01,
-                    3.587865112678535e-02]
+        expected = [0.76642059128493,  0.32621734202885,  0.03587865112678]
         assert_array_almost_equal(val, expected, 14)
+        assert_array_almost_equal(abs_error,  0.003897032262116954)
 
     def test_padua_val_grid(self):
         domain = [0, 1, 0, 1]
         a, b, c, d = domain
         points = padua_points(21, domain)
-        C0f, _abs_error = padua_fit(points, example_functions, 0)
+        C0f, abs_error = padua_fit(points, example_functions, 0)
         X1 = np.linspace(a, b, 2)
         X2 = np.linspace(c, d, 2)
         val = padua_val(X1, X2, C0f, domain, use_meshgrid=True)
-        expected = [[7.664205912849229e-01, 1.0757071952145181e-01],
-                    [2.703371615911344e-01, 3.5734971024838565e-02]]
+
+        expected = [[0.76642059128493,  0.10757071952145],
+                    [0.27033716159114,  0.03573497102484]]
         assert_array_almost_equal(val, expected, 14)
+        assert_array_almost_equal(abs_error,  0.0022486904061664046)
 
 
 if __name__ == "__main__":