From 5bfb2bf4b33f431e6d0730fd61eac34111f6ae23 Mon Sep 17 00:00:00 2001
From: Mahmoud Osman <mahmoud.ahmedemamosman@lpnhe.in2p3.fr>
Date: Thu, 17 Apr 2025 10:24:06 +0200
Subject: [PATCH] Calling gravity solver params from ini file, plotting slices
 of density field

---
 analysis.py  | 29 +++++++++++++++++++++++++++++
 example0.py  | 18 +++++++++++-------
 ini_file.ini | 18 +++++++++---------
 3 files changed, 49 insertions(+), 16 deletions(-)

diff --git a/analysis.py b/analysis.py
index 4d26a53..5c99f26 100644
--- a/analysis.py
+++ b/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)
+
diff --git a/example0.py b/example0.py
index dea70f8..756d547 100644
--- a/example0.py
+++ b/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'
-    ai = 0.05 # Initial scale factor
-    af = 1.0 # Final scale factor
-    supersampling = 2
-    forcesampling = 2
-    nsteps = 20 # Number of steps in the PM solver
+    which_model = config['gravity']['which_model']
+    ai = float(config['gravity']['ai']) # Initial scale factor
+    af = float(config['gravity']['af']) # Final scale factor
+    supersampling = int(config['gravity']['supersampling'])
+    forcesampling = int(config['gravity']['forcesampling'])
+    nsteps = int(config['gravity']['nsteps']) # Number of steps in the PM solver
 
     chain = borg.forward.ChainForwardModel(box)
 
diff --git a/ini_file.ini b/ini_file.ini
index 8259c0a..1adb910 100644
--- a/ini_file.ini
+++ b/ini_file.ini
@@ -1,9 +1,9 @@
 [system]
 console_output = borg_log
 VERBOSE_LEVEL = 2
-N0 = 32
-N1 = 32
-N2 = 32
+N0 = 64
+N1 = 64
+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
-af = 1.0 # Final scale factor
+ai = 0.05
+af = 1.0
 supersampling = 2
 forcesampling = 2
-nsteps = 20 # Number of steps in the PM solver
\ No newline at end of file
+nsteps = 20