From fbc86507b86f54a32d4ae9e6c2baf539c68596ad Mon Sep 17 00:00:00 2001 From: Christopher Wallis Date: Tue, 18 Jul 2017 17:40:32 +0100 Subject: [PATCH] unit tests added --- .gitignore | 1 + data/make_DES_cat.py | 2 +- src/python/cy_mass_mapping.pyx | 12 ++- src/python/mass_mapping.py | 2 +- tests/massmappy_test.py | 173 +++++++++++++++++++++++++++++++++ tests/massmappy_test.pyc | Bin 0 -> 7157 bytes 6 files changed, 184 insertions(+), 6 deletions(-) create mode 100644 tests/massmappy_test.py create mode 100644 tests/massmappy_test.pyc diff --git a/.gitignore b/.gitignore index 320df5f..6f3d058 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ src/python/*.c src/python/*.so data/DES.mat +build/ diff --git a/data/make_DES_cat.py b/data/make_DES_cat.py index d012520..981e068 100644 --- a/data/make_DES_cat.py +++ b/data/make_DES_cat.py @@ -5,7 +5,7 @@ # Download data from here https://des.ncsa.illinois.edu/releases/sva1/content # You will need sva1_gold_r1.0_wlinfo.fits.gz, sva1_gold_r1.1_im3shape.fits.gz and sva1_gold_r1.0_bpz_point.fits.gz -# We would like to thank Chihway Chang and Joe Zunts for their help with constructing tha catalogue +# We would like to thank Chihway Chang and Joe Zunts for their help with constructing the catalogue info_check = os.path.exists('data/sva1_gold_r1.0_wlinfo.fits.gz') im3shape_check = os.path.exists('data/sva1_gold_r1.1_im3shape.fits.gz') diff --git a/src/python/cy_mass_mapping.pyx b/src/python/cy_mass_mapping.pyx index b2efa70..54198e8 100644 --- a/src/python/cy_mass_mapping.pyx +++ b/src/python/cy_mass_mapping.pyx @@ -66,12 +66,14 @@ def generate_kappa_lm_hp(np.ndarray[double, ndim=1, mode="c"] Cl not None, int L cdef np.ndarray[complex, ndim=1] k_lm cdef int el, em, index - k_lm = np.empty(((L*L+L)/2,), dtype=complex) + k_lm = np.empty((L*(L+1)/2,), dtype=complex) if seed > -1: np.random.seed(seed) - k_lm[0] = 0.0; k_lm[1] = 0.0; k_lm[2] = 0.0; k_lm[3] = 0.0; + k_lm[mm.cy_healpy_lm2ind(0, 0, L)] = 0.0 + k_lm[mm.cy_healpy_lm2ind(1, 0, L)] = 0.0 + k_lm[mm.cy_healpy_lm2ind(1, 1, L)] = 0.0 for el in range(2,L): index = mm.cy_healpy_lm2ind(el, 0, L) @@ -139,8 +141,10 @@ def kappa_lm_to_gamma_lm_hp(np.ndarray[double complex, ndim=1, mode="c"] k_lm no gamma_E_lm = np.empty((L*(L+1)/2,), dtype=complex) gamma_B_lm = np.zeros((L*(L+1)/2,), dtype=complex) - gamma_E_lm[0] = 0.0; gamma_E_lm[1] = 0.0; gamma_E_lm[2] = 0.0; gamma_E_lm[3] = 0.0 - + gamma_E_lm[mm.cy_healpy_lm2ind(0, 0, L)] = 0.0; + gamma_E_lm[mm.cy_healpy_lm2ind(1, 0, L)] = 0.0; + gamma_E_lm[mm.cy_healpy_lm2ind(1, 1, L)] = 0.0; + for ell in range(2,L): D_ell = sqrt((ell+2.0)*(ell-1.0)/((ell+1.0)*ell)) for em in range(0,ell+1): diff --git a/src/python/mass_mapping.py b/src/python/mass_mapping.py index 3d87cb6..b4c570a 100644 --- a/src/python/mass_mapping.py +++ b/src/python/mass_mapping.py @@ -12,7 +12,7 @@ Method = "MW" Add_noise = False -Reduce_shear = False +Reduce_shear = True Cls = np.loadtxt("data/cls_ap.txt") diff --git a/tests/massmappy_test.py b/tests/massmappy_test.py new file mode 100644 index 0000000..7f84db6 --- /dev/null +++ b/tests/massmappy_test.py @@ -0,0 +1,173 @@ +import numpy as np +import cy_mass_mapping as mm +import cy_healpy_mass_mapping as hp_mm +import pyssht as ssht + + +def test_lm_conversion_functions(): + ''' To test the lm conversion functions + ''' + + L = 3 + + flm_mw = np.array([0., -1+1j, 1+0j, 1+1j, 2-2j, -2+1j, 2+0j, 2+1j, 2+2j]) + flm_hp = np.array([0., 1+0j, 2+0j, 1+1j, 2+1j, 2+2j]) + + flm_hp_test = mm.lm2lm_hp(flm_mw, 3) + flm_mw_test = mm.lm_hp2lm(flm_hp, 3) + + np.testing.assert_almost_equal(flm_hp_test, flm_hp, decimal=7, err_msg='mm.lm2lm_hp failed to convert properly', verbose=True) + np.testing.assert_almost_equal(flm_mw_test, flm_mw, decimal=7, err_msg='mm.lm2lm_mw failed to convert properly', verbose=True) + + +def test_hp_index2lm(): + ''' To test the MW index functions + ''' + + L = 3 + + + np.testing.assert_equal(mm.healpy_ind2lm(0,L)[0],0,err_msg='hp_mm.healpy_ind2lm failed to get the correct ell') + np.testing.assert_equal(mm.healpy_ind2lm(0,L)[1],0,err_msg='hp_mm.healpy_ind2lm failed to get the correct em') + + np.testing.assert_equal(mm.healpy_ind2lm(1,L)[0],1,err_msg='hp_mm.healpy_ind2lm failed to get the correct ell') + np.testing.assert_equal(mm.healpy_ind2lm(1,L)[1],0,err_msg='hp_mm.healpy_ind2lm failed to get the correct em') + + np.testing.assert_equal(mm.healpy_ind2lm(2,L)[0],2,err_msg='hp_mm.healpy_ind2lm failed to get the correct ell') + np.testing.assert_equal(mm.healpy_ind2lm(2,L)[1],0,err_msg='hp_mm.healpy_ind2lm failed to get the correct em') + + np.testing.assert_equal(mm.healpy_ind2lm(3,L)[0],1,err_msg='hp_mm.healpy_ind2lm failed to get the correct ell') + np.testing.assert_equal(mm.healpy_ind2lm(3,L)[1],1,err_msg='hp_mm.healpy_ind2lm failed to get the correct em') + + np.testing.assert_equal(mm.healpy_ind2lm(4,L)[0],2,err_msg='hp_mm.healpy_ind2lm failed to get the correct ell') + np.testing.assert_equal(mm.healpy_ind2lm(4,L)[1],1,err_msg='hp_mm.healpy_ind2lm failed to get the correct em') + + np.testing.assert_equal(mm.healpy_ind2lm(5,L)[0],2,err_msg='hp_mm.healpy_ind2lm failed to get the correct ell') + np.testing.assert_equal(mm.healpy_ind2lm(5,L)[1],2,err_msg='hp_mm.healpy_ind2lm failed to get the correct em') + +def test_hp_lm2index(): + + L=3 + + np.testing.assert_equal(mm.healpy_lm2ind(0,0,L),0,err_msg='hp_mm.healpy_lm2ind failed to get the correct index') + + np.testing.assert_equal(mm.healpy_lm2ind(1,0,L),1,err_msg='hp_mm.healpy_lm2ind failed to get the correct index') + + np.testing.assert_equal(mm.healpy_lm2ind(2,0,L),2,err_msg='hp_mm.healpy_lm2ind failed to get the correct index') + + np.testing.assert_equal(mm.healpy_lm2ind(1,1,L),3,err_msg='hp_mm.healpy_lm2ind failed to get the correct index') + + np.testing.assert_equal(mm.healpy_lm2ind(2,1,L),4,err_msg='hp_mm.healpy_lm2ind failed to get the correct index') + + np.testing.assert_equal(mm.healpy_lm2ind(2,2,L),5,err_msg='hp_mm.healpy_lm2ind failed to get the correct index') + +def test_generate_kappa_lm_hp(): + + L = 3 + Cl = np.array([1E2,1E4,1E6]) + seed = 1 + mm.generate_kappa_lm_hp(Cl, L, seed=seed) + + kappa_prep = np.array([0.0+0j, 0.0+0j,1624.34536366+0j, 0.0+0j,\ + -611.75641365-528.17175226j,-1072.96862216+865.40762932j]) + + kappa_gen = mm.generate_kappa_lm_hp(Cl, L, seed=seed) + + np.testing.assert_almost_equal(kappa_gen, kappa_prep, decimal=7, err_msg='mm.generate_kappa_lm_hp failed to generate standard output', verbose=True) + + +def test_generate_kappa_lm_mw(): + + L = 3 + Cl = np.array([1E2,1E4,1E6]) + seed = 1 + + kappa_gen = mm.generate_kappa_lm_mw(Cl, L, seed=seed) + + kappa_prep = np.array([0.0+0j, 0.0+0j, 0.0+0j, 0.0+0j,\ + -1072.96862216-865.40762932j, 611.75641365-528.17175226j,\ + 1624.34536366 +0j, -611.75641365-528.17175226j,\ + -1072.96862216+865.40762932j]) + + np.testing.assert_almost_equal(kappa_gen, kappa_prep, decimal=7, err_msg='mm.generate_kappa_lm_hp failed to generate standard output', verbose=True) + +def test_lm_transform_mw(): + + L = 20 + Cl = np.ones(L) + seed = 1 + + kappa_lm = mm.generate_kappa_lm_mw(Cl, L, seed=seed) + + gamma_lm = mm.kappa_lm_to_gamma_lm_mw(kappa_lm, L) + + kappa_lm_rec = mm.gamma_lm_to_kappa_lm_mw(gamma_lm, L) + + np.testing.assert_almost_equal(kappa_lm_rec, kappa_lm, decimal=7, err_msg="lm gamma and kappa MW transoforms are not invereses") + +def test_lm_transform_hp(): + + L = 20 + Cl = np.ones(L) + seed = 1 + + kappa_lm = mm.generate_kappa_lm_hp(Cl, L, seed=seed) + + gamma_lm_E, gamma_lm_B = mm.kappa_lm_to_gamma_lm_hp(kappa_lm, L) + + kappa_lm_E_rec, kappa_lm_B_rec = mm.gamma_lm_to_kappa_lm_hp(gamma_lm_E, gamma_lm_B, L) + + np.testing.assert_almost_equal(kappa_lm_E_rec, kappa_lm, decimal=7, err_msg="lm gamma and kappa HEALPix transoforms are not invereses (E)") + np.testing.assert_almost_equal(kappa_lm_B_rec+1.0, np.ones(kappa_lm.size), decimal=7, err_msg="lm gamma and kappa HEALPix transoforms are not invereses (B)") + +def test_gamma_transform_mw(): + + L = 20 + Cl = np.ones(L) + seed = 1 + Method = "MW" + + kappa_lm = mm.generate_kappa_lm_mw(Cl, L, seed=seed) + + k_mw = ssht.inverse(kappa_lm, L, Reality=True, Method=Method) + + gamma_lm = mm.kappa_lm_to_gamma_lm_mw(kappa_lm, L) + + gamma = ssht.inverse(gamma_lm, L, Method=Method, Spin=2) + + k_rec_mw = mm.gamma_to_kappa_mw(gamma, L, Method=Method) + + np.testing.assert_almost_equal(k_rec_mw, k_mw, decimal=7, err_msg="gamma and kappa transoforms are not invereses") + +def test_reduced_shear_transform_mw(): + + L = 20 + Cl = np.ones(L)*1E-4 + seed = 1 + Method = "MW" + + kappa_lm = mm.generate_kappa_lm_mw(Cl, L, seed=seed) + + k_mw = ssht.inverse(kappa_lm, L, Reality=True, Method=Method) + + gamma_lm = mm.kappa_lm_to_gamma_lm_mw(kappa_lm, L) + + gamma = ssht.inverse(gamma_lm, L, Method=Method, Spin=2) + + shear = gamma/(1.0-k_mw) + + k_rec_mw = mm.reduced_shear_to_kappa_mw(shear, L, Iterate=True, Method=Method, tol_error=1E-10) + + np.testing.assert_almost_equal(k_rec_mw, k_mw, decimal=7, err_msg="reduced shear and kappa transoforms are not invereses") + + +if __name__ == '__main__': + test_lm_conversion_functions() #1 + test_hp_index2lm() #2 + test_hp_lm2index() #3 + test_generate_kappa_lm_hp() #4 + test_generate_kappa_lm_mw() #5 + test_lm_transform_mw() #6 + test_lm_transform_hp() #7 + test_gamma_transform_mw() #8 + test_reduced_shear_transform_mw() #9 \ No newline at end of file diff --git a/tests/massmappy_test.pyc b/tests/massmappy_test.pyc new file mode 100644 index 0000000000000000000000000000000000000000..75d95bd8db8257a1047131f014eda503e880674c GIT binary patch literal 7157 zcmc&&TW=&s6|U~N*yDS=H+ygCL`yVA$r=_25(@0>E-Ml#Yod0dD65E8XS#Q$z0=b} zcNy;lTMF6;kdT1aNIb(Y03k(DEcQ9#ArA;45eeSHizOr;5DI+fR9CrY?3vwU<7Mow zuC6|P>QvQtzB;E$|D2t=a-sWLQ{+Dz_Xe))Llgo2EzuF;rlsq)xLMNmvbb5%^{Tiz zAp%>x4dF=S_DkO5LW_ZART@)oo7Q=GAS1Z(SbCdVL1tivEk=8jrbx$7x*GaAo&U zV8A;bW?KZXU{M6eL~vXLCq!^k1gAu>Bw8>7l2$=-Sp=s=a7F}YMQ~2Esw`j& zX^V}L*oHPs+_l9{iCfBI8_HEnnUz$;MpbM>6785T%oOGbvoj~;1U@WTWRnuY$53a4 zFkwb0SP<*=^E8;JP&h9qjtn#D$ad&-VyBsOe-Nfwlyse2o86{FC7Y^6C5X0g@sB5` zubHlj~3y)YV+88)NX??hDRhA>UNIBT(SaeT28U+TnOyXV~U zqfQt&GRaMpPA^S*VcO~QXE6GWBnt;Wi?_ynww7D9_ddGDZ||~hxU3qC-Ftb1kzo9j z(5@cYuKMJlgc*Bdltg1xx26}ZW27AiyKU4b(I0>%M7Jji!hV|ieQt;UVg_c>4o3l^gFVQ_LtMmuugg-uLTzwdDrSI=qp=oZuy-~l&yXn>-$D>XEjTkba^psLDB2MwpCJ> ztx^-U^jY1-USFPt9(w1Y*L^H|2FrxAYz~DeSM4uaPuPprqWuM{Ni8%I8WP!({R>Y; zPv3x_mc&MxTz*@SpQ5r%{tvb*4@ddo9z(TE1QEf9x)37O>5K!k&zx3%pjhS2{ za>s~pyT%Ub5?8zI{@*5E zVQIx$u@7Dc9rE80ALWp_PXguR{1ISOz6-f=2(A|7P#@_#(4ia>64-4=%0Y*b8;5hq zVSfj=KodtuIjx9!z7tOygL1}hD%B~>{8VA4LpLkT46U#eHnVUv90>Fiw?wI`n8iuzZ>{z;3S)} zw<+Dj!aN-@9@zzMkaF@^{uZ>GP zv`b)pv_24TeHVRZizsp*0Qwi<2~}%Gc>yIG!~*OEZ=zE40@^@1?*{n~Q$6R-h`{$# zY_xDNPb`i!Tv51csj!1o#Rg(zwDa~(aRat#O1Evr=oDzJMZ`Z$+iK?UROD9LrV$0U z{5bX<@TSA;;vk}wso%|#TS*#cj-Q53H<5@kk<5n3WjTjqrsK^Z>4q7zoqHDgBn451 zq-El2@i2B1OI zs*r<#Z&I96R0|Ak;EW2fq7la-kwy~jLia<`bC5X%NzY%u^5RR;o$+vT>ep9tzTLZ= zJ)EYiD;$ujZI}*|6fg;Qqk)M7fsw#$_x2e{g+9@np>hD)@AeVuQewtCwnnCReem(B zr2O7M^*Vudu<7a}#c#WZ^U?i1UEz1hK4U*_F@h=FFbMtw9>yRTNk-0HCJyCuLgF+q zm%u2MIQOHzIt zjO9K_x+phcU5Zu-b$E9K#L0d~r(U_`w5#OY`9s){5wX%<3hXg^o(FrL*D)pW4vmU; z(~nU(h(Bzl$8nNP=tF`HZXC*mK#!b44Py`xu4B zE&)#M!+OX3M*qaj&n!sd=HLGA@Ao~GjKRM3{2+t&5t v5VoP=D5jBOqUYiy*zAPg;=@I