#CCL cosmology
cosmo_ccl = ccl.Cosmology(Omega_c = 0.30711 - 0.048254, Omega_b = 0.048254, h = 0.677, sigma8 = 0.8822714165197718, n_s=0.96, Omega_k = 0, transfer_function='eisenstein_hu')
#ccl_cosmo_set_high_prec (cosmo_ccl)

cosmo_numcosmo, dist, ps_lin, ps_nln, hmfunc = create_nc_obj (cosmo_ccl)

psf = hmfunc.peek_psf ()

"""
Explanation: initialize the Cosmological models
End of explanation
"""


#CosmoSim_proxy model
#M_0, z_0
theta_pivot = [3e14/0.71, 0.6]
#\mu_0, a_\mu^z, a_\mu^M
theta_mu = [3.19, -0.7, 2]
#\sigma_0, a_\sigma^z, a_\sigma^M
theta_sigma = [0.33, 0.,-0.08]

area = (0.25)*4*np.pi / 100.0
lnRl = 1.0
lnRu = 2.0
zl = 0.25
zu = 1.0

#Numcosmo_proxy model
cluster_z = nc.ClusterRedshift.new_from_name("NcClusterRedshiftNodist{'z-min': <%20.15e>, 'z-max':<%20.15e>}" % (zl, zu))
cluster_m = nc.ClusterMass.new_from_name("NcClusterMassAscaso{'M0':<%20.15e>,'z0':<%20.15e>,'lnRichness-min':<%20.15e>, 'lnRichness-max':<%20.15e>}" % (3e14/(0.71),0.6, lnRl, lnRu))
cluster_m.param_set_by_name('mup0', 3.19)
cluster_m.param_set_by_name('mup1', 2/np.log(10))
cluster_m.param_set_by_name('mup2', -0.7/np.log(10))
cluster_m.param_set_by_name('sigmap0', 0.33)
cluster_m.param_set_by_name('sigmap1', -0.08/np.log(10))
cluster_m.param_set_by_name('sigmap2', 0/np.log(10))

"""
Explanation: Define proxy modelling
Use a mass proxy, define the probability for observing a proxy given a mass and redhsift
$$
P(\log\lambda|M,z) = N(\mu(M,z), \sigma^2(M,z))
$$
the mean is
$$
\mu(M,z) = \mu_0 + a_\mu^M\log_{10}\frac{M}{M_0} + a_\mu^z\log_{10}\frac{1+z}{1+z_0} 
$$
variance is
$$
\sigma(M,z) = \sigma_0 + a_\sigma^M\log_{10}\frac{M}{M_0} + a_\sigma ^z\log_{10}\frac{1+z}{1+z_0} 
$$
End of explanation
"""


#Numcosmo Cluster Abundance

#First we need to define the multiplicity function here we will use the tinker
mulf = nc.MultiplicityFuncTinker.new()
mulf.set_linear_interp (True)
mulf.set_mdef(nc.MultiplicityFuncMassDef.CRITICAL)
mulf.set_Delta(200)
#Second we need to construct a filtered power spectrum 

hmf = nc.HaloMassFunction.new(dist,psf,mulf)
hmf.set_area(area)


ca = nc.ClusterAbundance.new(hmf,None)
mset = ncm.MSet.new_array([cosmo_numcosmo,cluster_m,cluster_z])

"""
Explanation: initialize the ClusterAbundance object
End of explanation
"""


ncount = Nc.DataClusterNCount.new (ca, "NcClusterRedshiftNodist", "NcClusterMassAscaso")
ca.prepare(cosmo_numcosmo,cluster_z,cluster_m)
mset = Ncm.MSet.new_array ([cosmo_numcosmo, cluster_z, cluster_m])

rng = Ncm.RNG.seeded_new (None, 1234)

ncount.init_from_sampling (mset, area, rng)
print (ncount.desc)

"""
Explanation: Creating a DataClusterNCount and initializing it using resample
Here we are using the seed 1234.
End of explanation
"""


#Bins definition
z_edges = np.linspace(zl, zu, 5)
proxy_edges = np.linspace(lnRl, lnRu, 4)

ncount.del_bins ()
for rich_l, rich_u in zip(proxy_edges, proxy_edges[1:]):
    for zl, zu in zip(z_edges,z_edges[1:]):
        lnM_obs_lb = Ncm.Vector.new_array ([rich_l])
        lnM_obs_ub = Ncm.Vector.new_array ([rich_u])
        z_obs_lb = Ncm.Vector.new_array ([zl])
        z_obs_ub = Ncm.Vector.new_array ([zu])
        
        ncount.add_bin (lnM_obs_lb, lnM_obs_ub, z_obs_lb, z_obs_ub)

ncount.bin_data ()

ncount.catalog_save ("ncount_ascaso.fits", True)

ncount.catalog_load ("ncount_ascaso.fits")

"""
Explanation: Binning data
End of explanation
"""


