Commit 9d53de33 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

formatting

parent f4863980
Pipeline #76443 passed with stages
in 16 minutes and 1 second
......@@ -56,7 +56,7 @@ def measure_duccfft(a, nrepeat, nthr):
def measure_scipy_fftpack(a, nrepeat, nthr):
import scipy.fftpack
tmin = 1e38
if nthr!=1:
if nthr != 1:
raise NotImplementedError("scipy.fftpack does not support multiple threads")
for i in range(nrepeat):
t0 = time()
......@@ -78,9 +78,8 @@ def measure_scipy_fft(a, nrepeat, nthr):
def measure_numpy_fft(a, nrepeat, nthr):
import scipy.fft
tmin = 1e38
if nthr!=1:
if nthr != 1:
raise NotImplementedError("numpy.fft does not support multiple threads")
for i in range(nrepeat):
t0 = time()
......@@ -113,12 +112,12 @@ def bench_nd(ndim, nmax, nthr, ntry, tp, funcs, nrepeat, ttl="", filename="",
shp = np.array([duccfft.good_size(sz) for sz in shp])
print(" {0:4d}/{1}: shape={2} ...".format(n, ntry, shp), end=" ", flush=True)
a = (rng.random(shp)-0.5 + 1j*(rng.random(shp)-0.5)).astype(tp)
output=[]
output = []
for func, res in zip(funcs, results):
tmp = func(a, nrepeat, nthr)
res.append(tmp[0])
output.append(tmp[1])
print("{0:5.2e}/{1:5.2e} = {2:5.2f} L2 error={3}".format(results[0][n],results[1][n],results[0][n]/results[1][n],_l2error(output[0],output[1])))
print("{0:5.2e}/{1:5.2e} = {2:5.2f} L2 error={3}".format(results[0][n], results[1][n], results[0][n]/results[1][n], _l2error(output[0], output[1])))
results = np.array(results)
plt.title("{}: {}D, {}, max_extent={}".format(
ttl, ndim, str(tp), nmax))
......@@ -132,7 +131,7 @@ def bench_nd(ndim, nmax, nthr, ntry, tp, funcs, nrepeat, ttl="", filename="",
funcs = (measure_duccfft, measure_fftw)
ttl = "duccfft/FFTW()"
ntry=100
ntry = 100
nthr = 1
nice_sizes = True
bench_nd(1, 8192, nthr, ntry, "c16", funcs, 10, ttl, "1d.png", nice_sizes)
......
......@@ -6,7 +6,7 @@ rng = np.random.default_rng(42)
def _l2error(a, b, axes):
return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2))/np.log2(np.max([2,np.prod(np.take(a.shape,axes))]))
return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2))/np.log2(np.max([2, np.prod(np.take(a.shape, axes))]))
def fftn(a, axes=None, inorm=0, out=None, nthreads=1):
......
import time
import math
import numpy as np
import ducc0.healpix as ph
rng = np.random.default_rng(42)
def report (name,vlen,ntry,nside,isnest,perf):
print (name,": ",perf*1e-6,"MOps/s",sep="")
def report(name, vlen, ntry, nside, isnest, perf):
print(name, ": ", perf*1e-6, "MOps/s", sep="")
def random_ptg(vlen):
res=np.empty((vlen,2),dtype=np.float64)
res[:,0]=np.arccos((rng.random(vlen)-0.5)*2)
res[:,1]=rng.random(vlen)*2*math.pi
res = np.empty((vlen, 2), dtype=np.float64)
res[:, 0] = np.arccos((rng.random(vlen)-0.5)*2)
res[:, 1] = rng.random(vlen)*2*np.pi
return res
def random_pix(nside,vlen):
return rng.integers(low=0,high=12*nside*nside-1,size=vlen,dtype=np.int64)
def random_pix(nside, vlen):
return rng.integers(low=0, high=12*nside*nside-1, size=vlen, dtype=np.int64)
def dummy(vlen):
inp=np.zeros(vlen,dtype=np.int64)
_ = np.zeros(vlen, dtype=np.int64)
def genperf(func,fname,inp,vlen,ntry,nside,isnest):
cnt=0
t=time.time()
while (cnt<ntry):
def genperf(func, fname, inp, vlen, ntry, nside, isnest):
cnt = 0
t = time.time()
while cnt < ntry:
func(inp)
cnt+=1
t=time.time()-t
p=(vlen*ntry)/t
report (fname,vlen,ntry,nside,isnest,p)
def perf_pix2ang(vlen,ntry,nside,isnest):
inp=random_pix(nside,vlen)
base=ph.Healpix_Base (nside,"NEST" if isnest else "RING")
genperf(base.pix2ang,"pix2ang",inp,vlen,ntry,nside,isnest)
def perf_ang2pix(vlen,ntry,nside,isnest):
inp=random_ptg(vlen)
base=ph.Healpix_Base (nside,"NEST" if isnest else "RING")
genperf(base.ang2pix,"ang2pix",inp,vlen,ntry,nside,isnest)
def perf_pix2vec(vlen,ntry,nside,isnest):
inp=random_pix(nside,vlen)
base=ph.Healpix_Base (nside,"NEST" if isnest else "RING")
genperf(base.pix2vec,"pix2vec",inp,vlen,ntry,nside,isnest)
def perf_vec2pix(vlen,ntry,nside,isnest):
inp=ph.ang2vec(random_ptg(vlen))
base=ph.Healpix_Base (nside,"NEST" if isnest else "RING")
genperf(base.vec2pix,"vec2pix",inp,vlen,ntry,nside,isnest)
def perf_ring2nest(vlen,ntry,nside,isnest):
inp=random_pix(nside,vlen)
base=ph.Healpix_Base (nside,"NEST" if isnest else "RING")
genperf(base.ring2nest,"ring2nest",inp,vlen,ntry,nside,isnest)
def perf_nest2ring(vlen,ntry,nside,isnest):
inp=random_pix(nside,vlen)
base=ph.Healpix_Base (nside,"NEST" if isnest else "RING")
genperf(base.nest2ring,"nest2ring",inp,vlen,ntry,nside,isnest)
def perf_neighbors(vlen,ntry,nside,isnest):
inp=random_pix(nside,vlen)
base=ph.Healpix_Base (nside,"NEST" if isnest else "RING")
genperf(base.neighbors,"neighbors",inp,vlen,ntry,nside,isnest)
def suite (vlen,ntry,nside,isnest):
print ("vlen=",vlen,", ","NEST" if isnest else "RING",sep="")
cnt += 1
t = time.time()-t
p = (vlen*ntry)/t
report(fname, vlen, ntry, nside, isnest, p)
def perf_pix2ang(vlen, ntry, nside, isnest):
inp = random_pix(nside, vlen)
base = ph.Healpix_Base(nside, "NEST" if isnest else "RING")
genperf(base.pix2ang, "pix2ang", inp, vlen, ntry, nside, isnest)
def perf_ang2pix(vlen, ntry, nside, isnest):
inp = random_ptg(vlen)
base = ph.Healpix_Base(nside, "NEST" if isnest else "RING")
genperf(base.ang2pix, "ang2pix", inp, vlen, ntry, nside, isnest)
def perf_pix2vec(vlen, ntry, nside, isnest):
inp = random_pix(nside, vlen)
base = ph.Healpix_Base(nside, "NEST" if isnest else "RING")
genperf(base.pix2vec, "pix2vec", inp, vlen, ntry, nside, isnest)
def perf_vec2pix(vlen, ntry, nside, isnest):
inp = ph.ang2vec(random_ptg(vlen))
base = ph.Healpix_Base(nside, "NEST" if isnest else "RING")
genperf(base.vec2pix, "vec2pix", inp, vlen, ntry, nside, isnest)
def perf_ring2nest(vlen, ntry, nside, isnest):
inp = random_pix(nside, vlen)
base = ph.Healpix_Base(nside, "NEST" if isnest else "RING")
genperf(base.ring2nest, "ring2nest", inp, vlen, ntry, nside, isnest)
def perf_nest2ring(vlen, ntry, nside, isnest):
inp = random_pix(nside, vlen)
base = ph.Healpix_Base(nside, "NEST" if isnest else "RING")
genperf(base.nest2ring, "nest2ring", inp, vlen, ntry, nside, isnest)
def perf_neighbors(vlen, ntry, nside, isnest):
inp = random_pix(nside, vlen)
base = ph.Healpix_Base(nside, "NEST" if isnest else "RING")
genperf(base.neighbors, "neighbors", inp, vlen, ntry, nside, isnest)
def suite(vlen, ntry, nside, isnest):
print("vlen=", vlen, ", ", "NEST" if isnest else "RING", sep="")
dummy(vlen)
perf_pix2ang(vlen,ntry,nside,isnest)
perf_ang2pix(vlen,ntry,nside,isnest)
perf_pix2vec(vlen,ntry,nside,isnest)
perf_vec2pix(vlen,ntry,nside,isnest)
perf_neighbors(vlen,ntry,nside,isnest)
nside=512
ntry=1000
print ("nside=",nside,sep="")
for vlen in (1,10,100,1000,10000):
perf_pix2ang(vlen, ntry, nside, isnest)
perf_ang2pix(vlen, ntry, nside, isnest)
perf_pix2vec(vlen, ntry, nside, isnest)
perf_vec2pix(vlen, ntry, nside, isnest)
perf_neighbors(vlen, ntry, nside, isnest)
nside = 512
ntry = 1000
print("nside=", nside, sep="")
for vlen in (1, 10, 100, 1000, 10000):
for isnest in (True, False):
suite(vlen,ntry,nside,isnest)
perf_ring2nest(vlen,ntry,nside,isnest)
perf_nest2ring(vlen,ntry,nside,isnest)
suite(vlen, ntry, nside, isnest)
perf_ring2nest(vlen, ntry, nside, isnest)
perf_nest2ring(vlen, ntry, nside, isnest)
print()
import ducc0.healpix as ph
import numpy as np
import math
rng = np.random.default_rng(42)
def random_ptg(vlen):
res = np.empty((vlen, 2), dtype=np.float64)
res[:,0] = np.arccos((rng.random(vlen)-0.5)*2)
res[:,1] = rng.random(vlen)*2*math.pi
res[:, 0] = np.arccos((rng.random(vlen)-0.5)*2)
res[:, 1] = rng.random(vlen)*2*np.pi
return res
def check_pixangpix(vlen,ntry,nside,isnest):
base = ph.Healpix_Base (nside, "NEST" if isnest else "RING")
cnt = 0
while cnt < ntry:
cnt += 1
def check_pixangpix(vlen, ntry, nside, isnest):
base = ph.Healpix_Base(nside, "NEST" if isnest else "RING")
for cnt in range(ntry):
inp = rng.integers(low=0, high=12*nside*nside-1, size=vlen)
out = base.ang2pix(base.pix2ang(inp))
if not np.array_equal(inp, out):
raise ValueError("Test failed")
def check_vecpixvec(vlen, ntry, nside, isnest):
base = ph.Healpix_Base (nside, "NEST" if isnest else "RING")
cnt = 0
while cnt < ntry:
cnt += 1
base = ph.Healpix_Base(nside, "NEST" if isnest else "RING")
for cnt in range(ntry):
inp = ph.ang2vec(random_ptg(vlen))
out = base.pix2vec(base.vec2pix(inp))
if np.any(ph.v_angle(inp,out) > base.max_pixrad()):
if np.any(ph.v_angle(inp, out) > base.max_pixrad()):
raise ValueError("Test failed")
def check_pixangvecpix(vlen, ntry, nside, isnest):
base = ph.Healpix_Base (nside, "NEST" if isnest else "RING")
cnt = 0
while cnt < ntry:
cnt += 1
base = ph.Healpix_Base(nside, "NEST" if isnest else "RING")
for cnt in range(ntry):
inp = rng.integers(low=0, high=12*nside*nside-1, size=vlen)
out = base.vec2pix(ph.ang2vec(base.pix2ang(inp)))
if not np.array_equal(inp,out):
if not np.array_equal(inp, out):
raise ValueError("Test failed")
def check_pixvecangpix(vlen, ntry, nside, isnest):
base = ph.Healpix_Base (nside, "NEST" if isnest else "RING")
cnt = 0
while cnt < ntry:
cnt += 1
base = ph.Healpix_Base(nside, "NEST" if isnest else "RING")
for cnt in range(ntry):
inp = rng.integers(low=0, high=12*nside*nside-1, size=vlen)
out = base.ang2pix(ph.vec2ang(base.pix2vec(inp)))
if not np.array_equal(inp,out):
if not np.array_equal(inp, out):
raise ValueError("Test failed")
def check_pixvecpix(vlen,ntry,nside,isnest):
base=ph.Healpix_Base (nside,"NEST" if isnest else "RING")
cnt=0
while (cnt<ntry):
cnt+=1
inp=rng.integers(low=0,high=12*nside*nside-1,size=vlen)
out=base.vec2pix(base.pix2vec(inp))
if (np.array_equal(inp,out)==False):
def check_pixvecpix(vlen, ntry, nside, isnest):
base = ph.Healpix_Base(nside, "NEST" if isnest else "RING")
for cnt in range(ntry):
inp = rng.integers(low=0, high=12*nside*nside-1, size=vlen)
out = base.vec2pix(base.pix2vec(inp))
if not np.array_equal(inp, out):
raise ValueError("Test failed")
def check_ringnestring(vlen,ntry,nside):
base=ph.Healpix_Base (nside,"NEST")
cnt=0
while (cnt<ntry):
cnt+=1
inp=rng.integers(low=0,high=12*nside*nside-1,size=vlen)
out=base.nest2ring(base.ring2nest(inp))
if (np.array_equal(inp,out)==False):
def check_ringnestring(vlen, ntry, nside):
base = ph.Healpix_Base(nside, "NEST")
for cnt in range(ntry):
inp = rng.integers(low=0, high=12*nside*nside-1, size=vlen)
out = base.nest2ring(base.ring2nest(inp))
if not np.array_equal(inp, out):
raise ValueError("Test failed")
def check_pixxyfpix(vlen,ntry,nside,isnest):
base=ph.Healpix_Base (nside,"NEST" if isnest else "RING")
cnt=0
while (cnt<ntry):
cnt+=1
inp=rng.integers(low=0,high=12*nside*nside-1,size=vlen)
out=base.xyf2pix(base.pix2xyf(inp))
if (np.array_equal(inp,out)==False):
def check_pixxyfpix(vlen, ntry, nside, isnest):
base = ph.Healpix_Base(nside, "NEST" if isnest else "RING")
for cnt in range(ntry):
inp = rng.integers(low=0, high=12*nside*nside-1, size=vlen)
out = base.xyf2pix(base.pix2xyf(inp))
if not np.array_equal(inp, out):
raise ValueError("Test failed")
def check_vecangvec(vlen,ntry):
cnt=0
while (cnt<ntry):
cnt+=1
inp=random_ptg(vlen)
out=ph.vec2ang(ph.ang2vec(inp))
if (np.any(np.greater(np.abs(inp-out),1e-10))):
def check_vecangvec(vlen, ntry):
for cnt in range(ntry):
inp = random_ptg(vlen)
out = ph.vec2ang(ph.ang2vec(inp))
if np.any(np.greater(np.abs(inp-out), 1e-10)):
raise ValueError("Test failed")
check_vecangvec(1000,1000)
for nside in (1,32,512,8192,32768*8):
check_ringnestring(1000,1000,nside)
for isnest in (False,True):
check_vecpixvec(1000,1000,nside,isnest)
check_pixangpix(1000,1000,nside,isnest)
check_pixvecpix(1000,1000,nside,isnest)
check_pixxyfpix(1000,1000,nside,isnest)
check_pixangvecpix(1000,1000,nside,isnest)
check_pixvecangpix(1000,1000,nside,isnest)
isnest=False
for nside in (3,7,514,8167,32768*8+7):
check_vecpixvec(1000,1000,nside,isnest)
check_pixangpix(1000,1000,nside,isnest)
check_pixvecpix(1000,1000,nside,isnest)
check_pixxyfpix(1000,1000,nside,isnest)
check_pixangvecpix(1000,1000,nside,isnest)
check_pixvecangpix(1000,1000,nside,isnest)
check_vecangvec(1000, 1000)
for nside in (1, 32, 512, 8192, 32768*8):
check_ringnestring(1000, 1000, nside)
for isnest in (False, True):
check_vecpixvec(1000, 1000, nside, isnest)
check_pixangpix(1000, 1000, nside, isnest)
check_pixvecpix(1000, 1000, nside, isnest)
check_pixxyfpix(1000, 1000, nside, isnest)
check_pixangvecpix(1000, 1000, nside, isnest)
check_pixvecangpix(1000, 1000, nside, isnest)
isnest = False
for nside in (3, 7, 514, 8167, 32768*8+7):
check_vecpixvec(1000, 1000, nside, isnest)
check_pixangpix(1000, 1000, nside, isnest)
check_pixvecpix(1000, 1000, nside, isnest)
check_pixxyfpix(1000, 1000, nside, isnest)
check_pixangvecpix(1000, 1000, nside, isnest)
check_pixvecangpix(1000, 1000, nside, isnest)
......@@ -8,9 +8,11 @@ import ducc0.sht as sht
import numpy as np
from time import time
def _l2error(a, b):
return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2))
# set maximum multipole moment
lmax = 2047
# maximum m. For SHTOOLS this is alway equal to lmax, if I understand correctly.
......@@ -53,7 +55,7 @@ nlat = lmax+1
job.set_gauss_geometry(nlat, nlon)
# go from a_lm to map
t0=time()
t0 = time()
map = job.alm2map(alm)
print("time for map synthesis: {}s".format(time()-t0))
......@@ -63,12 +65,12 @@ print("time for map synthesis: {}s".format(time()-t0))
# number of pixels on each iso-latitude ring, which cannot be represented by 2D
# arrays (e.g. Healpix)
t0=time()
t0 = time()
alm2 = job.map2alm(map)
print("time for map analysis: {}s".format(time()-t0))
# make sure input was recovered accurately
print("L2 error: ", _l2error(alm,alm2))
print("L2 error: ", _l2error(alm, alm2))
print("testing Driscoll-Healy grid")
......@@ -80,7 +82,7 @@ nlat = 2*lmax+2
job.set_dh_geometry(nlat, nlon)
# go from a_lm to map
t0=time()
t0 = time()
map = job.alm2map(alm)
print("time for map synthesis: {}s".format(time()-t0))
......@@ -90,9 +92,9 @@ print("time for map synthesis: {}s".format(time()-t0))
# number of pixels on each iso-latitude ring, which cannot be represented by 2D
# arrays (e.g. Healpix)
t0=time()
t0 = time()
alm2 = job.map2alm(map)
print("time for map analysis: {}s".format(time()-t0))
# make sure input was recovered accurately
print("L2 error: ", _l2error(alm,alm2))
print("L2 error: ", _l2error(alm, alm2))
......@@ -3,9 +3,11 @@ import ducc0.misc as misc
import numpy as np
from time import time
def _l2error(a, b):
return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2))
lmax = 1023
mmax = lmax
......@@ -30,20 +32,20 @@ nlat = lmax+1
print("Converting them to Fejer1 grid with {}x{} points".format(nlat, nlon))
job.set_fejer1_geometry(nlat, nlon)
t0=time()
t0 = time()
map = job.alm2map(alm)
print("time for map synthesis: {}s".format(time()-t0))
nlat2 = 2*lmax+3
t0=time()
map2 = misc.upsample_to_cc(map.reshape((nlat,nlon)), nlat2, False, False)
t0 = time()
map2 = misc.upsample_to_cc(map.reshape((nlat, nlon)), nlat2, False, False)
print("time for upsampling: {}s".format(time()-t0))
job.set_cc_geometry(nlat2, nlon)
t0=time()
t0 = time()
alm2 = job.map2alm(map2.reshape((-1,)))
print("time for map analysis: {}s".format(time()-t0))
# make sure input was recovered accurately
print("L2 error: ", _l2error(alm,alm2))
print("L2 error: ", _l2error(alm, alm2))
nlon = 2*mmax+2
nlat = lmax+2
......@@ -52,20 +54,20 @@ nlat = lmax+2
print("Converting them to Clenshaw-Curtis grid with {}x{} points".format(nlat, nlon))
job.set_cc_geometry(nlat, nlon)
t0=time()
t0 = time()
map = job.alm2map(alm)
print("time for map synthesis: {}s".format(time()-t0))
nlat2 = 2*lmax+3
t0=time()
map2 = misc.upsample_to_cc(map.reshape((nlat,nlon)), nlat2, True, True)
t0 = time()
map2 = misc.upsample_to_cc(map.reshape((nlat, nlon)), nlat2, True, True)
print("time for upsampling: {}s".format(time()-t0))
job.set_cc_geometry(nlat2, nlon)
t0=time()
t0 = time()
alm2 = job.map2alm(map2.reshape((-1,)))
print("time for map analysis: {}s".format(time()-t0))
# make sure input was recovered accurately
print("L2 error: ", _l2error(alm,alm2))
print("L2 error: ", _l2error(alm, alm2))
nlon = 2*mmax+2
......@@ -75,17 +77,17 @@ nlat = lmax+1
print("Converting them to McEwen-Wiaux grid with {}x{} points".format(nlat, nlon))
job.set_mw_geometry(nlat, nlon)
t0=time()
t0 = time()
map = job.alm2map(alm)
print("time for map synthesis: {}s".format(time()-t0))
nlat2 = 2*lmax+3
t0=time()
map2 = misc.upsample_to_cc(map.reshape((nlat,nlon)), nlat2, False, True)
t0 = time()
map2 = misc.upsample_to_cc(map.reshape((nlat, nlon)), nlat2, False, True)
print("time for upsampling: {}s".format(time()-t0))
job.set_cc_geometry(nlat2, nlon)
t0=time()
t0 = time()
alm2 = job.map2alm(map2.reshape((-1,)))
print("time for map analysis: {}s".format(time()-t0))
# make sure input was recovered accurately
print("L2 error: ", _l2error(alm,alm2))
print("L2 error: ", _l2error(alm, alm2))
......@@ -7,6 +7,7 @@ import matplotlib.pyplot as plt
rng = np.random.default_rng(42)
def nalm(lmax, mmax):
return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
......@@ -15,11 +16,11 @@ def random_alm(lmax, mmax, ncomp):
res = rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \
+ 1j*rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp))
# make a_lm with m==0 real-valued
res[0:lmax+1,:].imag = 0.
res[0:lmax+1, :].imag = 0.
return res
def compress_alm(alm,lmax):
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
......@@ -27,8 +28,8 @@ def compress_alm(alm,lmax):
return res
def myalmdot(a1,a2,lmax,mmax,spin):
return np.vdot(compress_alm(a1,lmax),compress_alm(np.conj(a2),lmax))
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):
......@@ -36,14 +37,14 @@ def convolve(alm1, alm2, lmax):
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)
job.set_triangular_alm_info(0, 0)