From d4cc24b570d246d4e8d917a5140fa16786f5f78d Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Mon, 6 Dec 2021 17:19:48 -0600
Subject: [PATCH] Add, test basis_func_to_polynomial

---
 modepy/modes.py    | 40 ++++++++++++++++++++++++++++++++++++++++
 test/test_modes.py | 39 +++++++++++++++++++++++++++++++++++++++
 2 files changed, 79 insertions(+)

diff --git a/modepy/modes.py b/modepy/modes.py
index 20b6a057..d1a714b9 100644
--- a/modepy/modes.py
+++ b/modepy/modes.py
@@ -32,6 +32,8 @@
 from modepy.spaces import FunctionSpace, PN, QN
 from modepy.shapes import Shape, Simplex, Hypercube
 
+from pymbolic.mapper import IdentityMapper
+
 if TYPE_CHECKING:
     import pymbolic.primitives
 
@@ -796,6 +798,44 @@ def symbolicize_function(
     else:
         return result
 
+
+class _DeclutterMapper(IdentityMapper):
+    def map_common_subexpression(self, expr):
+        return self.rec(expr.child)
+
+    def map_if(self, expr):
+        # Assumes that the general-but-singular expression is in the "then_"
+        # branch.
+        return self.rec(expr.then)
+
+
+class _DivisionFoundError(Exception):
+    pass
+
+
+class _DivisionRaiser(IdentityMapper):
+    def map_quotient(self, expr):
+        raise _DivisionFoundError()
+
+
+def basis_func_to_polynomial(expr):
+    draise = _DivisionRaiser()
+    declutter = _DeclutterMapper()
+
+    from pymbolic.interop.sympy import PymbolicToSympyMapper, SympyToPymbolicMapper
+    p2s = PymbolicToSympyMapper()
+    s2p = SympyToPymbolicMapper()
+
+    import sympy as sp
+
+    clean_expr = declutter(expr)
+    sympy_expr = sp.sympify(p2s(clean_expr)).simplify()
+    poly_expr = s2p(sympy_expr)
+
+    draise(poly_expr)
+
+    return poly_expr
+
 # }}}
 
 
diff --git a/test/test_modes.py b/test/test_modes.py
index 308cc071..fca3f679 100644
--- a/test/test_modes.py
+++ b/test/test_modes.py
@@ -287,6 +287,45 @@ def test_symbolic_basis(shape, order, basis_getter):
 # }}}
 
 
+@pytest.mark.parametrize("shape", [
+    mp.Simplex(1),
+    mp.Simplex(2),
+    mp.Simplex(3),
+    ])
+@pytest.mark.parametrize("order", [3, 6])
+@pytest.mark.parametrize("basis_getter", [
+    (mp.orthonormal_basis_for_space),
+    ])
+def test_nonsing_onb(shape, order, basis_getter):
+    logging.basicConfig(level="INFO")
+    basis = basis_getter(mp.space_for_shape(shape, order), shape)
+    sym_basis = [mp.symbolicize_function(f, shape.dim) for f in basis.functions]
+
+    rng = np.random.Generator(np.random.PCG64(17))
+    r = mp.random_nodes_for_shape(shape, 10000, rng=rng)
+
+    from modepy.modes import basis_func_to_polynomial
+
+    for i, (sym_f, f) in enumerate(zip(sym_basis, basis.functions)):
+        poly_f = basis_func_to_polynomial(sym_f)
+
+        logger.info("BASISFUNC %d: %s", i, poly_f)
+
+        sym_val = MyEvaluationMapper({"r": r, "abs": abs})(poly_f)
+        ref_val = f(r)
+
+        ref_norm = la.norm(ref_val, np.inf)
+        err = la.norm(sym_val-ref_val, np.inf)
+        if ref_norm:
+            err = err/ref_norm
+
+        logger.info("ERROR: %.5g", err)
+
+        assert np.allclose(sym_val, ref_val)
+
+    # }}}
+
+
 # {{{ test_modal_coeffs_by_projection
 
 @pytest.mark.parametrize("dim", [2, 3])