test_totalconvolve.py 3.4 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
2
3
import numpy as np
import pytest
from numpy.testing import assert_
Martin Reinecke's avatar
Martin Reinecke committed
4
import ducc_0_1.totalconvolve as totalconvolve
Martin Reinecke's avatar
Martin Reinecke committed
5
import ducc_0_1.sht as sht
6
import ducc_0_1.misc as misc
Martin Reinecke's avatar
Martin Reinecke committed
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

pmp = pytest.mark.parametrize


def _l2error(a, b):
    return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2))


def _assert_close(a, b, epsilon):
    err = _l2error(a, b)
    if (err >= epsilon):
        print("Error: {} > {}".format(err, epsilon))
    assert_(err<epsilon)


def nalm(lmax, mmax):
    return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)


26
27
28
def random_alm(rng, lmax, mmax, ncomp):
    res = rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \
     + 1j*rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp))
Martin Reinecke's avatar
Martin Reinecke committed
29
    # make a_lm with m==0 real-valued
30
    res[0:lmax+1,:].imag = 0.
Martin Reinecke's avatar
Martin Reinecke committed
31
32
33
34
    return res


def convolve(alm1, alm2, lmax):
Martin Reinecke's avatar
Martin Reinecke committed
35
    job = sht.sharpjob_d()
Martin Reinecke's avatar
Martin Reinecke committed
36
37
38
39
40
41
42
    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)


Martin Reinecke's avatar
Martin Reinecke committed
43
44
45
46
47
48
49
50
51
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):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
52
    return np.vdot(compress_alm(a1,lmax),compress_alm(np.conj(a2),lmax))
Martin Reinecke's avatar
Martin Reinecke committed
53
54


55
56
57
58
@pmp("lkmax", [(13,13),(2,1),(30,15),(35,2)])
@pmp("ncomp", [1, 3])
@pmp("separate", [True, False])
def test_against_convolution(lkmax, ncomp, separate):
Martin Reinecke's avatar
Martin Reinecke committed
59
    lmax, kmax = lkmax
60
61
62
    rng = np.random.default_rng(42)
    slm = random_alm(rng, lmax, lmax, ncomp)
    blm = random_alm(rng, lmax, kmax, ncomp)
Martin Reinecke's avatar
Martin Reinecke committed
63

Martin Reinecke's avatar
Martin Reinecke committed
64
    inter = totalconvolve.PyInterpolator(slm, blm, separate, lmax, kmax,
65
                                         epsilon=1e-8, nthreads=2)
Martin Reinecke's avatar
Martin Reinecke committed
66
67
    nptg = 50
    ptg = np.zeros((nptg,3))
68
69
70
    ptg[:,0] = rng.uniform(0, np.pi, nptg)
    ptg[:,1] = rng.uniform(0, 2*np.pi, nptg)
    ptg[:,2] = rng.uniform(-np.pi, np.pi, nptg)
Martin Reinecke's avatar
Martin Reinecke committed
71
72
73

    res1 = inter.interpol(ptg)

74
75
76
77
78
    blm2 = np.zeros((nalm(lmax,lmax), ncomp))+0j
    blm2[0:blm.shape[0],:] = blm
    res2 = np.zeros((nptg, ncomp))
    for c in range(ncomp):
        for i in range(nptg):
79
            rbeam=misc.rotate_alm(blm2[:,c], lmax, ptg[i,2],ptg[i,0],ptg[i,1])
80
81
82
83
84
85
86
87
88
89
            res2[i,c] = convolve(slm[:,c], rbeam, lmax).real
    if separate:
        _assert_close(res1, res2, 1e-7)
    else:
        _assert_close(res1[:,0], np.sum(res2,axis=1), 1e-7)

@pmp("lkmax", [(13,13),(2,1),(30,15),(35,2)])
@pmp("ncomp", [1, 3])
@pmp("separate", [True, False])
def test_adjointness(lkmax, ncomp, separate):
Martin Reinecke's avatar
Martin Reinecke committed
90
    lmax, kmax = lkmax
91
92
93
    rng = np.random.default_rng(42)
    slm = random_alm(rng, lmax, lmax, ncomp)
    blm = random_alm(rng, lmax, kmax, ncomp)
Martin Reinecke's avatar
Martin Reinecke committed
94
    nptg=100000
95
    ptg=rng.uniform(0.,1.,nptg*3).reshape(nptg,3)
Martin Reinecke's avatar
Martin Reinecke committed
96
97
98
    ptg[:,0]*=np.pi
    ptg[:,1]*=2*np.pi
    ptg[:,2]*=2*np.pi
Martin Reinecke's avatar
Martin Reinecke committed
99
    foo = totalconvolve.PyInterpolator(slm,blm,separate,lmax, kmax, epsilon=1e-6, nthreads=2)
Martin Reinecke's avatar
Martin Reinecke committed
100
    inter1=foo.interpol(ptg)
101
    ncomp2 = inter1.shape[1]
102
    fake = rng.uniform(0.,1., (ptg.shape[0],ncomp2))
Martin Reinecke's avatar
Martin Reinecke committed
103
    foo2 = totalconvolve.PyInterpolator(lmax, kmax, ncomp2, epsilon=1e-6, nthreads=2)
Martin Reinecke's avatar
Martin Reinecke committed
104
    foo2.deinterpol(ptg.reshape((-1,3)), fake)
105
    bla=foo2.getSlm(blm)
106
107
108
    v1 = np.sum([myalmdot(slm[:,c], bla[:,c], lmax, lmax, 0) for c in range(ncomp)])
    v2 = np.sum([np.vdot(fake[:,c],inter1[:,c]) for c in range(ncomp2)])
    _assert_close(v1, v2, 1e-12)