Histograms of convergence mapsΒΆ

import sys

####################################################
#######LensTools functionality######################
####################################################

from lenstools import ConvergenceMap,Ensemble,GaussianNoiseGenerator
from lenstools.legacy.index import PDF,Indexer
from lenstools.simulations.igs1 import IGS1

#####################################################

import numpy as np
import astropy.units as u
import matplotlib.pyplot as plt

from lenstools.utils.mpi import MPIWhirlPool

import logging
import argparse

#########################################################################################
#########This function gets called on every map image and computes the histograms########
#########################################################################################

def compute_histograms(map_id,simulation_set,smoothing_scales,index,generator,bin_edges):

        assert len(index.descriptor_list) == len(smoothing_scales)

        z = 1.0

        #Get map name to analyze
        map_name = simulation_set.getNames(z=z,realizations=[map_id])[0]

        #Load the convergence map
        convergence_map = ConvergenceMap.load(map_name)

        #Generate the shape noise map
        noise_map = generator.getShapeNoise(z=z,ngal=15.0*u.arcmin**-2,seed=map_id)

        #Add the noise
        convergence_map += noise_map

        #Measure the features
        hist_output = np.zeros(index.size)
        for n,descriptor in enumerate(index.descriptor_list):

                logging.debug("Processing {0} x {1}".format(map_name,smoothing_scales[n]))

                smoothed_map = convergence_map.smooth(smoothing_scales[n])
                v,hist_output[descriptor.first:descriptor.last] = smoothed_map.pdf(bin_edges)

        #Return the histograms in array format
        return hist_output

############################################################
########################Main loop###########################
############################################################

if __name__=="__main__":

        #Initialize MPI pool
        try:
                pool = MPIWhirlPool()
        except ValueError:
                pool = None

        #Parse command line arguments
        parser = argparse.ArgumentParser()
        parser.add_argument("-v","--verbose",action="store_true",default=False,dest="verbose",help="Display debug info")
        parser.add_argument("-p","--path",action="store",dest="path",default="/Users/andreapetri/Documents/Columbia/spurious_shear/convergence_maps",help="Root path of IGS1 simulations")
        parser.add_argument("-n","--num_realizations",dest="num_realizations",action="store",type=int,default=3,help="How many realizations to process")

        cmd_args = parser.parse_args()

        #Set logging level
        if cmd_args.verbose:
                logging.basicConfig(level=logging.DEBUG)
        else:
                logging.basicConfig(level=logging.INFO)

        if (pool is not None) and not(pool.is_master()):

                pool.wait()
                sys.exit(0)

        #Root path of IGS1 maps
        root_path = cmd_args.path
        num_realizations = cmd_args.num_realizations

        #Smoothing scales in arcmin
        smoothing_scales = [ theta*u.arcmin for theta in [0.1,0.5,1.0,2.0] ]
        bin_edges = np.ogrid[-0.15:0.15:128j]
        bin_midpoints = 0.5*(bin_edges[1:] + bin_edges[:-1])

        #Create smoothing scale index for the histogram
        idx = Indexer.stack([PDF(bin_edges) for scale in smoothing_scales])

        #Create IGS1 simulation set object to look for the right simulations
        simulation_set = IGS1(root_path=root_path)

        #Look at a sample map
        sample_map = ConvergenceMap.load(simulation_set.getNames(z=1.0,realizations=[1])[0])

        #Initialize Gaussian shape noise generator
        generator = GaussianNoiseGenerator.forMap(sample_map)

        #Measure the histograms and load the data in the ensemble
        map_ensemble = Ensemble.compute(range(1,num_realizations+1),callback_loader=compute_histograms,pool=pool,simulation_set=simulation_set,smoothing_scales=smoothing_scales,index=idx,generator=generator,bin_edges=bin_edges)

        if pool is not None:
                pool.close()

        ##########################################################################################################################################
        ###############################Ensemble data available at this point for covariance, PCA, etc...##########################################
        ##########################################################################################################################################

        #Plot results to check
        fig,ax = plt.subplots(len(smoothing_scales),1)
        for i in range(len(smoothing_scales)):

                mean = map_ensemble.mean(0).values[idx[i].first:idx[i].last]
                error = np.sqrt(map_ensemble.cov().values.diagonal()[idx[i].first:idx[i].last])

                ax[i].errorbar(bin_midpoints,mean,yerr=error)

                ax[i].set_xlabel(r"$\kappa$")
                ax[i].set_ylabel(r"$P(\kappa)$")
                ax[i].set_title(r"${0:.1f}^\prime={1:.1f}$pix".format(smoothing_scales[i].value,(smoothing_scales[i] * sample_map.data.shape[0]/(sample_map.side_angle)).decompose().value))


        fig.tight_layout()
        fig.savefig("histograms.png")

You run this typing:

python histograms.py -p <path_to_your_simulations> -n <number_of_realizations>

Or, if you have a MPI installation and want to run on multiple processors:

mpiexec -n <number_of_processors> python histograms.py -p <path_to_your_simulations> -n <number_of_realizations>

This is how the result looks like

../_images/histograms.png