Check for NaN and inf when calculating the goodness-of-fit statistic …

…in CompareSpectra
tomasstolker · Jun 26, 2024 · 16206ae · 16206ae
1 parent c680303
commit 16206ae
Show file tree

Hide file tree

Showing 5 changed files with 21 additions and 13 deletions.
diff --git a/species/data/companion_data/companion_data.json b/species/data/companion_data/companion_data.json
@@ -712,7 +712,6 @@
             "Chauvin et al. 2004",
             "Gaia Early Data Release 3",
             "Mohanty et al. 2007",
-            "Mohanty et al. 200z",
             "Patience et al. 2010",
             "Song et al. 2006"
         ]

diff --git a/species/data/model_data/model_spectra.py b/species/data/model_data/model_spectra.py
@@ -192,8 +192,6 @@ def add_model_grid(
             flush=True,
         )
 
-        print(data_folder)
-        print(data_file)
         extract_tarfile(str(data_file), str(data_folder), member_list=member_list)
 
         print(" [DONE]")

diff --git a/species/fit/compare_spectra.py b/species/fit/compare_spectra.py
@@ -266,12 +266,15 @@ def spectral_type(
                                     g_k = 0.0
                                     c_k_spec = []
 
-                                c_k = np.sum(c_numer) / np.sum(c_denom)
+                                idx_select = np.isfinite(c_numer) & np.isfinite(c_denom)
+
+                                c_k = np.sum(c_numer[idx_select]) / np.sum(c_denom[idx_select])
                                 c_k_spec.append(c_k)
 
                                 chi_sq = (
-                                    spec_item[indices, 1] - c_k * flux_resample
-                                ) / spec_item[indices, 2]
+                                    spec_item[indices, 1][idx_select] - c_k * flux_resample[idx_select]
+                                ) / spec_item[indices, 2][idx_select]
+                                print(chi_sq)
 
                                 g_k += np.sum(w_i * chi_sq**2)
 

diff --git a/species/util/fit_util.py b/species/util/fit_util.py
@@ -303,7 +303,7 @@ def get_residuals(
 
         n_param = dset.attrs["n_param"]
 
-        if "fixed_param" in hdf5_file[f"results/fit/{tag}/"]:
+        if "n_fixed" in dset.attrs:
             n_fixed = dset.attrs["n_fixed"]
 
         else:

diff --git a/species/util/plot_util.py b/species/util/plot_util.py
@@ -325,12 +325,12 @@ def update_labels(param: List[str], object_type: str = "planet") -> List[str]:
         index = param.index("fsed")
         param[index] = r"$f_\mathrm{sed}$"
 
-    if "fsed_1" in param:
-        index = param.index("fsed_1")
+    if "fsed_0" in param:
+        index = param.index("fsed_0")
         param[index] = r"$f_\mathrm{sed,1}$"
 
-    if "fsed_2" in param:
-        index = param.index("fsed_2")
+    if "fsed_1" in param:
+        index = param.index("fsed_1")
         param[index] = r"$f_\mathrm{sed,2}$"
 
     if "f_clouds" in param:
@@ -1078,14 +1078,19 @@ def quantity_unit(
 
             elif item_split[0] == "logg":
                 quantity.append(f"logg_{param_index}")
-                unit.append("")
+                unit.append(None)
                 label.append(rf"$\log g_\mathrm{{{param_index+1}}}$")
 
             elif item_split[0] == "feh":
                 quantity.append(f"feh_{param_index}")
-                unit.append("")
+                unit.append(None)
                 label.append(rf"[Fe/H]$_\mathrm{{{param_index+1}}}$")
 
+            elif item_split[0] == "fsed":
+                quantity.append(f"fsed_{param_index}")
+                unit.append(None)
+                label.append(rf"$f_\mathrm{{sed,{param_index+1}}}$")
+
             elif item_split[0] == "radius":
                 quantity.append(f"radius_{param_index}")
 
@@ -1556,6 +1561,9 @@ def create_model_label(
         elif param_item[:3] == "feh":
             value = f"{model_param[param_item]:{param_fmt['feh']}}"
 
+        elif param_item[:4] == "fsed":
+            value = f"{model_param[param_item]:{param_fmt['fsed']}}"
+
         elif param_item[:6] == "radius":
             if object_type == "planet":
                 value = f"{model_param[param_item]:{param_fmt['radius']}}"