brinckmann / montepython_public

Public repository for the Monte Python Code
MIT License
93 stars 77 forks source link

Math domain error in 'compute_lkl' #312

Closed tyxie2003 closed 1 year ago

tyxie2003 commented 1 year ago

Hello everyone,

I am trying to run montepython(3.5.0) on mock data, with modification on the form of tensor power spectrum in class(3.0.2). Then I encounter the following error which cause the chain to stop immediately:

Traceback (most recent call last):
  File "/exports/home/tianyi/data/mpp/montepython/MontePython.py", line 40, in <module>
    sys.exit(run())
  File "/data/tianyi/mpp/montepython/run.py", line 45, in run
    sampler.run(cosmo, data, command_line)
  File "/data/tianyi/mpp/montepython/sampler.py", line 46, in run
    mcmc.chain(cosmo, data, command_line)
  File "/data/tianyi/mpp/montepython/mcmc.py", line 787, in chain
    newloglike = sampler.compute_lkl(cosmo, data)
  File "/data/tianyi/mpp/montepython/sampler.py", line 776, in compute_lkl
    value = likelihood.loglkl(cosmo, data)
  File "/data/tianyi/mpp/montepython/likelihood_class.py", line 1467, in loglkl
    lkl = self.compute_lkl(cl, cosmo, data)
  File "/data/tianyi/mpp/montepython/likelihood_class.py", line 1631, in compute_lkl
    (det_mix/det_the + math.log(det_the/det_obs) - num_modes)
ValueError: math domain error

After seeking into the file likelihood_class.py I find this is caused by negative value of the determine of the theoretical covariant matrix (Cov_the) defined in line 1554 (with B mode True and delensing False). As I output the problematic place, it seems that cl['tt'] or cl['bb'] went negative, but I don't know how exactly it becomes that. (Even the defination of cl['tt'] is not find by me in likelihood_class.py)

Usually the same fault appears when I didn't modify the tensor power spectrum in a correct manner. However, this time I second checked it to be reasonable, and tried a set of parameter and get normal results. But as I move one of the parameter, then generate another fiducial file, then do MCMC, the chain stop due to the above error. The fiducial file is generated successfully, the power spectrum can be calculate correctly using CLASS, but as it start sampling, the error appears to every chain and shut down them quickly.

My modification is about setting a peak on the power spectrum. The moved parameter is the position of the peak.

I'll appreciate it a lot if anyone could help me resolve this problem. Thanks in advance.

Best regards Reverseinfinity

tyxie2003 commented 1 year ago

After a few days' seeking and thinking I have the following a few conclusions.

And I think my problem is solved by this last step. I will close the issue soon, if I don't find more problems.

mtagliazucchi commented 1 year ago

Hi, I had exactly the same problem (but with a different implementation). When this error occurs, Montepython stops completely the chain . Is it possible to modify it so that it discards the set of parameters that caused negative Cl's and keep updating the chain with new proposed parameters?

Thank you so much.

tyxie2003 commented 1 year ago

Hi,

The problem I met is due to the interpolation process in CLASS. My modification has made Cl change fast w.r.t. $\ell$, so the spline interpolation gave oscillating features and some value became negative. I tried several ways to fix it like judge the negative point and modify it, but the most efficient and simplest of them is to change the spline interpolation to linear interpolation, and tighten the step length of calculating Cl.

Montepython use the function harmonic_cl_at_l or lensing_cl_at_l in CLASS. These functions interpolation is due to the function array_interpolate_spline. I changed the last function to be linear instead.

As to my circumstance, before I change the interpolation the chain stops almost the moment I start it, so I don't think discarding points is a good choice, because otherwise there's hardly points left. The negative value comes from bad interpolation and is avoidable. However, the circumstance may be different for you and I hope these can help you.