Comments (6)
spin0and2
and apin0and2_pure
cases now also work. The former has a speedup of around 6 with nspec=1
, and the latter around 3, but for the pure case there is still some remaining optimization potential.
# NOTE: this will not work with the official ducc0 release, but only with
# an experimental version from the "vector_wigners" branch.
# To install, first uninstall possibly existing ducc0 packages, then do
# pip install --no-binary ducc0 --user git+https://gitlab.mpcdf.mpg.de/mtr/ducc.git@vector_wigners
import numpy as np
from time import time
# We have to set nthreads before importing pspy, otherwise the
# environment variable change will have no effect.
nthreads=1
import os
os.environ["OMP_NUM_THREADS"]=str(nthreads)
# This must happen after setting OMP_NUM_THREADS!
from pspy.mcm_fortran.mcm_fortran import mcm_compute as mcm_fortran
import ducc0
# This routine is more complicated than mcm00_ducc, since a few multiplication
# steps are carried out in Python in pspy, and since the array indices are
# a bit different. Overall this should not have noticeable impact on performance
# at higher lmax.
def mcm00_pspy(spec, lmax):
nspec = spec.shape[0]
lrange_spec = np.arange(spec.shape[1])
res=np.zeros((nspec, lmax+1, lmax+1))
mcmtmp = np.empty((lmax+1, lmax+1))
for i in range(nspec):
wcl = spec[i]*(2*lrange_spec+1)
mcm_fortran.calc_coupling_spin0(wcl, lmax+1, lmax+1, lmax+1, mcmtmp.T)
mcm_fortran.fill_upper(mcmtmp.T)
mcmtmp *= (np.arange(2, lmax+3)*2+1.)/(4*np.pi)
res[i, 2:, 2:] = mcmtmp[:-2,:-2]
return res
def mcm02_pspy(spec, lmax):
nspec = spec.shape[0]
lrange_spec = np.arange(spec.shape[2])
res=np.zeros((nspec, 5, lmax+1, lmax+1))
mcmtmp = np.empty((5, lmax+1, lmax+1))
for i in range(nspec):
wcl = spec[i]*((2*lrange_spec+1).reshape((1,-1)))
mcm_fortran.calc_coupling_spin0and2(wcl[0], wcl[1], wcl[2], wcl[3], lmax+1, lmax+1, lmax+1, mcmtmp.T)
for j in range(5):
mcm_fortran.fill_upper(mcmtmp[j].T)
mcmtmp *= (np.arange(2, lmax+3)*2+1.)/(4*np.pi)
res[i, :, 2:, 2:] = mcmtmp[:,:-2,:-2]
return res
def mcm02_pure_pspy(spec, lmax):
nspec = spec.shape[0]
lrange_spec = np.arange(spec.shape[2])
res=np.zeros((nspec, 5, lmax+1, lmax+1))
mcmtmp = np.empty((5, lmax+1, lmax+1))
for i in range(nspec):
wcl = spec[i]*((2*lrange_spec+1).reshape((1,-1)))
mcm_fortran.calc_mcm_spin0and2_pure(wcl[0], wcl[1], wcl[2], wcl[3], mcmtmp.T)
mcmtmp *= (np.arange(2, lmax+3)*2+1.)/(4*np.pi)
res[i, :, 2:, 2:] = mcmtmp[:,:-2,:-2]
return res
def mcm00_ducc(spec, lmax):
return ducc0.misc.coupling_matrix_00(spec, lmax, nthreads=nthreads, algo=1)
def mcm02_ducc(spec, lmax):
return ducc0.misc.coupling_matrix_spin0and2(spec, lmax, nthreads=nthreads, algo=1)
def mcm02_pure_ducc(spec, lmax):
return ducc0.misc.coupling_matrix_spin0and2_pure(spec, lmax, nthreads=nthreads, algo=1)
# lmax up to which the MCM will be computed
lmax=1000
# number of spectra to process simultaneously
nspec=1
print()
print("Mode coupling matrix computation comparison")
print(f"nspec={nspec}, lmax={lmax}, nthreads={nthreads}")
print()
print("Spin 0 case:")
# we generate the spectra up to 2*lmax+1 to use all Wigner 3j symbols
# but this could also be lower.
spec = np.random.normal(size=(nspec, 2*lmax+1))
t0=time()
mcm_pspy = mcm00_pspy(spec, lmax)
print(f"pspy time: {time()-t0}s")
t0=time()
mcm_ducc = mcm00_ducc(spec, lmax)
print(f"ducc time: {time()-t0}s")
# pspy does not provide the first two rows and columns, so we zero them
# in the ducc result for comparison
mcm_ducc[:,0:2,:]=0
mcm_ducc[:,:,0:2]=0
# compare the results
print(f"L2 error between solutions: {ducc0.misc.l2error(mcm_pspy,mcm_ducc)}")
print()
print("Spin 0and2 case:")
# we generate the spectra up to 2*lmax+1 to use all Wigner 3j symbols
# but this could also be lower.
spec = np.random.normal(size=(nspec, 4, 2*lmax+1))
t0=time()
mcm_pspy = mcm02_pspy(spec, lmax)
print(f"pspy time: {time()-t0}s")
t0=time()
mcm_ducc = mcm02_ducc(spec, lmax)
print(f"ducc time: {time()-t0}s")
# pspy does not provide the first two rows and columns, so we zero them
# in the ducc result for comparison
mcm_ducc[:,:,0:2,:]=0
mcm_ducc[:,:,:,0:2]=0
# compare the results
print(f"L2 error between solutions: {ducc0.misc.l2error(mcm_pspy,mcm_ducc)}")
print()
print("Spin 0and2_pure case:")
# we generate the spectra up to 2*lmax+1 to use all Wigner 3j symbols
# but this could also be lower.
spec = np.random.normal(size=(nspec, 4, 2*lmax+1))
t0=time()
mcm_pspy = mcm02_pure_pspy(spec, lmax)
print(f"pspy time: {time()-t0}s")
t0=time()
mcm_ducc = mcm02_pure_ducc(spec, lmax)
print(f"ducc time: {time()-t0}s")
# pspy does not provide the first two rows and columns, so we zero them
# in the ducc result for comparison
mcm_ducc[:,:,0:2,:]=0
mcm_ducc[:,:,:,0:2]=0
# compare the results
print(f"L2 error between solutions: {ducc0.misc.l2error(mcm_pspy,mcm_ducc)}")
from pspy.
this is very impressive !
from pspy.
The spin0and2_pure
version is now vectorized and also exhibits a speedup factor of roughly 6.
During testing I noticed that the 5th array computed by calc_mcm_spin0and2_pure
always seems to be practically zero, no matter which power spectra I'm providing. My own code behaves in the same way. But if this array is in fact identically zero, why do we spend time computing it?
from pspy.
So Martin, this will be available in a future version of ducc ? we are about to finish the ACT dr6 analysis, so the code is kind of frozen at the moment, but really happy to make the change after that !
from pspy.
I'm absolutely planning to include this in a future ducc
release!
The reason for this early "unofficial" demo is that I'd like to get as many things right in the official version as possible, and that I want to include feedback from potential users.
A few questions that are still open for me:
- the
spin0
andspin_0and2
cases return matrices that are "almost" symmetric (up to a factor2*l1 + 1
), and this factor can be applied easily at a later point. Should I make use of this symmetry and only return a triangular matrix, saving a factor 2 in memory? pspy
does not compute values forel1<2
orel2<2
; should I compute those nevertheless, in case that other codes might want them?- what is the matter with the 5th component in
spin_0and2_pure
? Why is it always zero?
I realize now that I chose a fairly unconventional (and intrusive) way of asking for this kind of feedback, by opening an issue in a perfectly well-working code, and I'm sorry about that! If you prefer, I can certainly close the issue here and we can continue discussion on https://github.com/mreineck/ducc.
from pspy.
The code has now moved to the ducc0
branch (although it is still experimental), and a demo has been added to the ducc0
repository at python/demos/mcm_demo.py
. I'll continue to gather feedback before freezing the interface at some point.
from pspy.
Related Issues (8)
- possible python 2 to 3 bug leads to wrong type for taskrange HOT 1
- Flushes print command HOT 1
- Make sure to keep wcs attribute of enmap HOT 1
- pspy requiring ducc but not installing it as a dependency HOT 1
- Invalid memory accesses in calc_coupling_elem_spin0and2_pure? HOT 2
- Idea for speedup in calc_mcm_spin0and2_pure HOT 1
- Last component of `calc_mcm_spin0and2_pure` always seems to be (practically) zero
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from pspy.