21cmfast/21CMMC

CoreLightConeModule build_model_data always the same lightcone

Closed this issue · 1 comments

NOTE : copy from here

Hi,

I am trying to do a likelihood_NeutralFraction MCMC :

  1. CoreLightConeModule
core_LC = mcmc.CoreLightConeModule( 
    redshift = 5.5,                 
    max_redshift = 12.0,             
    user_params = dict(       
        HII_DIM = 50,         
        BOX_LEN = 125.0
    ),
    z_step_factor=1.04,   
    regenerate=False,
    nchunks=4,
    do_spin_temp=True,
    change_seed_every_iter=False,
) 
  1. likelihood
likelihood_NeutralFraction = mcmc.LikelihoodNeutralFraction()  
  1. MCMC:
chain = mcmc.run_mcmc(
    core_LC, 
    likelihood_NeutralFraction, #[likelihood_planck, likelihood_NeutralFraction], #, likelihood_greig],
    datadir='data',          # Directory for all outputs
    model_name=model_name,   # Filename of main chain output
    params=dict(             # Parameter dict as described above.
        F_STAR10   = [ -1.5, -2.5, -0.5, 1.0  ],
        F_ESC10    = [  -1.0, -2.0,  0.0, 1.0  ],
               ),
    
    walkersRatio=4,          # The number of walkers will be walkersRatio*nparams
    burninIterations=0,      # Number of iterations to save as burnin. Recommended to leave as zero.
    sampleIterations=10,  # Number of iterations to sample, per walker.
    threadCount=8,          # Number of processes to use in MCMC (best as a factor of walkersRatio)
    continue_sampling=False, # Whether to contine sampling from previous run *up to* sampleIterations.
) 

The problem is that the return lnL is always 0.
I looked at the core and likelihood module and found that the actual modeled neutral fraction is identical between walkers and iterations:

In CoreLightConeModule I print the params (they do change as expected), but after the computation of the lightcone lightcone.global_xHI is always the same !

    def build_model_data(self, ctx):
        # Update parameters
        astro_params, cosmo_params = self._update_params(ctx.getParams())
        
        print( 'IN BUILD astro_params :', astro_params.F_ESC10, astro_params.ALPHA_ESC )

        print( self.redshift[0], self.max_redshift, self.z_step_factor, self.initial_conditions_seed )
        
        # Call C-code
        lightcone = p21.run_lightcone(
            
            redshift=self.redshift[0],
            max_redshift=self.max_redshift,
            
            astro_params=astro_params, 
            cosmo_params=cosmo_params, 
            
            user_params=self.user_params,
            flag_options=self.flag_options,
            
            z_step_factor=self.z_step_factor,
            regenerate=self.regenerate,
            random_seed=self.initial_conditions_seed,
            write=self.io_options['cache_ionize'],
            direc=self.io_options['cache_dir'],
        )
        
        print( 'IN BUILD xHI :', lightcone.global_xHI )

        ctx.add('lightcone', lightcone)
>>> 
IN BUILD astro_params : -0.422045384281 -0.5
IN BUILD astro_params : -1.15392864038 -0.5
IN BUILD astro_params : -1.29628662729 -0.5
IN BUILD astro_params : -1.9993342796 -0.5
5.5 12.0 1.04 396696065585
5.5 12.0 1.04 396696065585
5.5 12.0 1.04 396696065585
5.5 12.0 1.04 396696065585
IN BUILD astro_params : -1.98190265686 -0.5
IN BUILD astro_params : -0.762196428892 -0.5
5.5 12.0 1.04 396696065585
5.5 12.0 1.04 396696065585
IN BUILD astro_params : -1.06946089891 -0.5
5.5 12.0 1.04 396696065585
IN BUILD astro_params : -0.565129474849 -0.5
5.5 12.0 1.04 396696065585
IN BUILD xHI : [ 0.96435517  0.95072514  0.93239552  0.90946776  0.8801595   0.84389597
  0.80032223  0.74773973  0.68238193  0.60641336  0.51849377  0.42002162
  0.31888756  0.22795989  0.15128654  0.09266903  0.04975107  0.02413185
  0.00956937]IN BUILD xHI : [ 0.96435517  0.95072514  0.93239552  0.90946776  0.8801595   0.84389597
  0.80032223  0.74773973  0.68238193  0.60641336  0.51849377  0.42002162
  0.31888756  0.22795989  0.15128654  0.09266903  0.04975107  0.02413185
  0.00956937]
lnprob :  0
lnprob :  0
IN BUILD xHI : [ 0.96435517  0.95072514  0.93239552  0.90946776  0.8801595   0.84389597
  0.80032223  0.74773973  0.68238193  0.60641336  0.51849377  0.42002162
  0.31888756  0.22795989  0.15128654  0.09266903  0.04975107  0.02413185
  0.00956937]
lnprob :  0
IN BUILD xHI : [ 0.96435517  0.95072514  0.93239552  0.90946776  0.8801595   0.84389597
  0.80032223  0.74773973  0.68238193  0.60641336  0.51849377  0.42002162
  0.31888756  0.22795989  0.15128654  0.09266903  0.04975107  0.02413185
  0.00956937]
lnprob :  0
IN BUILD xHI : [ 0.96435517  0.95072514  0.93239552  0.90946776  0.8801595   0.84389597
  0.80032223  0.74773973  0.68238193  0.60641336  0.51849377  0.42002162
  0.31888756  0.22795989  0.15128654  0.09266903  0.04975107  0.02413185
  0.00956937]

I am stuck here, I do not know why because the actual input params are different and should change the resulting xHI.

I try to switch of regenerate=False or to put the random_seed generator.
I get a warning for the first and I am not sure the second work too.

The defaults flag_options struct has USE_MASS_DEPENDENT_ZETA set to False. This should be set to True if you are using F_STAR and F_ESC.

It solves the problem and produces different light-cone global xHI (for neutral fraction and tau_e MCMC).