Check against old implementation

resolvelib/Makefile

-.PHONY: build
-
-CXXFLAGS_MINI:= -std=c++17 -Wfatal-errors -Wall -Wextra -shared -fPIC $(shell python3 -m pybind11 --includes) -Iducc/src -Ipybind11/include -fvisibility=hidden
-CXXFLAGS:= -O3 -march=native -ffast-math $(CXXFLAGS_MINI)
-#CXXFLAGS:= $(CXXFLAGS_MINI)
-
-
-OUT:= resolvelib$(shell python3-config --extension-suffix)
-
-$(OUT) : resolvelib.cpp
-	g++ $(CXXFLAGS) -o$(OUT) ducc/src/ducc0/infra/threading.cc $^
-
-build : $(OUT)

resolvelib/devel.py

@@ -17,6 +17,7 @@
 import resolve as rve
 import nifty8 as ift
 import resolvelib
+import numpy as np
 
 from cpp2py import Pybind11Operator
 from time import time
@@ -30,16 +31,23 @@ dom = {kk: dom[1:] for kk in pdom.labels}
 
 tgt = rve.default_sky_domain(pdom=pdom, sdom=sdom)
 
+opold = rve.polarization_matrix_exponential(tgt) @ rve.MultiFieldStacker(tgt, 0, tgt[0].labels)
+
 nthreads = 1
 loc = ift.from_random(dom)
+res0 = opold(loc)
 ntries = 100
 for nthreads  in [8]:
 #for nthreads  in range(1, 9):
     op = Pybind11Operator(dom, tgt, resolvelib.PolarizationMatrixExponential(nthreads))
+    assert op.domain is opold.domain
+    assert op.target is opold.target
     t0 = time()
     for _ in range(ntries):
-        op(loc)
+        res = op(loc)
     print(f"Nthreads {nthreads}: {(time()-t0)/ntries:.2f} s")
+
+    np.testing.assert_allclose(res0.val, res.val)
 exit()
 op(loc)
 exit()

resolvelib/resolvelib.cc

@@ -82,9 +82,9 @@ class PolarizationMatrixExponential {
                           auto &ouu,
                           auto &ovv
                           ){
-              auto pol{qq*qq + uu*uu + vv*vv};
-              oii = 0.5 * (exp(ii+pol) + exp(ii-pol));
-              auto tmp{0.5 * (exp(ii+pol) - exp(ii-pol)) / sqrt(pol)};
+              auto pol{sqrt(qq*qq + uu*uu + vv*vv)};
+              oii = exp(ii) * cosh(pol);
+              auto tmp{exp(ii) * sinh(pol) / pol};
               oqq = tmp * qq;
               ouu = tmp * uu;
               ovv = tmp * vv;