Getting Started¶
The rsdfit
Executable¶
The pyRSD package includes a rsdfit
executable that is capable of
running parameter estimation using the power spectrum models
GalaxySpectrum
and QuasarSpectrum
to fit power spectrum data input by the user.
After installing the pyRSD package, the rsdfit
command should be
located on your PATH
. We can then print the help message for the command
via
$ rsdfit -h
usage: rsdfit [-h] [--version] {mcmc,nlopt,restart,analyze} ...
From more help on each of the subcommands, type:
rsdfit mcmc -h
rsdfit nlopt -h
rsdfit restart -h
rsdfit analyze -h
fitting redshift space power spectrum observations with the `pyRSD` model
positional arguments:
{mcmc,nlopt,restart,analyze}
mcmc find the best-fitting parameters using emcee to run
MCMC chains
nlopt use nonlinear optimization to find the best-fitting
using the LBFGS algorithm
restart restart a parameter estimation fit from an existing
result
analyze analyze the results of an MCMC parameter fit
optional arguments:
-h, --help show this help message and exit
--version print the pyRSD version and exit
The rsdfit
command contains four sub-commands: nlopt, mcmc, restart,
and analyze. The most important of these sub-commands are nlopt
and mcmc
. This allows the user to run fresh parameter fits from scratch, which
is the main use case for the rsdfit
executable.
The calling sequence for the mcmc
command is
$ rsdfit mcmc -h
usage: rsdfit [-h] [--version] {mcmc,nlopt,restart,analyze} ...
From more help on each of the subcommands, type:
rsdfit mcmc -h
rsdfit nlopt -h
rsdfit restart -h
rsdfit analyze -h mcmc
[-h] [-m MODEL] [-p PARAMS] [--silent] -w WALKERS -i ITERATIONS
[-n NCHAINS] -o FOLDER [--debug] [--no-save-model]
optional arguments:
-h, --help show this help message and exit
-m MODEL, --model MODEL
file name holding the model path
-p PARAMS, --params PARAMS
file name holding the driver, theory, and data
parameters
--silent silence the standard output to the console
-w WALKERS number of emcee walkers to run the MCMC chain
(required)
-i ITERATIONS number of steps to run in the MCMC chain (required)
-n NCHAINS, --nchains NCHAINS
number of MCMC chains to run concurrently
-o FOLDER, --output FOLDER
the folder where the results will be written
(required)
--debug whether to print more info about the mpi4py.Pool
object
--no-save-model do not save the model instance
and for the nlopt
command is
$ rsdfit nlopt -h
usage: rsdfit [-h] [--version] {mcmc,nlopt,restart,analyze} ...
From more help on each of the subcommands, type:
rsdfit mcmc -h
rsdfit nlopt -h
rsdfit restart -h
rsdfit analyze -h nlopt
[-h] [-m MODEL] [-p PARAMS] [--silent] -i ITERATIONS -o FOLDER
[--debug] [--no-save-model]
optional arguments:
-h, --help show this help message and exit
-m MODEL, --model MODEL
file name holding the model path
-p PARAMS, --params PARAMS
file name holding the driver, theory, and data
parameters
--silent silence the standard output to the console
-i ITERATIONS the maximum number of iterations to run (required)
-o FOLDER, --output FOLDER
the folder where the results will be written
(required)
--debug whether to print more info about the mpi4py.Pool
object
--no-save-model do not save the model instance
The three most important options for these sub-commands are:
-p, —params
The name of the main parameter file to load; this holds all of the relevant information about the data, theory, and main driver parameters. We’ll discuss these parameter files in more detail in the next section.
-o, —output
This is the name of the directory where the results will be saved. Each time
rsdfit
is run, the parameter file is saved to a file “params.dat” in this directory as well. If this folder already exists and contains a “params.dat” file, those parameters will be used and results will be added to the existing directory.-m, —model
The name of a file holding a
GalaxySpectrum
orQuasarSpectrum
object to be loaded by the code. This argument is optional, but providing an already initialized model file will save a significant amount of time, since the code won’t need to initialize a model from scratch. See Model Initialization for more details on initializing and saving models.
Parameter Files¶
The rsdfit
command is configured through a single main parameter file.
This parameter file is responsible for not only configuring the main rsdfit
executable, but also for specifying the relevant theory parameters and
information about the data measurements to load. To separate the different
classes of parameters, the parameter names are prefixed with an additional
identifier, one of
- driver: general parameters that determine the
rsdfit
configuration - data: the parameters specifying the data and covariance matrix to load
- theory: the theoretical parameters for the desired pyRSD model, which determines the free pararameters, constrained parameters, etc, during the fitting procedure
- model: parameters specifying the
GalaxySpectrum
orQuasarSpectrum
configuration; these are parameters that are passed to the__init__()
function of these model classes
We’ll be exploring each of these parameter types in much more detail in the next few sections, but for now, let’s take a quick look at an example parameter file:
#-------------------------------------------------------------------------------
# driver_params
#-------------------------------------------------------------------------------
driver.burnin = 0
driver.epsilon = 0.02
driver.init_from = 'fiducial'
driver.init_scatter = 0.0
driver.lbfgs_epsilon = {'Nsat_mult': 0.01, 'f1h_cBs': 0.01}
driver.lbfgs_options = {'ftol': 1e-10, 'xtol': 1e-10, 'gtol': 1e-05}
driver.lbfgs_use_priors = True
driver.solver_type = 'nlopt'
driver.start_from = None
driver.test_convergence = False
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# data params
#-------------------------------------------------------------------------------
data.covariance = '$(PYRSD_DATA)/examples/runPB_poles_gaussian_cov.dat'
data.covariance_Nmocks = 0
data.covariance_rescaling = 1.0
data.data_file = '$(PYRSD_DATA)/examples/runPB_galaxy_poles.dat'
data.ells = [0, 2, 4]
data.fitting_range = [(0.02, 0.4), (0.02, 0.4), (0.02, 0.4)]
data.grid_file = '$(PYRSD_DATA)/examples/runPB_pkmu_grid.dat'
data.mode = 'poles'
data.mu_bounds = None
data.statistics = ['pole_0', 'pole_2', 'pole_4']
data.usedata = range(0, 3)
data.window_file = None
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# theory params
#-------------------------------------------------------------------------------
theory.F_AP = {'value': 1.0, 'vary': False, 'analytic': False, 'name': 'F_AP', 'expr': 'alpha_par/alpha_perp'}
theory.N = {'value': 0, 'vary': False, 'analytic': False, 'min': 0, 'name': 'N', 'lower': 0, 'prior': 'uniform', 'fiducial': 0, 'upper': 500}
theory.NcBs = {'value': 29608.657627046545, 'vary': False, 'analytic': False, 'name': 'NcBs', 'fiducial': 45000.0, 'expr': 'f1h_cBs / (fcB*(1 - fs)*nbar)'}
theory.NsBsB = {'value': 94583.2118641765, 'vary': False, 'analytic': False, 'name': 'NsBsB', 'fiducial': 94500.0, 'expr': 'f1h_sBsB / (fsB**2 * fs**2 * nbar) * (fcB*(1 - fs) - fs*(1-fsB))'}
theory.Nsat_mult = {'value': 2.4, 'vary': True, 'sigma': 0.2, 'min': 2.0, 'name': 'Nsat_mult', 'prior': 'normal', 'fiducial': 2.4, 'mu': 2.4, 'analytic': False}
theory.alpha = {'value': 1.0, 'vary': False, 'analytic': False, 'name': 'alpha', 'expr': '(alpha_perp**2 * alpha_par)**(1./3)'}
theory.alpha_par = {'value': 1.0, 'vary': True, 'analytic': False, 'name': 'alpha_par', 'prior': 'uniform', 'fiducial': 1.0, 'upper': 1.2, 'lower': 0.8}
theory.alpha_perp = {'value': 1.0, 'vary': True, 'analytic': False, 'name': 'alpha_perp', 'prior': 'uniform', 'fiducial': 1.0, 'upper': 1.2, 'lower': 0.8}
theory.b1 = {'value': 2.124276, 'vary': False, 'analytic': False, 'name': 'b1', 'expr': '(1 - fs)*b1_c + fs*b1_s'}
theory.b1_c = {'value': 1.9981383928571428, 'vary': False, 'analytic': False, 'name': 'b1_c', 'expr': '(1 - fcB)*b1_cA + fcB*b1_cB'}
theory.b1_cA = {'value': 1.9, 'vary': True, 'analytic': False, 'name': 'b1_cA', 'prior': 'uniform', 'fiducial': 1.9, 'upper': 2.5, 'lower': 1.2}
theory.b1_cB = {'value': 3.0028260869565218, 'vary': False, 'analytic': False, 'name': 'b1_cB', 'fiducial': 2.84, 'expr': '(1-fsB)/(1+fsB*(1./Nsat_mult - 1)) * b1_sA + (1 - (1-fsB)/(1+fsB*(1./Nsat_mult - 1))) * b1_sB'}
theory.b1_s = {'value': 3.2109999999999994, 'vary': False, 'analytic': False, 'name': 'b1_s', 'expr': '(1 - fsB)*b1_sA + fsB*b1_sB'}
theory.b1_sA = {'value': 2.755, 'vary': False, 'analytic': False, 'name': 'b1_sA', 'fiducial': 2.63, 'expr': 'gamma_b1sA*b1_cA'}
theory.b1_sB = {'value': 3.8949999999999996, 'vary': False, 'analytic': False, 'name': 'b1_sB', 'fiducial': 3.62, 'expr': 'gamma_b1sB*b1_cA'}
theory.b1sigma8 = {'value': 1.29580836, 'vary': False, 'analytic': False, 'name': 'b1sigma8', 'expr': 'b1*sigma8_z'}
theory.delta_sigsA = {'value': 1.0, 'vary': False, 'sigma': 0.2, 'min': 0.0, 'name': 'delta_sigsA', 'mu': 1.0, 'prior': 'normal', 'fiducial': 1.0, 'analytic': False}
theory.delta_sigsB = {'value': 1.0, 'vary': False, 'sigma': 0.2, 'min': 0.0, 'name': 'delta_sigsB', 'mu': 1.0, 'prior': 'normal', 'fiducial': 1.0, 'analytic': False}
theory.epsilon = {'value': 0.0, 'vary': False, 'analytic': False, 'name': 'epsilon', 'expr': '(alpha_perp/alpha_par)**(-1./3) - 1.0'}
theory.f = {'value': 0.78, 'vary': True, 'analytic': False, 'name': 'f', 'prior': 'uniform', 'fiducial': 0.78, 'upper': 1.0, 'lower': 0.6}
theory.f1h_cBs = {'value': 1.0, 'vary': True, 'sigma': 0.75, 'min': 0, 'name': 'f1h_cBs', 'mu': 1.0, 'prior': 'normal', 'fiducial': 1.0, 'analytic': False}
theory.f1h_sBsB = {'value': 4.0, 'vary': True, 'sigma': 1.0, 'min': 0.0, 'name': 'f1h_sBsB', 'prior': 'normal', 'fiducial': 4.0, 'mu': 4.0, 'analytic': False}
theory.f_so = {'value': 0.0, 'vary': False, 'sigma': 0.02, 'name': 'f_so', 'mu': 0.04, 'prior': 'normal', 'fiducial': 0.0, 'analytic': False}
theory.fcB = {'value': 0.08898809523809524, 'vary': False, 'analytic': False, 'min': 0, 'name': 'fcB', 'max': 1, 'fiducial': 0.089, 'expr': 'fs / (1 - fs) * (1 + fsB*(1./Nsat_mult - 1))'}
theory.fs = {'value': 0.104, 'vary': True, 'analytic': False, 'min': 0.0, 'name': 'fs', 'max': 1, 'prior': 'uniform', 'fiducial': 0.104, 'upper': 0.25, 'lower': 0.0}
theory.fsB = {'value': 0.4, 'vary': True, 'analytic': False, 'min': 0.0, 'name': 'fsB', 'max': 1, 'prior': 'uniform', 'fiducial': 0.4, 'upper': 1.0, 'lower': 0.0}
theory.fsigma8 = {'value': 0.4758, 'vary': False, 'analytic': False, 'name': 'fsigma8', 'expr': 'f*sigma8_z'}
theory.gamma_b1cB = {'value': 0.4, 'vary': False, 'sigma': 0.2, 'min': 0.0, 'name': 'gamma_b1cB', 'max': 1, 'prior': 'normal', 'fiducial': 0.4, 'mu': 0.4, 'analytic': False}
theory.gamma_b1sA = {'value': 1.45, 'vary': True, 'sigma': 0.3, 'min': 1.0, 'name': 'gamma_b1sA', 'prior': 'normal', 'fiducial': 1.45, 'mu': 1.45, 'analytic': False}
theory.gamma_b1sB = {'value': 2.05, 'vary': True, 'sigma': 0.3, 'min': 1.0, 'name': 'gamma_b1sB', 'prior': 'normal', 'fiducial': 2.05, 'mu': 2.05, 'analytic': False}
theory.nbar = {'value': 0.0004235857693396528, 'vary': False, 'analytic': False, 'name': 'nbar', 'fiducial': 0.0004235857693396528}
theory.sigma8_z = {'value': 0.61, 'vary': True, 'analytic': False, 'name': 'sigma8_z', 'prior': 'uniform', 'fiducial': 0.61, 'upper': 0.9, 'lower': 0.3}
theory.sigma_c = {'value': 1.0, 'vary': True, 'analytic': False, 'name': 'sigma_c', 'prior': 'uniform', 'fiducial': 1.0, 'upper': 3.0, 'lower': 0.0}
theory.sigma_sA = {'value': 3.5, 'vary': True, 'analytic': False, 'name': 'sigma_sA', 'prior': 'uniform', 'fiducial': 3.5, 'upper': 8.0, 'lower': 2.0}
theory.sigma_sB = {'value': 4.7072613620519554, 'vary': False, 'analytic': False, 'name': 'sigma_sB', 'prior': 'uniform', 'fiducial': 5, 'expr': 'sigma_sA * sigmav_from_bias(sigma8_z, b1_sB) / sigmav_from_bias(sigma8_z, b1_sA)', 'lower': 3.0, 'upper': 10.0}
theory.sigma_so = {'value': 0.0, 'vary': False, 'analytic': False, 'name': 'sigma_so', 'lower': 1.0, 'prior': 'uniform', 'fiducial': 0.0, 'upper': 7}
#-------------------------------------------------------------------------------
theory_extra.b2_00__0 = 'A0 for generic b2_00'
theory_extra.gamma_b1sA = 'the relative fraction of b1_sA to b1_cA'
theory_extra.delta_sigsA = 'the relative fraction of sigma_sA to sigma_s'
theory_extra.nbar = 'the mean number density in (h/Mpc)^3'
theory_extra.Nsat_mult = 'the mean number of satellites in halos with >1 sat'
theory_extra.gamma_b1cB = 'the relative fraction of b1_cB, varying linear between b1_sA and b1_s'
theory_extra.delta_sigsB = 'the relative fraction of sigma_sB to sigma_s'
theory_extra.f_nbar = 'fraction multiplying nbar'
theory_extra.f1h_sBsB = 'fraction multiplying 1-halo term, NsBsB'
theory_extra.b2_00__4 = 'A4 for generic b2_00'
theory_extra.f1h_cBs = 'fraction multiplying 1-halo term, NcBs'
theory_extra.gamma_b1sB = 'the relative fraction of b1_sB to b1_cA'
theory_extra.b2_00__2 = 'A2 for generic b2_00'
#-------------------------------------------------------------------------------
# model params
#-------------------------------------------------------------------------------
model.Pdv_model_type = 'jennings'
model.correct_mu2 = False
model.correct_mu4 = False
model.params = 'runPB.ini'
model.fog_model = 'modified_lorentzian'
model.include_2loop = False
model.interpolate = True
model.max_mu = 4
model.transfer_fit = 'CLASS'
model.use_P00_model = True
model.use_P01_model = True
model.use_P11_model = True
model.use_Pdv_model = True
model.use_Phm_model = True
model.use_mean_bias = False
model.use_so_correction = False
model.use_tidal_bias = False
model.use_vlah_biasing = True
model.vel_disp_from_sims = False
model.z = 0.55
#-------------------------------------------------------------------------------
The parameter file uses a simple key/value assignment syntax to define
the parameters. It is also possible to use environment variables when
defining parameters, was was done in the above example for $(PYRSD_DATA)
.
To learn more about each section of the parameter file, please see the next sections: