Aquila-Consortium/libsharp2
0
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

improve cython performance, declare types

Browse source
This commit is contained in:
Andrea Zonca 2018-01-16 08:30:54 -08:00
parent 41668b72c3
commit 158323aced
1 changed files with 15 additions and 7 deletions
Download patch file Download diff file Expand all files Collapse all files

18
python/libsharp/libsharp.pyx
View file

@ -247,6 +247,7 @@ cdef class alm_info:
sharp_destroy_alm_info(self.ainfo) sharp_destroy_alm_info(self.ainfo)
self.ainfo = NULL self.ainfo = NULL
@cython.boundscheck(False)
def almxfl(self, np.ndarray[double, ndim=3, mode='c'] alm, np.ndarray[double, ndim=2, mode='c'] fl): def almxfl(self, np.ndarray[double, ndim=3, mode='c'] alm, np.ndarray[double, ndim=2, mode='c'] fl):
"""Multiply Alm by a Ell based array """Multiply Alm by a Ell based array
@ -263,22 +264,29 @@ cdef class alm_info:
None, it modifies alms in-place None, it modifies alms in-place
""" """
mvstart = 0 cdef int mvstart = 0
has_multiple_beams = alm.shape[2] > 1 and fl.shape[1] > 1 cdef bint has_multiple_beams = alm.shape[2] > 1 and fl.shape[1] > 1
for m in self.mval(): cdef int f, i_m, m, num_ells, i_l, i_signal, i_pol, i_mv
for i_m in range(self.ainfo.nm):
m = self.ainfo.mval[i_m]
f = 1 if (m==0) else 2 f = 1 if (m==0) else 2
num_ells = self.ainfo.lmax + 1 - m num_ells = self.ainfo.lmax + 1 - m
if not has_multiple_beams: if not has_multiple_beams:
for i_signal in range(alm.shape[0]):
for i_pol in range(alm.shape[1]):
for i_l in range(num_ells): for i_l in range(num_ells):
l = m + i_l l = m + i_l
alm[:,:,mvstart + f*i_l:mvstart + f*i_l +f] *= fl[l] for i_mv in range(mvstart + f*i_l, mvstart + f*i_l +f):
alm[i_signal, i_pol, i_mv] *= fl[l, 0]
else: else:
for i_signal in range(alm.shape[0]): for i_signal in range(alm.shape[0]):
for i_pol in range(alm.shape[1]): for i_pol in range(alm.shape[1]):
for i_l in range(num_ells): for i_l in range(num_ells):
l = m + i_l l = m + i_l
alm[i_signal, i_pol, mvstart + f*i_l:mvstart + f*i_l +f] *= fl[l, i_pol] for i_mv in range(mvstart + f*i_l, mvstart + f*i_l +f):
alm[i_signal, i_pol, i_mv] *= fl[l, i_pol]
mvstart += f * num_ells mvstart += f * num_ells
cdef class triangular_order(alm_info): cdef class triangular_order(alm_info):