Multiply alms by function of ell

[msyriac]

Does enlib.lensing or enlib.curvedsky have the healpy equivalent of almxfl?
http://healpy.readthedocs.io/en/latest/generated/healpy.sphtfunc.almxfl.html

This would be useful. The specific use case I have in mind is getting kappa maps instead of phi maps from tenki. Because phi maps have very red spectra, it is hard to work with non-periodic part-sky cut-outs of them.

EDIT: See this PR for phi->kappa
#17

[msyriac]

I used lmul from sharp.pyx here:
37070d8#diff-085ebf46c806d5b76657cd7b7d8d9433R109

Maybe it's worth packaging that into a function in curvedsky.py?

[amaurea]

Yes, that might be a good idea. Though lmul is pretty simple to use already, as long as you have the alm_info.

[ajvanengelen]

I almost submitted a new issue when I saw this existing thread.

This is currently giving an out-of-memory error for me. The reason is that ainfo.lmul() calls np.eye() which generates an NxM matrix -- which is too big to store in memory. Here N and M are both huge (for lmax = 8000 they are both 32012001). I am looking into using the healpy tools to multiply by a function of ell instead. (Not sure whether the l,m orderings are the same in the two packages.)

[msyriac]

Hi @ajvanengelen , are you encountering memory issues in lensing.rand_map ? If yes, you can use the delta_theta argument that Sigurd added recently, that does the map-generation and lensing in chunks of declination.

https://github.com/amaurea/enlib/blob/master/lensing.py#L92

[ajvanengelen]

Hi Mat, It is just the phi to kappa conversion.

[ajvanengelen]

Here is my backtrace

MemoryError Traceback (most recent call last)
/global/u2/e/engelen/simsS1516/py/signalSims.py in ()
75 print 'rank', rank, 'doing iii' , iii, ', iMin', iMin, ', iMax', iMax
76
---> 77 uTquMap, lTquMap, pMap, kappaMap = lensing.rand_map((3,)+shape, wcs, ps, lmax=p['LMAX'], output="ulpk", verbose=True)
78 uTebMap, uTebCls[iii] = tqu2teb(uTquMap,p['LMAX_NYQ'], wantCl = True)
79 lTebMap, lTebCls[iii] = tqu2teb(lTquMap,p['LMAX_NYQ'], wantCl = True)

/global/homes/e/engelen/enlib/lensing.pyc in rand_map(shape, wcs, ps_lensinput, lmax, maplmax, dtype, seed, oversample, spin, output, geodesic, verbose, delta_theta)
116 kappa_map = enmap.empty(shape[-2:], wcs, dtype=dtype)
117 l = np.arange(ainfo.lmax+1.0)
--> 118 kappa_alm = ainfo.lmul(phi_alm, l*(l+1)/2)
119 for i1 in range(0, shape[-2], bsize):
120 curvedsky.alm2map(kappa_alm, kappa_map[...,i1:i1+bize,:])

/global/homes/e/engelen/enlib/sharp/sharp.pyx in enlib.sharp.sharp.alm_info.lmul (sharp.c:7333)()
258 if out is None: out = alm.copy()
259 if out.ndim == 1: out = out[None]
--> 260 if lmat.ndim == 1:
261 lmat = np.eye(out.shape[0])[:,:,None]*lmat
262 lmat = lmat.astype(out.real.dtype, copy=False)

/global/common/cori/software/python/2.7-anaconda/lib/python2.7/site-packages/numpy/lib/twodim_base.pyc in eye(N, M, k, dtype)
231 if M is None:
232 M = N
--> 233 m = zeros((N, M), dtype=dtype)
234 if k >= M:
235 return m

MemoryError:

[amaurea]

This is currently giving an out-of-memory error for me. The reason is that ainfo.lmul() calls np.eye() which generates an NxM matrix -- which is too big to store in memory. Here N and M are both huge (for lmax = 8000 they are both 32012001).

What is the shape of the arrays you are passing as arguments to lmul? The error you are getting here indicates that the wrong axis of one of the arrays is being interpreted as the number of stokes components.

[amaurea]

Hi @ajvanengelen , are you encountering memory issues in lensing.rand_map? If yes, you can use the delta_theta argument that Sigurd added recently, that does the map-generation and lensing in chunks of declination.

The error he is getting won't be solved by that. Looking at the code, for his error to happen he must be passing in a 1-d scaling function together with an alm array which is not 1d, but which has lm as the first axis. I don't see how that could happen. Maybe he has an old version?

[ajvanengelen]

Hi Sigurd, It is the current version. I am not calling lmul on my own, but rather am just asking for 'k' as one of the outputs for this routine. Cori is down right now so I can't answer your question about the shapes of the arrays that are getting passed down. Alex

…

On 12/19/2017 04:53 PM, amaurea wrote: > Hi @ajvanengelen , are you encountering memory issues in lensing.rand_map? >If yes, you can use the delta_theta argument that Sigurd added recently, that >does the map-generation and lensing in chunks of declination. The error he is getting won't be solved by that. Looking at the code, for his error to happen he must be passing in a 1-d scaling function together with an alm array which is not 1d, but which has lm as the first axis. I don't see how that could happen. Maybe he has an old version? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#16 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/APAXHLCXOpUyecbSIE4ZEXLlW-Vutikdks5tCDBsgaJpZM4PLW3->.

[ajvanengelen]

OK, in line 118 of lensing.py, the two arguments to lmul -- phi_alm and l(l+1)/2, have size 32012001 and 8001 respectively Alex

…

On 12/19/2017 05:26 PM, Alexander van Engelen wrote: Hi Sigurd, It is the current version. I am not calling lmul on my own, but rather am just asking for 'k' as one of the outputs for this routine. Cori is down right now so I can't answer your question about the shapes of the arrays that are getting passed down. Alex On 12/19/2017 04:53 PM, amaurea wrote: > > Hi @ajvanengelen , are you encountering memory issues in > lensing.rand_map? > >If yes, you can use the delta_theta argument that Sigurd added > recently, that > >does the map-generation and lensing in chunks of declination. > > The error he is getting won't be solved by that. Looking at the code, for > his error to happen he must be passing in a 1-d scaling function together > with an alm array which is not 1d, but which has lm as the first > axis. I don't > see how that could happen. Maybe he has an old version? > > — > You are receiving this because you were mentioned. > Reply to this email directly, view it on GitHub > <#16 (comment)>, > or mute the thread > <https://github.com/notifications/unsubscribe-auth/APAXHLCXOpUyecbSIE4ZEXLlW-Vutikdks5tCDBsgaJpZM4PLW3->. >