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.1 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.These are not tested yet.

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_SPECINDEX, 1.00);
 39   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_SIGMA8,    0.90);
 40   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_XCDM_W,   -1.10);
 41 
 42   /**************************************************************************** 
 43    * Printing the parameters used.
 44    ****************************************************************************/
 45   printf ("# Model parameters:\n"); 
 46   ncm_model_params_log_all (NCM_MODEL (cosmo));
 47 
 48   /**************************************************************************** 
 49    * Printing some distances up to redshift 1.0.
 50    ****************************************************************************/ 
 51   for (i = 0; i < 10; i++)
 52   {
 53     gdouble z = 1.0 / 9.0 * i;
 54     gdouble cd = ncm_c_hubble_radius () * nc_distance_comoving (dist, cosmo, z);
 55     printf ("% 10.8f % 20.15g\n", z, cd);
 56   }
 57 
 58   /**************************************************************************** 
 59    * Freeing objects.
 60    ****************************************************************************/ 
 61   nc_distance_free (dist);
 62   ncm_model_free (NCM_MODEL (cosmo));
 63 
 64   return 0;
 65 }

This example has the following output:

# Model parameters:
                    70                  0.25                   0.7                  2.72                  0.24                 3.046                  0.05                     1                   0.9                  -1.1
 0.00000000                    0
 0.11111111      326.28966196004
 0.22222222     637.556054363559
 0.33333333     932.377752746278
 0.44444444     1210.13776409848
 0.55555556     1470.85990963654
 0.66666667       1715.030465771
 0.77777778     1943.43693407981
 0.88888889     2157.03863437956
 1.00000000     2356.87151080258

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

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

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.ns      = 1.0
 39 cosmo.props.sigma8  = 0.9
 40 cosmo.props.w       = -1.0
 41 
 42 #
 43 #  Creating a new Modelset and set cosmo as the HICosmo model to be used.
 44 #
 45 mset = Ncm.MSet ()
 46 mset.set (cosmo)
 47 
 48 #
 49 #  Setting parameters Omega_c and w to be fitted and change parameter
 50 #  Omega_x -> Omega_k.
 51 #
 52 ##cosmo.de_omega_x2omega_k ()
 53 cosmo.props.Omegac_fit = True
 54 cosmo.props.w_fit = True
 55 
 56 #
 57 #  Creating a new Distance object optimized to redshift 2.
 58 #
 59 dist = Nc.Distance (zf = 2.0)
 60 
 61 #
 62 #  Creating a new Data object from distance modulus catalogs.
 63 #
 64 snia = Nc.DataDistMu.new_from_id (dist, Nc.DataSNIAId.SIMPLE_UNION2_1)
 65 
 66 #
 67 #  Creating a new Dataset and add snia to it.
 68 #
 69 dset = Ncm.Dataset ()
 70 dset.append_data (snia)
 71 
 72 #
 73 #  Creating a Likelihood from the Dataset.
 74 #
 75 lh = Ncm.Likelihood (dataset = dset)
 76 
 77 #
 78 #  Creating a Fit object of type NLOPT using the fitting algorithm ln-neldermead to
 79 #  fit the Modelset mset using the Likelihood lh and using a numerical differentiation
 80 #  algorithm (NUMDIFF_FORWARD) to obtain the gradient (if needed).
 81 #
 82 fit = Ncm.Fit.new (Ncm.FitType.NLOPT, "ln-neldermead", lh, mset, Ncm.FitGradType.NUMDIFF_FORWARD)
 83 
 84 #
 85 #  Running the fitter printing messages.
 86 #
 87 fit.run (Ncm.FitRunMsgs.SIMPLE)
 88 
 89 #
 90 #  Printing fitting informations.
 91 #
 92 fit.log_info ()
 93 
 94 #
 95 #  Calculating the parameters covariance using numerical differentiation.
 96 #
 97 fit.numdiff_m2lnL_covar ()
 98 
 99 #
100 #  Printing the covariance matrix.
101 # 
102 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.ns      = 1.0
 39 cosmo.props.sigma8  = 0.9
 40 cosmo.props.w       = -1.0
 41 
 42 #
 43 #  A new Modelset with cosmo as the HICosmo model to be used.
 44 #
 45 mset = Ncm.MSet ()
 46 mset.set (cosmo)
 47 
 48 #
 49 #  Setting parameters Omega_c, Omega_x and w to be fitted (and change 
 50 #  parameter Omega_x -> Omega_k).
 51 #
 52 #cosmo.de_omega_x2omega_k ()
 53 cosmo.props.Omegac_fit = True
 54 cosmo.props.Omegax_fit = True
 55 cosmo.props.w_fit = True
 56 
 57 #
 58 #  A new Distance object optimized to redshift 2.
 59 #
 60 dist = Nc.Distance (zf = 2.0)
 61 
 62 #
 63 #  A new Data object from distance modulus catalogs.
 64 #  A new Data object from BAO catalogs.
 65 #
 66 snia = Nc.DataDistMu.new_from_id (dist, Nc.DataSNIAId.SIMPLE_UNION2_1)
 67 bao = Nc.data_bao_create (dist, Nc.DataBaoId.A_EISENSTEIN2005)
 68 
 69 #
 70 #  A new Dataset with snia and bao set.
 71 #
 72 dset = Ncm.Dataset ()
 73 dset.append_data (snia)
 74 dset.append_data (bao)
 75 
 76 #
 77 #  Creating a Likelihood from the Dataset.
 78 #
 79 lh = Ncm.Likelihood (dataset = dset)
 80 
 81 #
 82 #  Creating a Fit object of type NLOPT using the fitting algorithm ln-neldermead to
 83 #  fit the Modelset mset using the Likelihood lh and using a numerical differentiation
 84 #  algorithm (NUMDIFF_FORWARD) to obtain the gradient (if needed).
 85 #
 86 fit = Ncm.Fit.new (Ncm.FitType.NLOPT, "ln-neldermead", lh, mset, Ncm.FitGradType.NUMDIFF_FORWARD)
 87 
 88 #
 89 #  Running the fitter printing messages.
 90 #
 91 fit.run (Ncm.FitRunMsgs.SIMPLE)
 92 
 93 #
 94 #  Printing fitting informations.
 95 #
 96 fit.log_info ()
 97 
 98 #
 99 #  Calculating the parameters covariance using numerical differentiation.
100 #
101 fit.numdiff_m2lnL_covar ()
102 
103 #
104 #  Printing the covariance matrix.
105 # 
106 fit.log_covar ()
107 
108 #
109 #  Creating a new Likelihood ratio test object.
110 #  First we create two PIndex indicating which parameter
111 #    we are going to study.
112 # 
113 p1 = Ncm.MSetPIndex.new (cosmo.id (), Nc.HICosmoDEParams.OMEGA_C)
114 p2 = Ncm.MSetPIndex.new (cosmo.id (), Nc.HICosmoDEXCDMParams.W)
115 
116 lhr2d = Ncm.LHRatio2d.new (fit, p1, p2, 1.0e-3)
117 
118 #
119 #  Calculating the confidence region using the Likelihood ratio test.
120 #  Also calculate using the Fisher matrix approach.
121 #
122 cr_rg = lhr2d.conf_region (0.6826, 300.0, Ncm.FitRunMsgs.SIMPLE)
123 fisher_rg = lhr2d.fisher_border (0.6826, 300.0, Ncm.FitRunMsgs.SIMPLE)
124 
125 cr_p1array = cr_rg.p1.dup_array ()
126 cr_p2array = cr_rg.p2.dup_array ()
127 
128 fisher_p1array = fisher_rg.p1.dup_array ()
129 fisher_p2array = fisher_rg.p2.dup_array ()
130 
131 #
132 #  Ploting the confidence regions obtained from both methods.
133 #
134 
135 plt.title ("Confidence regions (%f)" % (cr_rg.clevel * 100.0))
136 plt.plot (cr_p1array, cr_p2array, 'r', label="Likelihood Ratio")
137 plt.plot (fisher_p1array, fisher_p2array, 'b-', label="Fisher Matrix")
138 
139 plt.xlabel(r'$\Omega_c$')
140 plt.ylabel(r'$w$')
141 
142 plt.legend(loc=4)
143 
144 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   NcHIReion *reion;
  8   NcHICosmo *cosmo;
  9   NcDistance *dist;
 10   NcWindow *wp;
 11   NcTransferFunc *tf;
 12   NcMatterVar *vp;
 13   NcGrowthFunc *gf;
 14   NcMultiplicityFunc *mulf;
 15   NcMassFunction *mf;
 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 cosmological distance objects optimizied to perform calculations
 37    * up to redshift 2.0.
 38    ****************************************************************************/  
 39   dist = nc_distance_new (2.0);
 40  
 41   /**************************************************************************** 
 42    * New windown function 'NcWindowTophat'
 43    ****************************************************************************/  
 44   wp = nc_window_new_from_name ("NcWindowTophat");
 45 
 46   /**************************************************************************** 
 47    * New transfer function 'NcTransferFuncEH' using the Einsenstein, Hu
 48    * fitting formula.
 49    ****************************************************************************/  
 50   tf = nc_transfer_func_new_from_name ("NcTransferFuncEH");
 51 
 52   /**************************************************************************** 
 53    * New matter variance object using FFT method for internal calculations and
 54    * the window and transfer functions defined above.
 55    ****************************************************************************/  
 56   vp = nc_matter_var_new (NC_MATTER_VAR_FFT, wp, tf);
 57 
 58   /**************************************************************************** 
 59    * New growth function
 60    ****************************************************************************/  
 61   gf = nc_growth_func_new ();
 62 
 63   /**************************************************************************** 
 64    * New multiplicity function 'NcMultiplicityFuncTinkerMean'
 65    ****************************************************************************/  
 66   mulf = nc_multiplicity_func_new_from_name ("NcMultiplicityFuncTinkerMean");
 67 
 68   /**************************************************************************** 
 69    * New mass function object using the objects defined above.
 70    ****************************************************************************/  
 71   mf = nc_mass_function_new (dist, vp, gf, mulf);
 72 
 73   /**************************************************************************** 
 74    * Setting values for the cosmological model, those not set stay in the
 75    * default values. Remeber to use the _orig_ version to set the original
 76    * parameters in case when a reparametrization is used.
 77    ****************************************************************************/ 
 78   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_H0,        70.0);
 79   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_OMEGA_C,    0.25);
 80   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_OMEGA_X,    0.7);
 81   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_T_GAMMA0,   1.0);
 82   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_OMEGA_B,    0.05);
 83   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_SPECINDEX,  1.0);
 84   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_SIGMA8,     0.9);
 85   ncm_model_orig_param_set (NCM_MODEL (cosmo), NC_HICOSMO_DE_XCDM_W,    -1.1);
 86 
 87   /**************************************************************************** 
 88    * Printing the parameters used.
 89    ****************************************************************************/
 90   printf ("# Model parameters:\n#"); 
 91   ncm_model_params_log_all (NCM_MODEL (cosmo));
 92 
 93   /**************************************************************************** 
 94    * Printing the growth function and its derivative with respect to z 
 95    * up to redshift 1.
 96    ****************************************************************************/ 
 97   nc_growth_func_prepare (gf, cosmo);
 98 
 99   for (i = 0; i < np; i++)
100   {
101     gdouble z = 1.0 / (np - 1.0) * i;
102     gdouble gfz = nc_growth_func_eval (gf, cosmo, z);
103     gdouble dgfz = nc_growth_func_eval_deriv (gf, cosmo, z);
104     printf ("% 10.8f % 20.15g % 20.15g\n", z, gfz, dgfz);
105   }
106   printf ("\n\n");
107 
108   /**************************************************************************** 
109    * Printing the transfer function and the matter power spectrum in the 
110    * kh (in unities of h/Mpc) interval [1e-3, 1e3] 
111    ****************************************************************************/
112   nc_transfer_func_prepare (tf, reion, cosmo);
113 
114   for (i = 0; i < np; i++)
115   {
116     gdouble lnkh = log (1e-3) + log (1e6) / (np - 1.0) * i;
117     gdouble tfkh = nc_transfer_func_eval (tf, cosmo, exp (lnkh));
118     gdouble Pmkh = nc_transfer_func_matter_powerspectrum (tf, cosmo, exp (lnkh));
119     printf ("% 10.8f % 20.15g % 20.15g\n", exp (lnkh), tfkh, Pmkh);
120   }
121   printf ("\n\n");
122 
123   /**************************************************************************** 
124    * Printing the variance filtered with the tophat windown function using
125    * scales R in the interval [5, 50] at redshift 0.3.
126    * First calculates the growth function at z = 0.3 and then the spectrum
127    * amplitude from the sigma8 parameter.
128    ****************************************************************************/
129   nc_matter_var_prepare (vp, reion, cosmo);
130   {
131     gdouble Dz = nc_growth_func_eval (gf, cosmo, 0.3);
132     gdouble A = nc_matter_var_sigma8_sqrtvar0 (vp, cosmo);
133     gdouble prefac = Dz * Dz * A * A;
134     
135     for (i = 0; i < np; i++)
136     {
137       gdouble lnR = log (5.0) + log (10.0) / (np - 1.0) * i;
138       gdouble sigma2 = prefac * nc_matter_var_var0 (vp, cosmo, lnR);
139       gdouble dsigma2_dlnR = nc_matter_var_dlnvar0_dlnR (vp, cosmo, lnR);
140       printf ("% 10.8f % 20.15g % 20.15g\n", exp (lnR), sigma2, dsigma2_dlnR);
141     }
142     printf ("\n\n");
143   }
144 
145   /**************************************************************************** 
146    * Printing the mass function integrated in the mass interval [1e14, 1e16]
147    * for the redhshifts in the interval [0, 2.0] and area 200 squared degree.
148    ****************************************************************************/
149   nc_mass_function_set_area_sd (mf, 200.0);
150   nc_mass_function_set_eval_limits (mf, cosmo, log (1e14), log (1e16), 0.0, 2.0);
151   nc_mass_function_prepare (mf, reion, cosmo);
152 
153   for (i = 0; i < np; i++)
154   {
155     gdouble z = 2.0 / (np - 1.0) * i;
156     gdouble dndz = nc_mass_function_dn_dz (mf, cosmo, log(1e14), log(1e16), z, FALSE);
157     printf ("% 10.8f % 20.15g\n", z, dndz);
158   }
159   printf ("\n\n");
160 
161 
162   /**************************************************************************** 
163    * Freeing objects.
164    ****************************************************************************/ 
165   nc_distance_free (dist);
166   ncm_model_free (NCM_MODEL (cosmo));
167   nc_hireion_free (reion);
168   nc_window_free (wp);
169   nc_transfer_func_free (tf);
170   nc_matter_var_free (vp);
171   nc_growth_func_free (gf);
172   nc_multiplicity_func_free (mulf);
173   nc_mass_function_free (mf);
174 
175   return 0;
176 }

This example has the following output:

# Model parameters:
#                    70                  0.25                   0.7                     1                  0.24                 3.046                  0.05                     1                   0.9                  -1.1
 0.00000000                    1   -0.514645211187295
 0.02040816    0.989521827289464   -0.512155429309301
 0.04081633    0.979098119617743   -0.509310739491474
 0.06122449     0.96873589196847   -0.506133995937869
 0.08163265    0.958441692579252   -0.502647773351458
 0.10204082    0.948221611431494   -0.498874223297092
 0.12244898    0.938081287321928   -0.494835143213762
 0.14285714     0.92802591533946   -0.490552348617287
 0.16326531    0.918060256604688   -0.486046753489208
 0.18367347     0.90818865174697   -0.481338498434204
 0.20408163    0.898415034030169   -0.476447154136105
 0.22448980    0.888742942225605   -0.471391972970399
 0.24489796    0.879175535760649   -0.466191038842278
 0.26530612     0.86971561170594    -0.46086158890044
 0.28571429    0.860365620892682   -0.455420201147502
 0.30612245    0.851127683346253   -0.449882875589346
 0.32653061    0.842003605808585   -0.444264418958837
 0.34693878    0.832994899455343   -0.438578803761855
 0.36734694    0.824102796531267   -0.432839360624843
 0.38775510    0.815328266181285   -0.427058678812051
 0.40816327     0.80667203182617   -0.421248291402332
 0.42857143    0.798134587799352    -0.41541898537597
 0.44897959    0.789716214668303    -0.40958094757446
 0.46938776    0.781416993968273   -0.403743512371455
 0.48979592    0.773236824548903   -0.397915319647694
 0.51020408    0.765175435290975   -0.392104312509405
 0.53061224    0.757232400220416   -0.386317785454479
 0.55102041    0.749407150363481   -0.380562410233295
 0.57142857    0.741698987365771   -0.374844248405062
 0.59183673    0.734107094314609   -0.369168817507294
 0.61224490    0.726630547895928    -0.36354107332947
 0.63265306    0.719268328149603   -0.357965498840659
 0.65306122    0.712019329184082   -0.352446072148664
 0.67346939    0.704882367998991   -0.346986367842658
 0.69387755    0.697856193420111    -0.34158947062989
 0.71428571    0.690939494583063   -0.336258151869116
 0.73469388    0.684130908145143   -0.330994813726842
 0.75510204    0.677429025592986   -0.325801493511353
 0.77551020    0.670832399867168   -0.320679920128959
 0.79591837    0.664339551447366   -0.315631539169878
 0.81632653    0.657948973975885    -0.31065753690592
 0.83673469    0.651659139417519   -0.305758865804175
 0.85714286    0.645468502902642   -0.300936232632084
 0.87755102    0.639375506990631   -0.296190150623245
 0.89795918    0.633378585915658   -0.291520941518565
 0.91836735    0.627476169011182   -0.286928756593282
 0.93877551    0.621666684160382   -0.282413590835157
 0.95918367    0.615948560791005   -0.277975298428775
 0.97959184    0.610320232624505   -0.273613605766789
 1.00000000    0.604780140148379   -0.269328124307583


 0.00100000    0.999805328836874  0.00099961069557061
 0.00132571     0.99966036548697   0.0013248110038449
 0.00175751    0.999408916610617  0.00175543356821907
 0.00232995    0.998974629521377  0.00232517613258463
 0.00308884    0.998228668319402  0.00307791055504462
 0.00409492    0.996956115835978  0.00407002410849081
 0.00542868    0.994802993737763  0.00537239634122765
 0.00719686    0.991193154320994  0.00707065170884944
 0.00954095    0.985194564691879  0.00926053017245344
 0.01264855    0.975287700876948   0.0120311270012985
 0.01676833    0.958925135476369    0.015419106248539
 0.02222996    0.931708734586102   0.0192974137879651
 0.02947052    0.886491810088735   0.0231599282977442
 0.03906940    0.815930668834037   0.0260101735325903
 0.05179475    0.726956269760714   0.0273717325212831
 0.06866488    0.652751656188497   0.0292570584059584
 0.09102982    0.608453134180805    0.033700623703963
 0.12067926    0.553582842935875   0.0369826388441461
 0.15998587    0.477966307235999   0.0365490589605499
 0.21209509    0.407543542501762   0.0352272421382249
 0.28117687    0.334015235521443   0.0313698285817094
 0.37275937    0.264888077550762   0.0261549198951282
 0.49417134    0.202854253187564   0.0203350753857352
 0.65512856    0.150128256891569   0.0147656107317776
 0.86851137    0.107719872554174   0.0100778333402021
 1.15139540   0.0752208251004834  0.00651479381825275
 1.52641797   0.0513422137708357  0.00402367273917402
 2.02358965   0.0343913610414364  0.00239343245510617
 2.68269580   0.0226840078767555  0.00138041925156272
 3.55648031   0.0147706213896594  0.000775921776189658
 4.71486636  0.00951269449406509  0.000426654653121894
 6.25055193  0.00606782572434713  0.000230136002447234
 8.28642773  0.00383743951101831  0.000122025454123775
 10.98541142  0.00240823254484467  6.3710816194812e-05
 14.56348478  0.00150080180911841  3.28028815112764e-05
 19.30697729  0.000929403783999379  1.66772008206936e-05
 25.59547923  0.000572276426899725   8.382527350344e-06
 33.93221772  0.000350566272802228  4.17015797585198e-06
 44.98432669  0.000213756088939872  2.0554088104965e-06
 59.63623317  0.000129792925234918  1.00464411637065e-06
 79.06043211  7.85140104015962e-05   4.873640672225e-07
 104.81131342  4.73333022366793e-05  2.34823616311431e-07
 138.94954944  2.84479990994849e-05  1.12450293666329e-07
 184.20699693  1.70501130001047e-05  5.35501443336509e-08
 244.20530945  1.01930476529716e-05  2.5372497078206e-08
 323.74575428  6.07966647435448e-06  1.19664020805526e-08
 429.19342601  3.61860031217541e-06  5.61997383816258e-09
 568.98660290  2.14962619374807e-06  2.62922608119628e-09
 754.31200634  1.27471880566688e-06  1.22568763887568e-09
 1000.00000000  7.54666527284974e-07  5.69521567404362e-10


 5.00000000     1.82741749827407    -2.24929762765026
 5.24056567      1.6429318825971    -2.28009947787663
 5.49270571     1.47493241403446    -2.31095510706532
 5.75697700     1.32219257952954    -2.34183286681469
 6.03396320     1.18355153578356    -2.37270076477371
 6.32427608     1.05791325800824    -2.40352736917292
 6.62855683    0.944244617694199    -2.43428104986324
 6.94747747    0.841574098280314    -2.46492970566601
 7.28174239    0.748990090344704    -2.49544217653731
 7.63208984    0.665638924877999    -2.52578787846441
 7.99929360    0.590723206131626    -2.55593595388316
 8.38416468    0.523499595796996     -2.5858560655214
 8.78755312    0.463276880023895    -2.61551870146876
 9.21034985    0.409413777555954    -2.64489339217275
 9.65348864    0.361316639681813     -2.6739499764721
 10.11794824    0.318437327223747    -2.70266205737228
 10.60475444    0.280270728486028    -2.73100549380315
 11.11498241    0.246352524317859    -2.75895817056931
 11.64975905    0.216256917494735    -2.78649830826439
 12.21026547    0.189594331622082    -2.81360304121439
 12.79773961    0.166009225462142    -2.84024511306227
 13.41347898    0.145177827356924    -2.86639976168025
 14.05884349    0.126805984664047     -2.8920475942434
 14.73525851    0.110626993078777    -2.91717223761638
 15.44421798   0.0963995862965737    -2.94175681427814
 16.18728771   0.0839060393075224    -2.96577688699737
 16.96610886   0.0729503048055597    -2.98920061142061
 17.78240153   0.0633562382468015    -3.01199019960636
 18.63796860   0.0549658740629602     -3.0341280685365
 19.53469969   0.0476378265108445    -3.05559623589061
 20.47457531   0.0412458222256763    -3.07633382073925
 21.45967130    0.035677355172138    -3.09627714098131
 22.49216334   0.0308323789745028    -3.11537003140223
 23.57433182   0.0266220242410552     -3.1335430076925
 24.70856681   0.0229675229284545    -3.15069896513195
 25.89737340    0.019799263662057     -3.1666734130659
 27.14337720   0.0170557865410416    -3.18126467662513
 28.44933015   0.0146829227480062    -3.19422505388484
 29.81811658   0.0126330691989378    -3.20516121869808
 31.25275963   0.0108644371250924    -3.21348908975713
 32.75642784  0.00934051005264519    -3.21826737268197
 34.33244225  0.00802928574680967    -3.21889859337341
 35.98428365  0.00690264753207981    -3.21549443391641
 37.71560032  0.00593537466694092    -3.20976091658665
 39.53021605  0.00510509654493155    -3.20401614384522
 41.43213864  0.00439199268648945    -3.19993910489024
 43.42556869  0.00377897660424108    -3.19875396917298
 45.51490890  0.00325143287213659    -3.20121466693814
 47.70477382  0.00279696196289161    -3.20752079294866
 50.00000000  0.00240508505401569    -3.21755054030204


 0.00000000                    0
 0.04081633     166.601940542437
 0.08163265     615.627205223794
 0.12244898     1274.53402582782
 0.16326531     2076.81279650744
 0.20408163     2963.13202316309
 0.24489796     3882.06577491307
 0.28571429     4790.42546753145
 0.32653061     5653.23313226245
 0.36734694     6443.40382415774
 0.40816327     7141.19218857123
 0.44897959     7733.47500147512
 0.48979592     8212.93009331476
 0.53061224     8577.16648819165
 0.57142857     8827.85186700109
 0.61224490     8969.87358800412
 0.65306122     9010.55832646067
 0.69387755     8958.97018039727
 0.73469388     8825.29467468308
 0.77551020     8620.31472492615
 0.81632653     8354.97605095242
 0.85714286     8040.03846350063
 0.89795918     7685.80569937709
 0.93877551     7301.92566265266
 0.97959184     6897.25197014062
 1.02040816     6479.75768272437
 1.06122449     6056.49232327919
 1.10204082     5633.57396495167
 1.14285714     5216.20878785993
 1.18367347     4808.73150115181
 1.22448980      4414.6608384507
 1.26530612     4036.76516964334
 1.30612245     3677.13424934428
 1.34693878     3337.25364515969
 1.38775510     3018.07938812796
 1.42857143     2720.11049369967
 1.46938776     2443.45866438653
 1.51020408     2187.91194056493
 1.55102041     1952.99557536303
 1.59183673     1738.02534770985
 1.63265306     1542.15539694903
 1.67346939     1364.42108844583
 1.71428571     1203.77563259119
 1.75510204     1059.12164973301
 1.79591837     929.337973389153
 1.83673469     813.301394140972
 1.87755102     709.904895187837
 1.91836735      618.07176072769
 1.95918367     536.766739170608
 2.00000000     465.004159735957


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 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 #  New cosmological distance objects optimizied to perform calculations
 31 #  up to redshift 2.0.
 32 #
 33 dist = Nc.Distance.new (2.0)
 34 
 35 #
 36 # New windown function 'NcWindowTophat'
 37 #
 38 wp =  Nc.Window.new_from_name ("NcWindowTophat")
 39 
 40 #
 41 # New transfer function 'NcTransferFuncEH' using the Einsenstein, Hu
 42 # fitting formula.
 43 #
 44 tf = Nc.TransferFunc.new_from_name ("NcTransferFuncEH")
 45 
 46 #
 47 # New matter variance object using FFT method for internal calculations and
 48 # the window and transfer functions defined above.
 49 #
 50 vp = Nc.MatterVar.new (Nc.MatterVarStrategy.FFT, wp, tf)
 51 
 52 #
 53 # New growth function
 54 #
 55 gf = Nc.GrowthFunc.new ()
 56 
 57 #
 58 # New multiplicity function 'NcMultiplicityFuncTinkerMean'
 59 #
 60 mulf = Nc.MultiplicityFunc.new_from_name ("NcMultiplicityFuncTinkerMean")
 61 
 62 #
 63 # New mass function object using the objects defined above.
 64 #
 65 mf = Nc.MassFunction.new (dist, vp, gf, mulf)
 66 
 67 #
 68 #  Setting values for the cosmological model, those not set stay in the
 69 #  default values. Remember to use the _orig_ version to set the original
 70 #  parameters when a reparametrization is used.
 71 #
 72 cosmo.props.H0      = 70.0
 73 cosmo.props.Omegab  = 0.05
 74 cosmo.props.Omegac  = 0.25
 75 cosmo.props.Omegax  = 0.70
 76 cosmo.props.Tgamma0 = 2.72
 77 cosmo.props.ns      = 1.0
 78 cosmo.props.sigma8  = 0.9
 79 cosmo.props.w       = -1.0
 80 
 81 #
 82 #  Printing the parameters used.
 83 #
 84 print "# Model parameters: ", 
 85 cosmo.params_log_all ()
 86 
 87 #
 88 # Number of points to build the plots
 89 #
 90 np = 2000
 91 divfac = 1.0 / (np - 1.0)
 92 
 93 #
 94 #  Calculating growth and its derivative in the [0, 2] redshift
 95 #  range.
 96 #
 97 za = []
 98 Da = []
 99 dDa = []
100 
101 gf.prepare (cosmo)
102 
103 for i in range (0, np):
104   z = 2.0 * divfac * i
105   D = gf.eval (cosmo, z)
106   dD = gf.eval_deriv (cosmo, z)
107   za.append (z)
108   Da.append (D)
109   dDa.append (dD)
110   print z, D
111 
112 #
113 #  Ploting growth function.
114 #
115 
116 plt.title ("Growth Function")
117 plt.plot (za, Da, 'r', label="D")
118 plt.plot (za, dDa, 'b--', label="dD/dz")
119 plt.xlabel('$z$')
120 plt.legend(loc=2)
121 
122 plt.savefig ("growth_func.png")
123 plt.clf ()
124 
125 #
126 # Calculating the transfer function and the matter power spectrum in the
127 # kh (in unities of h/Mpc) interval [1e-3, 1e3]
128 #
129 
130 kha = []
131 Ta = []
132 Pma = []
133 
134 tf.prepare (reion, cosmo)
135 
136 for i in range (0, np):
137   lnkh = log (1e-4) +  log (1e7) * divfac * i
138   kh = exp (lnkh)
139   T = tf.eval (cosmo, kh)
140   Pm = 1.0e3 / 7.0 * tf.matter_powerspectrum (cosmo, kh)
141   kha.append (kh)
142   Ta.append (T)
143   Pma.append (Pm)
144 
145 #
146 #  Ploting transfer and matter power spectrum
147 #
148 
149 plt.title ("Transfer Function and Matter Power Spectrum")
150 plt.xscale('log')
151 plt.plot (kha, Ta, 'r', label="T(kh)")
152 plt.plot (kha, Pma, 'b--', label="P_m(kh)")
153 plt.xlabel('$k$')
154 plt.legend(loc=1)
155 
156 plt.savefig ("transfer_func.png")
157 plt.clf ()
158 
159 #
160 # Calculating the variance filtered with the tophat windown function using
161 # scales R in the interval [5, 50] at redshift 0.3.
162 # First calculates the growth function at z = 0.3 and then the spectrum
163 # amplitude from the sigma8 parameter.
164 #
165 
166 vp.prepare (reion, cosmo)
167 
168 Dz = gf.eval (cosmo, 0.3)
169 A  = vp.sigma8_sqrtvar0 (cosmo)
170 prefact = A * A * Dz * Dz
171 
172 Ra = []
173 sigma2a = []
174 dlnsigma2a = []
175 
176 for i in range (0, np):
177   lnR = log (5.0) +  log (10.0) * divfac * i
178   R = exp (lnR)
179   sigma2 = prefact * vp.var0 (cosmo, lnR)
180   dlnsigma2 = vp.dlnvar0_dlnR (cosmo, lnR) 
181   Ra.append (R)
182   sigma2a.append (sigma2)
183   dlnsigma2a.append (dlnsigma2)
184 
185 #
186 #  Ploting filtered matter variance
187 #
188 
189 plt.title ("Variance and Variance Derivative")
190 plt.plot (Ra, sigma2a, 'r', label='$\sigma^2(\ln(R))$')
191 plt.plot (Ra, dlnsigma2a, 'b--', label='$d\ln(\sigma^2)/d\ln(R)$')
192 plt.xlabel('$R$')
193 plt.legend(loc=1)
194 
195 plt.savefig ("matter_var.png")
196 plt.clf ()
197 
198 #
199 # Calculating the mass function integrated in the mass interval [1e14, 1e16]
200 # for the redhshifts in the interval [0, 2.0] and area 200 squared degree.
201 #
202 
203 mf.set_area_sd (200.0)
204 mf.set_eval_limits (cosmo, log (1e14), log(1e16), 0.0, 2.0)
205 mf.prepare (reion, cosmo)
206 
207 dndza = []
208 
209 for i in range (0, np):
210   dndz = mf.dn_dz (cosmo, log(1e14), log(1e16), za[i], True)
211   dndza.append (dndz)
212 
213 #
214 #  Ploting the mass function
215 #
216 
217 plt.title ("Mass Function")
218 plt.plot (za, dndza, 'r', label='$dn/dz$')
219 plt.xlabel('$z$')
220 plt.legend(loc=1)
221 
222 plt.savefig ("mass_function.png")
223 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 #reion.z_to_tau (cosmo)
 74 
 75 #
 76 # Preparing the Class backend object
 77 #
 78 Bcbe.prepare (prim, reion, cosmo)
 79 
 80 Cls1 = Ncm.Vector.new (lmax + 1)
 81 Cls2 = Ncm.Vector.new (lmax + 1)
 82 
 83 Bcbe.get_TT_Cls (Cls1)
 84 
 85 prim.props.a = 0
 86 Bcbe.prepare (prim, reion, cosmo)
 87 Bcbe.get_TT_Cls (Cls2)
 88 
 89 Cls1_a = Cls1.dup_array ()
 90 Cls2_a = Cls2.dup_array ()
 91 
 92 Cls1_a = np.array (Cls1_a[2:])
 93 Cls2_a = np.array (Cls2_a[2:])
 94 
 95 ell = np.array (range (2, lmax + 1))
 96 
 97 Cls1_a = ell * (ell + 1.0) * Cls1_a
 98 Cls2_a = ell * (ell + 1.0) * Cls2_a
 99 
100 #
101 #  Ploting ionization history.
102 #
103 
104 plt.title (r'Modified and non-modified $C_\ell$')
105 plt.xscale('log')
106 plt.plot (ell, Cls1_a, 'r', label="Modified")
107 plt.plot (ell, Cls2_a, 'b--', label="Non-modified")
108 
109 plt.xlabel(r'$\ell$')
110 plt.ylabel(r'$C_\ell$')
111 plt.legend(loc=2)
112 
113 plt.savefig ("hiprim_Cls.png")
114 
115 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)

Contact

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

Valid XHTML 1.0 Strict

Valid CSS!