draws more than 1 confidence interval

LSSTDESC · Jul 3, 2024 · 64d2f4a · 64d2f4a
1 parent 3512609
commit 64d2f4a
Show file tree

Hide file tree

Showing 2 changed files with 50 additions and 17 deletions.
diff --git a/augur/postprocess.py b/augur/postprocess.py
@@ -46,7 +46,9 @@ def postprocess(config):
     ls = pconfig["linestyle"] if "linestyle" in pconfig.keys() else "-"
     edgecolors = pconfig["linecolor"] if "linecolor" in pconfig.keys() else "k"
     facecolors = pconfig["facecolor"] if "facecolor" in pconfig.keys() else "none"
-    CL = pconfig["CL"] if "CL" in pconfig.keys() else 0.68
+    CL = pconfig["CL"] if "CL" in pconfig else [0.68,]
+    if not isinstance(CL, list):
+        CL = [CL,]
     f, ax = plt.subplots(npars, npars, figsize=(xsize, ysize))
 
     try:
@@ -64,9 +66,12 @@ def postprocess(config):
             inv_fisher[0, 1] = inv_cache[i, j]
             inv_fisher[1, 0] = inv_cache[j, i]
             inv_fisher[1, 1] = inv_cache[j, j]
-            sig0, sig1 = draw_fisher_ellipses(ax[j, i], inv_fisher, facecolors,
-                                              edgecolors, ls, lw,
-                                              mu=(centers[i], centers[j]), CL=CL)
+            alpha_count = 1.0
+            for cl in CL:
+                sig0, sig1 = draw_fisher_ellipses(ax[j, i], inv_fisher, facecolors,
+                                                  edgecolors, ls, lw,
+                                                  mu=(centers[i], centers[j]), CL=cl, alpha=alpha_count)
+                alpha_count -= 0.4
             # Create the 1D histogram for i, i
             if j == i+1:
                 xarr = np.linspace(centers[i]-5*sig0,
@@ -103,9 +108,13 @@ def postprocess(config):
             inv_fisher[0, 1] = inv_cache[ii, jj]
             inv_fisher[1, 0] = inv_cache[jj, ii]
             inv_fisher[1, 1] = inv_cache[jj, jj]
-            sig0, sig1 = draw_fisher_ellipses(ax, inv_fisher, facecolors,
-                                              edgecolors, ls, lw,
-                                              mu=(centers[ii], centers[jj]), CL=CL)
+            alpha_count = 1.0
+            for cl in CL:
+                sig0, sig1 = draw_fisher_ellipses(ax, inv_fisher, facecolors,
+                                                  edgecolors, ls, lw,
+                                                  mu=(centers[ii], centers[jj]), CL=cl, 
+                                                  alpha=alpha_count)
+                alpha_count -= 0.4
             ax.set_xlim(-sig0+centers[ii], sig0+centers[ii])
             ax.set_ylim(-sig1+centers[jj], sig1+centers[jj])
             ax.set_xlabel(pair0)
@@ -114,19 +123,40 @@ def postprocess(config):
             f.savefig(os.path.join(outdir, f"{pair0}--{pair1}.pdf"))
 
     # w0 -- wa plots are always made
+    # iw = np.where(keys == "w0")[0][0]
+    # iwa = np.where(keys == "wa")[0][0]
+    # sig_w0 = np.sqrt(inv_cache[iw, iw])
+    # sig_wa = np.sqrt(inv_cache[iwa, iwa])
+    # FOM, FOM2 = get_FoM_all(fisher, iw, iwa, CL)
+    # fisher_table = astropy.table.Table([[CL], [FOM], [FOM2], [sig_w0], [sig_wa]],
+    #                                    names=("CL", "FoM", "FoM (alt.)",
+    #                                           "sigma_w0", "sigma_wa"))
+    # fisher_table.write(pconfig["latex_table"], format="latex")
+    # Initialize table with column names only
+    column_names = ["CL", "FoM", "FoM (alt.)", "sigma_w0", "sigma_wa"]
+    fisher_table = astropy.table.Table(names=column_names)
+
+    # Find indices for w0 and wa
     iw = np.where(keys == "w0")[0][0]
     iwa = np.where(keys == "wa")[0][0]
-    sig_w0 = np.sqrt(inv_cache[iw, iw])
-    sig_wa = np.sqrt(inv_cache[iwa, iwa])
-    FOM, FOM2 = get_FoM_all(fisher, iw, iwa, CL)
-    fisher_table = astropy.table.Table([[CL], [FOM], [FOM2], [sig_w0], [sig_wa]],
-                                       names=("CL", "FoM", "FoM (alt.)",
-                                              "sigma_w0", "sigma_wa"))
-    fisher_table.write(pconfig["latex_table"], format="latex")
+
+    # Temporary list to hold row data
+    row_data = []
+
+    # Iterate over each CL value
+    for cl in CL:
+        sig_w0 = np.sqrt(inv_cache[iw, iw])
+        sig_wa = np.sqrt(inv_cache[iwa, iwa])
+        FOM, FOM2 = get_FoM_all(fisher, iw, iwa, cl)
+
+        fisher_table.add_row([cl, FOM, FOM2, sig_w0, sig_wa])
+
+    # Write the table to a LaTeX file
+    fisher_table.write(pconfig["latex_table"], format="latex", overwrite=True)
 
 
 def draw_fisher_ellipses(ax, inv_F, facecolors, edgecolors, linestyles, linewidth,
-                         mu=[0, 0], CL=0.95):
+                         mu=[0, 0], CL=0.95, alpha=1.0):
     """
     Draw uncertainty ellipses for a set of Fisher matrices.
 
@@ -172,7 +202,7 @@ def draw_fisher_ellipses(ax, inv_F, facecolors, edgecolors, linestyles, linewidt
     ellipses = EllipseCollection(
         widths, heights, angles, units="xy", offsets=mu,
         transOffset=ax.transData, facecolors=facecolors,
-        edgecolors=edgecolors, linestyles=linestyles, linewidth=linewidth)
+        edgecolors=edgecolors, linestyles=linestyles, linewidth=linewidth, alpha=alpha)
     ax.add_collection(ellipses)
 
     sig0 = np.sqrt(C[0, 0])

diff --git a/examples/config_test.yml b/examples/config_test.yml
@@ -112,6 +112,9 @@ fisher:
 postprocess:
   latex_table: output/latex_table.tex
   triangle_plot: output/triangle_plot.pdf
+  facecolor: blue
   pairplots: [(w0, wa), (omega_c, sigma8)]
   outdir: 'output/'
-  CL: 0.68
+  CL: 
+    - 0.68
+    - 0.95