In [8]:
############
#### Pangloss
############
from pangloss import BackgroundCatalog, ForegroundCatalog, \
TrueHaloMassDistribution, Kappamap, Shearmap
from pandas import DataFrame
# initialize background and foreground
import numpy as np
np.random.seed(0)
B = BackgroundCatalog(N=10.0, domain=[1.41, 1.4, -1.41, -1.4], field=[0, 0, 0, 0])
F = ForegroundCatalog.guo()
# initialize maps from Hilbert et al 2009
K = Kappamap.example()
S = Shearmap.example()
# drill lightcones and compute shears
B.drill_lightcones(radius=2.0, foreground=F, smooth_corr=False)
B.lens_by_map(K, S)
B.lens_by_halos(lookup_table=False, smooth_corr=False, relevance_lim=0)
print B.halo_mass_log_likelihood()
In [9]:
############
#### MassInference
############
from massinference.angle import Angle
from massinference.catalog import SourceCatalog, FastSampleHaloCatalogFactory, \
MutableHaloMassCatalog, HaloCatalog
from massinference.distribution import MassPrior
from massinference.inference import log_likelihood
from massinference.lenser import MapLenser
from massinference.lightcone import LightconeManager
from massinference.map import KappaMap, ShearMap
from massinference.plot import Limits
radius=2.0
sigma_e = 0.2
limits = Limits(Angle.from_radian(0.0246091424531), Angle.from_radian(0.0244346095279),
Angle.from_radian(-0.0246091424531), Angle.from_radian(-0.0244346095279))
limits_with_perimeter = limits.add_perimeter(Angle.from_arcmin(radius))
# get same source catalog that we used above
keymap = {
'ID': SourceCatalog.ID,
'RA': SourceCatalog.RA,
'Dec': SourceCatalog.DEC,
'z_obs': SourceCatalog.Z,
'e1_int': SourceCatalog.E1,
'e2_int': SourceCatalog.E2,
'eMod': SourceCatalog.E_MAG,
'ePhi': SourceCatalog.E_PHI,
}
source_catalog = SourceCatalog.from_file('/Users/user/Desktop/background.csv', keymap)
max_z = source_catalog.dataframe[SourceCatalog.Z].max()
base = MutableHaloMassCatalog.default(limits_with_perimeter, max_z)
halo_catalog_factory = FastSampleHaloCatalogFactory(base, MassPrior())
# lens by KappaMap and ShearMap
e1, e2 = MapLenser(KappaMap.default(), ShearMap.default()).lens(source_catalog)
# lens by halos and lightcones
lightcone_manager = LightconeManager(source_catalog, halo_catalog_factory, radius)
predictions = lightcone_manager.run(1)
print log_likelihood(predictions[0], e1, e2, sigma_e)