Measure the power spectrum of a convergence map (MPI)ΒΆ
When you need to measure the power spectrum of many different convergence maps (for example 1000 different realizations of a field of view), you are effectively measuring an Ensemble of quantities. lenstools provides Ensemble computing tools through the Ensemble
class. This class supports parallel operations too through MPI, like in this example. First define a function that measures the power spectrum out of a single map:
#The operations on convergence maps are handled with the ConvergenceMap class
from lenstools import ConvergenceMap
def measure_power_spectrum(filename,l_edges):
conv_map = ConvergenceMap.load(filename)
l,Pl = conv_map.powerSpectrum(l_edges)
return Pl
This is the actual code:
from lenstools.statistics.ensemble import Ensemble
from lenstools.utils.decorators import Parallelize
import logging
import numpy as np
import matplotlib.pyplot as plt
logging.basicConfig(level=logging.DEBUG)
@Parallelize.masterworker
def main(pool):
l_edges = np.arange(200.0,50000.0,200.0)
l = 0.5*(l_edges[:-1] + l_edges[1:])
conv_ensemble = Ensemble.compute(["Data/conv1.fit","Data/conv2.fit","Data/conv3.fit","Data/conv4.fit"],callback_loader=measure_power_spectrum,pool=pool,l_edges=l_edges)
fig,ax = plt.subplots()
for n in range(len(conv_ensemble)):
ax.plot(l,l*(l+1)*conv_ensemble.iloc[n]/(2.0*np.pi),label="Map {0}".format(n+1),linestyle="--")
mean = conv_ensemble.mean(0)
errors = np.sqrt(conv_ensemble.covariance().values.diagonal())
ax.errorbar(l,l*(l+1)*mean/(2.0*np.pi),yerr=l*(l+1)*errors/(2.0*np.pi),label="Mean")
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel(r"$l$")
ax.set_ylabel(r"$l(l+1)P_l/2\pi$")
ax.legend(loc="upper left")
plt.savefig("power_ensemble.png")
if __name__=="__main__":
main(None)
And this is the result
