remove kwarg stability_crit from hist_classified_stable_vs_hull_dist() and cumulative_clf_metric()

Author: janosh

matbench_discovery/plots.py

@@ -1,6 +1,6 @@
 from __future__ import annotations
 
-from typing import Any, Literal, get_args
+from typing import Any, Literal
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -15,7 +15,6 @@
 __author__ = "Janosh Riebesell"
 __date__ = "2022-08-05"
 
-StabilityCriterion = Literal["energy", "energy+std", "energy-std"]
 WhichEnergy = Literal["true", "pred"]
 AxLine = Literal["x", "y", "xy", ""]
 
@@ -72,14 +71,12 @@
 def hist_classified_stable_vs_hull_dist(
     e_above_hull_pred: pd.Series,
     e_above_hull_true: pd.Series,
-    std_pred: pd.Series = None,
     ax: plt.Axes = None,
     which_energy: WhichEnergy = "true",
-    stability_crit: StabilityCriterion = "energy",
     stability_threshold: float = 0,
     show_threshold: bool = True,
     x_lim: tuple[float | None, float | None] = (-0.4, 0.4),
-    rolling_accuracy: float = 0.02,
+    rolling_accuracy: float | None = 0.02,
 ) -> tuple[plt.Axes, dict[str, float]]:
     """
     Histogram of the energy difference (either according to DFT ground truth [default]
@@ -98,21 +95,16 @@ def hist_classified_stable_vs_hull_dist(
             energy.
         e_above_hull_true (pd.Series): energy diff to convex hull according to DFT
             ground truth.
-        std_pred (pd.Series, optional): standard deviation of the model's predicted
-            formation energy.
         ax (plt.Axes, optional): matplotlib axes to plot on.
         which_energy (WhichEnergy, optional): Whether to use the true formation energy
             or the model's predicted formation energy for the histogram.
-        stability_crit (StabilityCriterion, optional): Whether to add/subtract the
-            model's predicted uncertainty from its energy prediction when measuring
-            predicted stability.
         stability_threshold (float, optional): set stability threshold as distance to
             convex hull in eV/atom, usually 0 or 0.1 eV.
         show_threshold (bool, optional): Whether to plot stability threshold as dashed
             vertical line.
         x_lim (tuple[float | None, float | None]): x-axis limits.
-        rolling_accuracy (float): Rolling accuracy window size in eV / atom. Set to 0 to
-            disable. Defaults to 0.01.
+        rolling_accuracy (float): Rolling accuracy window size in eV / atom. Set to None
+            or 0 to disable. Defaults to 0.01.
 
     Returns:
         tuple[plt.Axes, dict[str, float]]: plot axes and classification metrics
@@ -122,17 +114,7 @@ def hist_classified_stable_vs_hull_dist(
     """
     ax = ax or plt.gca()
 
-    if stability_crit not in get_args(StabilityCriterion):
-        raise ValueError(
-            f"Invalid {stability_crit=} must be one of {get_args(StabilityCriterion)}"
-        )
-
     test = e_above_hull_pred + e_above_hull_true
-    if stability_crit == "energy+std":
-        test += std_pred
-    elif stability_crit == "energy-std":
-        test -= std_pred
-
     # --- histogram of DFT-computed distance to convex hull
     if which_energy == "true":
         actual_pos = e_above_hull_true <= stability_threshold
@@ -348,8 +330,6 @@ def cumulative_clf_metric(
     e_above_hull_error: pd.Series,
     e_above_hull_true: pd.Series,
     metric: Literal["precision", "recall"],
-    std_pred: pd.Series = None,
-    stability_crit: StabilityCriterion = "energy",
     stability_threshold: float = 0,  # set stability threshold as distance to convex
     # hull in eV / atom, usually 0 or 0.1 eV
     ax: plt.Axes = None,
@@ -370,9 +350,6 @@ def cumulative_clf_metric(
         e_above_hull_true (str, optional): Column name with convex hull distance values.
             Defaults to "e_above_hull".
         metric ('precision' | 'recall', optional): Metric to plot.
-        stability_crit ('energy' | 'energy+std' | 'energy-std', optional): Whether to
-            use energy+/-std as stability stability_crit where std is the model
-            predicted uncertainty for the energy it stipulated. Defaults to "energy".
         stability_threshold (float, optional): Max distance from convex hull before
             material is considered unstable. Defaults to 0.
         label (str, optional): Model name used to identify its liens in the legend.
@@ -391,15 +368,6 @@ def cumulative_clf_metric(
     e_above_hull_error = e_above_hull_error.sort_values()
     e_above_hull_true = e_above_hull_true.loc[e_above_hull_error.index]
 
-    if stability_crit not in get_args(StabilityCriterion):
-        raise ValueError(
-            f"Invalid {stability_crit=} must be one of {get_args(StabilityCriterion)}"
-        )
-    if stability_crit == "energy+std":
-        e_above_hull_error += std_pred
-    elif stability_crit == "energy-std":
-        e_above_hull_error -= std_pred
-
     true_pos_mask = (e_above_hull_true <= stability_threshold) & (
         e_above_hull_error <= stability_threshold
     )

scripts/hist_classified_stable_vs_hull_dist.py

@@ -1,13 +1,7 @@
 # %%
-import pandas as pd
-
-from matbench_discovery import ROOT, today
-from matbench_discovery.load_preds import df_wbm
-from matbench_discovery.plots import (
-    StabilityCriterion,
-    WhichEnergy,
-    hist_classified_stable_vs_hull_dist,
-)
+from matbench_discovery import today
+from matbench_discovery.load_preds import load_df_wbm_with_preds
+from matbench_discovery.plots import WhichEnergy, hist_classified_stable_vs_hull_dist
 
 __author__ = "Rhys Goodall, Janosh Riebesell"
 __date__ = "2022-06-18"
@@ -25,56 +19,43 @@
 
 
 # %%
-df = pd.read_csv(
-    # f"{ROOT}/data/2022-06-11-from-rhys/wren-mp-initial-structures.csv"
-    f"{ROOT}/models/wrenformer/2022-11-15-wrenformer-IS2RE-preds.csv"
-).set_index("material_id")
-
-df["e_above_hull"] = df_wbm.e_above_hull_mp2020_corrected_ppd_mp
+model_name = "Wrenformer"
+df_wbm = load_df_wbm_with_preds(models=[model_name]).round(3)
 
 
 # %%
-nan_counts = df.isna().sum()
-assert all(nan_counts == 0), f"df should not have missing values: {nan_counts}"
-
-# target_col = "e_form_target"
-target_col = "e_form_per_atom"
-stability_crit: StabilityCriterion = "energy"
+target_col = "e_form_per_atom_mp2020_corrected"
 which_energy: WhichEnergy = "true"
-
-if "std" in stability_crit:
-    # TODO column names to compute standard deviation from are currently hardcoded
-    # needs to be updated when adding non-aviary models with uncertainty estimation
-    var_aleatoric = (df.filter(like="_ale_") ** 2).mean(axis=1)
-    var_epistemic = df.filter(regex=r"_pred_\d").var(axis=1, ddof=0)
-    std_total = (var_epistemic + var_aleatoric) ** 0.5
-else:
-    std_total = None
-
-# make sure we average the expected number of ensemble member predictions
-pred_cols = df.filter(regex=r"_pred_\d").columns
-assert len(pred_cols) == 10
+# std_factor=0,+/-1,+/-2,... changes the criterion for material stability to
+# energy+std_factor*std. energy+std means predicted energy plus the model's uncertainty
+# in the prediction have to be on or below the convex hull to be considered stable. This
+# reduces the false positive rate, but increases the false negative rate. Vice versa for
+# energy-std. energy+std should be used for cautious exploration, energy-std for
+# exhaustive exploration.
+std_factor = 0
+
+# TODO column names to compute standard deviation from are currently hardcoded
+# needs to be updated when adding non-aviary models with uncertainty estimation
+var_aleatoric = (df_wbm.filter(like="_ale_") ** 2).mean(axis=1)
+var_epistemic = df_wbm.filter(regex=r"_pred_\d").var(axis=1, ddof=0)
+std_total = (var_epistemic + var_aleatoric) ** 0.5
+std_total = df_wbm[f"{model_name}_std"]
 
 ax, metrics = hist_classified_stable_vs_hull_dist(
-    e_above_hull_pred=df[pred_cols].mean(axis=1) - df[target_col],
-    e_above_hull_true=df.e_above_hull,
+    e_above_hull_pred=df_wbm[model_name] - std_factor * std_total - df_wbm[target_col],
+    e_above_hull_true=df_wbm.e_above_hull_mp2020_corrected_ppd_mp,
     which_energy=which_energy,
-    stability_crit=stability_crit,
-    std_pred=std_total,
     # stability_threshold=-0.05,
-    # rolling_accuracy=0,
+    rolling_accuracy=0,
 )
 
 fig = ax.figure
 fig.set_size_inches(10, 9)
 
-ax.legend(
-    loc="center left",
-    frameon=False,
-    title=f"Enrichment Factor = {metrics['enrichment']:.3}",
-)
+legend_title = f"Enrichment Factor = {metrics['enrichment']:.3}"
+ax.legend(loc="center left", frameon=False, title=legend_title)
 
 
 # %%
-fig_name = f"{today}-wren-wbm-hull-dist-hist-{which_energy=}-{stability_crit=}"
+fig_name = f"{today}-wren-wbm-hull-dist-hist-{which_energy=}"
 # fig.savefig(f"{ROOT}/figures/{fig_name}.pdf")

scripts/hist_classified_stable_vs_hull_dist_batches.py

@@ -2,7 +2,6 @@
 from matbench_discovery import ROOT, today
 from matbench_discovery.load_preds import load_df_wbm_with_preds
 from matbench_discovery.plots import (
-    StabilityCriterion,
     WhichEnergy,
     hist_classified_stable_vs_hull_dist,
     plt,
@@ -31,7 +30,6 @@
 
 # %%
 which_energy: WhichEnergy = "true"
-stability_crit: StabilityCriterion = "energy"
 fig, axs = plt.subplots(2, 3, figsize=(18, 9))
 
 model_name = "Wrenformer"
@@ -44,7 +42,6 @@
         e_above_hull_pred=batch_df[model_name] - batch_df[target_col],
         e_above_hull_true=batch_df[e_above_hull_col],
         which_energy=which_energy,
-        stability_crit=stability_crit,
         ax=ax,
     )
 
@@ -59,7 +56,6 @@
     e_above_hull_pred=df_wbm[model_name] - df_wbm[target_col],
     e_above_hull_true=df_wbm[e_above_hull_col],
     which_energy=which_energy,
-    stability_crit=stability_crit,
     ax=axs.flat[-1],
 )
 

scripts/metrics_table.py

@@ -135,7 +135,6 @@
     cmap="viridis",
     # gmap=np.log10(df_table) # for log scaled color map
 )
-df_styled
 
 
 # %%
@@ -145,4 +144,5 @@
 }
 df_styled.set_table_styles([dict(selector=sel, props=styles[sel]) for sel in styles])
 
-df_styled.to_html(f"{ROOT}/figures/{today}-metrics-table.html")
+html_path = f"{ROOT}/figures/{today}-metrics-table.html"
+# df_styled.to_html(html_path)

scripts/precision_recall.py

@@ -3,7 +3,7 @@
 
 from matbench_discovery import ROOT, today
 from matbench_discovery.load_preds import load_df_wbm_with_preds
-from matbench_discovery.plots import StabilityCriterion, cumulative_clf_metric, plt
+from matbench_discovery.plots import cumulative_clf_metric, plt
 
 __author__ = "Rhys Goodall, Janosh Riebesell"
 
@@ -18,10 +18,6 @@
 
 target_col = "e_form_per_atom_mp2020_corrected"
 e_above_hull_col = "e_above_hull_mp2020_corrected_ppd_mp"
-
-
-# %%
-stability_crit: StabilityCriterion = "energy"
 colors = "tab:blue tab:orange teal tab:pink black red turquoise tab:purple".split()
 
 
@@ -41,7 +37,6 @@
         color=color,
         label=f"{model_name}\n{F1=:.3}",
         project_end_point="xy",
-        stability_crit=stability_crit,
         ax=ax_prec,
         metric="precision",
     )
@@ -52,7 +47,6 @@
         color=color,
         label=f"{model_name}\n{F1=:.3}",
         project_end_point="xy",
-        stability_crit=stability_crit,
         ax=ax_recall,
         metric="recall",
     )
@@ -65,7 +59,7 @@
 # x-ticks every 10k materials
 # ax.set(xticks=range(0, int(ax.get_xlim()[1]), 10_000))
 
-fig.suptitle(f"{today} {stability_crit=}")
+fig.suptitle(f"{today} {model_name}")
 xlabel_cumulative = "Materials predicted stable sorted by hull distance"
 fig.text(0.5, -0.08, xlabel_cumulative, ha="center")