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

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