Numerical Cosmology - NumCosmo

The NumCosmo is a free software C library whose main purposes are to test cosmological models using observational data and to provide a set of tools to perform cosmological calculations. Particularly, the current version has implemented three different probes: cosmic microwave background (CMB),supernovae type Ia (SNeIa) and large scale structure (LSS) information, such as baryonic acoustic oscillations (BAO) and galaxy cluster abundance. The code supports a joint analysis of these data and the parameter space can include cosmological and phenomenological parameters. It is worth emphasizing that NumCosmo matter power spectrum and CMB codes were written independently of other implementations such as CMBFAST, CAMB, etc.

The library is structured in such way to simplify the inclusion of non-standard cosmological models. Besides the functions related to cosmological quantities, this library also implements mathematical and statistical tools. The former was developed to enable the inclusion of other probes and/or theoretical models and to optimize the codes. The statistical framework comprises algorithms which define likelihood functions, minimization, Monte Carlo, Fisher Matrix and profile likelihood methods.

Download and installation

Version 0.13.3 of NumCosmo is the latest version available from our web site:

NumCosmo uses autotools for configuration and installation. The standard procedure is

./configure && make && sudo make install

More about compilation and installation here.

There are also some pre-compiled packages in openSUSE Build Service.

Mailing Lists

Subscribe to the numcosmo-help mailing list for discussion of questions and problems about using NumCosmo.

Examples

A simple example in C. Download the source.

To compile the example use:

gcc -Wall example_simple.c -o example_simple `pkg-config numcosmo --libs --cflags`
  1 #include <glib.h>
  2 #include <numcosmo/numcosmo.h>
  3 
  4 gint
  5 main (gint argc, gchar *argv[])
  6 {
  7   NcHICosmo *cosmo;
  8   NcDistance *dist;
  9   gint i;
 10 
 11   /**************************************************************************** 
 12    * Initializing the library objects, this must be called before 
 13    * any other library function.
 14    ****************************************************************************/  
 15   ncm_cfg_init ();
 16   
 17   /**************************************************************************** 
 18    * New homogeneous and isotropic cosmological model NcHICosmoDEXcdm.
 19    ****************************************************************************/  
 20   cosmo = nc_hicosmo_new_from_name (NC_TYPE_HICOSMO, "NcHICosmoDEXcdm");
 21 
 22   /**************************************************************************** 
 23    * New cosmological distance objects optimizied to perform calculations
 24    * up to redshift 2.0.
 25    ****************************************************************************/  
 26   dist = nc_distance_new (2.0);
 27  
 28   /**************************************************************************** 
 29    * Setting values for the cosmological model, those not set stay in the
 30    * default values. Remeber to use the _orig_ version to set the original
 31    * parameters in case when a reparametrization is used.
 32    ****************************************************************************/ 
 33   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_H0,       70.00);
 34   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_OMEGA_C,   0.25);
 35   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_OMEGA_X,   0.70);
 36   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_T_GAMMA0,  2.72);
 37   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_OMEGA_B,   0.05);
 38   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_XCDM_W,   -1.10);
 39 
 40   /**************************************************************************** 
 41    * Printing the parameters used.
 42    ****************************************************************************/
 43   printf ("# Model parameters:\n"); 
 44   ncm_model_params_log_all (NCM_MODEL (cosmo));
 45 
 46   /**************************************************************************** 
 47    * Printing some distances up to redshift 1.0.
 48    ****************************************************************************/ 
 49   for (i = 0; i < 10; i++)
 50   {
 51     gdouble z = 1.0 / 9.0 * i;
 52     gdouble cd = nc_hicosmo_RH_Mpc (cosmo) * nc_distance_comoving (dist, cosmo, z);
 53     printf ("% 10.8f % 20.15g\n", z, cd);
 54   }
 55 
 56   /**************************************************************************** 
 57    * Freeing objects.
 58    ****************************************************************************/ 
 59   nc_distance_free (dist);
 60   ncm_model_free (NCM_MODEL (cosmo));
 61 
 62   return 0;
 63 }

This example has the following output:

# Model parameters:
                    70                  0.25                   0.7                  2.72                  0.24                 3.046                  0.05                  -1.1
 0.00000000                    0
 0.11111111     466.128088514343
 0.22222222     910.794363376513
 0.33333333     1331.96821820897
 0.44444444     1728.76823442641
 0.55555556     2101.22844233792
 0.66666667        2450.04352253
 0.77777778     2776.33847725688
 0.88888889     3081.48376339937
 1.00000000     3366.95930114654

The same example in Python. Download the source.

To try this example you must have PyGObject installed, and numcosmo built with --enable-introspection option.

  1 #!/usr/bin/python2
  2 
  3 import math
  4 import gi
  5 gi.require_version('NumCosmo', '1.0')
  6 gi.require_version('NumCosmoMath', '1.0')
  7 
  8 from gi.repository import NumCosmo as Nc
  9 from gi.repository import NumCosmoMath as Ncm
 10 
 11 #
 12 #  Initializing the library objects, this must be called before 
 13 #  any other library function.
 14 #
 15 Ncm.cfg_init ()
 16 
 17 #
 18 #  New homogeneous and isotropic cosmological model NcHICosmoDEXcdm 
 19 #
 20 cosmo = Nc.HICosmo.new_from_name (Nc.HICosmo, "NcHICosmoDEXcdm")
 21 cosmo.set_reparam (Nc.HICosmoDEReparamCMB.new (cosmo.len ()))
 22 
 23 #
 24 #  New cosmological distance objects optimizied to perform calculations
 25 #  up to redshift 2.0.
 26 #
 27 dist = Nc.Distance.new (2.0)
 28 
 29 #
 30 #  Setting values for the cosmological model, those not set stay in the
 31 #  default values. Remember to use the _orig_ version to set the original
 32 #  parameters when a reparametrization is used.
 33 #
 34 
 35 #
 36 # C-like
 37 #
 38 cosmo.orig_param_set (Nc.HICosmoDEParams.H0,       70.00)
 39 cosmo.orig_param_set (Nc.HICosmoDEParams.OMEGA_C,   0.25)
 40 cosmo.orig_param_set (Nc.HICosmoDEParams.OMEGA_X,   0.70)
 41 cosmo.orig_param_set (Nc.HICosmoDEParams.T_GAMMA0,  2.72)
 42 cosmo.orig_param_set (Nc.HICosmoDEParams.OMEGA_B,   0.05)
 43 cosmo.orig_param_set (Nc.HICosmoDEXCDMParams.W,    -1.10)
 44 
 45 #
 46 # OO-like
 47 #
 48 cosmo.props.H0      = 70.00
 49 cosmo.props.Omegab  =  0.04
 50 cosmo.props.Omegac  =  0.25
 51 cosmo.props.Omegax  =  0.70
 52 cosmo.props.Tgamma0 =  2.72
 53 cosmo.props.w       = -1.10
 54 
 55 #
 56 #  Printing the parameters used.
 57 #
 58 print "# Model parameters: ", 
 59 cosmo.params_log_all ()
 60 
 61 #
 62 #  Printing some distances up to redshift 1.0.
 63 #
 64 for i in range (0, 10):
 65   z  = 1.0 / 9.0 * i
 66   cd = cosmo.RH_Mpc () * dist.comoving (cosmo, z)
 67   print "% 10.8f % 20.15g" % (z, cd)
 68 
 69 dist = Nc.Distance.new (2000.0)
 70 
 71 dist.prepare (cosmo)
 72 
 73 print dist.theta100CMB (cosmo)
 74 print cosmo.zt (5.0)
 75 
 76 

An example of the recombination object in Python. Download the source.

To try this example you must have PyGObject installed, matplotlib matplotlib and numcosmo built with --enable-introspection option.

  1 #!/usr/bin/python2
  2 
  3 import gi
  4 gi.require_version('NumCosmo', '1.0')
  5 gi.require_version('NumCosmoMath', '1.0')
  6 
  7 from math import *
  8 from gi.repository import GObject
  9 import matplotlib.pyplot as plt
 10 from gi.repository import NumCosmo as Nc
 11 from gi.repository import NumCosmoMath as Ncm
 12 
 13 #
 14 #  Initializing the library objects, this must be called before 
 15 #  any other library function.
 16 #
 17 Ncm.cfg_init ()
 18 
 19 #
 20 #  New homogeneous and isotropic cosmological model NcHICosmoDEXcdm 
 21 #
 22 cosmo = Nc.HICosmo.new_from_name (Nc.HICosmo, "NcHICosmoDEXcdm")
 23 
 24 #
 25 #  New homogeneous and isotropic reionization object.
 26 #
 27 reion = Nc.HIReionCamb.new ()
 28 
 29 #
 30 # Adding submodels to the main cosmological model.
 31 #
 32 cosmo.add_submodel (reion)
 33 
 34 #
 35 #  New recombination object configured to calculate up to redshift 
 36 #  10^9 and precision 10^-7.
 37 #
 38 recomb = Nc.Recomb.new_from_name ("NcRecombSeager{'prec':<1.0e-7>, 'zi':<1e9>}")
 39 
 40 #
 41 #  Setting values for the cosmological model, those not set stay in the
 42 #  default values. Remeber to use the _orig_ version to set the original
 43 #  parameters in case when a reparametrization is used.
 44 #
 45 
 46 #
 47 # C-like
 48 #
 49 cosmo.orig_param_set (Nc.HICosmoDEParams.H0,        70.0)
 50 cosmo.orig_param_set (Nc.HICosmoDEParams.OMEGA_C,   0.25)
 51 cosmo.orig_param_set (Nc.HICosmoDEParams.OMEGA_X,   0.70)
 52 cosmo.orig_param_set (Nc.HICosmoDEParams.T_GAMMA0,  2.72)
 53 cosmo.orig_param_set (Nc.HICosmoDEParams.OMEGA_B,   0.05)
 54 cosmo.orig_param_set (Nc.HICosmoDEXCDMParams.W,    -1.00)
 55 
 56 #
 57 # OO-like
 58 #
 59 cosmo.props.H0      = 70.0
 60 cosmo.props.Omegab  = 0.05
 61 cosmo.props.Omegac  = 0.25
 62 cosmo.props.Omegax  = 0.70
 63 cosmo.props.Tgamma0 = 2.72
 64 cosmo.props.w       = -1.0
 65 
 66 #
 67 #  Preparing recomb with cosmo.
 68 #
 69 recomb.prepare (cosmo)
 70 
 71 #
 72 #  Calculating Xe, equilibrium Xe, v_tau and its derivatives.
 73 #
 74 x_a      = []
 75 x_dtau_a = []
 76 Xe_a     = []
 77 Xefi_a   = []
 78 XHI_a    = []
 79 XHII_a   = []
 80 XHeI_a   = []
 81 XHeII_a  = []
 82 XHeIII_a = []
 83 
 84 dtau_dlambda_a = []
 85 x_E_a = []
 86 x_Etaub_a = []
 87 v_tau_a = []
 88 dv_tau_dlambda_a = []
 89 d2v_tau_dlambda2_a = []
 90 
 91 for i in range (10000):
 92   alpha = -log (1.0e4) + (log (1.0e4) - log (1.0e2)) / 9999.0 * i
 93   Xe = recomb.Xe (cosmo, alpha)
 94   x = exp (-alpha)
 95   Xefi     = recomb.equilibrium_Xe (cosmo, x)
 96 
 97   XHI_i    = recomb.equilibrium_XHI (cosmo, x)
 98   XHII_i   = recomb.equilibrium_XHII (cosmo, x)
 99   XHeI_i   = recomb.equilibrium_XHeI (cosmo, x)
100   XHeII_i  = recomb.equilibrium_XHeII (cosmo, x)
101   XHeIII_i = recomb.equilibrium_XHeIII (cosmo, x)
102 
103   v_tau = recomb.v_tau (cosmo, alpha)
104   dv_tau_dlambda = recomb.dv_tau_dlambda (cosmo, alpha)
105   d2v_tau_dlambda2 = recomb.d2v_tau_dlambda2 (cosmo, alpha)
106   
107   x_a.append (x)
108   Xe_a.append (Xe)
109   Xefi_a.append (Xefi)
110 
111   XHI_a.append (XHI_i)
112   XHII_a.append (XHII_i)
113   XHeI_a.append (XHeI_i)
114   XHeII_a.append (XHeII_i)
115   XHeIII_a.append (XHeIII_i)
116 
117   v_tau_a.append (-v_tau)
118   dv_tau_dlambda_a.append (-dv_tau_dlambda / 10.0)
119   d2v_tau_dlambda2_a.append (-d2v_tau_dlambda2 / 200.0)
120 
121 for i in range (10000):
122   alpha = -log (1.0e10) + (log (1.0e10) - log (1.0e2)) / 9999.0 * i
123   x = exp (-alpha)
124   dtau_dlambda = abs(recomb.dtau_dlambda (cosmo, alpha))
125   x_E = x / sqrt (cosmo.E2 (x))
126   x_E_taub = x_E / dtau_dlambda
127   
128   x_dtau_a.append (x)
129   x_E_a.append (x_E)
130   x_Etaub_a.append (x_E_taub)
131   dtau_dlambda_a.append (dtau_dlambda)
132 
133 #
134 #  Ploting ionization history.
135 #
136 
137 plt.title ("Ionization History")
138 plt.xscale('log')
139 plt.plot (x_a, Xe_a, 'r', label="Recombination")
140 plt.plot (x_a, Xefi_a, 'b--', label="Equilibrium")
141 plt.xlabel('$x$')
142 plt.ylabel('$X_{e^-}$')
143 plt.legend(loc=2)
144 
145 plt.savefig ("recomb_Xe.png")
146 
147 plt.clf ()
148 
149 #
150 #  Ploting all components history.
151 #
152 
153 plt.title ("Fractions Equilibrium History")
154 plt.xscale('log')
155 plt.plot (x_a, XHI_a, 'r', label=r'$X_\mathrm{HI}$')
156 plt.plot (x_a, XHII_a, 'g', label=r'$X_\mathrm{HII}$')
157 plt.plot (x_a, XHeI_a, 'b', label=r'$X_\mathrm{HeI}$')
158 plt.plot (x_a, XHeII_a, 'y', label=r'$X_\mathrm{HeII}$')
159 plt.plot (x_a, XHeIII_a, 'm', label=r'$X_\mathrm{HeIII}$')
160 plt.xlabel(r'$x$')
161 plt.ylabel(r'$X_*$')
162 plt.legend(loc=6)
163 
164 plt.ylim([-0.1,1.1])
165 
166 plt.savefig ("recomb_eq_fractions.png")
167 
168 plt.clf ()
169 
170 #
171 #  Ploting visibility function and derivatives.
172 #
173 
174 (lambdam, lambdal, lambdau) = recomb.v_tau_lambda_features (cosmo, 2.0 * log (10.0))
175 
176 plt.title ("Visibility Function and Derivatives")
177 plt.xscale('log')
178 plt.plot (x_a, v_tau_a, 'r', label=r'$v_\tau$')
179 plt.plot (x_a, dv_tau_dlambda_a, 'b-', label=r'$\frac{1}{10}\frac{dv_\tau}{d\lambda}$')
180 plt.plot (x_a, d2v_tau_dlambda2_a, 'g--', label=r'$\frac{1}{200}\frac{d^2v_\tau}{d\lambda^2}$')
181 plt.legend()
182 plt.legend(loc=3)
183 
184 #
185 #  Annotating max and width.
186 #
187 
188 v_tau_max = -recomb.v_tau (cosmo, lambdam)
189 
190 plt.annotate (r'$v_\tau^{max}$, $z=%5.2f$' % (exp(-lambdam)-1), xy=(exp(-lambdam), v_tau_max),  xycoords='data',
191               xytext=(0.1, 0.95), textcoords='axes fraction',
192               arrowprops=dict(facecolor='black', shrink=0.05))
193 
194 v_tau_u = -recomb.v_tau (cosmo, lambdau)
195 
196 plt.annotate (r'$v_\tau=10^{-2}v_\tau^{max}$, $z=%5.2f$' % (exp(-lambdau)-1), xy=(exp(-lambdau), v_tau_u),  xycoords='data',
197               xytext=(0.02, 0.75), textcoords='axes fraction',
198               arrowprops=dict(facecolor='black', shrink=0.05))
199 
200 v_tau_l = -recomb.v_tau (cosmo, lambdal)
201 
202 plt.annotate (r'$v_\tau=10^{-2}v_\tau^{max}$, $z=%5.2f$' % (exp(-lambdal)-1), xy=(exp(-lambdal), v_tau_l),  xycoords='data',
203               xytext=(0.65, 0.25), textcoords='axes fraction',
204               arrowprops=dict(facecolor='black', shrink=0.05))
205 
206 #
207 #  Annotating value of zstar.
208 #
209 
210 lambdastar = recomb.tau_zstar (cosmo)
211 
212 plt.annotate (r'$z^\star=%5.2f$' % (exp(-lambdastar)-1), 
213               xy=(0.65, 0.95), xycoords='axes fraction')
214 
215 plt.savefig ("recomb_v_tau.png")
216 
217 plt.clf ()
218 
219 #
220 #  Ploting dtau_dlambda
221 #
222 
223 plt.title (r'$d\tau/d\lambda$')
224 plt.xscale('log')
225 plt.yscale('log')
226 plt.plot (x_dtau_a, dtau_dlambda_a, 'r', label=r'$\vert\mathrm{d}\tau/\mathrm{d}\lambda\vert$')
227 plt.plot (x_dtau_a, x_E_a, 'b', label=r'$x/E$')
228 plt.plot (x_dtau_a, (x_Etaub_a), 'g', label=r'$x/(E\bar{\tau})$')
229 
230 plt.legend()
231 plt.legend(loc=2)
232 
233 plt.savefig ("recomb_dtau_dlambda.png")
234 

This example produces:

Ionization Fraction History Visibility Function

An example of Supernova Ia model fitting in Python. Download the source.

To try this example you must have PyGObject installed and numcosmo built with --enable-introspection option.

  1 #!/usr/bin/python2
  2 
  3 import gi
  4 gi.require_version('NumCosmo', '1.0')
  5 gi.require_version('NumCosmoMath', '1.0')
  6 
  7 from math import *
  8 from gi.repository import GObject
  9 import matplotlib.pyplot as plt
 10 from gi.repository import NumCosmo as Nc
 11 from gi.repository import NumCosmoMath as Ncm
 12 
 13 #
 14 #  Initializing the library objects, this must be called before 
 15 #  any other library function.
 16 #
 17 Ncm.cfg_init ()
 18 
 19 #
 20 #  New homogeneous and isotropic cosmological model NcHICosmoDEXcdm 
 21 #
 22 cosmo = Nc.HICosmo.new_from_name (Nc.HICosmo, "NcHICosmoDEXcdm")
 23 
 24 #
 25 #  Setting values for the cosmological model, those not set stay in the
 26 #  default values. Remeber to use the _orig_ version to set the original
 27 #  parameters in case when a reparametrization is used.
 28 #
 29 
 30 #
 31 # OO-like
 32 #
 33 cosmo.props.H0      = 70.0
 34 cosmo.props.Omegab  = 0.05
 35 cosmo.props.Omegac  = 0.25
 36 cosmo.props.Omegax  = 0.70
 37 cosmo.props.Tgamma0 = 2.72
 38 cosmo.props.w       = -1.0
 39 
 40 #
 41 #  Creating a new Modelset and set cosmo as the HICosmo model to be used.
 42 #
 43 mset = Ncm.MSet ()
 44 mset.set (cosmo)
 45 
 46 #
 47 #  Setting parameters Omega_c and w to be fitted and change parameter
 48 #  Omega_x -> Omega_k.
 49 #
 50 ##cosmo.de_omega_x2omega_k ()
 51 cosmo.props.Omegac_fit = True
 52 cosmo.props.w_fit = True
 53 
 54 #
 55 #  Creating a new Distance object optimized to redshift 2.
 56 #
 57 dist = Nc.Distance (zf = 2.0)
 58 
 59 #
 60 #  Creating a new Data object from distance modulus catalogs.
 61 #
 62 snia = Nc.DataDistMu.new_from_id (dist, Nc.DataSNIAId.SIMPLE_UNION2_1)
 63 
 64 #
 65 #  Creating a new Dataset and add snia to it.
 66 #
 67 dset = Ncm.Dataset ()
 68 dset.append_data (snia)
 69 
 70 #
 71 #  Creating a Likelihood from the Dataset.
 72 #
 73 lh = Ncm.Likelihood (dataset = dset)
 74 
 75 #
 76 #  Creating a Fit object of type NLOPT using the fitting algorithm ln-neldermead to
 77 #  fit the Modelset mset using the Likelihood lh and using a numerical differentiation
 78 #  algorithm (NUMDIFF_FORWARD) to obtain the gradient (if needed).
 79 #
 80 fit = Ncm.Fit.new (Ncm.FitType.NLOPT, "ln-neldermead", lh, mset, Ncm.FitGradType.NUMDIFF_FORWARD)
 81 
 82 #
 83 #  Running the fitter printing messages.
 84 #
 85 fit.run (Ncm.FitRunMsgs.SIMPLE)
 86 
 87 #
 88 #  Printing fitting informations.
 89 #
 90 fit.log_info ()
 91 
 92 #
 93 #  Calculating the parameters covariance using numerical differentiation.
 94 #
 95 fit.numdiff_m2lnL_covar ()
 96 
 97 #
 98 #  Printing the covariance matrix.
 99 # 
100 fit.log_covar ()

This example produces:

#----------------------------------------------------------------------------------
# Model fitting. Interating using:
#  - solver:            NLOpt:ln-neldermead
#  - differentiation:   Numerical differentiantion (forward)
#.............
#  Minimum found with precision:  1.00000e-08 (|df|/f or |dx|)
#  Elapsed time: 00 days, 00:00:00.0360950
#  iteration            [000043]
#  function evaluations [000043]
#  gradient evaluations [000000]
#  degrees of freedom   [000577]
#  m2lnL     =     562.224308438599
#  Fit parameters:
#    [ 0.230882725715637  ] [-1.00940240383148   ] 
#----------------------------------------------------------------------------------
# Data used:
#   - Union2.1 sample
#----------------------------------------------------------------------------------
# Model[00]:
#   - XCDM - Constant EOS
#----------------------------------------------------------------------------------
# Model parameters
#   -      H0[00]:  70                  [FIXED]
#   -  Omegac[01]:  0.230882725715637   [FREE]
#   -  Omegak[02]:  0                   [FIXED]
#   - Tgamma0[03]:  2.72                [FIXED]
#   -  Omegab[04]:  0.05                [FIXED]
#   -      ns[05]:  1                   [FIXED]
#   -  sigma8[06]:  0.9                 [FIXED]
#   -       w[07]: -1.00940240383148    [FREE]
#----------------------------------------------------------------------------------
# Fitted parameters covariance matrix
#                                              -------------------------------
# Omegac[0001] =  0.2309      +/-  0.07425     |  1           | -0.9653      |
#      w[0007] = -1.009       +/-  0.2014      | -0.9653      |  1           |
#                                              -------------------------------

An example of Supernova Ia and BAO model fitting and confidence region calculation in Python. Download the source.

To try this example you must have PyGObject installed, matplotlib matplotlib and numcosmo built with --enable-introspection option.

  1 #!/usr/bin/python2
  2 
  3 import gi
  4 gi.require_version('NumCosmo', '1.0')
  5 gi.require_version('NumCosmoMath', '1.0')
  6 
  7 from math import *
  8 from gi.repository import GObject
  9 import matplotlib.pyplot as plt
 10 from gi.repository import NumCosmo as Nc
 11 from gi.repository import NumCosmoMath as Ncm
 12 
 13 #
 14 #  Initializing the library objects, this must be called before 
 15 #  any other library function.
 16 #
 17 Ncm.cfg_init ()
 18 
 19 #
 20 #  New homogeneous and isotropic cosmological model NcHICosmoDEXcdm 
 21 #
 22 cosmo = Nc.HICosmo.new_from_name (Nc.HICosmo, "NcHICosmoDEXcdm")
 23 
 24 #
 25 #  Setting values for the cosmological model, those not set stay in the
 26 #  default values. Remeber to use the _orig_ version to set the original
 27 #  parameters in case when a reparametrization is used.
 28 #
 29 
 30 #
 31 # OO-like
 32 #
 33 cosmo.props.H0      = 70.0
 34 cosmo.props.Omegab  = 0.05
 35 cosmo.props.Omegac  = 0.25
 36 cosmo.props.Omegax  = 0.70
 37 cosmo.props.Tgamma0 = 2.72
 38 cosmo.props.w       = -1.0
 39 
 40 #
 41 #  A new Modelset with cosmo as the HICosmo model to be used.
 42 #
 43 mset = Ncm.MSet ()
 44 mset.set (cosmo)
 45 
 46 #
 47 #  Setting parameters Omega_c, Omega_x and w to be fitted (and change 
 48 #  parameter Omega_x -> Omega_k).
 49 #
 50 #cosmo.de_omega_x2omega_k ()
 51 cosmo.props.Omegac_fit = True
 52 cosmo.props.Omegax_fit = True
 53 cosmo.props.w_fit = True
 54 
 55 #
 56 #  A new Distance object optimized to redshift 2.
 57 #
 58 dist = Nc.Distance (zf = 2.0)
 59 
 60 #
 61 #  A new Data object from distance modulus catalogs.
 62 #  A new Data object from BAO catalogs.
 63 #
 64 snia = Nc.DataDistMu.new_from_id (dist, Nc.DataSNIAId.SIMPLE_UNION2_1)
 65 bao = Nc.data_bao_create (dist, Nc.DataBaoId.A_EISENSTEIN2005)
 66 
 67 #
 68 #  A new Dataset with snia and bao set.
 69 #
 70 dset = Ncm.Dataset ()
 71 dset.append_data (snia)
 72 dset.append_data (bao)
 73 
 74 #
 75 #  Creating a Likelihood from the Dataset.
 76 #
 77 lh = Ncm.Likelihood (dataset = dset)
 78 
 79 #
 80 #  Creating a Fit object of type NLOPT using the fitting algorithm ln-neldermead to
 81 #  fit the Modelset mset using the Likelihood lh and using a numerical differentiation
 82 #  algorithm (NUMDIFF_FORWARD) to obtain the gradient (if needed).
 83 #
 84 fit = Ncm.Fit.new (Ncm.FitType.NLOPT, "ln-neldermead", lh, mset, Ncm.FitGradType.NUMDIFF_FORWARD)
 85 
 86 #
 87 #  Running the fitter printing messages.
 88 #
 89 fit.run (Ncm.FitRunMsgs.SIMPLE)
 90 
 91 #
 92 #  Printing fitting informations.
 93 #
 94 fit.log_info ()
 95 
 96 #
 97 #  Calculating the parameters covariance using numerical differentiation.
 98 #
 99 fit.numdiff_m2lnL_covar ()
100 
101 #
102 #  Printing the covariance matrix.
103 # 
104 fit.log_covar ()
105 
106 #
107 #  Creating a new Likelihood ratio test object.
108 #  First we create two PIndex indicating which parameter
109 #    we are going to study.
110 # 
111 p1 = Ncm.MSetPIndex.new (cosmo.id (), Nc.HICosmoDEParams.OMEGA_C)
112 p2 = Ncm.MSetPIndex.new (cosmo.id (), Nc.HICosmoDEXCDMParams.W)
113 
114 lhr2d = Ncm.LHRatio2d.new (fit, p1, p2, 1.0e-3)
115 
116 #
117 #  Calculating the confidence region using the Likelihood ratio test.
118 #  Also calculate using the Fisher matrix approach.
119 #
120 cr_rg = lhr2d.conf_region (0.6826, 300.0, Ncm.FitRunMsgs.SIMPLE)
121 fisher_rg = lhr2d.fisher_border (0.6826, 300.0, Ncm.FitRunMsgs.SIMPLE)
122 
123 cr_p1array = cr_rg.p1.dup_array ()
124 cr_p2array = cr_rg.p2.dup_array ()
125 
126 fisher_p1array = fisher_rg.p1.dup_array ()
127 fisher_p2array = fisher_rg.p2.dup_array ()
128 
129 #
130 #  Ploting the confidence regions obtained from both methods.
131 #
132 
133 plt.title ("Confidence regions (%f)" % (cr_rg.clevel * 100.0))
134 plt.plot (cr_p1array, cr_p2array, 'r', label="Likelihood Ratio")
135 plt.plot (fisher_p1array, fisher_p2array, 'b-', label="Fisher Matrix")
136 
137 plt.xlabel(r'$\Omega_c$')
138 plt.ylabel(r'$w$')
139 
140 plt.legend(loc=4)
141 
142 plt.savefig ("snia_bao_rg_omegac_w.png")

This example has the following output and figure:

Confidence Regions
#----------------------------------------------------------------------------------
# Model fitting. Interating using:
#  - solver:            NLOpt:ln-neldermead
#  - differentiation:   Numerical differentiantion (forward)
#......................
#  Minimum found with precision:  1.00000e-08 (|df|/f or |dx|)
#  Elapsed time: 00 days, 00:00:00.1035530
#  iteration            [000137]
#  function evaluations [000137]
#  gradient evaluations [000000]
#  degrees of freedom   [000577]
#  m2lnL     =     562.202808924736
#  Fit parameters:
#    [ 0.225227921129909  ] [ 0.584140145983326  ] [-1.18852980893615   ] 
#----------------------------------------------------------------------------------
# Data used:
#   - Union2.1 sample
#   - Eisenstein 2005, BAO Sample A
#----------------------------------------------------------------------------------
# Model[00]:
#   - XCDM - Constant EOS
#----------------------------------------------------------------------------------
# Model parameters
#   -      H0[00]:  70                  [FIXED]
#   -  Omegac[01]:  0.225227921129909   [FREE]
#   -  Omegax[02]:  0.584140145983326   [FREE]
#   - Tgamma0[03]:  2.72                [FIXED]
#   -  Omegab[04]:  0.05                [FIXED]
#   -      ns[05]:  1                   [FIXED]
#   -  sigma8[06]:  0.9                 [FIXED]
#   -       w[07]: -1.18852980893615    [FREE]
#----------------------------------------------------------------------------------
# Fitted parameters covariance matrix
#                                              ----------------------------------------------
# Omegac[0001] =  0.2252      +/-  0.01881     |  1           |  0.1347      |  0.1098      |
# Omegax[0002] =  0.5841      +/-  0.4915      |  0.1347      |  1           |  0.9963      |
#      w[0007] = -1.189       +/-  0.9275      |  0.1098      |  0.9963      |  1           |
#                                              ----------------------------------------------
#----------------------------------------------------------------------------------
# Likelihood ratio confidence region at 68.260%, bestfit [  0.22522792   -1.1885298]:
#
#  looking root in interval [           0    1.5149868]:
#............
#  root found at   0.51764498 with precision 1.00000000e-05.
#  border found at   0.51764498.
#.
#  looking root in interval [   7.0669608    8.6377572]:
#.......
#  root found at    7.9006802 with precision 1.00000000e-05.
#  looking root in interval [   7.1152821    8.6860784]:
#......
.
.
. Several messages like the above.
.
.
#  root found at    14.179015 with precision 1.00000000e-05.
#  looking root in interval [   13.393617    14.964413]:
#.....
#  root found at    14.180544 with precision 1.00000000e-05.
#  looking root in interval [   13.395146    14.965943]:
#.......
#  root found at    14.182127 with precision 1.00000000e-05.
#  Start found at [270], ending...

Cluster abundance example in C. Download the source.

To compile the example use:

gcc -Wall example_ca.c -o example_ca `pkg-config numcosmo --libs --cflags`
  1 #include <glib.h>
  2 #include <numcosmo/numcosmo.h>
  3 
  4 gint
  5 main (gint argc, gchar *argv[])
  6 {
  7   NcHIPrim  *prim;
  8   NcHIReion *reion;
  9   NcHICosmo *cosmo;
 10   NcDistance *dist;
 11   NcTransferFunc *tf;
 12   NcMultiplicityFunc *mulf;
 13   NcHaloMassFunction *mf;
 14   NcPowspecML *psml;
 15   NcmPowspecFilter *psf;
 16   guint np = 50;
 17   gint i;
 18 
 19   /**************************************************************************** 
 20    * Initializing the library objects, this must be called before 
 21    * any other library function.
 22    ****************************************************************************/  
 23   ncm_cfg_init ();
 24   
 25   /**************************************************************************** 
 26    * New homogeneous and isotropic cosmological model NcHICosmoDEXcdm.
 27    ****************************************************************************/  
 28   cosmo = nc_hicosmo_new_from_name (NC_TYPE_HICOSMO, "NcHICosmoDEXcdm");
 29 
 30   /**************************************************************************** 
 31    * New homogeneous and isotropic reionization object.
 32    ****************************************************************************/  
 33   reion = NC_HIREION (nc_hireion_camb_new ());
 34 
 35   /**************************************************************************** 
 36    * New homogeneous and isotropic primordial object.
 37    ****************************************************************************/  
 38   prim  = NC_HIPRIM (nc_hiprim_power_law_new ());
 39 
 40   /**************************************************************************** 
 41    * Adding the submodels to the main cosmology model.
 42    ****************************************************************************/  
 43   ncm_model_add_submodel (NCM_MODEL (cosmo), NCM_MODEL (reion));
 44   ncm_model_add_submodel (NCM_MODEL (cosmo), NCM_MODEL (prim));
 45 
 46   /**************************************************************************** 
 47    * New cosmological distance objects optimizied to perform calculations
 48    * up to redshift 2.0.
 49    ****************************************************************************/  
 50   dist = nc_distance_new (2.0);
 51  
 52   /**************************************************************************** 
 53    * New transfer function 'NcTransferFuncEH' using the Einsenstein and Hu
 54    * fitting formula.
 55    ****************************************************************************/  
 56   tf = nc_transfer_func_new_from_name ("NcTransferFuncEH");
 57 
 58   /**************************************************************************** 
 59    * New linear matter power spectrum object based of the EH transfer function.
 60    ****************************************************************************/  
 61   psml = NC_POWSPEC_ML (nc_powspec_ml_transfer_new (tf));
 62   ncm_powspec_require_kmin (NCM_POWSPEC (psml), 1.0e-3);
 63   ncm_powspec_require_kmax (NCM_POWSPEC (psml), 1.0e3);
 64 
 65   /**************************************************************************** 
 66    * Apply a tophat filter to the psml object, set best output interval.
 67    ****************************************************************************/     
 68   psf = ncm_powspec_filter_new (NCM_POWSPEC (psml), NCM_POWSPEC_FILTER_TYPE_TOPHAT);
 69   ncm_powspec_filter_set_best_lnr0 (psf);
 70 
 71   /**************************************************************************** 
 72    * New multiplicity function 'NcMultiplicityFuncTinkerMean'
 73    ****************************************************************************/  
 74   mulf = nc_multiplicity_func_new_from_name ("NcMultiplicityFuncTinkerMean");
 75 
 76   /**************************************************************************** 
 77    * New mass function object using the objects defined above.
 78    ****************************************************************************/  
 79   mf = nc_halo_mass_function_new (dist, psf, mulf);
 80 
 81   /**************************************************************************** 
 82    * Setting values for the cosmological model, those not set stay in the
 83    * default values. Remeber to use the _orig_ version to set the original
 84    * parameters in case when a reparametrization is used.
 85    ****************************************************************************/ 
 86   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_H0,        70.0);
 87   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_OMEGA_C,    0.25);
 88   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_OMEGA_X,    0.7);
 89   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_T_GAMMA0,   1.0);
 90   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_OMEGA_B,    0.05);
 91   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_XCDM_W,    -1.1);
 92 
 93   /**************************************************************************** 
 94    * Printing the parameters used.
 95    ****************************************************************************/
 96   printf ("# Model parameters:\n#"); 
 97   ncm_model_params_log_all (NCM_MODEL (cosmo));
 98 
 99   /**************************************************************************** 
100    * Printing the growth function and its derivative with respect to z 
101    * up to redshift 1.
102    ****************************************************************************/ 
103   ncm_powspec_prepare (NCM_POWSPEC (psml), NCM_MODEL (cosmo));
104 
105   {
106     NcGrowthFunc *gf = nc_powspec_ml_transfer_peek_gf (NC_POWSPEC_ML_TRANSFER (psml));
107     for (i = 0; i < np; i++)
108     {
109       gdouble z = 1.0 / (np - 1.0) * i;
110       gdouble gfz  = nc_growth_func_eval (gf, cosmo, z);
111       gdouble dgfz = nc_growth_func_eval_deriv (gf, cosmo, z);
112       printf ("% 10.8f % 21.15g % 21.15g\n", z, gfz, dgfz);
113     }
114     printf ("\n\n");
115   }
116 
117   /**************************************************************************** 
118    * Printing the matter power spectrum in the k (in unities of Mpc^-1) 
119    * interval [1e-3, 1e3] 
120    ****************************************************************************/
121 
122   for (i = 0; i < np; i++)
123   {
124     gdouble lnk = log (1e-3) + log (1e6) / (np - 1.0) * i;
125     gdouble Pmk = ncm_powspec_eval (NCM_POWSPEC (psml), NCM_MODEL (cosmo), 0.0, exp (lnk));
126     printf ("% 10.8f % 21.15g\n", exp (lnk), Pmk);
127   }
128   printf ("\n\n");
129 
130   /**************************************************************************** 
131    * Printing the variance filtered with the tophat windown function using
132    * scales R in the interval [5, 50] at redshift 0.3.
133    * First calculates the growth function at z = 0.3 and then the spectrum
134    * amplitude from the sigma8 parameter.
135    ****************************************************************************/
136   ncm_powspec_filter_prepare (psf, NCM_MODEL (cosmo));
137   {
138     for (i = 0; i < np; i++)
139     {
140       gdouble lnR = log (5.0) + log (10.0) / (np - 1.0) * i;
141       gdouble sigma2       = ncm_powspec_filter_eval_var_lnr (psf, 0.0, lnR);
142       gdouble dsigma2_dlnR = ncm_powspec_filter_eval_dvar_dlnr (psf, 0.0, lnR);
143       printf ("% 10.8f % 21.15g % 21.15g\n", exp (lnR), sigma2, dsigma2_dlnR);
144     }
145     printf ("\n\n");
146   }
147 
148   /**************************************************************************** 
149    * Printing the mass function integrated in the mass interval [1e14, 1e16]
150    * for the redhshifts in the interval [0, 2.0] and area 200 squared degree.
151    ****************************************************************************/
152   nc_halo_mass_function_set_area_sd (mf, 200.0);
153   nc_halo_mass_function_set_eval_limits (mf, cosmo, log (1e14), log (1e16), 0.0, 2.0);
154   nc_halo_mass_function_prepare (mf, cosmo);
155 
156   for (i = 0; i < np; i++)
157   {
158     gdouble z = 2.0 / (np - 1.0) * i;
159     gdouble dndz = nc_halo_mass_function_dn_dz (mf, cosmo, log (1.0e14), log (1.0e16), z, FALSE);
160     printf ("% 10.8f % 21.15g\n", z, dndz);
161   }
162   printf ("\n\n");
163 
164 
165   /**************************************************************************** 
166    * Freeing objects.
167    ****************************************************************************/ 
168   nc_distance_free (dist);
169   ncm_model_free (NCM_MODEL (cosmo));
170   nc_hireion_free (reion);
171   nc_transfer_func_free (tf);
172   nc_powspec_ml_free (psml);
173   ncm_powspec_filter_free (psf);
174   nc_multiplicity_func_free (mulf);
175   nc_halo_mass_function_free (mf);
176 
177   return 0;
178 }

This example has the following output:

# Model parameters:
#                    70                  0.25                   0.7                     1                  0.24                 3.046                  0.05                  -1.1
 0.00000000                     1    -0.514644654419114
 0.02040816     0.989521826725434    -0.512155244764271
 0.04081633     0.979098119606862    -0.509310781871472
 0.06122449     0.968735891665614     -0.50613407471797
 0.08163265     0.958441692132124    -0.502647748249237
 0.10204082     0.948221611280675    -0.498874137680349
 0.12244898      0.93808128729867     -0.49483519697729
 0.14285714     0.928025915007738    -0.490552408995221
 0.16326531     0.918060256240339    -0.486046711990544
 0.18367347     0.908188651681923    -0.481338429687734
 0.20408163     0.898415033977879    -0.476447217261769
 0.22448980     0.888742941916688    -0.471392011077301
 0.24489796     0.879175535498096    -0.466190989192638
 0.26530612     0.869715611706689    -0.460861545417226
 0.28571429     0.860365620820335     -0.45542026048715
 0.30612245       0.8511276830911    -0.449882893307777
 0.32653061     0.842003605649148    -0.444264370106883
 0.34693878     0.832994899495015    -0.438578786587639
 0.36734694     0.824102796456841    -0.432839407186525
 0.38775510     0.815328265994592    -0.427058681439474
 0.40816327       0.8066720317663     -0.42124825145089
 0.42857143     0.798134587860193    -0.415418975090372
 0.44897959      0.78971621462385     -0.40958094479373
 0.46938776      0.78141699403074    -0.403743513854987
 0.48979592     0.773236824534143    -0.397915319528739
 0.51020408     0.765175435359684    -0.392104311802892
 0.53061224        0.757232400225    -0.386317785507923
 0.55102041     0.749407150438857    -0.380562408230182
 0.57142857      0.74169898738897     -0.37484424964717
 0.59183673     0.734107094396005    -0.369168814847993
 0.61224490     0.726630547936096    -0.363541074461111
 0.63265306     0.719268328236883    -0.357965495849794
 0.65306122     0.712019329254614    -0.352446072213635
 0.67346939     0.704882368040839    -0.346986350728436
 0.69387755     0.697856193515343    -0.341589465472344
 0.71428571     0.690939494671464    -0.336258160472153
 0.73469388     0.684130908210908    -0.330994820634155
 0.75510204     0.677429025658355     -0.32580149243272
 0.77551020     0.670832399953898    -0.320679913900392
 0.79591837     0.664339551557803     -0.31563153477332
 0.81632653     0.657948974090976    -0.310657539214771
 0.83673469     0.651659139531841     -0.30575886723651
 0.85714286     0.645468502994453    -0.300936234717449
 0.87755102     0.639375507111642    -0.296190152122611
 0.89795918     0.633378586046076    -0.291520942016937
 0.91836735     0.627476169151325    -0.286928755459564
 0.93877551     0.621666684304661    -0.282413588008501
 0.95918367     0.615948560936719    -0.277975295225614
 0.97959184     0.610320232770929    -0.273613602896915
 1.00000000     0.604780140295806    -0.269328123759946


 0.00100000      20992.0105375543
 0.00132571      27611.5721478528
 0.00175751      36303.1990910879
 0.00232995      47696.4966245018
 0.00308884      62589.4327595911
 0.00409492      81966.0388032499
 0.00542868      106980.466407995
 0.00719686      138854.132747529
 0.00954095      178564.635905189
 0.01264855      226046.990158811
 0.01676833      278354.766962786
 0.02222996      326614.943405558
 0.02947052      356954.039913845
 0.03906940      370362.141215367
 0.05179475      401316.928256473
 0.06866488      460851.390548054
 0.09102982      485265.279663216
 0.12067926      472704.972426439
 0.15998587      442810.807939275
 0.21209509      385951.661667605
 0.28117687      313352.060997751
 0.37275937      237218.138630065
 0.49417134      168163.694698828
 0.65512856      112232.568462686
 0.86851137      71103.1506006262
 1.15139540      43131.9201085223
 1.52641797      25248.5095973815
 2.02358965      14351.9517177465
 2.68269580      7959.55929932173
 3.55648031      4321.88249317831
 4.71486636       2303.4007149377
 6.25055193       1207.3488223963
 8.28642773      623.404516903949
 10.98541142      317.542004464884
 14.56348478      159.767751713327
 19.30697729      79.4969439754696
 25.59547923      39.1614592604859
 33.93221772      19.1181395096212
 44.98432669      9.25762916023201
 59.63623317       4.4500984263147
 79.06043211       2.1250297017269
 104.81131342      1.00869633890308
 138.94954944     0.476209967563408
 184.20699693     0.223714012871928
 244.20530945      0.10462453118995
 323.74575428    0.0487288831369188
 429.19342601    0.0226099588574535
 568.98660290     0.010454531479576
 754.31200634   0.00481856229723604
 1000.00000000   0.00221433080313483


 5.00000000       352.22812124856     -703.519625517058
 5.24056567      320.453436307077     -649.479977995808
 5.49270571       291.14071333352     -598.710101386115
 5.75697700      264.139232089674     -551.085439702178
 6.03396320        239.3041250489     -506.481486640893
 6.32427608      216.496370885828     -464.773965222956
 6.62855683      195.582779060954      -425.83902696738
 6.94747747      176.435964620594     -389.553468888724
 7.28174239      158.934312424672     -355.794966123377
 7.63208984      142.961930128137     -324.442317558735
 7.99929360      128.408589374548     -295.375701450075
 8.38416468      115.169654809688      -268.47693769857
 8.78755312       103.14600068544     -243.629753217917
 9.21034985      92.2439149961179     -220.720046658131
 9.65348864      82.3749912664377     -199.636148697025
 10.11794824      73.4560082884816     -180.269074124235
 10.60475444       65.408798279596     -162.512762071774
 11.11498241      58.1601041001827     -146.264300946898
 11.64975905      51.6414263255456     -131.424134914323
 12.21026547      45.7888611058059     -117.896249140305
 12.79773961      40.5429298691392     -105.588331434523
 13.41347898      35.8484020236741     -94.4119083997072
 14.05884349        31.65411189039     -84.2824547052855
 14.73525851      27.9127711520905     -75.1194746254199
 15.44421798      24.5807781315379     -66.8465555074631
 16.18728771      21.6180252154133      -59.391393348745
 16.96610886      18.9877057208739     -52.6857911439787
 17.78240153      16.6561214597213     -46.6656311117018
 18.63796860      14.5924921936746       -41.27082229924
 19.53469969      12.7687680955594     -36.4452254047737
 20.47457531      11.1594462383292     -32.1365569241757
 21.45967130        9.741392029621     -28.2962749361501
 22.49216334      8.49366639739007     -24.8794489761155
 23.57433182      7.39735941510409     -21.8446165202251
 24.70856681      6.43543093606257     -19.1536286097571
 25.89737340      5.59255868845342     -16.7714870972825
 27.14337720      4.85499416834021     -14.6661758951481
 28.44933015      4.21042655917633     -12.8084884614094
 29.81811658      3.64785480567749     -11.1718535707401
 31.25275963      3.15746787873541     -9.73216119695467
 32.75642784      2.73053318813328     -8.46759007726422
 34.33244225      2.35929303262744      -7.3584382437539
 35.98428365      2.03686892365992     -6.38695750543312
 37.71560032      1.75717357927453     -5.53719260359912
 39.53021605      1.51483035309689     -4.79482573582747
 41.43213864      1.30509981817025     -4.14702786813802
 43.42556869      1.12381310758768     -3.58232090478772
 45.51490890     0.967311278784347     -3.09046212171767
 47.70477382     0.832388996160069     -2.66238609917912
 50.00000000     0.716237123697593     -2.29035718089057


 0.00000000                     0
 0.04081633      300.885914134605
 0.08163265      1166.22916177237
 0.12244898      2538.51906692107
 0.16326531       4359.3030876979
 0.20408163      6570.44656463484
 0.24489796      9115.21215512547
 0.28571429      11939.1998962209
 0.32653061      14991.0391813295
 0.36734694      18222.9174169036
 0.40816327      21590.9168667584
 0.44897959      25055.1967446792
 0.48979592      28580.0488877436
 0.53061224      32133.8402842676
 0.57142857      35688.8605789788
 0.61224490        39221.12371093
 0.65306122       42710.114597382
 0.69387755      46138.5078011454
 0.73469388      49491.8945964343
 0.77551020      52758.4764686533
 0.81632653      55928.7961852755
 0.85714286      58995.4539628103
 0.89795918      61952.8295067753
 0.93877551      64796.9019761646
 0.97959184      67525.0572590576
 1.02040816      70135.9059438574
 1.06122449      72629.1297306759
 1.10204082      75005.3428843361
 1.14285714      77265.9686247836
 1.18367347      79413.1387769998
 1.22448980      81449.5981945698
 1.26530612      83378.6295873068
 1.30612245      85203.9805479686
 1.34693878      86929.8101485711
 1.38775510      88560.6404107272
 1.42857143      90101.3153282961
 1.46938776       91556.968105198
 1.51020408      92932.9800735537
 1.55102041      94234.9684157802
 1.59183673      95468.7877746998
 1.63265306      96640.4768753469
 1.67346939      97756.2665631033
 1.71428571      98822.5523412273
 1.75510204      99845.8611765124
 1.79591837      100832.811842252
 1.83673469      101789.980199225
 1.87755102      102723.727712041
 1.91836735      103639.806279308
 1.95918367      104542.669319772
 2.00000000      105434.099964863


An example of the cluster abundance object in Python. Download the source.

To try this example you must have PyGObject installed, matplotlib matplotlib and numcosmo built with --enable-introspection option.

  1 #!/usr/bin/python2
  2 
  3 import gi
  4 gi.require_version('NumCosmo', '1.0')
  5 gi.require_version('NumCosmoMath', '1.0')
  6 
  7 import numpy
  8 from math import *
  9 from gi.repository import GObject
 10 import matplotlib.pyplot as plt
 11 from gi.repository import NumCosmo as Nc
 12 from gi.repository import NumCosmoMath as Ncm
 13 
 14 #
 15 #  Initializing the library objects, this must be called before 
 16 #  any other library function.
 17 #
 18 Ncm.cfg_init ()
 19 
 20 #
 21 #  New homogeneous and isotropic cosmological model NcHICosmoDEXcdm 
 22 #
 23 cosmo = Nc.HICosmo.new_from_name (Nc.HICosmo, "NcHICosmoDEXcdm")
 24 
 25 #
 26 #  New homogeneous and isotropic reionization object.
 27 #
 28 reion = Nc.HIReionCamb.new () 
 29 
 30 #
 31 #  New homogeneous and isotropic primordial object.
 32 #
 33 prim = Nc.HIPrimPowerLaw.new () 
 34 
 35 #
 36 # Adding submodels to the main cosmological model.
 37 #
 38 cosmo.add_submodel (reion)
 39 cosmo.add_submodel (prim)
 40 
 41 #
 42 #  New cosmological distance objects optimizied to perform calculations
 43 #  up to redshift 2.0.
 44 #
 45 dist = Nc.Distance.new (2.0)
 46 
 47 #
 48 # New transfer function 'NcTransferFuncEH' using the Einsenstein, Hu
 49 # fitting formula.
 50 #
 51 tf = Nc.TransferFunc.new_from_name ("NcTransferFuncEH")
 52 
 53 #
 54 # New linear matter power spectrum object based of the EH transfer function.
 55 # 
 56 psml = Nc.PowspecMLTransfer.new (tf)
 57 psml.require_kmin (1.0e-3)
 58 psml.require_kmax (1.0e3)
 59 
 60 #
 61 # Apply a tophat filter to the psml object, set best output interval.
 62 #
 63 psf = Ncm.PowspecFilter.new (psml, Ncm.PowspecFilterType.TOPHAT)
 64 psf.set_best_lnr0 ()
 65 
 66 #
 67 # New multiplicity function 'NcMultiplicityFuncTinkerMean'
 68 #
 69 mulf = Nc.MultiplicityFunc.new_from_name ("NcMultiplicityFuncTinkerMean")
 70 
 71 #
 72 # New mass function object using the objects defined above.
 73 #
 74 mf = Nc.HaloMassFunction.new (dist, psf, mulf)
 75 
 76 #
 77 #  Setting values for the cosmological model, those not set stay in the
 78 #  default values. Remember to use the _orig_ version to set the original
 79 #  parameters when a reparametrization is used.
 80 #
 81 cosmo.props.H0      = 70.0
 82 cosmo.props.Omegab  = 0.05
 83 cosmo.props.Omegac  = 0.25
 84 cosmo.props.Omegax  = 0.70
 85 cosmo.props.Tgamma0 = 2.72
 86 cosmo.props.w       = -1.0
 87 
 88 #
 89 #  Printing the parameters used.
 90 #
 91 print "# Model parameters: ", 
 92 cosmo.params_log_all ()
 93 
 94 #
 95 # Number of points to build the plots
 96 #
 97 np = 2000
 98 divfac = 1.0 / (np - 1.0)
 99 
100 #
101 #  Calculating growth and its derivative in the [0, 2] redshift
102 #  range.
103 #
104 za = []
105 Da = []
106 dDa = []
107 
108 gf = psml.peek_gf ()
109 gf.prepare (cosmo)
110 
111 for i in range (0, np):
112   z = 2.0 * divfac * i
113   D = gf.eval (cosmo, z)
114   dD = gf.eval_deriv (cosmo, z)
115   za.append (z)
116   Da.append (D)
117   dDa.append (dD)
118   print z, D
119 
120 #
121 #  Ploting growth function.
122 #
123 
124 plt.title ("Growth Function")
125 plt.plot (za, Da, 'r', label="D")
126 plt.plot (za, dDa, 'b--', label="dD/dz")
127 plt.xlabel('$z$')
128 plt.legend(loc=2)
129 #plt.yscale ('log')
130 
131 plt.savefig ("growth_func.png")
132 plt.clf ()
133 
134 #
135 # Calculating the transfer function and the matter power spectrum in the
136 # kh (in unities of h/Mpc) interval [1e-3, 1e3]
137 #
138 
139 ka = []
140 Ta = []
141 Pma = []
142 
143 psml.prepare (cosmo)
144 
145 for i in range (0, np):
146   lnk = log (1e-3) +  log (1e6) * divfac * i
147   k   = exp (lnk)
148   T   = tf.eval (cosmo, k)
149   Pm  = psml.eval (cosmo, 0.0, k)
150   ka.append (k)
151   Ta.append (T)
152   Pma.append (Pm)
153 
154 #
155 #  Ploting transfer and matter power spectrum
156 #
157 
158 plt.title ("Transfer Function and Matter Power Spectrum")
159 plt.yscale('log')
160 plt.xscale('log')
161 plt.plot (ka, Ta, 'r', label="T(kh)")
162 plt.plot (ka, Pma, 'b--', label="P_m(kh)")
163 plt.xlabel('$k$')
164 plt.legend(loc=1)
165 
166 plt.savefig ("transfer_func.png")
167 plt.clf ()
168 
169 #
170 # Calculating the variance filtered with the tophat windown function using
171 # scales R in the interval [5, 50] at redshift 0.3.
172 # First calculates the growth function at z = 0.3 and then the spectrum
173 # amplitude from the sigma8 parameter.
174 #
175 psf.prepare (cosmo)
176 
177 Ra = []
178 sigma2a = []
179 dlnsigma2a = []
180 
181 for i in range (0, np):
182   lnR       = log (5.0) +  log (10.0) * divfac * i
183   R         = exp (lnR)
184   sigma2    = psf.eval_var_lnr (0.0, lnR)
185   dlnsigma2 = psf.eval_dlnvar_dlnr (0.0, lnR) 
186   Ra.append (R)
187   sigma2a.append (sigma2)
188   dlnsigma2a.append (dlnsigma2)
189 
190 #
191 #  Ploting filtered matter variance
192 #
193 plt.title ("Variance and Variance Derivative")
194 plt.plot (Ra, sigma2a, 'r', label='$\sigma^2(\ln(R))$')
195 plt.plot (Ra, dlnsigma2a, 'b--', label='$d\ln(\sigma^2)/d\ln(R)$')
196 plt.xlabel('$R$')
197 plt.legend(loc=1)
198 
199 plt.savefig ("matter_var.png")
200 plt.clf ()
201 
202 #
203 # Calculating the mass function integrated in the mass interval [1e14, 1e16]
204 # for the redhshifts in the interval [0, 2.0] and area 200 squared degree.
205 #
206 
207 mf.set_area_sd (200.0)
208 mf.set_eval_limits (cosmo, log (1e14), log(1e16), 0.0, 2.0)
209 mf.prepare (cosmo)
210 
211 dndza = []
212 
213 for i in range (0, np):
214   dndz = mf.dn_dz (cosmo, log (1.0e14), log (1.0e16), za[i], True)
215   dndza.append (dndz)
216 
217 #
218 #  Ploting the mass function
219 #
220 
221 plt.title ("Mass Function")
222 plt.plot (za, dndza, 'r', label='$dn/dz$')
223 plt.xlabel('$z$')
224 plt.legend(loc=1)
225 
226 plt.savefig ("mass_function.png")
227 plt.clf ()

This example produces:

Growth Function Transfer Function and Matter Power Spectrum Variance and Variance Derivative Mass Function

An example of the NcHIPrim (primordial cosmology object) usage in Python. Download the source.

To try this example you must have PyGObject installed, matplotlib matplotlib and numcosmo built with --enable-introspection option.

This example uses Class as backend, to more information about Class see the links in the documentation: NcHIPertBoltzmannCBE

  1 #!/usr/bin/python2
  2 
  3 import gi
  4 import math
  5 import numpy as np
  6 import matplotlib.pyplot as plt
  7 
  8 gi.require_version('NumCosmo', '1.0')
  9 gi.require_version('NumCosmoMath', '1.0')
 10 
 11 from gi.repository import GObject
 12 from gi.repository import NumCosmo as Nc
 13 from gi.repository import NumCosmoMath as Ncm
 14 
 15 from py_hiprim_example import PyHIPrimExample
 16 
 17 #
 18 # Script params
 19 #
 20 lmax = 2500
 21 
 22 #
 23 # Creating a new instance of PyHIPrimExample
 24 #
 25 prim = PyHIPrimExample ()
 26 
 27 print "# As        = ", prim.props.As
 28 print "# P (k = 1) = ", prim.SA_powspec_k (1.0)
 29 print "# (a, b, c) = ( ", prim.props.a, ", ", prim.props.b, ", ", prim.props.c, " )"
 30 
 31 #
 32 #  New CLASS backend precision object
 33 #  Lets also increase k_per_decade_primordial since we are
 34 #  dealing with a modified spectrum.
 35 #
 36 cbe_prec = Nc.CBEPrecision.new ()
 37 cbe_prec.props.k_per_decade_primordial = 50.0
 38 
 39 #
 40 #  New CLASS backend object
 41 #
 42 cbe = Nc.CBE.prec_new (cbe_prec)
 43 
 44 #
 45 #  Set precision parameters
 46 #
 47 #ser.set_property_from_key_file (cbe_prec, "cl_ref.pre")
 48 
 49 #
 50 #  New CLASS backend object
 51 #
 52 Bcbe = Nc.HIPertBoltzmannCBE.full_new (cbe)
 53 Bcbe.set_TT_lmax (lmax)
 54 Bcbe.set_target_Cls (Nc.DataCMBDataType.TT)
 55 Bcbe.set_lensed_Cls (True)
 56 
 57 #
 58 # Makes sure that Bcbe contain the default values
 59 #
 60 #Bcbe.props.precision.assert_default ()
 61 
 62 #
 63 #  New homogeneous and isotropic cosmological model NcHICosmoDEXcdm
 64 #
 65 cosmo = Nc.HICosmo.new_from_name (Nc.HICosmo, "NcHICosmoDEXcdm")
 66 cosmo.omega_x2omega_k ()
 67 cosmo.param_set_by_name ("Omegak", 0.0)
 68 
 69 #
 70 #  New homogeneous and isotropic reionization object
 71 #
 72 reion = Nc.HIReionCamb.new ()
 73 
 74 #
 75 # Adding submodels to the main cosmological model.
 76 #
 77 cosmo.add_submodel (reion)
 78 cosmo.add_submodel (prim)
 79 
 80 #
 81 # Preparing the Class backend object
 82 #
 83 #Bcbe.prepare (prim, reion, cosmo)
 84 Bcbe.prepare (cosmo)
 85 
 86 Cls1 = Ncm.Vector.new (lmax + 1)
 87 Cls2 = Ncm.Vector.new (lmax + 1)
 88 
 89 Bcbe.get_TT_Cls (Cls1)
 90 
 91 prim.props.a = 0
 92 Bcbe.prepare (cosmo)
 93 Bcbe.get_TT_Cls (Cls2)
 94 
 95 Cls1_a = Cls1.dup_array ()
 96 Cls2_a = Cls2.dup_array ()
 97 
 98 Cls1_a = np.array (Cls1_a[2:])
 99 Cls2_a = np.array (Cls2_a[2:])
100 
101 ell = np.array (range (2, lmax + 1))
102 
103 Cls1_a = ell * (ell + 1.0) * Cls1_a
104 Cls2_a = ell * (ell + 1.0) * Cls2_a
105 
106 #
107 #  Ploting ionization history.
108 #
109 
110 plt.title (r'Modified and non-modified $C_\ell$')
111 plt.xscale('log')
112 plt.plot (ell, Cls1_a, 'r', label="Modified")
113 plt.plot (ell, Cls2_a, 'b--', label="Non-modified")
114 
115 plt.xlabel(r'$\ell$')
116 plt.ylabel(r'$C_\ell$')
117 plt.legend(loc=2)
118 
119 plt.savefig ("hiprim_Cls.png")
120 
121 plt.clf ()

This example produces:

Cls

In this example we used a implementation of NcHIPrim done in python. The primordial power spectrum is modified by multiplying the usual power law by an oscillatory factor, i.e., $P(k) = A_s\left(\frac{k}{k_\text{pivot}}\right)^{n_s - 1}\left[1 + a^2\cos\left(b\frac{k}{k_\text{pivot}} + c\right)^2\right]$.

  1 import gi
  2 import math
  3 
  4 gi.require_version('NumCosmo', '1.0')
  5 gi.require_version('NumCosmoMath', '1.0')
  6 
  7 from gi.repository import GObject
  8 from gi.repository import NumCosmo as Nc
  9 from gi.repository import NumCosmoMath as Ncm
 10 
 11 #
 12 #  Initializing the library objects, this must be called before
 13 #  any other library function.
 14 #
 15 Ncm.cfg_init ()
 16 
 17 #
 18 # New ModelBuilder object, defines a new model NcHIPrimExample implementing
 19 # the Nc.HIPrim abstract class.
 20 # 
 21 mb = Ncm.ModelBuilder.new (Nc.HIPrim, "NcHIPrimExample", "A example primordial model")
 22 
 23 #
 24 # New parameter A_s to describe the spectrum amplitud (it is usually better
 25 # to work with ln(10^10As))
 26 #
 27 mb.add_sparam ("A_s", "As", 0.0, 1.0, 0.1, 0.0, 1.0e-9, Ncm.ParamType.FREE)
 28 
 29 #
 30 # New parameter n_s to describe the spectral index
 31 #
 32 mb.add_sparam ("n_s", "ns", 0.5, 1.5, 0.1, 0.0, 0.96,   Ncm.ParamType.FREE)
 33 
 34 #
 35 # New parameters a, b and c to describe the spectrum modification
 36 #
 37 mb.add_sparam ("a", "a", 0.0,   1.0, 0.01, 0.0,   0.5, Ncm.ParamType.FREE)
 38 mb.add_sparam ("b", "b", 0.0, 1.0e4, 0.10, 0.0, 100.0, Ncm.ParamType.FREE)
 39 mb.add_sparam ("c", "c", 0.0,   6.0, 0.10, 0.0,   0.0, Ncm.ParamType.FREE)
 40 
 41 #
 42 # Creates a new GObject, it is not a Python object yet!
 43 #
 44 GNcHIPrimExample = mb.create ()
 45 
 46 #
 47 # Workaraound to make the pygobject "see" the new object
 48 #
 49 GObject.new (GNcHIPrimExample)
 50 
 51 #
 52 # Gets the Python version of the object (.pytype) and register
 53 # it as a PyGObject object.
 54 #
 55 NcHIPrimExample = GNcHIPrimExample.pytype
 56 GObject.type_register (NcHIPrimExample)
 57 
 58 #
 59 # Creating a new class implementing our object NcHIPrimExample
 60 #
 61 class PyHIPrimExample (NcHIPrimExample):
 62   #
 63   # Defining some property which is not part of the model paramers.
 64   # All model parameters must be defined by the ModelBuilder.
 65   #
 66   some_property = GObject.Property (type = str)
 67   
 68   #
 69   # Implementing the adiabatic spectrum function
 70   #
 71   def do_lnSA_powspec_lnk (self, lnk):
 72     lnk0 = self.get_lnk_pivot ()
 73     lnka = lnk - lnk0
 74     ka = math.exp (lnka)
 75     As = self.props.As
 76     ns = self.props.ns
 77     a  = self.props.a
 78     b  = self.props.b
 79     c  = self.props.c
 80     a2 = a * a
 81     
 82     return (ns - 1.0) * lnka + math.log (As) + math.log1p (a2 * math.cos (b * ka + c)**2) 
 83 
 84 #
 85 # Register our new Python class PyHIPrimExample
 86 #
 87 GObject.type_register(PyHIPrimExample)

An example of power spectrum in Python. Download the source.

To try this example you must have PyGObject installed, matplotlib matplotlib and numcosmo built with --enable-introspection option.

  1 #!/usr/bin/python2
  2 
  3 import gi
  4 gi.require_version('NumCosmo', '1.0')
  5 gi.require_version('NumCosmoMath', '1.0')
  6 
  7 import sys
  8 import time
  9 from math import *
 10 import numpy as np
 11 from gi.repository import GObject
 12 import matplotlib
 13 import matplotlib.pyplot as plt
 14 from gi.repository import NumCosmo as Nc
 15 from gi.repository import NumCosmoMath as Ncm
 16 
 17 matplotlib.rcParams.update({'font.size': 11})
 18 
 19 #
 20 #  Initializing the library objects, this must be called before 
 21 #  any other library function.
 22 #
 23 Ncm.cfg_init ()
 24 
 25 #
 26 #  New homogeneous and isotropic cosmological model NcHICosmoDEXcdm 
 27 #
 28 cosmo = Nc.HICosmo.new_from_name (Nc.HICosmo, "NcHICosmoDEXcdm")
 29 cosmo.omega_x2omega_k ()
 30 cosmo.param_set_by_name ("Omegak", 0.0)
 31 cosmo.param_set_by_name ("w", -1.0)
 32 cosmo.param_set_by_name ("Omegab", 0.04909244421)
 33 cosmo.param_set_by_name ("Omegac", 0.26580755578)
 34 
 35 reion = Nc.HIReionCamb.new ()
 36 prim  = Nc.HIPrimPowerLaw.new ()
 37 
 38 cosmo.param_set_by_name ("H0", 67.31)
 39 
 40 prim.param_set_by_name ("n_SA", 0.9658)
 41 prim.param_set_by_name ("ln10e10ASA", 3.0904)
 42 
 43 reion.param_set_by_name ("z_re", 9.9999)
 44 
 45 cosmo.add_submodel (reion)
 46 cosmo.add_submodel (prim)
 47 
 48 ps_cbe  = Nc.PowspecMLCBE.new ()
 49 ps_eh   = Nc.PowspecMLTransfer.new (Nc.TransferFuncEH.new ())
 50 
 51 z_min = 0.0
 52 z_max = 2.0
 53 zdiv  = 0.49999999999
 54 
 55 k_min = 1.0e-5
 56 k_max = 1.0e3
 57 
 58 nk = 2000
 59 nR = 2000
 60 Rh8 = 8.0 / cosmo.h ()
 61 
 62 ps_cbe.set_kmin (k_min)
 63 ps_eh.set_kmin (k_min)
 64 
 65 ps_cbe.set_kmax (k_max)
 66 ps_eh.set_kmax (k_max)
 67 
 68 ps_cbe.require_zi (z_min)
 69 ps_cbe.require_zf (z_max)
 70 
 71 ps_eh.require_zi (z_min)
 72 ps_eh.require_zf (z_max)
 73 
 74 ps_eh.prepare (cosmo)
 75 ps_cbe.prepare (cosmo)
 76 
 77 for z in np.arange (z_min, z_max, (z_max - z_min) * zdiv):
 78   k_a = []
 79   Pk_eh_a = []
 80   Pk_cbe_a = []
 81   for lnk in np.arange (log (ps_cbe.props.kmin), log (ps_cbe.props.kmax), log (k_max / k_min) / nk):
 82     k = exp (lnk)
 83     k2 = k * k
 84     k3 = k2 * k
 85     k_a.append (k)
 86     Pk_cbe_a.append (k3 * ps_cbe.eval (cosmo, z, k))
 87     Pk_eh_a.append (k3 * ps_eh.eval (cosmo, z, k))
 88     #print k, ps_eh.eval (cosmo, z, k) / ps_cbe.eval (cosmo, z, k)
 89   
 90   plt.plot (k_a, Pk_cbe_a, label = r'CLASS $z = %.2f$' % (z))
 91   plt.plot (k_a, Pk_eh_a, label  = r'EH    $z = %.2f$' % (z))
 92 
 93 plt.legend (loc="lower right")
 94 plt.xscale('log')
 95 plt.yscale('log')
 96 plt.savefig ("ps_cbe_eh.png")
 97 plt.clf ()
 98 
 99 for z in np.arange (z_min, z_max, (z_max - z_min) * zdiv * 0.5):
100   k_a = []
101   Pk_eh_a = []
102   Pk_cbe_a = []
103   for lnk in np.arange (log (ps_cbe.props.kmin), log (ps_cbe.props.kmax), log (k_max / k_min) / nk):
104     k = exp (lnk)
105     k2 = k * k
106     k3 = k2 * k
107     k_a.append (k)
108     Pk_cbe_a.append (k3 * ps_cbe.eval (cosmo, z, k))
109     Pk_eh_a.append (k3 * ps_eh.eval (cosmo, z, k))
110     #print k, ps_eh.eval (cosmo, z, k) / ps_cbe.eval (cosmo, z, k)
111   
112   plt.plot (k_a, np.abs (1.0 - np.array (Pk_eh_a) / np.array (Pk_cbe_a)), label = r'CLASS EH cmp $z = %.2f$' % (z))
113 
114 plt.legend (loc="lower right")
115 plt.xscale('log')
116 plt.yscale('log')
117 plt.savefig ("ps_diff_cbe_eh.png")
118 plt.clf ()
119 
120 #
121 # Filtering 
122 #
123 psf_cbe = Ncm.PowspecFilter.new (ps_cbe, Ncm.PowspecFilterType.TOPHAT)
124 psf_eh  = Ncm.PowspecFilter.new (ps_eh, Ncm.PowspecFilterType.TOPHAT)
125 
126 psf_cbe.set_best_lnr0 ()
127 psf_eh.set_best_lnr0 ()
128 
129 psf_cbe.prepare (cosmo)
130 psf_eh.prepare (cosmo)
131 
132 
133 print "# CBE sigma8 = % 20.15g, EH sigma8 = % 20.15g" % (psf_cbe.eval_sigma (0.0, Rh8), psf_eh.eval_sigma (0.0, Rh8))
134 print "# kmin % 20.15g kmax % 20.15g" % (ps_cbe.props.kmin, ps_cbe.props.kmax)
135 print "# Rmin % 20.15g Rmax % 20.15g" % (psf_cbe.get_r_min (), psf_cbe.get_r_max ())
136 
137 lnRmin = log (psf_cbe.get_r_min ())
138 lnRmax = log (psf_cbe.get_r_max ())
139 
140 for z in np.arange (z_min, z_max, (z_max - z_min) * zdiv):
141   Rh_a = []
142   sigma_eh_a  = []
143   sigma_cbe_a = []
144   for lnR in np.arange (lnRmin, lnRmax, (lnRmax - lnRmin) / nR):
145     R  = exp (lnR)
146     Rh = R * cosmo.h ()
147     
148     Rh_a.append (Rh)
149     sigma_cbe_a.append (psf_cbe.eval_sigma_lnr (z, lnR))
150     sigma_eh_a.append (psf_eh.eval_sigma_lnr (z, lnR))
151 
152   plt.plot (Rh_a, sigma_cbe_a, label = r'$\sigma$ CLASS $z = %.2f$' % (z))
153   plt.plot (Rh_a, sigma_cbe_a, label = r'$\sigma$ EH    $z = %.2f$' % (z))
154 
155 plt.legend (loc="lower left")
156 plt.xscale('log')
157 plt.yscale('log')
158 plt.savefig ("ps_var_cbe_eh.png")
159 plt.clf ()
160 
161 for z in np.arange (z_min, z_max, (z_max - z_min) * zdiv):
162   Rh_a = []
163   dvar_eh_a   = []
164   dvar_cbe_a  = []
165   for lnR in np.arange (lnRmin, lnRmax, (lnRmax - lnRmin) / nR):
166     R  = exp (lnR)
167     Rh = R * cosmo.h ()
168     
169     Rh_a.append (Rh)
170     dvar_cbe_a.append (psf_cbe.eval_dlnvar_dlnr (z, lnR))
171     dvar_eh_a.append (psf_eh.eval_dlnvar_dlnr (z, lnR))
172 
173   plt.plot (Rh_a, dvar_cbe_a, label  = r'$\frac{\mathrm{d}\ln\sigma^2}{\mathrm{d}\ln{}R}$ CLASS $z = %.2f$' % (z))
174   plt.plot (Rh_a, dvar_eh_a, label   = r'$\frac{\mathrm{d}\ln\sigma^2}{\mathrm{d}\ln{}R}$ EH    $z = %.2f$' % (z))
175 
176 plt.legend (loc="lower left")
177 plt.xscale('log')
178 plt.savefig ("ps_dvar_cbe_eh.png")
179 plt.clf ()
180 
181 zmaxnl = 10.0
182 z_max  = zmaxnl
183 pshf   = Nc.PowspecMNLHaloFit.new (ps_cbe, zmaxnl, 1.0e-5)
184 
185 pshf.set_kmin (k_min)
186 pshf.set_kmax (k_max)
187 pshf.require_zi (z_min)
188 pshf.require_zf (z_max)
189 
190 pshf.prepare (cosmo)
191 
192 for z in np.arange (z_min, z_max, (z_max - z_min) * zdiv):
193   k_a = []
194   Pk_cbe_a = []
195   Pknl_cbe_a = []
196   for lnk in np.arange (log (ps_cbe.props.kmin), log (ps_cbe.props.kmax), log (k_max / k_min) / nk):
197     k = exp (lnk)
198     k2 = k * k
199     k3 = k2 * k
200     k_a.append (k)
201     Pk_cbe_a.append (k3 * ps_cbe.eval (cosmo, z, k))
202     Pknl_cbe_a.append (k3 * pshf.eval (cosmo, z, k))
203   
204   plt.plot (k_a, Pk_cbe_a, label = r'CLASS $z = %.2f$' % (z))
205   plt.plot (k_a, Pknl_cbe_a, label  = r'CLASS+HaloFit $z = %.2f$' % (z))
206 
207 plt.legend (loc="lower right")
208 plt.xscale('log')
209 plt.yscale('log')
210 plt.savefig ("ps_cbe_halofit.png")
211 plt.clf ()
212 
213 psf_cbenl = Ncm.PowspecFilter.new (pshf, Ncm.PowspecFilterType.TOPHAT)
214 psf_cbenl.set_best_lnr0 ()
215 psf_cbenl.prepare (cosmo)
216 
217 print "# CBE+HaloFit sigma8 = % 20.15g" % (psf_cbenl.eval_sigma (0.0, Rh8))
218 
# CBE sigma8 =     0.84625009116749, EH sigma8 =    0.852830701094344
# kmin                1e-05 kmax                 1000
# Rmin  0.000999999999999999 Rmax               100000
# CBE+HaloFit sigma8 =    0.956987437165177

This example has the following output and figure:

CLASS and EH power spectra CLASS and EH power spectra comparison CLASS and EH power spectra + Halofit CLASS and EH power spectra filtered CLASS and EH power spectra filtered derivative

Contact

Please report bugs, ask questions and send comments/suggestions at numcosmo-help@nongnu.org

Valid XHTML 1.0 Strict

Valid CSS!