Calling gravity solver params from ini file, plotting slices of density field

analysis.py

@@ -7,6 +7,8 @@ from contextlib import contextmanager
 import os
 import Pk_library as PKL
 import glob
+import pylab as pl
+pl.ion()
 
 @contextmanager
 def suppress_stdout():
@@ -128,3 +130,30 @@ def get_spectra(ini_name, dirname, mcmc_steps, which_field='BORG_final_density',
         all_r[i,:] = Pk.XPk[:,0,0] / np.sqrt(Pk.Pk[:,0,1] * Pk.Pk[:,0,0])
         
     return k, Pk_true, all_pk, all_r
+
+def plot_dens_fields(mean_field, std_field, true_field, slice):
+    """
+    """
+    sides = mean_field.shape[0]
+    x = np.arange(sides)
+    y = x
+    X, Y = np.meshgrid(x, y)
+    fig, axes = pl.subplots(figsize=(16,7), ncols=3) 
+    _, _, _, im0 = axes[0].hist2d(X.ravel(), Y.ravel(), bins=(sides, sides), weights=true_field[slice].ravel())
+    axes[0].set_xlabel('x')
+    axes[0].set_ylabel('y')
+    axes[0].set_title('Generated mock density')
+    pl.colorbar(im0)
+
+    _, _, _, im1 = axes[1].hist2d(X.ravel(), Y.ravel(), bins=(sides, sides), weights=mean_field[slice].ravel())
+    axes[1].set_xlabel('x')
+    axes[1].set_ylabel('y') 
+    axes[1].set_title('Mean reconstructed density')	
+    pl.colorbar(im1)
+
+    _, _, _, im2 = axes[2].hist2d(X.ravel(), Y.ravel(), bins=(sides, sides), weights=std_field[slice].ravel())
+    axes[2].set_xlabel('x')
+    axes[2].set_ylabel('y') 
+    axes[2].set_title('Std of reconstructed density')	
+    pl.colorbar(im2)
+

example0.py

@@ -238,7 +238,7 @@ class MyLikelihood(borg.likelihood.BaseLikelihood):
 
 
 @borg.registerGravityBuilder
-def build_gravity_model(state: borg.likelihood.MarkovState, box: borg.forward.BoxModel) -> borg.forward.BaseForwardModel:
+def build_gravity_model(state: borg.likelihood.MarkovState, box: borg.forward.BoxModel, ini_fname=None) -> borg.forward.BaseForwardModel:
     """
     Builds the gravity model and returns the forward model chain.
 
@@ -255,13 +255,17 @@ def build_gravity_model(state: borg.likelihood.MarkovState, box: borg.forward.Bo
     global chain
     myprint("Building gravity model")
 
+    if ini_fname is None:
+        ini_fname=borg.getIniConfigurationFilename()
+    config = configparser.ConfigParser()
+    config.read(ini_fname)
     # READ FROM INI FILE
-    which_model = 'lpt'
+    which_model = config['gravity']['which_model']
-    ai = 0.05 # Initial scale factor
+    ai = float(config['gravity']['ai']) # Initial scale factor
-    af = 1.0 # Final scale factor
+    af = float(config['gravity']['af']) # Final scale factor
-    supersampling = 2
+    supersampling = int(config['gravity']['supersampling'])
-    forcesampling = 2
+    forcesampling = int(config['gravity']['forcesampling'])
-    nsteps = 20 # Number of steps in the PM solver
+    nsteps = int(config['gravity']['nsteps']) # Number of steps in the PM solver
 
     chain = borg.forward.ChainForwardModel(box)
 

ini_file.ini

@@ -1,9 +1,9 @@
 [system]
 console_output = borg_log
 VERBOSE_LEVEL = 2
-N0 = 32
+N0 = 64
-N1 = 32
+N1 = 64
-N2 = 32
+N2 = 64
 L0 = 500.0
 L1 = 500.0
 L2 = 500.0
@@ -20,7 +20,7 @@ bias_sampler_blocked= true
 ares_heat = 1.0
 
 [mcmc]
-number_to_generate = 20
+number_to_generate = 50
 random_ic = false
 init_random_scaling = 0.1
 
@@ -31,7 +31,7 @@ max_timesteps = 50
 mixing = 1
 
 [python]
-likelihood_path = /home/bartlett/borg_examples/example0.py
+likelihood_path = /home/mosman/borg_examples/example0.py
 
 [run]
 run_type = mock
@@ -53,12 +53,12 @@ beta = 1.5
 z0 = 0
 
 [mock]
-sigma_dens = 1
+sigma_dens = 1.
 
 [gravity]
 which_model = lpt
-ai = 0.05 # Initial scale factor
+ai = 0.05
-af = 1.0 # Final scale factor
+af = 1.0
 supersampling = 2
 forcesampling = 2
-nsteps = 20 # Number of steps in the PM solver
+nsteps = 20