info --want-covmat gives error with modified data.py

[shuishenshui]

Hi all,

I am running an extended neutrino model with free T_ncdm and degree of degeneracy g_f. Sampling {T_ncdm,g_f} are difficult to converge since they are very degenerated, thus I sample {N_eff, T_ncdm} instead, and calculate g_f inside data.py.

The chains runs well and I can plot contours with getdist. The acc rate is 0.28, Total number of accepted steps: 339255. I ran 10 chains, but some chains stop much earlier than N.

However, it always give this error when I try to run info --want-covmat:

--> Computing covariance matrix Traceback (most recent call last): File "montepython/MontePython.py", line 40, in sys.exit(run()) File "/home/moon/wuyunzhi/montepython_public-3.5/montepython/run.py", line 32, in run custom_command) File "/home/moon/wuyunzhi/montepython_public-3.5/montepython/run.py", line 191, in safe_initialisation cosmo, data, command_line, success = initialise(custom_command) File "/home/moon/wuyunzhi/montepython_public-3.5/montepython/initialise.py", line 59, in initialise analyze(command_line) File "/home/moon/wuyunzhi/montepython_public-3.5/montepython/analyze.py", line 151, in analyze enumerate(info.chain[info.sorted_indices[0], 2:])] IndexError: index 38 is out of bounds for axis 0 with size 38

This is my .param:

# Cosmological parameters list

data.parameters['omega_b']      = [  2.2377,   None, None,      0.015, 0.01, 'cosmo']
data.parameters['omega_cdm']    = [ 0.12010,   None, None,     0.0013,    1, 'cosmo']
data.parameters['100*theta_s']  = [ 1.04110,   None, None,    0.00030,    1, 'cosmo']
data.parameters['ln10^{10}A_s'] = [  3.0447,   None, None,      0.015,    1, 'cosmo']
data.parameters['n_s']          = [  0.9659,   None, None,     0.0042,    1, 'cosmo']
data.parameters['tau_reio']     = [  0.0543,  0.004, None,      0.008,    1, 'cosmo']

data.parameters['deg_ncdm'] = [1, 0, None, 0.1, 1, 'cosmo'] # calculated inside data.py, do not occur in log.param
data.parameters['m_ncdm'] = [0.06, 0, 0.8, 0.1, 1, 'cosmo'] 
data.parameters['T_ncdm'] = [0.7, 0, None, 0.1, 1, 'cosmo']
data.parameters['Neff_ncdm'] = [1, 0, 2.0 , 0.1, 1, 'cosmo']  

# Derived parameters

data.parameters['z_reio']          = [1, None, None, 0,     1,   'derived']
data.parameters['Omega_Lambda']    = [1, None, None, 0,     1,   'derived']
data.parameters['YHe']             = [1, None, None, 0,     1,   'derived']
data.parameters['H0']              = [0, None, None, 0,     1,   'derived']
data.parameters['A_s']             = [0, None, None, 0,  1e-9,   'derived']
data.parameters['sigma8']          = [0, None, None, 0,     1,   'derived']
data.parameters['z_rec']          = [0, None, None, 0,     1,   'derived']
data.parameters['ra_rec']          = [0, None, None, 0,     1,   'derived']
data.parameters['Neff']          = [0, None, None, 0,     1,   'derived']
# Other cosmo parameters (fixed parameters, precision parameters, etc.)

data.cosmo_arguments['sBBN file'] = data.path['cosmo']+'/external/bbn/sBBN.dat'
# BBN file path is automatically set to match CLASS version if 'sBBN file' is requested
# You can force the code to use the exact BBN file passed above with flag
#data.custom_bbn_file = True

data.cosmo_arguments['k_pivot'] = 0.05

data.cosmo_arguments['N_ur'] = 2.0328
data.cosmo_arguments['N_ncdm'] = 1

# These two are required to get sigma8 as a derived parameter
# (class must compute the P(k) until sufficient k)
data.cosmo_arguments['output'] = 'mPk'
data.cosmo_arguments['P_k_max_h/Mpc'] = 1.

# The Planck Lensing likelihood is more precise when the non-linear effects are taken
# into consideration. For this you can use halofit (default) or hmcode.
# If you are running an exotic model for which the non-linearities cannot be
# computed with either of these codes, you are advised to comment out the following line.
data.cosmo_arguments['non linear'] = 'halofit'

#nuissance
...

Modified data.py:

      elif elem == 'Neff_ncdm':
          self.cosmo_arguments['deg_ncdm'] = self.cosmo_arguments[elem]/self.cosmo_arguments['T_ncdm']**4 * (4./11.)**(4./3.)
          del self.cosmo_arguments[elem]

Thank you in advance!

[shuishenshui]

update: for {g_f, T_ncdm, m_ncdm} chains, this error also happens when calculating covmat.

[brinckmann]

Hi,

MontePython does not tolerate added comments at the end of input lines, i.e. this will cause problems (the comment needs to be on a separate line): data.parameters['deg_ncdm'] = [1, 0, None, 0.1, 1, 'cosmo'] # calculated inside data.py, do not occur in log.param I would not be very confident that the results are correct either, but maybe you can do a shorter run fixing this problem and checking if they seem to be converging to the same result.

Best, Thejs

[brinckmann]

Also, once this issue is fixed, I suspect you will have problems having both deg_ncdm and Neff_ncdm as varying parameters, should you not only have Neff_ncdm, since it defines deg_ncdm? It seems to me you're effectively sampling the same parameter twice.

Best, Thejs

[shuishenshui]

Hi Thejs,

Sorry for the late response because I ran the chains to verify the problem is solved. I removed the comment in .param and the error's gone and it convergences well now. Thank you so much!

Yes, you are right, the Neff_ncdm and deg_ncdm are degenerated, thus I set Neff_ncdm as input and calculate deg_ncdm in the data.py, which is then put into CLASS as the parameters.

Cheers, Yunzhi

[brinckmann]

Great! Thanks for reporting back.

Best, Thejs