GSL error for new parameterisation and USE_TS_FLUCT

[steven-murray]

This is mostly from @catherinewatkinson and Tom Binnie. This is an MWE script:

from py21cmmc import mcmc
import logging

logging.getLogger("21CMMC").setLevel(logging.DEBUG)

# ======== USER ADJUSTABLE VARIABLES 
THREADS = 8
WALK_RATIO = 2
ITER = 1
OLD_PARAMETERISATION = False
# ===================================

core = mcmc.CoreCoevalModule(
    redshift=[9],
    user_params=dict(HII_DIM=50, BOX_LEN=125.0),
    flag_options=dict(USE_MASS_DEPENDENT_ZETA=not OLD_PARAMETERISATION),
    do_spin_temp=True,
    regenerate=True,   # ensure each run is started fresh
    initial_conditions_seed=1234  # ensure each run is exactly the same.
)

# Now the likelihood...
datafiles = ["data/power_mcmc_data_%s.npz" % z for z in core.redshift]
power_spec = mcmc.Likelihood1DPowerCoeval(
    datafile=datafiles,
    noisefile=None,
    logk=False, min_k=0.1, max_k=1.0,
    simulate=True
)

# OLD PARAMATRISATION - for code development and proof of concepts
if OLD_PARAMETERISATION:
    params = dict(
        HII_EFF_FACTOR=[30.0, 10.0, 50.0, 3.0],
        ION_Tvir_MIN=[4.7, 4, 6, 0.1])
else:
    params = dict(
        F_STAR10=[-1.301029996, -3, 0, 0.1],
        ALPHA_STAR=[0.5, -0.5, 1, 0.05],
        F_ESC10=[-1, -3, 0, 0.1],
        ALPHA_ESC=[-0.5, -1, 0.5, 0.05],
        M_TURN=[8.698970004, 8, 10, 0.1],
        t_STAR=[0.5, 0, 1, 0.05],
        L_X=[40.5, 38, 42, 0.15],
        NU_X_THRESH=[500, 100, 1500, 50]
    )

chain = mcmc.run_mcmc(
    core, power_spec, datadir='data', model_name="power_only",
    params=params,
    walkersRatio=WALK_RATIO, 
    sampleIterations=ITER, threadCount=THREADS, continue_sampling=False
)

Using the develop-steven branch (commit 6954656379661c17eb672bfb007e1160fe8d20f9) and compiling with

LOG_LEVEL=SUPER_DEBUG pip install -e .

then running with python cat_test.py &> output.txt &, seemingly gives different errors on different runs (even for the same seed!). This is the contents of output.txt for one of the runs I did:

INFO:21CMMC:Initializing init and perturb boxes for the entire chain.
INFO:21CMMC:Initialization done.
DEBUG:21CMMC:PID=31884 Updating parameters: 
DEBUG:21CMMC:PID=31884 AstroParams: AstroParams(ALPHA_ESC:-0.5, ALPHA_STAR:0.5, F_ESC10:0.1, F_STAR10:0.05011872336272722, HII_EFF_FACTOR:30.0, ION_Tvir_MIN:49999.9995007974, L_X:1e+40, M_TURN:501187233.6272715, NU_X_THRESH:500.0, N_RSD_STEPS:20, R_BUBBLE_MAX:15.0, X_RAY_SPEC_INDEX:1.0, X_RAY_Tvir_MIN:49999.9995007974, t_STAR:0.5)
INFO:21CMMC:Existing init_boxes found and read in (seed=1234).
INFO:21CMMC:Existing z=9 perturb_field boxes found and read in (seed=1234).
DEBUG:21CMMC:redshifts: [35.225231081424155, 34.51493243276878, 33.81856120859684, 33.135844322153766, 32.46651404132722, 31.81030788365414, 31.166968513386408, 30.53624364057491, 29.917885922132264, 29.31165286483555, 28.71730673023093, 28.134614441402874, 27.563347491571445, 27.00328185448181, 26.45419789655079, 25.91588029073607, 25.388117932094186, 24.8707038549943, 24.363435151955194, 23.86611289407372, 23.37854205301345, 22.900531424522992, 22.431893553453914, 21.972444660248936, 21.522004568871505, 21.080396636148535, 20.647447682498562, 20.2229879240182, 19.806850905900195, 19.398873437157054, 18.998895526624562, 18.606760320220157, 18.222314039431527, 17.845405921011302, 17.475888157854218, 17.113615841033546, 16.758446902974065, 16.41024206173928, 16.068864766411057, 15.734181143540251, 15.406059944647303, 15.084372494752259, 14.768992641913979, 14.459796707758802, 14.156663438979217, 13.859473959783546, 13.568111725277987, 13.282462475762733, 13.002414191924247, 12.727857050906124, 12.458683383241297, 12.194787630628722, 11.936066304537963, 11.682417945625454, 11.433743083946524, 11.189944199947572, 10.95092568622311, 10.716593810022657, 10.4868566764928, 10.261624192640001, 10.040808032000001, 9.824321600000001, 9.61208, 9.404, 9.2, 9]
DEBUG:21CMMC:PID=31884 doing spin temp for z=35.225231081424155
DEBUG:21CMMC:In spin temp, redshift is 35.225231081424155 and prev_redshift is 0
2019-02-06 11:34:38 | DEBUG   | SpinTemperatureBox.c | ComputeTsBox:37 [pid=31884] | input values:
2019-02-06 11:34:38 | DEBUG   | SpinTemperatureBox.c | ComputeTsBox:38 [pid=31884] | redshift=35.225231, prev_redshift=0.000000 perturbed_field_redshift=35.225231
2019-02-06 11:34:38 | INFO    | UsefulFunctions.c | writeAstroParams:548 [pid=31884] | AstroParams: [HII_EFF_FACTOR=30.000000, ALPHA_STAR=0.500000, F_ESC10=0.100000, ALPHA_ESC=-0.500000, M_TURN=501187232.000000, R_BUBBLE_MAX=15.000000, L_X=10000000000000000303786028427003666890752.000000, NU_X_THRESH=500.000000, X_RAY_SPEC_INDEX=1.000000, F_STAR10=0.050119, F_STAR=0.500000, N_RSD_STEPS=0.000000]
2019-02-06 11:34:38 | SUPER-DEBUG | SpinTemperatureBox.c | ComputeTsBox:107 [pid=31884] | initialized
2019-02-06 11:34:38 | SUPER-DEBUG | SpinTemperatureBox.c | ComputeTsBox:110 [pid=31884] | initalising Ts Interp Arrays
2019-02-06 11:34:38 | SUPER-DEBUG | SpinTemperatureBox.c | ComputeTsBox:240 [pid=31884] | initalised Ts Interp Arrays
2019-02-06 11:34:38 | SUPER-DEBUG | SpinTemperatureBox.c | ComputeTsBox:275 [pid=31884] | Initialised PS
2019-02-06 11:34:38 | SUPER-DEBUG | SpinTemperatureBox.c | ComputeTsBox:286 [pid=31884] | About to initialise heat
2019-02-06 11:34:38 | SUPER-DEBUG | SpinTemperatureBox.c | ComputeTsBox:288 [pid=31884] | Initialised heat
2019-02-06 11:34:39 | SUPER-DEBUG | SpinTemperatureBox.c | ComputeTsBox:293 [pid=31884] | Initialised sigmaM interp table
2019-02-06 11:34:39 | SUPER-DEBUG | SpinTemperatureBox.c | ComputeTsBox:298 [pid=31884] | redshift greater than Z_HEAT_MAX
gsl: qag.c:261: ERROR: could not integrate function
Default GSL error handler invoked.

[steven-murray]

@BradGreig the issue in this case seems to be from the following lines in SpinTemperatureBox.c:

        if(!flag_options->M_MIN_in_Mass) {
            M_MIN = TtoM(redshift, astro_params->X_RAY_Tvir_MIN, mu_for_Ts);
            initialiseSigmaMInterpTable(M_MIN,1e20);
        }

When printing out M_MIN, it sometimes comes out as a VERY big number, and sometimes as nan.

[BradGreig]

@steven-murray actually, the issue is in the MWE not the code itself. For flag_options, you are passing a dictionary rather than actually using the FlagOptions class. As a result, M_MIN_in_Mass is being set to False, rather than True which it should automatically switch to when USE_MASS_DEPENDENT_ZETA = True is passed.

Alternatively, if you are going to stick to using a dict(), you must set M_MIN_in_Mass = True when passing the dict() directly to the function.

[steven-murray]

Ah, this is almost right. Actually, the code takes that dict and makes a FlagOptions class out of it, which should automatically set the M_MIN_in_Mass. However, there was a bug preventing this from happening. This is fixed with 1513fbbf9f66550bf008e5d66f62ada8219055be.