diff --git a/src/py21cmfast/inputs.py b/src/py21cmfast/inputs.py index 409bcaf29..4bb582131 100644 --- a/src/py21cmfast/inputs.py +++ b/src/py21cmfast/inputs.py @@ -290,31 +290,6 @@ class GlobalParams(StructInstanceWrapper): The peak gas temperatures behind the supersonic ionization fronts during reionization. VAVG: Avg value of the DM-b relative velocity [im km/s], ~0.9*SIGMAVCB (=25.86 km/s) normally. - STOC_MASS_TOL: - Mass tolerance for the stochastic halo sampling, a sample will be rejected if its mass sum falls - outside of the expected mass * (1 +- STOC_MASS_TOL) - SAMPLER_MIN_MASS: - Sets the minimum mass used by the halo sampler, halos below this mass will have the average contriubution - based on the integral of the given CMF. Greatly affects performance and memory usage. - MAXHALO_FACTOR: - Safety factor to multiply halo array sizes, total (all threads) array size will be UMF integral * MAXHALO_FACTOR - N_MASS_INTERP: - Number of mass bins for the sampler interpolation tables. - N_DELTA_INTERP: - Number of delta bins for the sampler interpolation tables. - N_PROB_INTERP: - Number of probability bins for the sampler interpolation tables. - MIN_LOGPROB: - Lower limit (in log probability) of the inverse CDF table. - SAMPLE_METHOD: - Sampling method to use for stochastic halos: - 0: Mass sampling from CMF - 1: N_halo sampling from CMF - 2: Sheth & Lemson 99 Fcoll sampling - 3: Darkforest (Qiu+20) binary splitting - AVG_BELOW_SAMPLER: - int flag to add the average halo proeprties in each cell between the source mass limit - and the mass sampled by the halo sampler """ def __init__(self, wrapped, ffi): @@ -510,11 +485,6 @@ class UserParams(StructWithDefaults): 0: GSL QAG adaptive integration, 1: Gauss-Legendre integration, previously forced in the interpolation tables, 2: Approximate integration, assuming sharp cutoffs and a triple power-law for sigma(M) based on EPS - INTEGRATION_METHOD_HALOS: int, optional - The integration method to use for all conditional MF integrals for halos from the sampler catalogues: - 0: GSL QAG adaptive integration, - 1: Gauss-Legendre integration, previously forced in the interpolation tables, - 2: Approximate integration, assuming sharp cutoffs and a triple power-law for sigma(M) based on EPS USE_2LPT: bool, optional Whether to use second-order Lagrangian perturbation theory (2LPT). Set this to True if the density field or the halo positions are extrapolated to @@ -532,6 +502,40 @@ class UserParams(StructWithDefaults): If STOC_MINIMUM_Z is not provided, we simply sample at the given redshift, unless USE_TS_FLUCT is given or we want a lightcone, in which case the minimum redshift is set to the minimum redshift of those fields. + SAMPLER_MIN_MASS: float, optional + The minimum mass to sample in the halo sampler when USE_HALO_FIELD and HALO_STOCHASTICITY are true, + decreasing this can drastically increase both compute time and memory usage. + SAMPLER_BUFFER_FACTOR: float, optional + The arrays for the halo sampler will have size of SAMPLER_BUFFER_FACTOR multiplied by the expected + number of halos in the box. Ideally this should be close to unity but one may wish to increase it to + test alternative scenarios + N_COND_INTERP: int, optional + The number of condition bins in the inverse CMF tables. + N_PROB_INTERP: int, optional + The number of probability bins in the inverse CMF tables. + MIN_LOGPROB: float, optional + The minimum log-probability of the inverse CMF tables. + SAMPLE_METHOD: int, optional + The sampling method to use in the halo sampler when calculating progenitor populations: + 0: Mass-limited CMF sampling, where samples are drawn until the expected mass is reached + 1: Number-limited CMF sampling, where we select a number of halos from the Poisson distribution + and then sample the CMF that many times + 2: Sheth et al 1999 Partition sampling, where the EPS collapsed fraction is sampled (gaussian tail) + and then the condition is updated using the conservation of mass. + 3: Parkinsson et al 2008 Binary split model as in DarkForest (Qiu et al 2021) where the EPS merger rate + is sampled on small internal timesteps such that only binary splits can occur. + NOTE: Sampling from the density grid will ALWAYS use number-limited sampling (method 1) + AVG_BELOW_SAMPLER: bool, optional + When switched on, an integral is performed in each cell between the minimum source mass and SAMPLER_MIN_MASS, + effectively placing the average halo population in each HaloBox cell below the sampler resolution. + When switched off, all halos below SAMPLER_MIN_MASS are ignored. This flag saves memory for larger boxes, + while still including the effects of smaller sources, albeit without stochasticity. + HALOMASS_CORRECTION: float, optional + This provides a corrective factor to the mass-limited (SAMPLE_METHOD==0) sampling, which multiplies the + expected mass from a condition by this number. The default value of 0.9 is calibrated to the mass-limited + sampling on a timestep of ZPRIME_STEP_FACTOR=1.02. + If ZPRIME_STEP_FACTOR is increased, this value should be set closer to 1. + This factor is also used in the partition (SAMPLE_METHOD==2) sampler, dividing nu(M) of each sample drawn. """ _ffi = ffi @@ -551,15 +555,25 @@ class UserParams(StructWithDefaults): "USE_INTERPOLATION_TABLES": None, "INTEGRATION_METHOD_ATOMIC": 1, "INTEGRATION_METHOD_MINI": 1, - "INTEGRATION_METHOD_HALOS": 0, "USE_2LPT": True, "MINIMIZE_MEMORY": False, "STOC_MINIMUM_Z": None, "KEEP_3D_VELOCITIES": False, + "SAMPLER_MIN_MASS": 5e7, + "SAMPLER_BUFFER_FACTOR": 2.0, + "MAXHALO_FACTOR": 2.0, + "N_COND_INTERP": 200, + "N_PROB_INTERP": 400, + "MIN_LOGPROB": -12, + "SAMPLE_METHOD": 0, + "AVG_BELOW_SAMPLER": False, + "HALOMASS_CORRECTION": 0.9, } _hmf_models = ["PS", "ST", "WATSON", "WATSON-Z", "DELOS"] _power_models = ["EH", "BBKS", "EFSTATHIOU", "PEEBLES", "WHITE", "CLASS"] + _sample_methods = ["MASS-LIMITED", "NUMBER-LIMITED", "PARTITION", "BINARY-SPLIT"] + _integral_methods = ["GSL-QAG", "GAUSS-LEGENDRE", "GAMMA-APPROX"] @property def USE_INTERPOLATION_TABLES(self): @@ -603,6 +617,38 @@ def HII_tot_num_pixels(self): """Total number of pixels in the low-res box.""" return self.NON_CUBIC_FACTOR * self.HII_DIM**3 + def _get_enum_property(self, prop, enum_list, propname=""): + """ + Retrieve a value for a property with defined enum list (see UserParams._power_models etc.). + + Arguments + --------- + prop: the hidden attribute to find in the enum list + enum_list: the list of parameter value strings + propname: the name of the property (for error messages) + + Returns + ------- + The index of prop within the parameter list, corresponding to it's integer value in + the C backend + """ + # if it's a string we grab the index of the list + if isinstance(prop, str): + val = enum_list.index(prop.upper()) + # otherwise it's a number so we leave it alone + else: + val = prop + + try: + val = int(val) + except (ValueError, TypeError) as e: + raise ValueError(f"{val} is an invalid value for {propname}") from e + + if not 0 <= val < len(enum_list): + raise ValueError(f"HMF must be an int between 0 and {len(enum_list) - 1}") + + return val + @property def POWER_SPECTRUM(self): """ @@ -622,17 +668,9 @@ def POWER_SPECTRUM(self): self._POWER_SPECTRUM = 5 return 5 else: - if isinstance(self._POWER_SPECTRUM, str): - val = self._power_models.index(self._POWER_SPECTRUM.upper()) - else: - val = self._POWER_SPECTRUM - - if not 0 <= val < len(self._power_models): - raise ValueError( - f"Power spectrum must be between 0 and {len(self._power_models) - 1}" - ) - - return val + return self._get_enum_property( + self._POWER_SPECTRUM, self._power_models, propname="POWER_SPECTRUM" + ) @property def HMF(self): @@ -640,22 +678,41 @@ def HMF(self): See hmf_model for a string representation. """ - if isinstance(self._HMF, str): - val = self._hmf_models.index(self._HMF.upper()) - else: - val = self._HMF + return self._get_enum_property(self._HMF, self._hmf_models, propname="HMF") - try: - val = int(val) - except (ValueError, TypeError) as e: - raise ValueError("Invalid value for HMF") from e + @property + def SAMPLE_METHOD(self): + """The SAMPLE_METHOD to use (an int, mapping to a given form). - if not 0 <= val < len(self._hmf_models): - raise ValueError( - f"HMF must be an int between 0 and {len(self._hmf_models) - 1}" - ) + See sample_method for a string representation. + """ + return self._get_enum_property( + self._SAMPLE_METHOD, self._sample_methods, propname="SAMPLE_METHOD" + ) - return val + @property + def INTEGRATION_METHOD_ATOMIC(self): + """The Integration method to use for atomic halos (an int, mapping to a given form). + + See integral_method_atomic for a string representation. + """ + return self._get_enum_property( + self._INTEGRATION_METHOD_ATOMIC, + self._integral_methods, + propname="INTEGRATION_METHOD_ATOMIC", + ) + + @property + def INTEGRATION_METHOD_MINI(self): + """The Integration method to use for atomic halos (an int, mapping to a given form). + + See integral_method_atomic for a string representation. + """ + return self._get_enum_property( + self._INTEGRATION_METHOD_MINI, + self._integral_methods, + propname="INTEGRATION_METHOD_MINI", + ) @property def hmf_model(self): @@ -668,13 +725,19 @@ def power_spectrum_model(self): return self._power_models[self.POWER_SPECTRUM] @property - def INTEGRATION_METHOD_HALOS(self): - """The integration methods other than QAG do not yet work for halos.""" - if self._INTEGRATION_METHOD_HALOS != 0: - warnings.warn( - "Only the QAG integrator currently works for the halo sampler, setting to 1 or 2 is for testing only" - ) - return self._INTEGRATION_METHOD_HALOS + def sample_method(self): + """String representation of the halo sampler method used.""" + return self._sample_methods[self.SAMPLE_METHOD] + + @property + def integration_method_atomic(self): + """String representation of the halo sampler method used.""" + return self._integral_methods[self.INTEGRATION_METHOD_ATOMIC] + + @property + def integration_method_mini(self): + """String representation of the halo sampler method used.""" + return self._intregal_methods[self.INTEGRATION_METHOD_MINI] @property def cell_size(self) -> un.Quantity[un.Mpc]: diff --git a/src/py21cmfast/src/21cmFAST.h b/src/py21cmfast/src/21cmFAST.h index f771a778e..e5eb55e90 100644 --- a/src/py21cmfast/src/21cmFAST.h +++ b/src/py21cmfast/src/21cmFAST.h @@ -31,10 +31,20 @@ struct UserParams{ bool USE_INTERPOLATION_TABLES; int INTEGRATION_METHOD_ATOMIC; int INTEGRATION_METHOD_MINI; - int INTEGRATION_METHOD_HALOS; bool USE_2LPT; bool MINIMIZE_MEMORY; bool KEEP_3D_VELOCITIES; + + //Halo Sampler Options + float SAMPLER_MIN_MASS; + double SAMPLER_BUFFER_FACTOR; + float MAXHALO_FACTOR; + int N_COND_INTERP; + int N_PROB_INTERP; + double MIN_LOGPROB; + int SAMPLE_METHOD; + bool AVG_BELOW_SAMPLER; + double HALOMASS_CORRECTION; }; struct AstroParams{ @@ -132,7 +142,7 @@ struct HaloBox{ float *n_ion; //weighted by F_ESC*PopN_ion float *halo_sfr; //for x-rays and Ts stuff float *halo_sfr_mini; //for x-rays and Ts stuff - float *whalo_sfr; + float *whalo_sfr; //SFR weighted by PopN_ion and F_ESC, used for Gamma12 double log10_Mcrit_LW_ave; }; @@ -301,6 +311,7 @@ void initialise_SFRD_Conditional_table(double min_density, double max_density, d float Fstar10, float Fstar7_MINI, int method, int method_mini, bool minihalos); void initialise_dNdM_tables(double xmin, double xmax, double ymin, double ymax, double growth1, double param, bool from_catalog); +void initialise_dNdM_inverse_table(double xmin, double xmax, double lnM_min, double growth1, double param, bool from_catalog); //evaluation of tables double EvaluateNionTs(double redshift, double Mlim_Fstar, double Mlim_Fesc); diff --git a/src/py21cmfast/src/FindHaloes.c b/src/py21cmfast/src/FindHaloes.c index 59ac48596..d6cf80873 100644 --- a/src/py21cmfast/src/FindHaloes.c +++ b/src/py21cmfast/src/FindHaloes.c @@ -106,11 +106,7 @@ LOG_DEBUG("redshift=%f", redshift); .HMF = user_params->HMF, }; - if(user_params->INTEGRATION_METHOD_HALOS == 1){ - initialise_GL(NGL_INT,log(M_MIN),log(Mmax_debug)); - } - - double nhalo_debug = IntegratedNdM(log(M_MIN), log(Mmax_debug), params, 1, user_params->INTEGRATION_METHOD_HALOS) * VOLUME; + double nhalo_debug = IntegratedNdM(log(M_MIN), log(Mmax_debug), params, 1, 0) * VOLUME; //expected halos above minimum filter mass LOG_DEBUG("DexM: We expect %.2f Halos between Masses [%.2e,%.2e]",nhalo_debug,M_MIN,Mmax_debug); } @@ -166,7 +162,7 @@ LOG_DEBUG("redshift=%f", redshift); //Pending a serious deep-dive into this algorithm, I will force DexM to use the fitted parameters to the // Sheth-Tormen mass function (as of right now, We do not even reproduce EPS results) - delta_crit = growth_factor*sheth_delc(Deltac/growth_factor, sigma_z0(M)); + delta_crit = growth_factor*sheth_delc_dexm(Deltac/growth_factor, sigma_z0(M)); // if(global_params.DELTA_CRIT_MODE == 1){ // //This algorithm does not use the sheth tormen OR Jenkins parameters, diff --git a/src/py21cmfast/src/Globals.h b/src/py21cmfast/src/Globals.h index 0508a6907..4517f4a12 100644 --- a/src/py21cmfast/src/Globals.h +++ b/src/py21cmfast/src/Globals.h @@ -81,17 +81,6 @@ struct GlobalParams{ float VAVG; - //Stochasticity Options - float STOC_MASS_TOL; - float SAMPLER_MIN_MASS; - float MAXHALO_FACTOR; - int N_MASS_INTERP; - int N_COND_INTERP; - int N_PROB_INTERP; - double MIN_LOGPROB; - int SAMPLE_METHOD; - int AVG_BELOW_SAMPLER; - bool USE_ADIABATIC_FLUCTUATIONS; }; @@ -153,10 +142,8 @@ extern struct GlobalParams global_params = { .OMtot = 1.0, .Y_He = 0.245, .wl = -1.0, - // .SHETH_b = 0.485, //In the literature this is 0.485 (RM08) or 0.5 (SMT01) Master 21cmFAST currently has 0.15 - // .SHETH_c = 0.615, //In the literature this is 0.615 (RM08) or 0.6 (SMT01) Master 21cmFAST currently has 0.05 - .SHETH_b = 0.15, - .SHETH_c = 0.05, + .SHETH_b = 0.15, //In the literature this is 0.485 (RM08) or 0.5 (SMT01) or 0.34 (Barkana+01) Master 21cmFAST currently has 0.15 + .SHETH_c = 0.05, //In the literature this is 0.615 (RM08) or 0.6 (SMT01) or 0.81 (Barkana+01) Master 21cmFAST currently has 0.05 .Zreion_HeII = 3.0, .FILTER = 0, .R_BUBBLE_MIN = 0.620350491, @@ -167,15 +154,5 @@ extern struct GlobalParams global_params = { .VAVG=25.86, - .STOC_MASS_TOL = 5.0, //effectively infinite, mass tolerance semi-deprecated - .SAMPLER_MIN_MASS = 5e7, - .MAXHALO_FACTOR = 2, - .N_MASS_INTERP = 200, - .N_COND_INTERP = 200, - .N_PROB_INTERP = 200, - .MIN_LOGPROB = -16, - .SAMPLE_METHOD = 0, - .AVG_BELOW_SAMPLER = 0, - .USE_ADIABATIC_FLUCTUATIONS = 1, }; diff --git a/src/py21cmfast/src/HaloBox.c b/src/py21cmfast/src/HaloBox.c index d6f9b6389..676f28efa 100644 --- a/src/py21cmfast/src/HaloBox.c +++ b/src/py21cmfast/src/HaloBox.c @@ -172,7 +172,6 @@ int set_fixed_grids(double redshift, double norm_esc, double alpha_esc, double M user_params_stoc->INTEGRATION_METHOD_MINI, flag_options_stoc->USE_MINI_HALOS, false); - //TODO: disable inverse table generation here with a flag or split up the functions initialise_dNdM_tables(min_density, max_density, lnMmin, lnMmax, growth_z, lnMcell, false); } @@ -438,8 +437,8 @@ int ComputeHaloBox(double redshift, struct UserParams *user_params, struct Cosmo } else{ //set below-resolution properties - if(global_params.AVG_BELOW_SAMPLER && M_min < global_params.SAMPLER_MIN_MASS){ - set_fixed_grids(redshift, norm_esc, alpha_esc, M_min, global_params.SAMPLER_MIN_MASS, ini_boxes, + if(user_params->AVG_BELOW_SAMPLER && M_min < user_params->SAMPLER_MIN_MASS){ + set_fixed_grids(redshift, norm_esc, alpha_esc, M_min, user_params->SAMPLER_MIN_MASS, ini_boxes, perturbed_field, previous_spin_temp, previous_ionize_box, grids, averages_box); //This is pretty redundant, but since the fixed grids have density units (X Mpc-3) I have to re-multiply before adding the halos. // I should instead have a flag to output the summed values in cell. (2*N_pixel > N_halo so generally i don't want to do it in the halo loop) @@ -450,7 +449,7 @@ int ComputeHaloBox(double redshift, struct UserParams *user_params, struct Cosmo grids->halo_sfr_mini[idx] *= cell_volume; grids->whalo_sfr[idx] *= cell_volume; } - LOG_DEBUG("finished subsampler M[%.2e %.2e]",M_min,global_params.SAMPLER_MIN_MASS); + LOG_DEBUG("finished subsampler M[%.2e %.2e]",M_min,user_params->SAMPLER_MIN_MASS); } #pragma omp parallel num_threads(user_params->N_THREADS) private(idx) { @@ -589,10 +588,10 @@ int ComputeHaloBox(double redshift, struct UserParams *user_params, struct Cosmo M_turn_a_avg /= total_n_halos; } - get_box_averages(redshift, norm_esc, alpha_esc, global_params.SAMPLER_MIN_MASS, global_params.M_MAX_INTEGRAL, M_turn_a_avg, M_turn_m_avg, averages_global); - get_box_averages(redshift, norm_esc, alpha_esc, M_min, global_params.SAMPLER_MIN_MASS, M_turn_a_avg, M_turn_m_avg, averages_subsampler); - get_box_averages(redshift, norm_esc, alpha_esc, global_params.SAMPLER_MIN_MASS, global_params.M_MAX_INTEGRAL, M_turn_a_nofb, M_turn_m_nofb_avg, averages_nofb); - get_box_averages(redshift, norm_esc, alpha_esc, M_min, global_params.SAMPLER_MIN_MASS, M_turn_a_nofb, M_turn_m_nofb_avg, averages_nofb_sub); + get_box_averages(redshift, norm_esc, alpha_esc, user_params->SAMPLER_MIN_MASS, global_params.M_MAX_INTEGRAL, M_turn_a_avg, M_turn_m_avg, averages_global); + get_box_averages(redshift, norm_esc, alpha_esc, M_min, user_params->SAMPLER_MIN_MASS, M_turn_a_avg, M_turn_m_avg, averages_subsampler); + get_box_averages(redshift, norm_esc, alpha_esc, user_params->SAMPLER_MIN_MASS, global_params.M_MAX_INTEGRAL, M_turn_a_nofb, M_turn_m_nofb_avg, averages_nofb); + get_box_averages(redshift, norm_esc, alpha_esc, M_min, user_params->SAMPLER_MIN_MASS, M_turn_a_nofb, M_turn_m_nofb_avg, averages_nofb_sub); } } @@ -609,7 +608,7 @@ int ComputeHaloBox(double redshift, struct UserParams *user_params, struct Cosmo LOG_DEBUG("Ratio: (HM %11.3e, NION %11.3e, SFR %11.3e, SFR_MINI %11.3e, WSFR %11.3e)",hm_avg/averages_global[0],nion_avg/averages_global[3], sfr_avg/averages_global[1],sfr_avg_mini/averages_global[2], wsfr_avg/averages_global[4]); - if(global_params.AVG_BELOW_SAMPLER && M_min < global_params.SAMPLER_MIN_MASS){ + if(user_params->AVG_BELOW_SAMPLER && M_min < user_params->SAMPLER_MIN_MASS){ LOG_DEBUG("SUB-SAMPLER",redshift); LOG_DEBUG("Exp. averages: (HM %11.3e, NION %11.3e, SFR %11.3e, SFR_MINI %11.3e, WSFR %11.3e)",averages_subsampler[0],averages_subsampler[3],averages_subsampler[1], averages_subsampler[2],averages_subsampler[4]); diff --git a/src/py21cmfast/src/SpinTemperatureBox.c b/src/py21cmfast/src/SpinTemperatureBox.c index 78f77f87f..fabdb1472 100644 --- a/src/py21cmfast/src/SpinTemperatureBox.c +++ b/src/py21cmfast/src/SpinTemperatureBox.c @@ -1258,7 +1258,7 @@ void ts_main(float redshift, float prev_redshift, struct UserParams *user_params zpp_for_evolve_list[R_ct],M_min_R[R_ct],M_max_R[R_ct],sigma_min[R_ct],sigma_max[R_ct]); } - //As far as I can tell, the only thing used from this is the X_e array + //Initialize heating interpolation arrays init_heat(); if(redshift >= global_params.Z_HEAT_MAX){ LOG_DEBUG("redshift greater than Z_HEAT_MAX"); diff --git a/src/py21cmfast/src/Stochasticity.c b/src/py21cmfast/src/Stochasticity.c index 32c8405e8..e6d11f4a1 100644 --- a/src/py21cmfast/src/Stochasticity.c +++ b/src/py21cmfast/src/Stochasticity.c @@ -2,13 +2,6 @@ * i.e sampling the halo mass function and * other halo relations.*/ -//TODO: Don't have every error be a ValueError - -//max number of attempts for mass tolerance before failure -#define MAX_ITERATIONS 1e2 -#define MAX_ITER_N 1e2 //for stoc_halo_sample (select N halos) how many tries for one N, this should be large to enforce a near-possion p(N) -#define MMAX_TABLES 1e14 - //buffer size (per cell of arbitrary size) in the sampling function #define MAX_HALO_CELL (int)1e5 @@ -61,22 +54,6 @@ void print_hs_consts(struct HaloSamplingConstants * c){ return; } -//The sigma interp table is regular in log mass, not sigma so we need to loop ONLY FOR METHOD=3 -double EvaluateSigmaInverse(double sigma, struct RGTable1D_f *s_table){ - int idx; - for(idx=0;idxy_arr[idx]) break; - } - if(idx == NMass){ - LOG_ERROR("sigma inverse out of bounds."); - Throw(TableEvaluationError); - } - double table_val_0 = s_table->x_min + idx*s_table->x_width; - double table_val_1 = s_table->x_min + (idx-1)*s_table->x_width; - double interp_point = (sigma - table_val_0)/(table_val_1-table_val_0); - - return table_val_0*(1-interp_point) + table_val_1*(interp_point); -} void Broadcast_struct_global_STOC(struct UserParams *user_params, struct CosmoParams *cosmo_params,struct AstroParams *astro_params, struct FlagOptions *flag_options){ cosmo_params_stoc = cosmo_params; @@ -90,7 +67,7 @@ double expected_nhalo(double redshift, struct UserParams *user_params, struct Co //minimum sampled mass Broadcast_struct_global_UF(user_params,cosmo_params); Broadcast_struct_global_PS(user_params,cosmo_params); - double M_min = global_params.SAMPLER_MIN_MASS; + double M_min = user_params->SAMPLER_MIN_MASS; //maximum sampled mass double M_max = RHOcrit * cosmo_params->OMm * VOLUME / HII_TOT_NUM_PIXELS; double growthf = dicke(redshift); @@ -98,7 +75,7 @@ double expected_nhalo(double redshift, struct UserParams *user_params, struct Co init_ps(); if(user_params->USE_INTERPOLATION_TABLES) - initialiseSigmaMInterpTable(M_min/2,M_max); + initialiseSigmaMInterpTable(M_min,M_max); struct parameters_gsl_MF_integrals params = { .redshift = redshift, @@ -106,10 +83,7 @@ double expected_nhalo(double redshift, struct UserParams *user_params, struct Co .HMF = user_params->HMF, }; - if(user_params->INTEGRATION_METHOD_HALOS == 1) - initialise_GL(NGL_INT, log(M_min), log(M_max)); - - result = IntegratedNdM(log(M_min), log(M_max), params, 1, user_params->INTEGRATION_METHOD_HALOS) * VOLUME; + result = IntegratedNdM(log(M_min), log(M_max), params, 1, 0) * VOLUME; LOG_DEBUG("Expected %.2e Halos in the box from masses %.2e to %.2e at z=%.2f",result,M_min,M_max,redshift); if(user_params->USE_INTERPOLATION_TABLES) @@ -118,16 +92,13 @@ double expected_nhalo(double redshift, struct UserParams *user_params, struct Co return result; } -//NOTE: if p(x) is uniform, p(log(1-x)) follows the exponential distribution -// But the gsl_ran_exponential function does the exact same thing but adds a mean double sample_dndM_inverse(double condition, struct HaloSamplingConstants * hs_constants, gsl_rng * rng){ - double p_in, min_prob; + double p_in, min_prob, result; p_in = gsl_rng_uniform(rng); - if(!hs_constants->from_catalog){ - p_in = log(p_in); - if(p_in < global_params.MIN_LOGPROB) p_in = global_params.MIN_LOGPROB; - } - return EvaluateRGTable2D(condition,p_in,&Nhalo_inv_table); + result = EvaluateNhaloInv(condition,p_in); + result = fmin(1,fmax(0,result)); //clip in case of extrapolation + result = result * hs_constants->M_cond; + return result; } //Set the constants that are calculated once per snapshot @@ -137,15 +108,17 @@ void stoc_set_consts_z(struct HaloSamplingConstants *const_struct, double redshi const_struct->z_out = redshift; const_struct->z_in = redshift_desc; - const_struct->M_min = global_params.SAMPLER_MIN_MASS; + const_struct->M_min = user_params_stoc->SAMPLER_MIN_MASS / user_params_stoc->SAMPLER_BUFFER_FACTOR; const_struct->lnM_min = log(const_struct->M_min); const_struct->M_max_tables = global_params.M_MAX_INTEGRAL; const_struct->lnM_max_tb = log(const_struct->M_max_tables); - init_ps(); if(user_params_stoc->USE_INTERPOLATION_TABLES){ - initialiseSigmaMInterpTable(const_struct->M_min / 2,const_struct->M_max_tables); + if(user_params_stoc->SAMPLE_METHOD == 3) + initialiseSigmaMInterpTable(const_struct->M_min/2,const_struct->M_max_tables); //the binary split needs to go below the resolution + else + initialiseSigmaMInterpTable(const_struct->M_min,const_struct->M_max_tables); } const_struct->sigma_min = EvaluateSigma(const_struct->lnM_min); @@ -162,8 +135,15 @@ void stoc_set_consts_z(struct HaloSamplingConstants *const_struct, double redshi const_struct->corr_star = 0; const_struct->from_catalog = 1; - initialise_dNdM_tables(const_struct->lnM_min, const_struct->lnM_max_tb,const_struct->lnM_min, const_struct->lnM_max_tb, + initialise_dNdM_tables(log(user_params_stoc->SAMPLER_MIN_MASS), const_struct->lnM_max_tb, const_struct->lnM_min, const_struct->lnM_max_tb, const_struct->growth_out, const_struct->growth_in, true); + if(user_params_stoc->SAMPLE_METHOD == 0 || user_params_stoc->SAMPLE_METHOD == 1){ + initialise_dNdM_inverse_table(log(user_params_stoc->SAMPLER_MIN_MASS), const_struct->lnM_max_tb, const_struct->lnM_min, + const_struct->growth_out, const_struct->growth_in, true); + } + if(user_params_stoc->SAMPLE_METHOD == 3){ + initialise_J_split_table(200,1e-4,20.,0.2); + } } else { double M_cond = RHOcrit * cosmo_params_stoc->OMm * VOLUME / HII_TOT_NUM_PIXELS; @@ -173,7 +153,10 @@ void stoc_set_consts_z(struct HaloSamplingConstants *const_struct, double redshi //for the table limits double delta_crit = get_delta_crit(user_params_stoc->HMF,const_struct->sigma_cond,const_struct->growth_out); const_struct->from_catalog = 0; - initialise_dNdM_tables(DELTA_MIN, MAX_DELTAC_FRAC*delta_crit, const_struct->lnM_min, const_struct->lnM_max_tb, const_struct->growth_out, const_struct->lnM_cond, false); + initialise_dNdM_tables(DELTA_MIN, MAX_DELTAC_FRAC*delta_crit, const_struct->lnM_min, const_struct->lnM_max_tb, + const_struct->growth_out, const_struct->lnM_cond, false); + initialise_dNdM_inverse_table(DELTA_MIN, MAX_DELTAC_FRAC*delta_crit, const_struct->lnM_min, + const_struct->growth_out, const_struct->lnM_cond, false); } LOG_DEBUG("Done."); @@ -198,32 +181,31 @@ void stoc_set_consts_cond(struct HaloSamplingConstants *const_struct, double con } //Here the condition is a cell of a given density, the volume/mass is given by the grid parameters else{ + //since the condition mass/sigma is already set all we need is delta const_struct->delta = cond_val; const_struct->cond_val = cond_val; - - //the splines don't work well for cells above Deltac, but there CAN be cells above deltac, since this calculation happens - //before the overlap, and since the smallest dexm mass is M_cell*(1.01^3) there *could* be a cell above Deltac not in a halo - //NOTE: all this does is prevent integration errors below since these cases are also dealt with in stoc_sample - if(cond_val > MAX_DELTAC_FRAC*get_delta_crit(user_params_stoc->HMF,const_struct->sigma_cond,const_struct->growth_out)){ - const_struct->expected_M = const_struct->M_cond; - const_struct->expected_N = 1; - return; - } - if(cond_val <= DELTA_MIN){ - const_struct->expected_M = 0; - const_struct->expected_N = 0; - return; - } } //Get expected N and M from interptables - // double EvaluateNhalo(double condition, double growthf, double lnMmin, double lnMmax, double sigma, double delta) - n_exp = EvaluateNhalo(const_struct->cond_val,const_struct->growth_out,const_struct->lnM_min - ,const_struct->lnM_max_tb,const_struct->M_cond,const_struct->sigma_cond,const_struct->delta); - m_exp = EvaluateMcoll(const_struct->cond_val,const_struct->growth_out,const_struct->lnM_min - ,const_struct->lnM_max_tb,const_struct->M_cond,const_struct->sigma_cond,const_struct->delta); - const_struct->expected_N = n_exp * const_struct->M_cond; - const_struct->expected_M = m_exp * const_struct->M_cond; + //the splines don't work well for cells above Deltac, but there CAN be cells above deltac, since this calculation happens + //before the overlap, and since the smallest dexm mass is M_cell*(1.01^3) there *could* be a cell above Deltac not in a halo + //NOTE: all this does is prevent integration errors below since these cases are also dealt with in stoc_sample + if(const_struct->delta > MAX_DELTAC_FRAC*get_delta_crit(user_params_stoc->HMF,const_struct->sigma_cond,const_struct->growth_out)){ + const_struct->expected_M = const_struct->M_cond; + const_struct->expected_N = 1; + } + else if(const_struct->delta <= DELTA_MIN){ + const_struct->expected_M = 0; + const_struct->expected_N = 0; + } + else{ + n_exp = EvaluateNhalo(const_struct->cond_val,const_struct->growth_out,const_struct->lnM_min, + const_struct->lnM_max_tb,const_struct->M_cond,const_struct->sigma_cond,const_struct->delta); + m_exp = EvaluateMcoll(const_struct->cond_val,const_struct->growth_out,const_struct->lnM_min, + const_struct->lnM_max_tb,const_struct->M_cond,const_struct->sigma_cond,const_struct->delta); + const_struct->expected_N = n_exp * const_struct->M_cond; + const_struct->expected_M = m_exp * const_struct->M_cond; + } return; } @@ -302,53 +284,20 @@ int add_properties_cat(struct UserParams *user_params, struct CosmoParams *cosmo /* Creates a realisation of halo properties by sampling the halo mass function and * conditional property PDFs, the number of halos is poisson sampled from the integrated CMF*/ int stoc_halo_sample(struct HaloSamplingConstants *hs_constants, gsl_rng * rng, int *n_halo_out, float *M_out){ - double M_min = hs_constants->M_min; - double mass_tol = global_params.STOC_MASS_TOL; double exp_N = hs_constants->expected_N; - double exp_M = hs_constants->expected_M; - double M_cond = hs_constants->M_cond; - double hm_sample, M_prog; + double hm_sample; int ii, nh; - int n_attempts=0, n_failures=0; int halo_count=0; double tbl_arg = hs_constants->cond_val; - for(n_failures=0;n_failures exp_M*(1+mass_tol))){ - nh = gsl_ran_poisson(rng,exp_N); - } - for(n_attempts=0;n_attempts exp_M*(1-mass_tol))){ - //using goto to break double loop - goto found_halo_sample; - } - } + nh = gsl_ran_poisson(rng,exp_N); + for(ii=0;ii= MAX_ITERATIONS){ - LOG_ERROR("passed max iter in sample, last attempt M=%.3e [%.3e, %.3e] Me %.3e Mt %.3e cond %.3e",tbl_arg,M_prog,exp_M*(1-mass_tol),exp_M*(1+mass_tol),exp_M); - Throw(ValueError); - } - - found_halo_sample: *n_halo_out = halo_count; - // LOG_ULTRA_DEBUG("Got %d (exp. %.2e) halos mass %.2e (exp. %.2e) %.2f | (%d,%d) att.", - // nh,exp_N,M_prog,exp_M,M_prog/exp_M - 1, n_failures, n_attempts); + *n_halo_out = halo_count; return 0; } @@ -367,18 +316,19 @@ double remove_random_halo(gsl_rng * rng, int n_halo, int *idx, double *M_prog, f return last_M_del; } -/*Convenience function to store all of the mass-based sample "corrections" i.e what we decide to do when the mass samples go above the trigger*/ +/*Function which "corrects" a mass sample after it exceeds the expected mass*/ //CURRENT IMPLEMENTATION: half the time I keep/throw away the last halo based on which sample is closer to the expected mass. // However this introduces a bias since the last halo is likely larger than average So the other half the time, -// I throw away random halos until we are again below exp_M, effectively the same process in reverse. which has the (exact?) opposite bias -int fix_mass_sample(gsl_rng * rng, double exp_M, int *n_halo_pt, double *M_tot_pt, float *M_out){ +// I throw away random halos until we are again below exp_M, effectively the same process in reverse. which has the opposite bias +void fix_mass_sample(gsl_rng * rng, double exp_M, int *n_halo_pt, double *M_tot_pt, float *M_out){ //Keep the last halo if it brings us closer to the expected mass //This is done by addition or subtraction over the limit to balance //the bias of the last halo being larger int random_idx; - int n_removed = 0; double last_M_del; - if(gsl_rng_uniform_int(rng,2)){ + bool sel = gsl_rng_uniform_int(rng,2); + // int sel = 1; + if(sel){ //LOG_ULTRA_DEBUG("Deciding to keep last halo M %.3e tot %.3e exp %.3e",M_out[*n_halo_pt-1],*M_tot_pt,exp_M); if(fabs(*M_tot_pt - M_out[*n_halo_pt-1] - exp_M) < fabs(*M_tot_pt - exp_M)){ //LOG_ULTRA_DEBUG("removed"); @@ -392,7 +342,6 @@ int fix_mass_sample(gsl_rng * rng, double exp_M, int *n_halo_pt, double *M_tot_p //here we remove by setting halo mass to zero, skipping it during the consolidation last_M_del = remove_random_halo(rng,*n_halo_pt,&random_idx,M_tot_pt,M_out); //LOG_ULTRA_DEBUG("Removed halo %d M %.3e tot %.3e",random_idx,last_M_del,*M_tot_pt); - n_removed++; } // if the sample with the last subtracted halo is closer to the expected mass, keep it @@ -401,121 +350,92 @@ int fix_mass_sample(gsl_rng * rng, double exp_M, int *n_halo_pt, double *M_tot_p M_out[random_idx] = last_M_del; // LOG_ULTRA_DEBUG("kept."); *M_tot_pt += last_M_del; - n_removed--; } } - return n_removed; - //Possible future improvements - // Try different target / trigger, which could overcome the M > exp_M issues better than a delta - // e.g: set trigger for the check at exp_M - (M_Max - exp_M), target at exp_M - // This provides a range of trigger points which can be absorbed into the below-resolution values - // Obviously breaks for exp_M < 0.5 M, think of other thresholds. - } /* Creates a realisation of halo properties by sampling the halo mass function and * conditional property PDFs, Sampling is done until there is no more mass in the condition - * Stochasticity is ignored below a certain mass threshold*/ + * Stochasticity is ignored below a certain mass threshold */ int stoc_mass_sample(struct HaloSamplingConstants * hs_constants, gsl_rng * rng, int *n_halo_out, float *M_out){ - //lnMmin only used for sampling, apply factor here - double mass_tol = global_params.STOC_MASS_TOL; double exp_M = hs_constants->expected_M; - if(hs_constants->from_catalog)exp_M *= 1.; //fudge factor for assuming that internal lagrangian volumes are independent - int n_halo_sampled, n_failures=0; + //The mass-limited sampling as-is has a slight bias to producing too many halos, + // which is independent of density or halo mass, + // this factor reduces the total expected mass to bring it into line with the CMF + exp_M *= user_params_stoc->HALOMASS_CORRECTION; + + int n_halo_sampled=0; double M_prog=0; double M_sample; - int n_removed; double tbl_arg = hs_constants->cond_val; - for(n_failures=0;n_failures= exp_M*(1-mass_tol))){ - //using goto to break double loop - goto found_halo_sample; - } - } - - //technically I don't need the if statement but it might be confusing otherwise - if(n_failures >= MAX_ITERATIONS){ - LOG_ERROR("passed max iter in sample, last attempt M=%.3e [%.3e, %.3e] cond %.3e",M_prog,exp_M*(1-mass_tol),exp_M*(1+mass_tol),tbl_arg); - Throw(ValueError); + M_prog += M_sample; + M_out[n_halo_sampled++] = M_sample; + // LOG_ULTRA_DEBUG("Sampled %.3e | %.3e %d",M_sample,M_prog,n_halo_sampled); } + // LOG_ULTRA_DEBUG("Before fix: %d %.3e",n_halo_sampled,M_prog); + //The above sample is above the expected mass, by up to 100%. I wish to make the average mass equal to exp_M + fix_mass_sample(rng,exp_M,&n_halo_sampled,&M_prog,M_out); + // LOG_ULTRA_DEBUG("After fix: %d %.3e",n_halo_sampled,M_prog); - found_halo_sample: *n_halo_out = n_halo_sampled; - // LOG_ULTRA_DEBUG("Got %d (exp.%.2e) halos mass %.2e (exp. %.2e) %.2f | %d att", - // n_halo_sampled,hs_constants->expected_N,M_prog,exp_M,M_prog/exp_M - 1,n_failures); + // LOG_ULTRA_DEBUG("Got %d (exp.%.2e) halos mass %.2e (exp. %.2e) %.2f", + // n_halo_sampled,hs_constants->expected_N,M_prog,exp_M,M_prog/exp_M - 1); + *n_halo_out = n_halo_sampled; return 0; } -//Sheth & Lemson partition model, Changes delta,M in the CMF as you sample +//Sheth & Lemson 1999 partition model, Changes delta,M in the CMF as you sample //This avoids discretization issues in the mass sampling, however..... //it has been noted to overproduce small halos (McQuinn+ 2007) //I don't know why sampling from Fcoll(M) is correct? // Do we not sample the same `particle` multiple times? (i.e 10x more samples for a single 10x mass halo) // How does the reduction of mass after each sample *exactly* cancel this 1/M effect //If you want a non-barrier-based CMF, I don't know how to implement it here -int stoc_sheth_sample(struct HaloSamplingConstants * hs_constants, gsl_rng * rng, int *n_halo_out, float *M_out){ +int stoc_partition_sample(struct HaloSamplingConstants * hs_constants, gsl_rng * rng, int *n_halo_out, float *M_out){ //lnMmin only used for sampling, apply factor here double exp_M = hs_constants->expected_M; double M_cond = hs_constants->M_cond; double d_cond = hs_constants->delta; double growthf = hs_constants->growth_out; - double M_min = hs_constants->M_min; double sigma_min = hs_constants->sigma_min; - double lnM_min = hs_constants->lnM_min; - double lnM_max_tb = hs_constants->lnM_max_tb; int n_halo_sampled; - double x_sample, sigma_sample, M_sample, M_remaining, delta_current; + double nu_sample, sigma_sample, M_sample, M_remaining, delta_current; double lnM_remaining, sigma_r, del_term; double tbl_arg = hs_constants->cond_val; n_halo_sampled = 0; - // LOG_ULTRA_DEBUG("Start: M %.2e (%.2e) d %.2e",M_cond,exp_M,d_cond); + LOG_ULTRA_DEBUG("Start: M %.2e (%.2e) d %.2e",M_cond,exp_M,d_cond); + double sigma_fudge_factor = user_params_stoc->HALOMASS_CORRECTION; //set initial amount (subtracted unresolved Mass) - M_remaining = exp_M; + M_remaining = M_cond; lnM_remaining = log(M_remaining); - double x_min; - - while(M_remaining > global_params.SAMPLER_MIN_MASS){ + double nu_min; + while(M_remaining > user_params_stoc->SAMPLER_MIN_MASS){ sigma_r = EvaluateSigma(lnM_remaining); + delta_current = (get_delta_crit(user_params_stoc->HMF,sigma_r,growthf) - d_cond)/(M_remaining/M_cond); del_term = delta_current*delta_current/growthf/growthf; - //Low x --> high sigma --> low mass, high x --> low sigma --> high mass - //|x| required to sample the smallest halo on the tables - x_min = sqrt(del_term/(sigma_min*sigma_min - sigma_r*sigma_r)); + nu_min = sqrt(del_term/(sigma_min*sigma_min - sigma_r*sigma_r)); //nu at minimum progenitor - //LOG_ULTRA_DEBUG("M_rem %.2e d %.2e sigma %.2e min %.2e xmin %.2f",M_remaining,delta_current,sigma_r,sigma_min,x_min); + LOG_ULTRA_DEBUG("M_rem %.2e d %.2e sigma %.2e min %.2e numin %.2e",M_remaining,delta_current,sigma_r,sigma_min,nu_min); //we use the gaussian tail distribution to enforce our Mmin limit from the sigma tables - x_sample = gsl_ran_ugaussian_tail(rng,x_min); - sigma_sample = sqrt(del_term/(x_sample*x_sample) + sigma_r*sigma_r); - M_sample = EvaluateSigmaInverse(sigma_sample,&Sigma_InterpTable); + nu_sample = gsl_ran_ugaussian_tail(rng,nu_min)/sigma_fudge_factor; + sigma_sample = sqrt(del_term/(nu_sample*nu_sample) + sigma_r*sigma_r); + // sigma_sample *= sigma_fudge_factor; + M_sample = EvaluateSigmaInverse(sigma_sample); M_sample = exp(M_sample); - //LOG_ULTRA_DEBUG("found Mass %d %.2e",n_halo_sampled,M_sample); + LOG_ULTRA_DEBUG("found Mass %d %.2e sigma %.2e nu %.2e dt %.2e",n_halo_sampled,M_sample,sigma_sample,nu_sample,del_term); M_out[n_halo_sampled++] = M_sample; M_remaining -= M_sample; lnM_remaining = log(M_remaining); @@ -525,14 +445,23 @@ int stoc_sheth_sample(struct HaloSamplingConstants * hs_constants, gsl_rng * rng return 0; } +double ComputeFraction_split( + double sigma_start, double sigmasq_start, double sigmasq_res, + double G1, double dd, double gamma1 +) { + double u_res = sigma_start*pow(sigmasq_res - sigmasq_start, -.5); + // LOG_ULTRA_DEBUG("Frac: u_res = %.2e, dd=%.2e, J=%.2e",u_res,dd,EvaluateJ(u_res,gamma1)); + return sqrt(2./PI)*EvaluateJ(u_res,gamma1)*G1/sigma_start*dd; +} + //binary splitting with small internal steps based on Parkinson+08, Bensen+16, Qiu+20 (Darkforest) //This code was mostly taken from Darkforest (Qiu+20) //NOTE: some unused variables here //Only works with adjusted EPS int stoc_split_sample(struct HaloSamplingConstants * hs_constants, gsl_rng * rng, int *n_halo_out, float *M_out){ double G0 = 1; - double gamma1 = 0.2; - double gamma2 = -0.4; + double gamma1 = 0.4; + double gamma2 = -0.2; double m_res = hs_constants->M_min; double lnm_res = hs_constants->lnM_min; double eps1 = 0.1; @@ -554,25 +483,20 @@ int stoc_split_sample(struct HaloSamplingConstants * hs_constants, gsl_rng * rng double m_prog1, m_prog2; int save; int n_out = 0; - int idx_first = -1; int idx = 0; - double growthf = hs_constants->growth_out; - //Can Accelerate (i.e use the J function to speed up) - double growth_d = 0; float d_points[MAX_HALO_CELL], m_points[MAX_HALO_CELL]; int n_points; //set initial points - d_points[0] = hs_constants->delta / growthf; + d_points[0] = Deltac / hs_constants->growth_in; m_points[0] = hs_constants->M_cond; - d_target = Deltac / growthf; + d_target = Deltac / hs_constants->growth_out; n_points = 1; + double M_total = 0; sigma_res = EvaluateSigma(lnm_res); sigmasq_res = sigma_res*sigma_res; - // LOG_DEBUG("Starting split %.2e %.2e",d_points[0],m_points[0]); - while(idx < n_points) { d_start = d_points[idx]; m_start = m_points[idx]; @@ -600,8 +524,7 @@ int stoc_split_sample(struct HaloSamplingConstants * hs_constants, gsl_rng * rng dd = dd_target; save = 1; } - growth_d = Deltac/(d_start + dd); - F = 1 - FgtrM_bias_fast(growth_d,d_start,sigma_res,sigma_start); + F = ComputeFraction_split(sigma_start, sigmasq_start, sigmasq_res, G1, dd, gamma1); } else { // Compute B and beta @@ -629,8 +552,7 @@ int stoc_split_sample(struct HaloSamplingConstants * hs_constants, gsl_rng * rng } N_upper = dN_dd*dd; // Compute F - growth_d = Deltac/(d_start + dd); - F = 1 - FgtrM_bias_fast(growth_d,d_start,sigma_res,sigma_start); + F = ComputeFraction_split(sigma_start, sigmasq_start, sigmasq_res, G1, dd, gamma1); // Generate random numbers and split the tree if (gsl_rng_uniform(rng) < N_upper) { q = pow(pow(q_res, eta) + pow_diff*gsl_rng_uniform(rng), 1./eta); @@ -647,9 +569,9 @@ int stoc_split_sample(struct HaloSamplingConstants * hs_constants, gsl_rng * rng q = 0.; // No split } } - // LOG_DEBUG("split i %d n %d m %.2e d %.2e",idx,n_points,m_start,d_start); - // LOG_DEBUG("q %.2e dd %.2e (%.2e %.2e) of %.2e",q,dd,dd1,dd2,dd_target); - // LOG_DEBUG("dNdd %.2e B %.2e pow %.2e eta %.2e ah %.2e G2 %.2e b %.2e",dN_dd,B,pow_diff,eta,alpha_half,G2,beta); + // LOG_ULTRA_DEBUG("split i %d n %d m %.2e d %.2e",idx,n_points,m_start,d_start); + // LOG_ULTRA_DEBUG("q %.2e F %.2e dd %.2e (%.2e %.2e) of %.2e",q,F,dd,dd1,dd2,dd_target); + // LOG_ULTRA_DEBUG("dNdd %.2e B %.2e pow %.2e eta %.2e ah %.2e G2 %.2e b %.2e",dN_dd,B,pow_diff,eta,alpha_half,G2,beta); // Compute progenitor mass m_prog1 = (1 - F - q)*m_start; m_prog2 = q*m_start; @@ -657,9 +579,11 @@ int stoc_split_sample(struct HaloSamplingConstants * hs_constants, gsl_rng * rng if (save) { if (m_prog1 > m_res) { M_out[n_out++] = m_prog1; + M_total += m_prog1; } if (m_prog2 > m_res) { M_out[n_out++] = m_prog2; + M_total += m_prog2; } } //if not finished yet, add them to the internal arrays @@ -670,11 +594,14 @@ int stoc_split_sample(struct HaloSamplingConstants * hs_constants, gsl_rng * rng // keep the active halo at zero until the end, but that's minor as it should only drift a few dozen else { if (m_prog1 > m_res){ + //replace current halo with the larger progenitor d_points[idx] = dd + d_start; m_points[idx] = m_prog1; + //since we replaced, do not advance the index idx--; } if (m_prog2 > m_res){ + //add the smaller progenitor to the end d_points[n_points] = dd + d_start; m_points[n_points++] = m_prog2; } @@ -682,6 +609,7 @@ int stoc_split_sample(struct HaloSamplingConstants * hs_constants, gsl_rng * rng idx++; } *n_halo_out = n_out; + // LOG_ULTRA_DEBUG("Total M = %.4e (%.4e) N = %d",M_total,M_total/hs_constants->M_cond,n_out); return 0; } @@ -696,12 +624,12 @@ int stoc_sample(struct HaloSamplingConstants * hs_constants, gsl_rng * rng, int int err; //If the expected mass is below our minimum saved mass, don't bother calculating //NOTE: some of these conditions are redundant with set_consts_cond() - if(hs_constants->delta <= DELTA_MIN || hs_constants->expected_M < global_params.SAMPLER_MIN_MASS){ + if(hs_constants->delta <= DELTA_MIN || hs_constants->expected_M < user_params_stoc->SAMPLER_MIN_MASS){ *n_halo_out = 0; return 0; } //if delta is above critical, form one big halo - if(hs_constants->delta > MAX_DELTAC_FRAC*get_delta_crit(user_params_stoc->HMF,hs_constants->sigma_cond,hs_constants->growth_out)){ + if(hs_constants->delta >= MAX_DELTAC_FRAC*get_delta_crit(user_params_stoc->HMF,hs_constants->sigma_cond,hs_constants->growth_out)){ *n_halo_out = 1; //Expected mass takes into account potential dexm overlap @@ -709,16 +637,17 @@ int stoc_sample(struct HaloSamplingConstants * hs_constants, gsl_rng * rng, int return 0; } - if(global_params.SAMPLE_METHOD == 0 || (global_params.SAMPLE_METHOD == 3 && !hs_constants->from_catalog)){ - err = stoc_mass_sample(hs_constants, rng, n_halo_out, M_out); - } - else if(global_params.SAMPLE_METHOD == 1){ + //We always use Number-Limited sampling for grid-based cases + if(user_params_stoc->SAMPLE_METHOD == 1 || !hs_constants->from_catalog){ err = stoc_halo_sample(hs_constants, rng, n_halo_out, M_out); } - else if(global_params.SAMPLE_METHOD == 2){ - err = stoc_sheth_sample(hs_constants, rng, n_halo_out, M_out); + else if(user_params_stoc->SAMPLE_METHOD == 0){ + err = stoc_mass_sample(hs_constants, rng, n_halo_out, M_out); + } + else if(user_params_stoc->SAMPLE_METHOD == 2){ + err = stoc_partition_sample(hs_constants, rng, n_halo_out, M_out); } - else if(global_params.SAMPLE_METHOD == 3){ + else if(user_params_stoc->SAMPLE_METHOD == 3){ err = stoc_split_sample(hs_constants, rng, n_halo_out, M_out); } else{ @@ -735,13 +664,9 @@ int stoc_sample(struct HaloSamplingConstants * hs_constants, gsl_rng * rng, int // will have to add properties here and output grids, instead of in perturbed int sample_halo_grids(gsl_rng **rng_arr, double redshift, float *dens_field, float *halo_overlap_box, struct HaloField *halofield_large, struct HaloField *halofield_out, struct HaloSamplingConstants *hs_constants){ int lo_dim = user_params_stoc->HII_DIM; - int hi_dim = user_params_stoc->DIM; double boxlen = user_params_stoc->BOX_LEN; double Mcell = hs_constants->M_cond; - double lnMcell = hs_constants->lnM_cond; - double Mmax_tb = hs_constants->M_max_tables; - double lnMmax_tb = hs_constants->lnM_max_tb; double Mmin = hs_constants->M_min; double lnMmin = hs_constants->lnM_min; double growthf = hs_constants->growth_out; @@ -757,24 +682,6 @@ int sample_halo_grids(gsl_rng **rng_arr, double redshift, float *dens_field, flo LOG_DEBUG("z = %f, Mmin = %e, Mmax = %e,volume = %.3e, D = %.3e",redshift,Mmin,Mcell,Mcell/RHOcrit/cosmo_params_stoc->OMm,growthf); LOG_DEBUG("Total Array Size %llu, array size per thread %llu (~%.3e GB total)",arraysize_total,arraysize_local,6.*arraysize_total*sizeof(int)/1e9); - //Since the conditional MF is extended press-schecter, we rescale by a factor equal to the ratio of the collapsed fractions (n_order == 1) of the UMF - - double ps_ratio = 1.; - if(user_params_stoc->HMF>1 && user_params_stoc->HMF<4){ - struct parameters_gsl_MF_integrals params = { - .redshift = hs_constants->z_out, - .growthf = growthf, - .HMF = 0, - }; - - if(user_params_stoc->INTEGRATION_METHOD_HALOS == 1) - initialise_GL(NGL_INT, lnMmin, lnMcell); - - ps_ratio = IntegratedNdM(lnMmin,lnMcell,params,2,user_params_stoc->INTEGRATION_METHOD_HALOS); - params.HMF = user_params_stoc->HMF; - ps_ratio /= IntegratedNdM(lnMmin,lnMcell,params,2,user_params_stoc->INTEGRATION_METHOD_HALOS); - } - double total_volume_excluded=0.; double total_volume_dexm=0.; double vol_conversion = pow((double)user_params_stoc->HII_DIM / (double)user_params_stoc->DIM,3); @@ -799,7 +706,7 @@ int sample_halo_grids(gsl_rng **rng_arr, double redshift, float *dens_field, flo unsigned long long int count=0; unsigned long long int istart = threadnum * arraysize_local; //debug total - double M_cell=0.; + double M_tot_cell=0.; //we need a private version struct HaloSamplingConstants hs_constants_priv; @@ -819,9 +726,6 @@ int sample_halo_grids(gsl_rng **rng_arr, double redshift, float *dens_field, flo count++; } -//I need the full radii list before the loop starts -#pragma omp barrier - #pragma omp for reduction(+:total_volume_excluded) for (x=0; xSAMPLER_MIN_MASS) continue; if(count >= arraysize_local){ - LOG_ERROR("Too many halos to fit in buffer, try raising global_params.MAXHALO_FACTOR"); + LOG_ERROR("More than %d halos (expected %.1f) with buffer size factor %d", + arraysize_local,arraysize_local/user_params_stoc->MAXHALO_FACTOR,user_params_stoc->MAXHALO_FACTOR); + LOG_ERROR("If you expected to have an above average halo number try raising user_params->MAXHALO_FACTOR"); Throw(ValueError); } @@ -863,25 +769,20 @@ int sample_halo_grids(gsl_rng **rng_arr, double redshift, float *dens_field, flo halofield_out->sfr_rng[istart + count] = prop_buf[1]; count++; - M_cell += hm_buf[i]; + M_tot_cell += hm_buf[i]; if((x+y+z) == 0){ LOG_ULTRA_DEBUG("Halo %d Mass %.2e Stellar %.2e SFR %.2e",i,hm_buf[i],prop_buf[0],prop_buf[1]); } } if((x+y+z) == 0){ - LOG_SUPER_DEBUG("Cell (%d %d %d): delta %.2f | N %d (exp. %.2e) | Total M %.2e (exp. %.2e) overlap defc %.2f ps_ratio %.2f",x,y,z, - delta,nh_buf,hs_constants_priv.expected_N,M_cell,hs_constants_priv.expected_M,mass_defc,ps_ratio); + LOG_SUPER_DEBUG("Cell (%d %d %d): delta %.2f | N %d (exp. %.2e) | Total M %.2e (exp. %.2e) overlap defc %.2f",x,y,z, + delta,nh_buf,hs_constants_priv.expected_N,M_tot_cell,hs_constants_priv.expected_M,mass_defc); } } } } LOG_SUPER_DEBUG("Thread %llu found %llu halos",threadnum,count); - if(count >= arraysize_local){ - LOG_ERROR("Ran out of memory, with %llu halos and local size %llu",count,arraysize_local); - Throw(ValueError); - } - istart_threads[threadnum] = istart; nhalo_threads[threadnum] = count; } @@ -985,7 +886,15 @@ int sample_halo_progenitors(gsl_rng ** rng_arr, double z_in, double z_out, struc for(jj=0;jjSAMPLER_MIN_MASS) continue; + + if(count >= arraysize_local){ + LOG_ERROR("More than %d halos (expected %d) with buffer size factor %d", + arraysize_local,arraysize_local/user_params_stoc->MAXHALO_FACTOR,user_params_stoc->MAXHALO_FACTOR); + LOG_ERROR("If you expected to have an above average halo number try raising user_params_stoc->MAXHALO_FACTOR"); + Throw(ValueError); + } + set_prop_rng(rng_arr[threadnum], 1, corr_arr, propbuf_in, propbuf_out); halofield_out->halo_masses[istart + count] = prog_buf[jj]; @@ -1013,11 +922,6 @@ int sample_halo_progenitors(gsl_rng ** rng_arr, double z_in, double z_out, struc M2,n_prog,hs_constants_priv.expected_N,M_prog,hs_constants_priv.expected_M); } } - if(count >= arraysize_local){ - LOG_ERROR("Ran out of memory, with %llu halos and local size %llu",count,arraysize_local); - Throw(ValueError); - } - istart_threads[threadnum] = istart; nhalo_threads[threadnum] = count; } @@ -1179,21 +1083,13 @@ int single_test_sample(struct UserParams *user_params, struct CosmoParams *cosmo LOG_DEBUG("SINGLE SAMPLE: z = (%.2f,%.2f), Mmin = %.3e, cond(%d)=[%.2e,%.2e,%.2e...]",z_out,z_in,hs_constants->M_min, n_condition,conditions[0],conditions[1],conditions[2]); - //Since the conditional MF is press-schecter, we rescale by a factor equal to the ratio of the collapsed fractions (n_order == 1) of the UMF - double ps_ratio = 1.; + struct parameters_gsl_MF_integrals integral_params = { .redshift = z_out, .growthf = hs_constants->growth_out, .HMF = user_params->HMF, }; - struct parameters_gsl_MF_integrals params_second = integral_params; - if(!hs_constants->from_catalog && user_params_stoc->HMF>1 && user_params_stoc->HMF<4){ - params_second.HMF = 0; - ps_ratio = IntegratedNdM(hs_constants->lnM_min,hs_constants->lnM_max_tb,params_second,2,user_params_stoc->INTEGRATION_METHOD_HALOS); - ps_ratio /= IntegratedNdM(hs_constants->lnM_min,hs_constants->lnM_max_tb,integral_params,2,user_params_stoc->INTEGRATION_METHOD_HALOS); - } - //halo catalogues + cell sums from multiple conditions, given M as cell descendant halos/cells //the result mapping is n_halo_total (1) (exp_n,exp_m,n_prog,m_prog) (n_desc) M_cat (n_prog_total) int n_halo_tot=0; @@ -1217,7 +1113,7 @@ int single_test_sample(struct UserParams *user_params, struct CosmoParams *cosmo n_halo_cond = 0; M_prog = 0; for(i=0;iSAMPLER_MIN_MASS) continue; M_prog += out_hm[i]; n_halo_cond++; @@ -1240,13 +1136,36 @@ int single_test_sample(struct UserParams *user_params, struct CosmoParams *cosmo } } //output descendant statistics - out_n_exp[j] = hs_constants_priv.expected_N; - out_m_exp[j] = hs_constants_priv.expected_M; out_n_cell[j] = n_halo_cond; out_m_cell[j] = M_prog; } + } - out_n_tot[0] = n_halo_tot; + out_n_tot[0] = n_halo_tot; + + //get expected values from the saved mass range + if(hs_constants->from_catalog){ + initialise_dNdM_tables(log(user_params->SAMPLER_MIN_MASS), hs_constants->lnM_max_tb, log(user_params->SAMPLER_MIN_MASS), + hs_constants->lnM_max_tb, hs_constants->growth_out, hs_constants->growth_in, true); + } + else{ + double delta_crit = get_delta_crit(user_params_stoc->HMF,hs_constants->sigma_cond,hs_constants->growth_out); + initialise_dNdM_tables(DELTA_MIN, MAX_DELTAC_FRAC*delta_crit, log(user_params->SAMPLER_MIN_MASS), hs_constants->lnM_max_tb, + hs_constants->growth_out, hs_constants->lnM_cond, false); + } + #pragma omp parallel + { + struct HaloSamplingConstants hs_constants_exp; + hs_constants_exp = *hs_constants; + double cond; + + #pragma omp for + for(j=0;jUSE_INTERPOLATION_TABLES){ diff --git a/src/py21cmfast/src/UsefulFunctions.c b/src/py21cmfast/src/UsefulFunctions.c index 25b9bb617..5fc6f3ffb 100644 --- a/src/py21cmfast/src/UsefulFunctions.c +++ b/src/py21cmfast/src/UsefulFunctions.c @@ -27,7 +27,8 @@ void seed_rng_threads(gsl_rng * rng_arr[], int seed){ unsigned int seeds[user_params_ufunc->N_THREADS]; // For multithreading, seeds for the RNGs are generated from an initial RNG (based on the input random_seed) and then shuffled (Author: Fred Davies) - int num_int = INT_MAX/16; //JD: this was taking a few seconds per snapshot so i reduced the number TODO: init the RNG once + // int num_int = INT_MAX/16; + int num_int = INT_MAX/256; //JD: this was taking a few seconds per snapshot so i reduced the number TODO: init the RNG once int i, thread_num; unsigned int *many_ints = (unsigned int *)malloc((size_t)(num_int*sizeof(unsigned int))); // Some large number of possible integers for (i=0; iM | M_in), the CDF of dNdM_conditional -//NOTE: Assumes you give it ymin as the minimum mass -void initialise_dNdM_tables(double xmin, double xmax, double ymin, double ymax, double growth1, double param, bool from_catalog){ - int nx,ny,np; - double lnM_cond,delta_crit; - int k_lim = from_catalog ? 1 : 0; +void initialise_dNdM_tables(double xmin, double xmax, double ymin, double ymax, double growth_out, double param, bool from_catalog){ + int nx; + double lnM_cond; double sigma_cond; - LOG_DEBUG("Initialising dNdM Table from [[%.2e,%.2e],[%.2e,%.2e]]",xmin,xmax,ymin,ymax); - LOG_DEBUG("D_out %.2e P %.2e up %d",growth1,param,from_catalog); + LOG_DEBUG("Initialising dNdM Tables from [%.2e,%.2e] (Intg. Limits %.2e %.2e)",xmin,xmax,ymin,ymax); + LOG_DEBUG("D_out %.2e P %.2e from_cat %d",growth_out,param,from_catalog); if(!from_catalog){ lnM_cond = param; sigma_cond = EvaluateSigma(lnM_cond); - //current barrier at the condition for bounds checking - delta_crit = get_delta_crit(user_params_it->HMF,sigma_cond,growth1); } - nx = global_params.N_COND_INTERP; - ny = global_params.N_MASS_INTERP; - np = global_params.N_PROB_INTERP; - double xa[nx], ya[ny], pa[np]; + nx = user_params_it->N_COND_INTERP; - int i,j,k; + double xa[nx]; + int i; //set up coordinate grids for(i=0;iN_THREADS) private(i) firstprivate(sigma_cond,lnM_cond) + { + double x; + double delta; + double M_cond; + + #pragma omp for + for(i=0;iHMF,sigma_cond,param)/param*growth_out; + } + else{ + delta = x; + } + + M_cond = exp(lnM_cond); + Nhalo_table.y_arr[i] = Nhalo_Conditional(growth_out,ymin,ymax,M_cond,sigma_cond,delta,0); + Mcoll_table.y_arr[i] = Mcoll_Conditional(growth_out,ymin,ymax,M_cond,sigma_cond,delta,0); + } + } + LOG_DEBUG("Done."); + +} + +struct rf_inv_params{ + double growthf; + double lnM_cond; + double M_cond; + double delta; + double sigma; + + double rf_norm; + double rf_target; +}; + +double dndm_inv_f(double lnM_min, void * params){ + struct rf_inv_params *p = (struct rf_inv_params *)params; + double integral = Nhalo_Conditional(p->growthf,lnM_min,p->lnM_cond,p->M_cond,p->sigma,p->delta,0); + //This ensures that we never find the root if the ratio is zero, since that will set to M_cond + double result = integral == 0 ? 2*user_params_it->MIN_LOGPROB : log(integral / p->rf_norm); + + return result - p->rf_target; +} + +//This table is N(>M | M_in), the CDF of dNdM_conditional +//NOTE: Assumes you give it ymin as the minimum lower-integral limit, and ymax as the maximum +// `param` is either the constant log condition mass for the grid case (!from_catalog) OR the descendant growth factor with from_catalog +void initialise_dNdM_inverse_table(double xmin, double xmax, double lnM_min, double growth_out, double param, bool from_catalog){ + LOG_DEBUG("Initialising dNdM Tables from [%.2e,%.2e] (Intg. Min. %.2e)",xmin,xmax,lnM_min); + LOG_DEBUG("D_out %.2e P %.2e up %d",growth_out,param,from_catalog); + + int nx = user_params_it->N_COND_INTERP; + int np = user_params_it->N_PROB_INTERP; + double xa[nx], pa[np]; + double rf_tol_abs = 1e-4; + double rf_tol_rel = 0.; + + double lnM_cond; + double sigma_cond; + double min_lp = user_params_it->MIN_LOGPROB; + if(!from_catalog){ + lnM_cond = param; + sigma_cond = EvaluateSigma(lnM_cond); + } + + int i,j,k; + //set up coordinate grids + for(i=0;iN_THREADS) private(i,j,k) firstprivate(delta_crit,sigma_cond,lnM_cond) + #pragma omp parallel num_threads(user_params_it->N_THREADS) private(i,j,k) firstprivate(sigma_cond,lnM_cond) { - double x,y,buf; - double norm,fcoll; - double lnM_prev,lnM_p; - double prob; - double p_prev,p_target; + double x; + double norm; + double lnM_lo,lnM_hi,lnM_guess; double delta; + double M_cond; + double lnM_init; + + //RF stuff + int status, iter; + const gsl_root_fsolver_type *T; + gsl_root_fsolver *solver; + gsl_function F; + struct rf_inv_params params_rf; + params_rf.growthf = growth_out; + + F.function = &dndm_inv_f; + F.params = ¶ms_rf; + + T = gsl_root_fsolver_brent; + solver = gsl_root_fsolver_alloc(T); #pragma omp for for(i=0;iHMF,sigma_cond,param)/param*growth1; - delta_crit = get_delta_crit(user_params_it->HMF,sigma_cond,growth1); + //Barrier at descendant mass scaled to progenitor redshift + delta = get_delta_crit(user_params_it->HMF,sigma_cond,param)/param*growth_out; } else{ delta = x; } - lnM_prev = ymin; - p_prev = 0; - - if(ymin >= lnM_cond){ - Nhalo_table.y_arr[i] = 0.; - Mcoll_table.y_arr[i] = 0.; - for(k=1;k= MAX_DELTAC_FRAC*delta_crit){ - Nhalo_table.y_arr[i] = 1/(exp(lnM_cond)); //one halo - Mcoll_table.y_arr[i] = 1.;//both this and Nhalo_table are multiplied by condition mass - for(k=1;kINTEGRATION_METHOD_HALOS == 1) - initialise_GL(NGL_INT, ymin, lnM_cond); + params_rf.M_cond = M_cond; + params_rf.lnM_cond = lnM_cond; + params_rf.delta = delta; + params_rf.sigma = sigma_cond; - norm = Nhalo_Conditional(growth1,ymin,ymax,exp(lnM_cond),sigma_cond,delta,user_params_it->INTEGRATION_METHOD_HALOS); - fcoll = Mcoll_Conditional(growth1,ymin,ymax,exp(lnM_cond),sigma_cond,delta,user_params_it->INTEGRATION_METHOD_HALOS); - Nhalo_table.y_arr[i] = norm; - Mcoll_table.y_arr[i] = fcoll; - // LOG_DEBUG("cond x: %.2e M [%.2e,%.2e] %.2e d %.2f D %.2f n %d ==> %.8e / %.8e",x,exp(ymin),exp(ymax),exp(lnM_cond),delta,growth1,i,norm,fcoll); + //NOTE: The total number density and collapsed fraction must be + norm = Nhalo_Conditional(growth_out,lnM_min,lnM_cond,M_cond,sigma_cond,delta,0); + // LOG_ULTRA_DEBUG("cond x: %.2e M_min %.2e M_cond %.2e d %.4f D %.2f n %d ==> %.8e",x,exp(lnM_min),exp(lnM_cond),delta,growth_out,i,norm); + params_rf.rf_norm = norm; - //if the condition has no halos set the dndm table directly since norm==0 breaks things + //if the condition has no halos set the dndm table directly to avoid integration and divide by zero if(norm==0){ for(k=1;kINTEGRATION_METHOD_HALOS == 1) - initialise_GL(NGL_INT, y, lnM_cond); - - y = ya[j]; - if(lnM_cond <= y){ - //setting to one guarantees samples at lower mass - //This fixes upper mass limits for the conditions - buf = 0.; - } - else{ - buf = Nhalo_Conditional(growth1,y,ymax,exp(lnM_cond),sigma_cond,delta,user_params_it->INTEGRATION_METHOD_HALOS); //Number density between ymin and y - } - prob = buf / norm; - //catch some norm errors - if(prob != prob){ - LOG_ERROR("Normalisation error in inverse halo table generation"); - Throw(TableGenerationError); - } - - //There are times where we have gone over the probability (machine precision) limit before reaching the mass limit - if(!from_catalog){ - if(prob == 0.){ - prob = global_params.MIN_LOGPROB; //to make sure we go over the limit we extrapolate to here - if(y > lnM_cond) y = lnM_cond; + Nhalo_inv_table.z_arr[i][np-1] = exp(lnM_min)/M_cond; + lnM_init = lnM_min; + for(k=np-2;k>=0;k--){ + iter = 0; + lnM_hi = lnM_cond; + params_rf.rf_target = pa[k]; + // LOG_ULTRA_DEBUG("Target %.6e",pa[k]); + gsl_root_fsolver_set(solver, &F, lnM_init, lnM_cond); + do{ + iter++; + status = gsl_root_fsolver_iterate(solver); + lnM_guess = gsl_root_fsolver_root(solver); + lnM_lo = gsl_root_fsolver_x_lower(solver); + lnM_hi = gsl_root_fsolver_x_upper(solver); + status = gsl_root_test_interval(lnM_lo, lnM_hi, rf_tol_abs, rf_tol_rel); + + // LOG_ULTRA_DEBUG("Current step %d | [%.6e,%.6e] - %.6e",iter,lnM_lo,lnM_hi,lnM_guess); + + if (status == GSL_SUCCESS){ + lnM_init = lnM_lo; + Nhalo_inv_table.z_arr[i][k] = exp(lnM_guess)/M_cond; + // LOG_ULTRA_DEBUG("Found %.6e --> %.6e",lnM_guess,exp(lnM_guess)/M_cond); } - else prob = log(prob); - } - // LOG_ULTRA_DEBUG("Int || x: %.2e (%d) y: %.2e (%d) ==> %.8e / %.8e",from_catalog ? exp(x) : x,i,exp(y),j,prob,p_prev); - //loop through the remaining spaces in the inverse table and fill them - while(prob <= pa[k] && k >= k_lim){ - //since we go ascending in y, prob > prob_prev - //NOTE: linear interpolation in (lnM,log(p)|p) - lnM_p = (p_prev-pa[k])*(y - lnM_prev)/(p_prev-prob) + lnM_prev; - Nhalo_inv_table.z_arr[i][k] = lnM_p; - - // LOG_ULTRA_DEBUG("Found c: %.2e p: (%.2e,%.2e,%.2e) (c %d, m %d, p %d) z: %.5e",from_catalog ? exp(x) : x,p_prev,pa[k],prob,i,j,k,exp(lnM_p)); - k--; + }while((status == GSL_CONTINUE) && (iter < MAX_ITER_RF)); + if(status!=GSL_SUCCESS) { + LOG_ERROR("gsl RF error occured! %d",status); + GSL_ERROR(status); } - //keep the value at the previous mass bin for interpolation - p_prev = prob; - lnM_prev = y; } } + gsl_root_fsolver_free(solver); } LOG_DEBUG("Done."); } +double J_integrand(double u, void *params){ + double gamma1 = *(double *) params; + return pow((1. + 1./u/u),gamma1*0.5); +} + +double integrate_J(double u_res, double gamma1){ + double result, error, lower_limit, upper_limit; + gsl_function F; + double rel_tol = 1e-4; //<- relative tolerance + int w_size = 1000; + gsl_integration_workspace * w + = gsl_integration_workspace_alloc (w_size); + + int status; + F.function = &J_integrand; + F.params = &gamma1; + lower_limit = 0.; + upper_limit = u_res; + + gsl_set_error_handler_off(); + status = gsl_integration_qag (&F, lower_limit, upper_limit, 0, rel_tol, + w_size, GSL_INTEG_GAUSS61, w, &result, &error); + + if(status!=0) { + LOG_ERROR("gsl integration error occured!"); + LOG_ERROR("J: gamma1 = %.4e u_res = %.4e",gamma1,u_res); + GSL_ERROR(status); + } + + gsl_integration_workspace_free (w); + + return result; +} + +void initialise_J_split_table(int Nbin, double umin, double umax, double gamma1){ + int i; + if(!J_split_table.allocated) + allocate_RGTable1D(Nbin,&J_split_table); + + J_split_table.x_min = 1e-3; + J_split_table.x_width = (umax - umin)/((double)Nbin-1); + + for(i=0;iUSE_INTERPOLATION_TABLES) return EvaluateRGTable1D(condition,&Nhalo_table); - return Nhalo_Conditional(growthf, lnMmin, lnMmax, M_cond, sigma, delta, user_params_it->INTEGRATION_METHOD_HALOS); + return Nhalo_Conditional(growthf, lnMmin, lnMmax, M_cond, sigma, delta, 0); } double EvaluateMcoll(double condition, double growthf, double lnMmin, double lnMmax, double M_cond, double sigma, double delta){ if(user_params_it->USE_INTERPOLATION_TABLES) return EvaluateRGTable1D(condition,&Mcoll_table); - return Mcoll_Conditional(growthf, lnMmax, lnMmax, M_cond, sigma, delta, user_params_it->INTEGRATION_METHOD_HALOS); + return Mcoll_Conditional(growthf, lnMmin, lnMmax, M_cond, sigma, delta, 0); +} + +//extrapolation function for log-probability based tables +//NOTE: this is very similar to the EvaluateRGTableX function, +// it may be worth allowing extrapolation there by simply setting the indices to the min/max +double extrapolate_dNdM_inverse(double condition, double lnp){ + double x_min = Nhalo_inv_table.x_min; + double x_width = Nhalo_inv_table.x_width; + int x_idx = (int)floor((condition - x_min)/x_width); + double x_table = x_min + x_idx*x_width; + double interp_point_x = (condition - x_table)/x_width; + + double extrap_point_y = (lnp - user_params_it->MIN_LOGPROB)/Nhalo_inv_table.y_width; + + //find the log-mass at the edge of the table for this condition + double xlimit = Nhalo_inv_table.z_arr[x_idx][0]*(interp_point_x) + + Nhalo_inv_table.z_arr[x_idx+1][0]*(1-interp_point_x); + double xlimit_m1 = Nhalo_inv_table.z_arr[x_idx][1]*(interp_point_x) + + Nhalo_inv_table.z_arr[x_idx+1][1]*(1-interp_point_x); + + double result = xlimit + (xlimit_m1 - xlimit)*(extrap_point_y); + + return result; } //This one is always a table double EvaluateNhaloInv(double condition, double prob){ - return EvaluateRGTable2D(condition,prob,&Nhalo_inv_table); + if(prob == 0.) + return 1.; //q == 1 -> condition mass + double lnp = log(prob); + if(lnp < user_params_it->MIN_LOGPROB) + return extrapolate_dNdM_inverse(condition,lnp); + return EvaluateRGTable2D(condition,lnp,&Nhalo_inv_table); +} + +double EvaluateJ(double u_res,double gamma1){ + if(fabs(gamma1) < FRACT_FLOAT_ERR) + return u_res; + //small u approximation + if(u_res < J_split_table.x_min) + return pow(u_res, 1. - gamma1)/(1. - gamma1); + double u_max = J_split_table.x_min + (J_split_table.n_bin-1)*J_split_table.x_width; + //asymptotic expansion + if(u_res > u_max) + return J_split_table.y_arr[J_split_table.n_bin - 1] + u_res - 0.5*gamma1*(1./u_res - 1./u_max); + return EvaluateRGTable1D(u_res,&J_split_table); +} + +//The sigma interp table is regular in log mass, not sigma so we need to loop ONLY FOR SAMPLE_METHOD==2 +//NOTE: This should be improved with its own RGTable but we do not often use this method +double EvaluateSigmaInverse(double sigma){ + if(!user_params_it->USE_INTERPOLATION_TABLES){ + LOG_ERROR("Cannot currently do sigma inverse without USE_INTERPOLATION_TABLES"); + Throw(ValueError); + } + int idx; + for(idx=0;idx Sigma_InterpTable.y_arr[idx]) break; + } + if(idx == NMass){ + LOG_ERROR("sigma inverse out of bounds."); + Throw(TableEvaluationError); + } + double sigma_left = Sigma_InterpTable.y_arr[idx-1]; + double sigma_right = Sigma_InterpTable.y_arr[idx]; + double table_val_left = Sigma_InterpTable.x_min + (idx-1)*Sigma_InterpTable.x_width; //upper lnM + double table_val_right = Sigma_InterpTable.x_min + (idx)*Sigma_InterpTable.x_width; //upper lnM + + double interp_point = (sigma - sigma_left)/(sigma_right - sigma_left); + + return table_val_left*(1-interp_point) + table_val_right*(interp_point); } diff --git a/src/py21cmfast/src/interpolation.c b/src/py21cmfast/src/interpolation.c index 792418561..b6aa96cac 100644 --- a/src/py21cmfast/src/interpolation.c +++ b/src/py21cmfast/src/interpolation.c @@ -144,9 +144,9 @@ double EvaluateRGTable2D(double x, double y, struct RGTable2D *table){ double left_edge, right_edge, result; - // LOG_DEBUG("2D Interp: val (%.2e,%.2e) min (%.2e,%.2e) wid (%.2e,%.2e)",x,y,x_min,y_min,x_width,y_width); - // LOG_DEBUG("2D Interp: idx (%d,%d) tbl (%.2e,%.2e) itp (%.2e,%.2e)",x_idx,y_idx,x_table,y_table,interp_point_x,interp_point_y); - // LOG_DEBUG("2D Interp: table corners (%.2e,%.2e,%.2e,%.2e)",table->z_arr[x_idx][y_idx],table->z_arr[x_idx][y_idx+1],table->z_arr[x_idx+1][y_idx],table->z_arr[x_idx+1][y_idx+1]); + // LOG_ULTRA_DEBUG("2D Interp: val (%.2e,%.2e) min (%.2e,%.2e) wid (%.2e,%.2e)",x,y,x_min,y_min,x_width,y_width); + // LOG_ULTRA_DEBUG("2D Interp: idx (%d,%d) tbl (%.2e,%.2e) itp (%.2e,%.2e)",x_idx,y_idx,x_table,y_table,interp_point_x,interp_point_y); + // LOG_ULTRA_DEBUG("2D Interp: table cell corners (%.2e,%.2e,%.2e,%.2e)",table->z_arr[x_idx][y_idx],table->z_arr[x_idx][y_idx+1],table->z_arr[x_idx+1][y_idx],table->z_arr[x_idx+1][y_idx+1]); left_edge = table->z_arr[x_idx][y_idx]*(1-interp_point_y) + table->z_arr[x_idx][y_idx+1]*(interp_point_y); right_edge = table->z_arr[x_idx+1][y_idx]*(1-interp_point_y) + table->z_arr[x_idx+1][y_idx+1]*(interp_point_y); diff --git a/src/py21cmfast/src/ps.c b/src/py21cmfast/src/ps.c index 92d3049ba..b0a386b37 100644 --- a/src/py21cmfast/src/ps.c +++ b/src/py21cmfast/src/ps.c @@ -38,6 +38,9 @@ static gsl_spline *erfc_spline; #define Mhalo_max (double)(1e16) #define MAX_DELTAC_FRAC (float)0.99 //max delta/deltac for the mass function integrals +#define JENKINS_a (0.73) //Jenkins+01, SMT has 0.707 +#define JENKINS_b (0.34) //Jenkins+01 fit from Barkana+01, SMT has 0.5 +#define JENKINS_c (0.81) //Jenkins+01 from from Barkana+01, SMT has 0.6 bool initialised_ComputeLF = false; @@ -738,16 +741,14 @@ double dsigmasqdm_z0(double M){ ////MASS FUNCTIONS BELOW////// /* sheth correction to delta crit */ -double sheth_delc(double del, double sig){ +double sheth_delc_dexm(double del, double sig){ return sqrt(SHETH_a)*del*(1. + global_params.SHETH_b*pow(sig*sig/(SHETH_a*del*del), global_params.SHETH_c)); } /*DexM uses a fit to this barrier to acheive MF similar to ST, Here I use the fixed version for the sampler*/ +//NOTE: if I made this a table it would save a pow call per condition in the sampler double sheth_delc_fixed(double del, double sig){ - double a = 0.707; - double sheth_b = 0.485; - double sheth_c = 0.615; - return sqrt(a)*del*(1. + sheth_b*pow(sig*sig/(a*del*del), sheth_c)); + return sqrt(JENKINS_a)*del*(1. + JENKINS_b*pow(sig*sig/(JENKINS_a*del*del), JENKINS_c)); } //Get the relevant excursion set barrier density given the user-specified HMF @@ -775,16 +776,14 @@ double dNdlnM_Delos(double growthf, double lnM){ const double exp_factor = -0.469; sigma = EvaluateSigma(lnM); - dsigmadm = EvaluatedSigmasqdm(lnM); - sigma_inv = 1/(sigma); - dsigmadm = dsigmadm * 0.5 * sigma_inv; + sigma_inv = 1/sigma; + dsigmadm = EvaluatedSigmasqdm(lnM) * (0.5*sigma_inv); //d(s^2)/dm z0 to dsdm nu = DELTAC_DELOS*sigma_inv/growthf; dfdnu = coeff_nu*pow(nu,index_nu)*exp(exp_factor*nu*nu); dfdM = dfdnu * fabs(dsigmadm) * sigma_inv; - //NOTE: unlike the other UMFs this is dNdlogM //NOTE: dfdM == constants*dNdlnM return dfdM*cosmo_params_ps->OMm*RHOcrit; } @@ -798,17 +797,14 @@ double dNdlnM_conditional_Delos(double growthf, double lnM, double delta_cond, d const double exp_factor = -0.469; sigma = EvaluateSigma(lnM); - dsigmadm = EvaluatedSigmasqdm(lnM); if(sigma < sigma_cond) return 0.; + dsigmadm = EvaluatedSigmasqdm(lnM) * 0.5; //d(s^2)/dm to s*dsdm sigdiff_inv = sigma == sigma_cond ? 1e6 : 1/(sigma*sigma - sigma_cond*sigma_cond); - sigma_inv = 1/sigma; - dsigmadm = dsigmadm * (0.5*sigma_inv); //d(s^2)/dm z0 to dsdm - nu = (DELTAC_DELOS - delta_cond)*sqrt(sigdiff_inv)/growthf; dfdnu = coeff_nu*pow(nu,index_nu)*exp(exp_factor*nu*nu); - dfdM = dfdnu * fabs(dsigmadm) * sigma_inv * sigma * sigma * sigdiff_inv; + dfdM = dfdnu * fabs(dsigmadm) * sigdiff_inv; //NOTE: like the other CMFs this is dNdlogM and leaves out // the (cosmo_params_ps->OMm)*RHOcrit @@ -819,48 +815,41 @@ double dNdlnM_conditional_Delos(double growthf, double lnM, double delta_cond, d //Sheth Tormen 2002 fit for the CMF, while the moving barrier does not allow for a simple rescaling, it has been found //That a taylor expansion of the barrier shape around the point of interest well approximates the simulations -double st_taylor_factor(double sig, double sig_cond, double delta_cond, double growthf){ +double st_taylor_factor(double sig, double sig_cond, double growthf, double *zeroth_order){ int i; - // double a = SHETH_a; - double a = 0.707; - double alpha = 0.615; //fixed instead of global_params.SHETH_c bc of DexM corrections - double beta = 0.485; //fixed instead of global_params.SHETH_b + double a = JENKINS_a; + double alpha = JENKINS_c; //fixed instead of global_params.SHETH_c bc of DexM corrections + double beta = JENKINS_b; //fixed instead of global_params.SHETH_b double del = Deltac/growthf; double sigsq = sig*sig; + double sigsq_inv = 1./sigsq; double sigcsq = sig_cond*sig_cond; - //See note below double sigdiff = sig == sig_cond ? 1e-6 : sigsq - sigcsq; - double dn_const = sqrt(a)*del*beta*pow(a*del*del,-alpha); - //define arrays of factors to save time and math calls - int n_fac[6] = {1,1,2,6,24,120}; - double a_fac[6]; - double s_fac[6]; - a_fac[0] = 1; - s_fac[0] = 1; - for(i=1;i<=5;i++){ - a_fac[i] = a_fac[i-1] * (alpha-i+1); - s_fac[i] = s_fac[i-1] * (-sigdiff); - } + // This array cumulatively builds the taylor series terms + // sigdiff^n / n! * df/dsigma (polynomial w alpha) + double t_array[6]; + t_array[0] = 1.; + for(i=1;i<6;i++) + t_array[i] = t_array[i-1] * (-sigdiff) / i * (alpha-i+1) * sigsq_inv; + //Sum small to large double result = 0.; - //Taylor expansion of the x^a part around (sigsq - sigcondsq) (summing small to large) - for(i=5;i>=1;i--){ - result += s_fac[i]/n_fac[i] * pow(sigsq,alpha-i)*a_fac[i]; - // LOG_ULTRA_DEBUG("%d term %.2e",i,result); + for(i=5;i>=0;i--){ + result += t_array[i]; } - result *= dn_const; - //add the constant terms from the 0th derivative of the barrier (condition delta independent of halo sigma) - // result += sqrt(a)*delta_crit - delta_cond; + double prefactor_1 = sqrt(a)*del; + double prefactor_2 = beta*pow(sigsq_inv*(a*del*del),-alpha); + + result = prefactor_1*(1 + prefactor_2*result); + *zeroth_order = prefactor_1*(1+prefactor_2); //0th order term gives the barrier for efficiency return result; } -//CMF Corresponding to the Sheth Mo Tormen HMF, here we assume that we are passing the correct delta2, -// which is the condition delta, the barrier delta1 is set by the mass, so it is passed usually as Deltac -//NOTE: Currently broken and I don't know why +//CMF Corresponding to the Sheth Mo Tormen HMF (Sheth+ 2002) double dNdM_conditional_ST(double growthf, double lnM, double delta_cond, double sigma_cond){ double sigma1, dsigmasqdm, Barrier, factor, sigdiff_inv, result; double delta_0 = delta_cond / growthf; @@ -868,8 +857,7 @@ double dNdM_conditional_ST(double growthf, double lnM, double delta_cond, double dsigmasqdm = EvaluatedSigmasqdm(lnM); if(sigma1 < sigma_cond) return 0.; - Barrier = sheth_delc_fixed(Deltac/growthf,sigma1); - factor = st_taylor_factor(sigma1,sigma_cond,delta_0,growthf) + (Barrier - delta_0); + factor = st_taylor_factor(sigma1,sigma_cond,growthf,&Barrier) - delta_0; sigdiff_inv = sigma1 == sigma_cond ? 1e6 : 1/(sigma1*sigma1 - sigma_cond*sigma_cond); @@ -1200,7 +1188,7 @@ double IntegratedNdM_QAG(double lnM_lo, double lnM_hi, struct parameters_gsl_MF_ double result, error, lower_limit, upper_limit; gsl_function F; // double rel_tol = FRACT_FLOAT_ERR*128; //<- relative tolerance - double rel_tol = 1e-3; //<- relative tolerance + double rel_tol = 1e-4; //<- relative tolerance int w_size = 1000; gsl_integration_workspace * w = gsl_integration_workspace_alloc (w_size); @@ -1550,6 +1538,8 @@ double Nhalo_Conditional(double growthf, double lnM1, double lnM2, double M_cond .delta = delta, }; + if(delta <= -1. || lnM1 >= log(M_cond)) + return 0.; //return 1 halo AT THE CONDITION MASS if delta is exceeded if(delta > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma,growthf)){ if(M_cond*(1-FRACT_FLOAT_ERR) <= exp(lnM2)) //this limit is not ideal, but covers floating point errors when we set lnM2 == log(M_cond) @@ -1557,8 +1547,6 @@ double Nhalo_Conditional(double growthf, double lnM1, double lnM2, double M_cond else return 0.; } - if(delta <= -1.) - return 0.; return IntegratedNdM(lnM1,lnM2,params,-1, method); } @@ -1571,6 +1559,8 @@ double Mcoll_Conditional(double growthf, double lnM1, double lnM2, double M_cond .delta = delta, }; + if(delta <= -1. || lnM1 >= log(M_cond)) + return 0.; //return 100% of mass AT THE CONDITION MASS if delta is exceeded if(delta > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma,growthf)){ if(M_cond*(1-FRACT_FLOAT_ERR) <= exp(lnM2)) //this limit is not ideal, but covers floating point errors when we set lnM2 == log(M_cond) @@ -1578,9 +1568,6 @@ double Mcoll_Conditional(double growthf, double lnM1, double lnM2, double M_cond else return 0.; } - if(delta <= -1.) - return 0.; - return IntegratedNdM(lnM1,lnM2,params,-2, method); } @@ -1597,12 +1584,13 @@ double Nion_ConditionalM_MINI(double growthf, double lnM1, double lnM2, double M .f_esc_norm = Fesc7, .Mlim_star = Mlim_Fstar, .Mlim_esc = Mlim_Fesc, - // .HMF = user_params_ps->HMF, - .HMF = 0, //FORCE EPS UNTIL THE OTHERS WORK + .HMF = user_params_ps->HMF, .sigma_cond = sigma2, .delta = delta2, }; + if(delta2 <= -1. || lnM1 >= log(M_cond)) + return 0.; //return 1 halo at the condition mass if delta is exceeded //NOTE: this will almost always be zero, due to the upper turover, // however this replaces an integral so it won't be slow @@ -1612,8 +1600,11 @@ double Nion_ConditionalM_MINI(double growthf, double lnM1, double lnM2, double M else return 0.; } - if(delta2 <= -1.) - return 0.; + + //If we don't have a corresponding CMF, use EPS and normalise + //NOTE: it's possible we may want to use another default + if(params.HMF != 0 && params.HMF != 1 && params.HMF != 4) + params.HMF = 0; // LOG_ULTRA_DEBUG("params: D=%.2e Mtl=%.2e Mtu=%.2e as=%.2e ae=%.2e fs=%.2e fe=%.2e Ms=%.2e Me=%.2e hmf=%d sig=%.2e del=%.2e", // growthf,MassTurnover,MassTurnover_upper,Alpha_star,Alpha_esc,Fstar7,Fesc7,Mlim_Fstar,Mlim_Fesc,0,sigma2,delta2); @@ -1632,12 +1623,13 @@ double Nion_ConditionalM(double growthf, double lnM1, double lnM2, double M_cond .f_esc_norm = Fesc10, .Mlim_star = Mlim_Fstar, .Mlim_esc = Mlim_Fesc, - // .HMF = user_params_ps->HMF, - .HMF = 0, //FORCE EPS UNTIL THE OTHERS WORK + .HMF = user_params_ps->HMF, .sigma_cond = sigma2, .delta = delta2, }; + if(delta2 <= -1. || lnM1 >= log(M_cond)) + return 0.; //return 1 halo at the condition mass if delta is exceeded if(delta2 > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma2,growthf)){ if(M_cond*(1-FRACT_FLOAT_ERR) <= exp(lnM2)) @@ -1645,8 +1637,11 @@ double Nion_ConditionalM(double growthf, double lnM1, double lnM2, double M_cond else return 0.; } - if(delta2 <= -1.) - return 0.; + + //If we don't have a corresponding CMF, use EPS and normalise + //NOTE: it's possible we may want to use another default + if(params.HMF != 0 && params.HMF != 1 && params.HMF != 4) + params.HMF = 0; // LOG_ULTRA_DEBUG("params: D=%.2e Mtl=%.2e as=%.2e ae=%.2e fs=%.2e fe=%.2e Ms=%.2e Me=%.2e sig=%.2e del=%.2e", // growthf,MassTurnover,Alpha_star,Alpha_esc,Fstar10,Fesc10,Mlim_Fstar,Mlim_Fesc,sigma2,delta2); diff --git a/src/py21cmfast/wrapper.py b/src/py21cmfast/wrapper.py index 12fc17da2..020d583c7 100644 --- a/src/py21cmfast/wrapper.py +++ b/src/py21cmfast/wrapper.py @@ -1120,7 +1120,7 @@ def determine_halo_list( hbuffer_size = lib.expected_nhalo( redshift, user_params(), cosmo_params(), astro_params(), flag_options() ) - hbuffer_size = int((hbuffer_size + 1) * global_params.MAXHALO_FACTOR) + hbuffer_size = int((hbuffer_size + 1) * user_params.MAXHALO_FACTOR) # set a minimum in case of fluctuation at high z hbuffer_size = int(max(hbuffer_size, 1e6)) @@ -1272,7 +1272,7 @@ def perturb_halo_list( hbuffer_size = lib.expected_nhalo( redshift, user_params(), cosmo_params(), astro_params(), flag_options() ) - hbuffer_size = int((hbuffer_size + 1) * global_params.MAXHALO_FACTOR) + hbuffer_size = int((hbuffer_size + 1) * user_params.MAXHALO_FACTOR) # set a minimum in case of fluctuation at high z hbuffer_size = int(max(hbuffer_size, 1e3)) else: @@ -1458,7 +1458,7 @@ def halo_box( not isinstance(perturbed_field, PerturbedField) or not perturbed_field.is_computed ): - if flag_options.FIXED_HALO_GRIDS or global_params.AVG_BELOW_SAMPLER: + if flag_options.FIXED_HALO_GRIDS or user_params.AVG_BELOW_SAMPLER: perturbed_field = perturb_field( user_params=user_params, cosmo_params=cosmo_params, diff --git a/tests/test_c_interpolation_tables.py b/tests/test_c_interpolation_tables.py index 1d7c968b3..f04bb8ac4 100644 --- a/tests/test_c_interpolation_tables.py +++ b/tests/test_c_interpolation_tables.py @@ -108,7 +108,7 @@ def test_sigma_table(name, plt): # for now this is acceptable # @pytest.mark.xfail @pytest.mark.parametrize("name", options_hmf) -def test_Massfunc_conditional_tables(name): +def test_Massfunc_conditional_tables(name, plt): redshift, kwargs = OPTIONS_HMF[name] opts = prd.get_all_options(redshift, **kwargs) @@ -174,9 +174,12 @@ def test_Massfunc_conditional_tables(name): N_cmfi_cell = ( N_cmfi_cell / max_in_d ) # to get P(>M) since the y-axis is the lower integral limit - mask_cell = N_cmfi_cell > np.exp( - global_params.MIN_LOGPROB - ) # we don't want to compare below the min prob + mask_cell = ( + (N_cmfi_cell > np.exp(global_params.MIN_LOGPROB)) + & (arg_list_inv_d[0] > -1) + & (arg_list_inv_d[0] < delta_crit) + ) + # we don't want to compare below the min prob or outside our delta range N_cmfi_cell = np.clip(N_cmfi_cell, np.exp(global_params.MIN_LOGPROB), 1) N_cmfi_cell = np.log(N_cmfi_cell) @@ -215,6 +218,14 @@ def test_Massfunc_conditional_tables(name): np.log(cell_mass), False, ) + lib.initialise_dNdM_inverse_table( + edges_d[0], + edges_d[-1], + edges_ln[0], + growth_out, + np.log(cell_mass), + False, + ) M_exp_cell = ( np.vectorize(lib.EvaluateMcoll)(edges_d[:-1], 0.0, 0.0, 0.0, 0.0, 0.0, 0.0) @@ -245,9 +256,12 @@ def test_Massfunc_conditional_tables(name): N_cmfi_halo = N_cmfi_halo / ( N_cmfi_halo[:, :1] + np.all(N_cmfi_halo == 0, axis=1)[:, None] ) # if all entries are zero, do not nan the row, just divide by 1 - mask_halo = ( - arg_list_inv_m[0] > arg_list_inv_m[1] - ) # we don't want to compare above the conditinos + mask_halo = (arg_list_inv_m[0] > arg_list_inv_m[1]) & ( + N_cmfi_cell > np.exp(global_params.MIN_LOGPROB) + ) # we don't want to compare above the conditions or below min probability + + N_cmfi_halo = np.clip(N_cmfi_halo, np.exp(global_params.MIN_LOGPROB), 1) + N_cmfi_halo = np.log(N_cmfi_halo) M_cmf_halo = ( np.vectorize(lib.Mcoll_Conditional)( @@ -284,6 +298,15 @@ def test_Massfunc_conditional_tables(name): growth_in, True, ) + lib.initialise_dNdM_inverse_table( + edges_ln[0], + edges_ln[-1], + edges_ln[0], + growth_out, + growth_in, + True, + ) + M_exp_halo = ( np.vectorize(lib.EvaluateMcoll)(edges_ln[:-1], 0.0, 0.0, 0.0, 0.0, 0.0, 0.0) * edges[:-1] @@ -300,6 +323,77 @@ def test_Massfunc_conditional_tables(name): # NOTE: The tables get inaccurate in the smallest halo bin where the condition mass approaches the minimum # We set the absolute tolerance to be insiginificant in sampler terms (~1% of the smallest halo) abs_tol_halo = 1e-2 + + if plt == mpl.pyplot: + xl = edges_d[:-1].size + sel = (xl * np.arange(6) / 6).astype(int) + massfunc_table_comparison_plot( + edges[:-1], + edges[sel], + np.exp(N_cmfi_halo[sel, :]), + np.exp(N_inverse_halo[sel, :]), + edges_d[sel], + np.exp(N_cmfi_cell[sel, :]), + np.exp(N_inverse_cell[sel, :]), + plt, + ) + + # mask out relevant quantities + print(mask_halo.shape, N_inverse_halo.shape, arg_list_inv_m[1].shape) + N_inverse_halo[~mask_halo] = arg_list_inv_m[1][~mask_halo] + N_inverse_cell[~mask_cell] = arg_list_inv_d[1][~mask_cell] + + print_failure_stats( + N_cmf_halo, + N_exp_halo, + [edges[:-1]], + abs_tol_halo, + RELATIVE_TOLERANCE, + "expected N halo", + ) + print_failure_stats( + M_cmf_halo, + M_exp_halo, + [edges[:-1]], + abs_tol_halo, + RELATIVE_TOLERANCE, + "expected M halo", + ) + + print_failure_stats( + N_cmf_cell, + N_exp_cell, + [edges_d[:-1]], + abs_tol_halo, + RELATIVE_TOLERANCE, + "expected N cell", + ) + print_failure_stats( + M_cmf_cell, + M_exp_cell, + [edges_d[:-1]], + abs_tol_halo, + RELATIVE_TOLERANCE, + "expected M cell", + ) + + print_failure_stats( + np.exp(arg_list_inv_m[1]), + np.exp(N_inverse_halo), + arg_list_inv_d, + 0.0, + RELATIVE_TOLERANCE, + "Inverse Halo", + ) + print_failure_stats( + np.exp(arg_list_inv_d[1]), + np.exp(N_inverse_cell), + arg_list_inv_m, + 0.0, + RELATIVE_TOLERANCE, + "Inverse cell", + ) + np.testing.assert_allclose( N_cmf_halo, N_exp_halo, atol=abs_tol_halo, rtol=RELATIVE_TOLERANCE ) @@ -1320,15 +1414,83 @@ def print_failure_stats(test, truth, input_arr, abs_tol, rel_tol, name): sel_failed = np.fabs(truth - test) > (abs_tol + np.fabs(truth) * rel_tol) if sel_failed.sum() > 0: print( - f"{name}: atol {abs_tol} rtol {rel_tol} subcube of failures [min] [max] {np.argwhere(sel_failed).min(axis=0)} {np.argwhere(sel_failed).max(axis=0)}" + f"{name}: atol {abs_tol} rtol {rel_tol} failed {sel_failed.sum()} of {sel_failed.size} {sel_failed.sum() / sel_failed.size * 100:.4f}%" + ) + print( + f"subcube of failures [min] [max] {np.argwhere(sel_failed).min(axis=0)} {np.argwhere(sel_failed).max(axis=0)}" ) for i, row in enumerate(input_arr): print( - f"failure range of inputs axis {i} {row[sel_failed].min()} {row[sel_failed].max()}" + f"failure range of inputs axis {i} {row[sel_failed].min():.2e} {row[sel_failed].max():.2e}" ) print( - f"failure range truth ({truth[sel_failed].min()},{truth[sel_failed].max()}) test ({test[sel_failed].min()},{test[sel_failed].max()})" + f"failure range truth ({truth[sel_failed].min():.3e},{truth[sel_failed].max():.3e}) test ({test[sel_failed].min():.3e},{test[sel_failed].max():.3e})" ) print( - f"max abs diff of failures {np.fabs(truth - test)[sel_failed].max()} relative {(np.fabs(truth - test) / truth)[sel_failed].max()}" + f"max abs diff of failures {np.fabs(truth - test)[sel_failed].max():.4e} relative {(np.fabs(truth - test) / truth)[sel_failed].max():.4e}" + ) + + +def massfunc_table_comparison_plot( + massbins, + conds_m, + N_cmfi_halo, + M_inverse_halo, + conds_d, + N_cmfi_cell, + M_inverse_cell, + plt, +): + fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(16, 8)) + for ax in axs: + ax.grid() + ax.set_xlabel("M_halo") + ax.set_xlim([1e7, 5e11]) + ax.set_xscale("log") + ax.set_ylabel("P(>M)") + ax.set_ylim([1e-8, 1.2]) + ax.set_yscale("log") + + [ + axs[0].plot( + massbins, + N_cmfi_halo[i, :], + color=f"C{i:d}", + linestyle="-", + label=f"M = {m:.3e}", + ) + for i, m in enumerate(conds_m) + ] + [ + axs[0].plot( + M_inverse_halo[i, :], + N_cmfi_halo[i, :], + color=f"C{i:d}", + linestyle=":", + linewidth=3, + ) + for i, m in enumerate(conds_m) + ] + axs[0].legend() + + [ + axs[1].plot( + massbins, + N_cmfi_cell[i, :], + color=f"C{i:d}", + linestyle="-", + label=f"d = {d:.3f}", + ) + for i, d in enumerate(conds_d) + ] + [ + axs[1].plot( + M_inverse_cell[i, :], + N_cmfi_cell[i, :], + color=f"C{i:d}", + linestyle=":", + linewidth=3, ) + for i, m in enumerate(conds_d) + ] + axs[1].legend() diff --git a/tests/test_halo_sampler.py b/tests/test_halo_sampler.py index 8a3068a46..f1e906e8e 100644 --- a/tests/test_halo_sampler.py +++ b/tests/test_halo_sampler.py @@ -1,8 +1,7 @@ import pytest +import matplotlib as mpl import numpy as np -from astropy import constants as c -from astropy import units as u from py21cmfast import AstroParams, CosmoParams, FlagOptions, UserParams, global_params from py21cmfast.c_21cmfast import ffi, lib @@ -14,10 +13,13 @@ options_hmf = list(cint.OPTIONS_HMF.keys()) +options_delta = [-0.9, 0, 1, 1.6] # cell densities to draw samples from +options_mass = [1e8, 1e9, 1e10, 1e11, 1e12] # halo masses to draw samples from + -# @pytest.mark.xfail @pytest.mark.parametrize("name", options_hmf) -def test_sampler(name): +@pytest.mark.parametrize("mass", options_mass) +def test_sampler_from_catalog(name, mass, plt): redshift, kwargs = cint.OPTIONS_HMF[name] opts = prd.get_all_options(redshift, **kwargs) @@ -31,159 +33,267 @@ def test_sampler(name): lib.Broadcast_struct_global_IT(up(), cp(), ap(), fo()) lib.Broadcast_struct_global_STOC(up(), cp(), ap(), fo()) - l10min = 8 - l10max = 11 - edges = np.logspace(l10min, l10max, num=10 * (l10max - l10min)) + l10min = np.log10(global_params.SAMPLER_MIN_MASS) + l10max = np.log10(mass) + edges = np.logspace(l10min, l10max, num=int(10 * (l10max - l10min))) bin_minima = edges[:-1] bin_maxima = edges[1:] dlnm = np.log(edges[1:]) - np.log(edges[:-1]) - # centres = (edges[:-1] * np.exp(dlnm / 2)).astype("f4") lib.init_ps() lib.initialiseSigmaMInterpTable( global_params.M_MIN_INTEGRAL, global_params.M_MAX_INTEGRAL ) - n_cond = 1000 - - conditions_d = np.array([-0.9, 0, 0.5, 1, 1.4]) - conditions_m = np.array([1e8, 1e9, 1e10, 1e11, 1e12]) + n_cond = 30000 z = 6.0 z_prev = 5.8 growth_prev = lib.dicke(z_prev) growthf = lib.dicke(z) - sigma_cond_m = np.vectorize(lib.sigma_z0)(conditions_m) + sigma_cond_m = lib.sigma_z0(mass) delta_cond_m = ( - np.vectorize(lib.get_delta_crit)(up.HMF, sigma_cond_m, growth_prev) - * growthf - / growth_prev + lib.get_delta_crit(up.HMF, sigma_cond_m, growth_prev) * growthf / growth_prev + ) + mass_dens = cp.cosmo.Om0 * cp.cosmo.critical_density(0).to("Mpc-3 M_sun").value + volume_total_m = mass * n_cond / mass_dens + + crd_in = np.zeros(3 * n_cond).astype("i4") + # HALO MASS CONDITIONS WITH FIXED z-step + cond_in = np.full(n_cond, fill_value=mass).astype("f4") # mass at z6 + + nhalo_out = np.zeros(1).astype("i4") + N_out = np.zeros(n_cond).astype("i4") + M_out = np.zeros(n_cond).astype("f8") + exp_M = np.zeros(n_cond).astype("f8") + exp_N = np.zeros(n_cond).astype("f8") + halomass_out = np.zeros(int(1e8)).astype("f4") + halocrd_out = np.zeros(int(3e8)).astype("i4") + + lib.single_test_sample( + up(), + cp(), + ap(), + fo(), + 12345, + n_cond, + ffi.cast("float *", cond_in.ctypes.data), + ffi.cast("int *", crd_in.ctypes.data), + z, + z_prev, + ffi.cast("int *", nhalo_out.ctypes.data), + ffi.cast("int *", N_out.ctypes.data), + ffi.cast("double *", exp_N.ctypes.data), + ffi.cast("double *", M_out.ctypes.data), + ffi.cast("double *", exp_M.ctypes.data), + ffi.cast("float *", halomass_out.ctypes.data), + ffi.cast("int *", halocrd_out.ctypes.data), + ) + + # since the tables are reallocated in the test sample function, we redo them here + lib.initialiseSigmaMInterpTable(edges[0] / 2, edges[-1]) + + # get CMF integrals in the same bins + bin_minima = edges[:-1] + bin_maxima = edges[1:] + binned_cmf = np.vectorize(lib.Nhalo_Conditional)( + growthf, + np.log(bin_minima), + np.log(bin_maxima), + mass, + sigma_cond_m, + delta_cond_m, + up.INTEGRATION_METHOD_HALOS, + ) + + hist, _ = np.histogram(halomass_out, edges) + mf_out = hist / volume_total_m / dlnm + binned_cmf = binned_cmf * n_cond / volume_total_m / dlnm * mass + + if plt == mpl.pyplot: + plot_sampler_comparison( + edges, + exp_N, + exp_M, + N_out, + M_out, + binned_cmf, + mf_out, + f"mass = {mass:.2e}", + plt, + ) + + np.testing.assert_allclose(N_out.mean(), exp_N[0], rtol=RELATIVE_TOLERANCE) + np.testing.assert_allclose(M_out.mean(), exp_M[0], rtol=RELATIVE_TOLERANCE) + np.testing.assert_allclose(mf_out, binned_cmf, rtol=RELATIVE_TOLERANCE) + + +@pytest.mark.parametrize("name", options_hmf) +@pytest.mark.parametrize("delta", options_delta) +def test_sampler_from_grid(name, delta, plt): + redshift, kwargs = cint.OPTIONS_HMF[name] + opts = prd.get_all_options(redshift, **kwargs) + + up = UserParams(opts["user_params"]) + cp = CosmoParams(opts["cosmo_params"]) + ap = AstroParams(opts["astro_params"]) + fo = FlagOptions(opts["flag_options"]) + up.update(USE_INTERPOLATION_TABLES=True) + lib.Broadcast_struct_global_PS(up(), cp()) + lib.Broadcast_struct_global_UF(up(), cp()) + lib.Broadcast_struct_global_IT(up(), cp(), ap(), fo()) + lib.Broadcast_struct_global_STOC(up(), cp(), ap(), fo()) + + lib.init_ps() + lib.initialiseSigmaMInterpTable( + global_params.M_MIN_INTEGRAL, global_params.M_MAX_INTEGRAL ) + n_cond = 30000 + + z = 6.0 + growthf = lib.dicke(z) + mass_dens = cp.cosmo.Om0 * cp.cosmo.critical_density(0).to("Mpc-3 M_sun").value cellvol = (up.BOX_LEN / up.HII_DIM) ** 3 cell_mass = cellvol * mass_dens - volume_total_m = conditions_d * n_cond / mass_dens - cond_mass_m = conditions_d + l10min = np.log10(global_params.SAMPLER_MIN_MASS) + l10max = np.log10(cell_mass) + edges = np.logspace(l10min, l10max, num=int(10 * (l10max - l10min))) + bin_minima = edges[:-1] + bin_maxima = edges[1:] + dlnm = np.log(edges[1:]) - np.log(edges[:-1]) - delta_cond_d = conditions_d - sigma_cond_d = np.full_like(delta_cond_d, lib.sigma_z0(cell_mass)) - volume_total_d = np.full_like(delta_cond_d, cellvol * n_cond) - cond_mass_d = cell_mass + sigma_cond = lib.sigma_z0(cell_mass) + volume_total = cellvol * n_cond crd_in = np.zeros(3 * n_cond).astype("i4") - # CELL DENSITY CONDITIONS WITH FIXED SIZE - for i, d in enumerate(conditions_d): - cond_in = np.full(n_cond, fill_value=d).astype("f4") # mass at z6 - - nhalo_out = np.zeros(1).astype("i4") - N_out = np.zeros(n_cond).astype("i4") - M_out = np.zeros(n_cond).astype("f8") - exp_M = np.zeros(n_cond).astype("f8") - exp_N = np.zeros(n_cond).astype("f8") - halomass_out = np.zeros(int(1e8)).astype("f4") - halocrd_out = np.zeros(int(3e8)).astype("i4") - - print(f"starting {i} d={d}", flush=True) - lib.single_test_sample( - up(), - cp(), - ap(), - fo(), - 12345, # TODO: homogenize - n_cond, - ffi.cast("float *", cond_in.ctypes.data), - ffi.cast("int *", crd_in.ctypes.data), - z, - -1, - ffi.cast("int *", nhalo_out.ctypes.data), - ffi.cast("int *", N_out.ctypes.data), - ffi.cast("double *", exp_N.ctypes.data), - ffi.cast("double *", M_out.ctypes.data), - ffi.cast("double *", exp_M.ctypes.data), - ffi.cast("float *", halomass_out.ctypes.data), - ffi.cast("int *", halocrd_out.ctypes.data), - ) - - # since the tables are reallocated in the test sample function, we redo them here - lib.initialiseSigmaMInterpTable(edges[0] / 2, edges[-1]) - - # get CMF integrals in the same bins - bin_minima = edges[:-1] - bin_maxima = edges[1:] - binned_cmf = np.vectorize(lib.Nhalo_Conditional)( - growthf, - np.log(bin_minima), - np.log(bin_maxima), - cell_mass, - sigma_cond_d[i], - delta_cond_d[i], - up.INTEGRATION_METHOD_HALOS, - ) + cond_in = np.full(n_cond, fill_value=delta).astype("f4") # mass at z6 + + nhalo_out = np.zeros(1).astype("i4") + N_out = np.zeros(n_cond).astype("i4") + M_out = np.zeros(n_cond).astype("f8") + exp_M = np.zeros(n_cond).astype("f8") + exp_N = np.zeros(n_cond).astype("f8") + halomass_out = np.zeros(int(1e8)).astype("f4") + halocrd_out = np.zeros(int(3e8)).astype("i4") + + lib.single_test_sample( + up(), + cp(), + ap(), + fo(), + 12345, # TODO: homogenize + n_cond, + ffi.cast("float *", cond_in.ctypes.data), + ffi.cast("int *", crd_in.ctypes.data), + z, + -1, + ffi.cast("int *", nhalo_out.ctypes.data), + ffi.cast("int *", N_out.ctypes.data), + ffi.cast("double *", exp_N.ctypes.data), + ffi.cast("double *", M_out.ctypes.data), + ffi.cast("double *", exp_M.ctypes.data), + ffi.cast("float *", halomass_out.ctypes.data), + ffi.cast("int *", halocrd_out.ctypes.data), + ) - hist, _ = np.histogram(halomass_out, edges) - mf_out = hist / volume_total_d[i] / dlnm - binned_cmf = binned_cmf * n_cond / volume_total_d[i] / dlnm * cond_mass_d + # since the tables are reallocated in the test sample function, we redo them here + lib.initialiseSigmaMInterpTable(edges[0] / 2, edges[-1]) - np.testing.assert_allclose(N_out.mean(), exp_N[0], rtol=RELATIVE_TOLERANCE) - np.testing.assert_allclose(M_out.mean(), exp_M[0], rtol=RELATIVE_TOLERANCE) - np.testing.assert_allclose(mf_out, binned_cmf, rtol=RELATIVE_TOLERANCE) + # get CMF integrals in the same bins + bin_minima = edges[:-1] + bin_maxima = edges[1:] + binned_cmf = np.vectorize(lib.Nhalo_Conditional)( + growthf, + np.log(bin_minima), + np.log(bin_maxima), + cell_mass, + sigma_cond, + delta, + up.INTEGRATION_METHOD_HALOS, + ) - # HALO MASS CONDITIONS WITH FIXED z-step - for i, m in enumerate(conditions_m): - cond_in = np.full(n_cond, fill_value=m).astype("f4") # mass at z6 - - nhalo_out = np.zeros(1).astype("i4") - N_out = np.zeros(n_cond).astype("i4") - M_out = np.zeros(n_cond).astype("f8") - exp_M = np.zeros(n_cond).astype("f8") - exp_N = np.zeros(n_cond).astype("f8") - halomass_out = np.zeros(int(1e8)).astype("f4") - halocrd_out = np.zeros(int(3e8)).astype("i4") - - print(f"starting {i} m={m:.2e}", flush=True) - lib.single_test_sample( - up(), - cp(), - ap(), - fo(), - 12345, - n_cond, - ffi.cast("float *", cond_in.ctypes.data), - ffi.cast("int *", crd_in.ctypes.data), - z, - z_prev, - ffi.cast("int *", nhalo_out.ctypes.data), - ffi.cast("int *", N_out.ctypes.data), - ffi.cast("double *", exp_N.ctypes.data), - ffi.cast("double *", M_out.ctypes.data), - ffi.cast("double *", exp_M.ctypes.data), - ffi.cast("float *", halomass_out.ctypes.data), - ffi.cast("int *", halocrd_out.ctypes.data), + hist, _ = np.histogram(halomass_out, edges) + mf_out = hist / volume_total / dlnm + binned_cmf = binned_cmf * n_cond / volume_total / dlnm * cell_mass + + if plt == mpl.pyplot: + plot_sampler_comparison( + edges, + exp_N, + exp_M, + N_out, + M_out, + binned_cmf, + mf_out, + f"delta = {delta:.2f}", + plt, ) - # since the tables are reallocated in the test sample function, we redo them here - lib.initialiseSigmaMInterpTable(edges[0] / 2, edges[-1]) - - # get CMF integrals in the same bins - bin_minima = edges[:-1] - bin_maxima = edges[1:] - binned_cmf = np.vectorize(lib.Nhalo_Conditional)( - growthf, - np.log(bin_minima), - np.log(bin_maxima), - sigma_cond_m[i], - delta_cond_m[i], - up.INTEGRATION_METHOD_HALOS, - ) + np.testing.assert_allclose(N_out.mean(), exp_N[0], rtol=RELATIVE_TOLERANCE) + np.testing.assert_allclose(M_out.mean(), exp_M[0], rtol=RELATIVE_TOLERANCE) + np.testing.assert_allclose(mf_out, binned_cmf, rtol=RELATIVE_TOLERANCE) + + +def plot_sampler_comparison( + bin_edges, exp_N, exp_M, N_array, M_array, exp_mf, mf_out, title, plt +): + fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(16, 8)) + + axst = [axs[0].twinx(), axs[1].twiny()] + + # total M axis + axs[1].set_xlabel("M") + # total N axis + axst[1].set_ylabel("N") + + # mass function axis + axs[0].set_title(title) + axs[0].set_ylim([1e-2, 1e4]) + axs[0].set_xlim([bin_edges[0], bin_edges[-1]]) + axs[0].set_yscale("log") + axs[0].set_ylabel("dn/dlnM") + + # ratio axis + axst[0].set_ylim([1e-1, 1e1]) + axst[0].set_yscale("log") + axst[0].set_ylabel("ratio") + + for ax in axs: + ax.grid() + ax.set_xscale("log") + + # log-spaced bins + dlnm = np.log(bin_edges[1:]) - np.log(bin_edges[:-1]) + bin_centres = (bin_edges[:-1] * np.exp(dlnm / 2)).astype("f4") + edges_n = np.arange(N_array.max() + 1) - 0.5 + centres_n = (edges_n[:-1] + edges_n[1:]) / 2 + + hist_n, _ = np.histogram(N_array, edges_n) + p_n = hist_n / N_array.size + hist_m, _ = np.histogram(M_array, bin_edges) + p_m = hist_m / M_array.size / dlnm # p(lnM) + + axst[0].loglog(bin_centres, mf_out / exp_mf, color="r", linewidth=2) + axst[0].loglog( + bin_centres, np.ones_like(exp_mf), color="r", linestyle=":", linewidth=1 + ) + + axs[0].loglog(bin_centres, mf_out, color="k", linewidth=2, label="Sample") + axs[0].loglog( + bin_centres, exp_mf, color="k", linestyle=":", linewidth=1, label="Expected" + ) - hist, _ = np.histogram(halomass_out, edges) - mf_out = hist / volume_total_m[i] / dlnm - binned_cmf = binned_cmf * n_cond / volume_total_m[i] / dlnm * cond_mass_m + axs[1].semilogx(bin_centres, p_m / p_m.max(), color="k", linewidth=2) + axs[1].axvline(exp_M[0], color="k", linestyle=":", linewidth=2) + axs[1].axvline(M_array.mean(), color="k", linestyle="-", linewidth=1) + axs[1].set_xlim(M_array.min() / 10, M_array.max() * 10) - np.testing.assert_allclose(N_out.mean(), exp_N[0], rtol=RELATIVE_TOLERANCE) - np.testing.assert_allclose(M_out.mean(), exp_M[0], rtol=RELATIVE_TOLERANCE) - np.testing.assert_allclose(mf_out, binned_cmf, rtol=RELATIVE_TOLERANCE) + axst[1].plot(centres_n, p_n / p_n.max(), color="r", linewidth=2) + axst[1].axvline(exp_N[0], color="r", linestyle=":", linewidth=2) + axst[1].axvline(N_array.mean(), color="r", linestyle="-", linewidth=1) + axst[1].set_xlim(N_array.min() - 1, N_array.max() + 1)