From bd30feabab42d417879a71ac32430b741d2b76c5 Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Fri, 5 Jun 2020 16:33:13 +0200
Subject: [PATCH] cleanup

---
 python/demos/totalconvolve_accuracy.py        | 90 +++++++++++++++++++
 .../not_yet_integrated}/compress_utils.h      |  0
 .../mr_util/not_yet_integrated}/crangeset.h   |  0
 .../mr_util/not_yet_integrated}/datatypes.h   |  0
 .../gridder_accuracy_tools.py                 |  0
 .../not_yet_integrated}/healpix_map.cc        |  0
 .../mr_util/not_yet_integrated}/healpix_map.h |  0
 .../mr_util/not_yet_integrated}/math_utils.h  |  0
 .../mr_util/not_yet_integrated}/misc_utils.h  |  0
 .../mr_util/not_yet_integrated}/moc.h         |  0
 .../mr_util/not_yet_integrated}/stacktrace.h  |  0
 .../mr_util/not_yet_integrated}/tofrac.c      |  0
 .../mr_util/not_yet_integrated}/vec.h         |  0
 13 files changed, 90 insertions(+)
 create mode 100644 python/demos/totalconvolve_accuracy.py
 rename {not_yet_integrated => src/mr_util/not_yet_integrated}/compress_utils.h (100%)
 rename {not_yet_integrated => src/mr_util/not_yet_integrated}/crangeset.h (100%)
 rename {not_yet_integrated => src/mr_util/not_yet_integrated}/datatypes.h (100%)
 rename {not_yet_integrated => src/mr_util/not_yet_integrated}/gridder_accuracy_tools.py (100%)
 rename {not_yet_integrated => src/mr_util/not_yet_integrated}/healpix_map.cc (100%)
 rename {not_yet_integrated => src/mr_util/not_yet_integrated}/healpix_map.h (100%)
 rename {not_yet_integrated => src/mr_util/not_yet_integrated}/math_utils.h (100%)
 rename {not_yet_integrated => src/mr_util/not_yet_integrated}/misc_utils.h (100%)
 rename {not_yet_integrated => src/mr_util/not_yet_integrated}/moc.h (100%)
 rename {not_yet_integrated => src/mr_util/not_yet_integrated}/stacktrace.h (100%)
 rename {not_yet_integrated => src/mr_util/not_yet_integrated}/tofrac.c (100%)
 rename {not_yet_integrated => src/mr_util/not_yet_integrated}/vec.h (100%)

diff --git a/python/demos/totalconvolve_accuracy.py b/python/demos/totalconvolve_accuracy.py
new file mode 100644
index 0000000..b9a7864
--- /dev/null
+++ b/python/demos/totalconvolve_accuracy.py
@@ -0,0 +1,90 @@
+import ducc_0_1.totalconvolve as totalconvolve
+import numpy as np
+import ducc_0_1.sht as sht
+import ducc_0_1.misc as misc
+import time
+import matplotlib.pyplot as plt
+
+np.random.seed(48)
+
+def nalm(lmax, mmax):
+    return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
+
+
+def random_alm(lmax, mmax, ncomp):
+    res = np.random.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \
+     + 1j*np.random.uniform(-1., 1., (nalm(lmax, mmax), ncomp))
+    # make a_lm with m==0 real-valued
+    res[0:lmax+1,:].imag = 0.
+    return res
+
+
+def compress_alm(alm,lmax):
+    res = np.empty(2*len(alm)-lmax-1, dtype=np.float64)
+    res[0:lmax+1] = alm[0:lmax+1].real
+    res[lmax+1::2] = np.sqrt(2)*alm[lmax+1:].real
+    res[lmax+2::2] = np.sqrt(2)*alm[lmax+1:].imag
+    return res
+
+
+def myalmdot(a1,a2,lmax,mmax,spin):
+    return np.vdot(compress_alm(a1,lmax),compress_alm(np.conj(a2),lmax))
+
+
+def convolve(alm1, alm2, lmax):
+    job = sht.sharpjob_d()
+    job.set_triangular_alm_info(lmax, lmax)
+    job.set_gauss_geometry(lmax+1, 2*lmax+1)
+    map = job.alm2map(alm1)*job.alm2map(alm2)
+    job.set_triangular_alm_info(0,0)
+    return job.map2alm(map)[0]*np.sqrt(4*np.pi)
+
+
+lmax=60
+kmax=13
+ncomp=1
+separate=True
+ncomp2 = ncomp if separate else 1
+
+# get random sky a_lm
+# the a_lm arrays follow the same conventions as those in healpy
+slm = random_alm(lmax, lmax, ncomp)
+
+# build beam a_lm
+blm = random_alm(lmax, kmax, ncomp)
+
+
+t0=time.time()
+# build interpolator object for slm and blm
+foo = totalconvolve.PyInterpolator(slm,blm,separate,lmax, kmax, epsilon=1e-4, nthreads=2)
+print("setup time: ",time.time()-t0)
+nth = lmax+1
+nph = 2*lmax+1
+
+
+# compute a convolved map at a fixed psi and compare it to a map convolved
+# "by hand"
+
+ptg = np.zeros((nth,nph,3))
+ptg[:,:,0] = (np.pi*(0.5+np.arange(nth))/nth).reshape((-1,1))
+ptg[:,:,1] = (2*np.pi*(0.5+np.arange(nph))/nph).reshape((1,-1))
+ptg[:,:,2] = np.pi*0.2
+t0=time.time()
+# do the actual interpolation
+bar=foo.interpol(ptg.reshape((-1,3))).reshape((nth,nph,ncomp2))
+print("interpolation time: ", time.time()-t0)
+plt.subplot(2,2,1)
+plt.imshow(bar[:,:,0])
+bar2 = np.zeros((nth,nph))
+blmfull = np.zeros(slm.shape)+0j
+blmfull[0:blm.shape[0],:] = blm
+for ith in range(nth):
+    rbeamth=misc.rotate_alm(blmfull[:,0], lmax, ptg[ith,0,2],ptg[ith,0,0],0)
+    for iph in range(nph):
+        rbeam=misc.rotate_alm(rbeamth, lmax, 0, 0, ptg[ith,iph,1])
+        bar2[ith,iph] = convolve(slm[:,0], rbeam, lmax).real
+plt.subplot(2,2,2)
+plt.imshow(bar2)
+plt.subplot(2,2,3)
+plt.imshow(bar2-bar[:,:,0])
+plt.show()
diff --git a/not_yet_integrated/compress_utils.h b/src/mr_util/not_yet_integrated/compress_utils.h
similarity index 100%
rename from not_yet_integrated/compress_utils.h
rename to src/mr_util/not_yet_integrated/compress_utils.h
diff --git a/not_yet_integrated/crangeset.h b/src/mr_util/not_yet_integrated/crangeset.h
similarity index 100%
rename from not_yet_integrated/crangeset.h
rename to src/mr_util/not_yet_integrated/crangeset.h
diff --git a/not_yet_integrated/datatypes.h b/src/mr_util/not_yet_integrated/datatypes.h
similarity index 100%
rename from not_yet_integrated/datatypes.h
rename to src/mr_util/not_yet_integrated/datatypes.h
diff --git a/not_yet_integrated/gridder_accuracy_tools.py b/src/mr_util/not_yet_integrated/gridder_accuracy_tools.py
similarity index 100%
rename from not_yet_integrated/gridder_accuracy_tools.py
rename to src/mr_util/not_yet_integrated/gridder_accuracy_tools.py
diff --git a/not_yet_integrated/healpix_map.cc b/src/mr_util/not_yet_integrated/healpix_map.cc
similarity index 100%
rename from not_yet_integrated/healpix_map.cc
rename to src/mr_util/not_yet_integrated/healpix_map.cc
diff --git a/not_yet_integrated/healpix_map.h b/src/mr_util/not_yet_integrated/healpix_map.h
similarity index 100%
rename from not_yet_integrated/healpix_map.h
rename to src/mr_util/not_yet_integrated/healpix_map.h
diff --git a/not_yet_integrated/math_utils.h b/src/mr_util/not_yet_integrated/math_utils.h
similarity index 100%
rename from not_yet_integrated/math_utils.h
rename to src/mr_util/not_yet_integrated/math_utils.h
diff --git a/not_yet_integrated/misc_utils.h b/src/mr_util/not_yet_integrated/misc_utils.h
similarity index 100%
rename from not_yet_integrated/misc_utils.h
rename to src/mr_util/not_yet_integrated/misc_utils.h
diff --git a/not_yet_integrated/moc.h b/src/mr_util/not_yet_integrated/moc.h
similarity index 100%
rename from not_yet_integrated/moc.h
rename to src/mr_util/not_yet_integrated/moc.h
diff --git a/not_yet_integrated/stacktrace.h b/src/mr_util/not_yet_integrated/stacktrace.h
similarity index 100%
rename from not_yet_integrated/stacktrace.h
rename to src/mr_util/not_yet_integrated/stacktrace.h
diff --git a/not_yet_integrated/tofrac.c b/src/mr_util/not_yet_integrated/tofrac.c
similarity index 100%
rename from not_yet_integrated/tofrac.c
rename to src/mr_util/not_yet_integrated/tofrac.c
diff --git a/not_yet_integrated/vec.h b/src/mr_util/not_yet_integrated/vec.h
similarity index 100%
rename from not_yet_integrated/vec.h
rename to src/mr_util/not_yet_integrated/vec.h
-- 
GitLab