Connection between chosen parameters by CosmoMC and execution of CAMB

[Vanessa Smer]

Hi everyone

I write about a question regarding the step between choosing values for the parameters that CosmoMC varies and the connection with CAMB. Does someone knows if there's a file where the program can be forced to give other sets of parameters before sending them to CMB_Cls_simple.f90 to be sent to CAMB to create the power spectrum?

The reason for this is because while using the Quint module, there are some combinations of values that will make it crash, and one way to avoid that is to make the program keep looking for a set of parameters that at least will be accepted by CAMB without interfering with the final parameter space found.

Any suggestions will be much appreciated.

[Jason Dossett]

What part of CAMB is crashing?

I dont know what version of CosmoMC you are using. Search your ini files for

Code: Select all

stop_on_error=

and set that to F.

That may help.

If it is the reionization module that is failing (this should be rarer). There is a workaround for that as well. Still do the same as above, but change reionization.f90 according to the diff below

Code: Select all

diff --git a/camb/reionization.f90 b/camb/reionization.f90
index 3c84981..8987356 100644
--- a/camb/reionization.f90
+++ b/camb/reionization.f90
@@ -1,5 +1,6 @@
 
 module Reionization
+ use Errors  !JD for ISiTGR
  use Precision
  use AMLutils
  implicit none  
@@ -298,16 +299,22 @@ end function Reionization_GetOptDepth
                   try_b = Reion%redshift
        end if
        if &#40;abs&#40;try_b - try_t&#41; < 2e-3/Reionization_AccuracyBoost&#41; exit
-       if &#40;i>100&#41; call mpiStop&#40;'Reionization_zreFromOptDepth&#58; failed to converge'&#41;
+       !JD Modified so that bad TGR parameters won't kill CosmoMC
+       if &#40;i>100&#41; then
+           call GlobalError&#40;'Reionization_zreFromOptDepth&#58; failed to converge',error_reionization&#41;
+           return
+           !call mpiStop&#40;'Reionization_zreFromOptDepth&#58; failed to converge'&#41; 
+       end if
   end do
-  
-  
+   
+   !JD Modified so that bad TGR parameters won't kill CosmoMC
    if &#40;abs&#40;tau - Reion%optical_depth&#41; > 0.002&#41; then
-    write &#40;*,*&#41; 'Reionization_zreFromOptDepth&#58; Did not converge to optical depth'
-    write &#40;*,*&#41; 'tau =',tau, 'optical_depth = ', Reion%optical_depth
-    write &#40;*,*&#41; try_t, try_b
-    call mpiStop&#40;&#41;
-  end if
+    !write &#40;*,*&#41; 'Reionization_zreFromOptDepth&#58; Did not converge to optical depth'
+    !write &#40;*,*&#41; 'tau =',tau, 'optical_depth = ', Reion%optical_depth
+    !write &#40;*,*&#41; try_t, try_b
+    !call mpiStop&#40;&#41;
+    call GlobalError&#40;'Reionization_zreFromOptDepth&#58; Did not converge to optical depth',error_reionization&#41;
+   end if
     
  end subroutine Reionization_zreFromOptDepth 

[Vanessa Smer]

Hi Jason

Thank you very much for replying. The part of CAMB that crashes is in the Quint module, after calculating the variable om1:

if (abs(om1-CP%omegav) > 1d-3) then
!if not, do binary search in the interval
OK=.false.
initial_phidot = astart*Quint_phidot_start(initial_phi2)
om2= Quint_GetOmegaFromInitial(astart,initial_phi2,initial_phidot, atol)
if (om1 > CP%omegav .or. om2 < CP%omegav) then
write (*,*) 'initial phi values must bracket required value '
write (*,*) 'om1, om2 = ', real(om1), real(om2)
stop
end if


There is good reason for it, as this part of the code makes sure the variables om1 and om2 are enclosing CP%omegav. What I am trying to do is to give CAMB values that will calculate a CP%omegav that is always between the range om1 < CP%omegav < om2, for the program not to stop and the chains of CosmoMC to reach convergence.

So, in case the values given by CosmoMC give a "bad" CP%omegav, then ask the program to calculate others, before calling CAMB, to avoid this issue.

Once again, thank you for your answer.

[Jason Dossett]

Hi Vanessa,

I don't have much experience with the Quint module, so dont know when exactly this routine is going to be called, but hopefully what I have below will get you started.

First, as in my previous post you will have to turn stop_on_error to F.

Next in the module where the below routine is called, insert

Code: Select all

use Errors

Then you can modify the code to:

Code: Select all

if (abs(om1-CP%omegav) > 1d−3) then 
!if not, do binary search in the interval 
    OK=.false. 
    initial_phidot = astart*Quint_phidot_start(initial_phi2) 
    om2= Quint_GetOmegaFromInitial(astart,initial_phi2,initial_phidot, atol) 
    if &#40;om1 > CP%omegav .or. om2 < CP%omegav&#41; then 
    write &#40;*,*&#41; 'initial phi values must bracket required value ' 
    write &#40;*,*&#41; 'om1, om2 = ', real&#40;om1&#41;, real&#40;om2&#41; 
    call GlobalError&#40;&#41;
end if 

You may have to insert a return statement so the to stop execution in that particular routine and maybe propagate it out using statements like:

Code: Select all

if(global_error_flag /=0) return

Hope that helps!

[Paul Anthony Fronteri]

Hey Jason.

Just a follow up question about this discussion. How would one go about to make it such that:

The value sent by COSMOMC, to CAMB, that makes it "crash" (then not accepted, or does give unphysical physics, or as in this case, will not give a suitable value for this module) is "send back" to COSMOMC, but then gives a zero likelihood for this value.

Or to put in other word, values that does not work, gives back zero likelihood in COSMOMC, and continues to search for another value, that will/or will not work.

Thanks for any reply.

[Jason Dossett]

Hi Paul,

For the most part the method above should work to do just that. Most of the places CAMB is called in CosmoMC, (mostly in CMB_Cls_simple.f90) have a return if error/=0 from CAMB. These return statements for non-zero error flag usually work to cause LogLike=LogZero and cosmoMC will reject that parameter point.

If you are implementing a new module and making a lot of changes to CosmoMC or CAMB you would just have to check that the error flags are propagated correctly.

You can see an example of this in ISiTGR (http://isit.gr) by looking at my changes to TranferGetMatterPower in camb/modules, GetCls in CMB_Cls_simple, and NonLinear_GetNonLinRatios in camb/halofit.f90. Given my current release is not for the latest version of CosmoMC, but the same principles should apply.