Merge Juno mcmc (#143)

- Bind juno mwr to python
- Can run emcee with juno mwr

---------

Co-authored-by: jihenghu <[email protected]>

.github/workflows/mac.yml

@@ -358,5 +358,5 @@ jobs:
         run: ./combine.py -o test
 
       - name: Test
-        working-directory: ${{github.workspace}}/examples/2024-CLi-juno-mwr
+        working-directory: ${{github.workspace}}/examples/2024-JHu-juno-mwr
         run: echo "cheers"

.github/workflows/main.yml

@@ -705,5 +705,5 @@ jobs:
         run: python3 combine.py -o test
 
       - name: Test
-        working-directory: ${{github.workspace}}/examples/2024-CLi-juno-mwr
+        working-directory: ${{github.workspace}}/examples/2024-JHu-juno-mwr
         run: echo "cheers"

CODEOWNERS

@@ -18,3 +18,4 @@ src/harp @happysky19
 src/opacity @happysky19
 src/nbody @astrophy6
 patches/??.cs_*.patch @cshsgy
+examples/2024-JHu-juno-mwr/ @jihenghu

cmake/examples/juno.cmake

@@ -1,11 +1,5 @@
 # configuration for Jupiter Juno forward and inversion calculation
 
-macro(SET_IF_EMPTY _variable)
-  if("${${_variable}}" STREQUAL "")
-    set(${_variable} ${ARGN})
-  endif()
-endmacro()
-
 # athena variables
 set_if_empty(NUMBER_GHOST_CELLS 2)
 set_if_empty(NVAPOR 2)
@@ -14,11 +8,8 @@ set_if_empty(NVAPOR 2)
 set_if_empty(NCLOUD 4)
 set_if_empty(NTRACER 2)
 
-# canoe task set(TASKLIST InversionTasks)
-
 # canoe configure
 set(HYDROSTATIC ON)
 set(NETCDF ON)
-# set(FITS ON)
 set(RadianceSolver lambert)
 set(PYTHON_BINDINGS ON)

cmake/examples/jupiter_dry.cmake

@@ -5,8 +5,7 @@ set_if_empty(NUMBER_GHOST_CELLS 3)
 
 # canoe configure
 set(PNETCDF ON)
-set(DISORT ON)
-set(PYTHON_BINDINGS ON)
 set(MPI ON)
 set(TASKLIST ImplicitHydroTasks)
 set(RSOLVER lmars)
+# set(DISORT ON) set(PYTHON_BINDINGS ON)

examples/2024-JHu-juno-mwr/inspect_emcee_chain.py

@@ -0,0 +1,67 @@
+# import emcee
+import corner
+
+# %config InlineBackend.figure_format = "retina"
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.stats import norm
+
+# reader = emcee.backends.HDFBackend("run_mcmc.h5")
+
+# tau = reader.get_autocorr_time()
+# burnin = int(2 * np.max(tau))
+# thin = int(0.5 * np.min(tau))
+# samples = reader.get_chain(discard=burnin, flat=True, thin=thin)
+# log_prob_samples = reader.get_log_prob(discard=burnin, flat=True, thin=thin)
+
+import h5py
+
+h5 = h5py.File("run_mcmc_5000.h5", "r")
+chain = h5["mcmc"]["chain"][180:]
+h5.close()
+
+flattened_chain = chain.reshape(-1, 2)
+
+# Create labels for the parameters
+labels = ["qNH3 [ppmv]", "Temperature [K]"]
+
+# Create the corner plot
+fig, ax = plt.subplots(2, 2, figsize=(10, 10))
+corner.hist2d(flattened_chain[:, 0], flattened_chain[:, 1], ax=ax[1, 0])
+ax[1, 0].set_xlabel("qNH3 [ppmv]")
+ax[1, 0].set_ylabel("Temperature [K]")
+ax[1, 0].set_xlim([0, 700])
+
+# Plot histograms for each parameter
+x = np.linspace(1, 700, 700)
+for i in range(2):
+    ax[i, i].hist(flattened_chain[:, i], bins=30, color="blue", alpha=0.7, density=True)
+    ax[i, i].set_xlabel(labels[i])
+    ax[i, i].set_ylabel("PDF")
+    means = np.mean(flattened_chain[:, i])
+    stdev = np.std(flattened_chain[:, i])
+    ax[i, i].axvline(
+        means,
+        color="b",
+        linestyle="--",
+        label=f"posterior: ({means:6.2f}, {stdev:5.2f}",
+    )
+
+    # Plot prior distribution
+    mean, stddev = [(300, 100), (169, 10)][i]
+
+    prior = norm.pdf(x, mean, stddev)
+    ax[i, i].plot(
+        x, prior, color="red", linestyle="--", label=f"Prior: norm({mean}, {stddev})"
+    )
+    ax[i, i].legend()
+    x = np.linspace(131, 200, 70)
+
+ax[0, 0].set_xlim([0, 700])
+
+ax[0, 1].axis("off")  # Disable the axis
+ax[0, 1].set_visible(False)  # Hide the axis
+# Show the plot
+plt.tight_layout()
+# Show the plot
+plt.savefig("emcee_cornerplot_5000.png")

examples/2024-JHu-juno-mwr/inspect_emcee_step.py

@@ -0,0 +1,42 @@
+# import emcee
+
+# %config InlineBackend.figure_format = "retina"
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.stats import norm
+
+# reader = emcee.backends.HDFBackend("run_mcmc.h5")
+
+# tau = reader.get_autocorr_time()
+# burnin = int(2 * np.max(tau))
+# thin = int(0.5 * np.min(tau))
+# samples = reader.get_chain(discard=burnin, flat=True, thin=thin)
+# log_prob_samples = reader.get_log_prob(discard=burnin, flat=True, thin=thin)
+
+import h5py
+
+h5 = h5py.File("run_mcmc_5000.h5", "r")
+chain = h5["mcmc"]["chain"][:]
+h5.close()
+
+flattened_chain = chain.reshape(-1, 2)
+
+# Create labels for the parameters
+labels = ["qNH3 [ppmv]", "Temperature [K]"]
+
+# Create the corner plot
+fig, ax = plt.subplots(2, 1, figsize=(17, 10))
+for iw in range(5):
+    ax[0].plot(range(5000), chain[:, iw, 0])
+
+ax[0].set_ylabel("qNH3 [ppmv]")
+ax[0].set_xlim([0, 5000])
+
+for iw in range(5):
+    ax[1].plot(range(5000), chain[:, iw, 1])
+ax[1].set_ylabel("Temperature [K]")
+ax[1].set_xlim([0, 5000])
+ax[1].set_xlabel("step")
+plt.tight_layout()
+# Show the plot
+plt.savefig("emcee_inspect_step_5000.png")

examples/2024-CLi-juno-mwr/juno_mwr.inp → examples/2024-JHu-juno-mwr/juno_mwr.inp

@@ -65,7 +65,7 @@ PressureScaleHeight = 30.E3
 
 <hydro>
 gamma       = 1.42     # gamma = C_p/C_v
-grav_acc1   = -23.3
+grav_acc1   = -27.01 #-23.3
 sfloor      = 0.
 
 <species>