From d8da71884f2c0996cf6e0df16dd8bc4c5e3dd152 Mon Sep 17 00:00:00 2001
From: Aaron Virshup <avirshup@gmail.com>
Date: Thu, 24 Aug 2017 09:09:34 -0700
Subject: [PATCH] Adaptively increase resolution for numerical integrals in
 tests

---
 moldesign/_tests/helpers.py            |  4 ++--
 moldesign/_tests/test_gaussian_math.py | 31 +++++++++++++++++++++-----
 2 files changed, 28 insertions(+), 7 deletions(-)

diff --git a/moldesign/_tests/helpers.py b/moldesign/_tests/helpers.py
index 07fc108..2105188 100644
--- a/moldesign/_tests/helpers.py
+++ b/moldesign/_tests/helpers.py
@@ -230,7 +230,7 @@ def assert_almost_equal(actual, desired, **kwargs):
 """ Decorator to disable tests that need internet if we can't connect to the network """
 
 
-def generate_grid(g, g2=None):
+def generate_grid(g, g2=None, npoints=64):
     from moldesign import mathutils
 
     if g2 is None: g2 = g
@@ -243,6 +243,6 @@ def generate_grid(g, g2=None):
     ranges = np.ones((3, 2))*5.0*width
     ranges[:, 0] *= -1
     ranges += ((g.center + g2.center)/2.0)[:, None]
-    grid = mathutils.VolumetricGrid(*ranges, npoints=64)
+    grid = mathutils.VolumetricGrid(*ranges, npoints=npoints)
     allpoints = grid.allpoints()
     return allpoints, grid
diff --git a/moldesign/_tests/test_gaussian_math.py b/moldesign/_tests/test_gaussian_math.py
index d2fc817..309b1a7 100644
--- a/moldesign/_tests/test_gaussian_math.py
+++ b/moldesign/_tests/test_gaussian_math.py
@@ -173,11 +173,32 @@ def _assert_numerical_analytical_overlaps_match(g1, g2):
     else:
         helpers.assert_almost_equal(prod.integral, olap)
 
-    allpoints, grid = helpers.generate_grid(g1, g2)
-    with np.errstate(under='ignore'):
-        prodvals = g1(allpoints) * g2(allpoints)
-    numsum = prodvals.sum() * grid.dx * grid.dy * grid.dz
-    helpers.assert_almost_equal(numsum, olap, decimal=4)
+    def assert_with_resolution(npoints):
+        allpoints, grid = helpers.generate_grid(g1, g2, npoints)
+        with np.errstate(under='ignore'):
+            prodvals = g1(allpoints) * g2(allpoints)
+        numsum = prodvals.sum() * grid.dx * grid.dy * grid.dz
+        helpers.assert_almost_equal(numsum, olap, decimal=4)
+
+    # If numerical isn't equal to analytical, try again with higher resolution
+    # to make sure the failure isn't due to a sparse grid:
+    try:
+        assert_with_resolution(64)
+    except AssertionError:
+        pass
+    else:
+        return
+
+    try:
+        assert_with_resolution(128)
+    except AssertionError:
+        pass
+    else:
+        return
+
+    assert_with_resolution(256)
+
+
 
 
 @pytest.mark.parametrize('withunits', [False, True])