diff --git a/pls/pls-model-inversion-null-space.png b/pls/pls-model-inversion-null-space.png new file mode 100644 index 0000000..f29a19f Binary files /dev/null and b/pls/pls-model-inversion-null-space.png differ diff --git a/pls/pls-model-inversion-null-space.py b/pls/pls-model-inversion-null-space.py new file mode 100644 index 0000000..1cb483c --- /dev/null +++ b/pls/pls-model-inversion-null-space.py @@ -0,0 +1,154 @@ +"""Generate the committed PNG for the PLS model-inversion / null-space figure. + +Mirrors the worked example in the pid-book page +``latent-variable-modelling/projection-to-latent-structures/pls-model-inversion-and-the-orthogonal-space.rst``: +a two-component PLS model of the cheddar-cheese data (trained on cheeses 5 to +30), inverted toward a target taste of 20.9. The score plot shows the +calibration cheeses, the direct-inversion solution, the null space (the line of +scores that all predict the target taste), and the O-PLS orthogonal space +projected into the PLS score space. The null space and the projected orthogonal +space coincide, which is the point of the page and of Garcia-Carrion et al. +(2025). + +The chapter shows the equivalent Plotly code; this committed figure is the +matplotlib rendering. + +Requires the ``process_improve`` package (``pip install process-improve``) for +``PLS`` and ``OPLS``. + +Usage:: + + python pls/pls-model-inversion-null-space.py [output_dir] + +``output_dir`` defaults to this script's own directory (``pls/``). +""" + +from __future__ import annotations + +import sys +from pathlib import Path + +import matplotlib as mpl + +mpl.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +from process_improve.multivariate import OPLS, PLS + +DATA_URL = "https://openmv.net/file/cheddar-cheese.csv" +X_COLUMNS = ["Acetic", "H2S", "Lactic"] +TARGET_TASTE = 20.9 # the main inversion, drawn as the orange null space +SECOND_TASTE = 47.9 # a second inversion (held-out cheese 4), drawn red dashed +THIRD_TASTE = 12.3 # the lowest held-out taste (cheese 1), drawn purple dotted + +DARK_BLUE = "#1f3d7a" # calibration cheeses +ORANGE = "#e6820a" # the null space (PLS) for the main target +RED = "#d62728" # the null space for the second target +PURPLE = "#7b3fa0" # the null space for the lowest target (a cool contrast) +GREEN = "#2e6f3e" # the orthogonal space (O-PLS), projected into PLS space +BLACK = "#111111" # the direct-inversion solution +plt.rcParams.update( + { + "font.size": 11, + "axes.grid": True, + "grid.alpha": 0.3, + "figure.dpi": 140, + } +) + + +def build_figure(out_dir: Path) -> None: + """Fit the models, invert, and render the score plot.""" + cheese = pd.read_csv(DATA_URL) + train = cheese.iloc[4:] # cheeses 5 to 30 + x_train = train[X_COLUMNS] + y_train = train[["Taste"]] + + pls = PLS(n_components=2).fit(x_train, y_train) + opls = OPLS(n_orthogonal_components=1).fit(x_train, y_train) + + result = pls.invert(y_desired=TARGET_TASTE) + tau = result.scores.to_numpy() + g_ns = result.null_space_basis.to_numpy().ravel() + + # Null-space line in the PLS score space. + steps = np.linspace(-4.0, 4.0, 60) + ns_line = np.array([tau + s * g_ns for s in steps]) + + # A second inversion, toward a different target taste. In a single-response + # model every target shares the same null-space direction, so this line is + # parallel to the first, just shifted to pass through its own solution. + result2 = pls.invert(y_desired=SECOND_TASTE) + tau2 = result2.scores.to_numpy() + ns_line2 = np.array([tau2 + s * g_ns for s in steps]) + + result3 = pls.invert(y_desired=THIRD_TASTE) + tau3 = result3.scores.to_numpy() + ns_line3 = np.array([tau3 + s * g_ns for s in steps]) + + # O-PLS orthogonal space: walk along the orthogonal loading in the input + # space (scaled), starting from the O-PLS design, then project into the PLS + # score space with the direct weights. These points should land on the NS. + opls_result = opls.invert(y_desired=TARGET_TASTE) + x_scaler = pls._x_scaler # noqa: SLF001 - reuse the fitted scaling for projection + x_design_scaled = x_scaler.transform(opls_result.x_new.to_frame().T).to_numpy().ravel() + os_direction = opls_result.orthogonal_space_basis.to_numpy().ravel() + os_direction = os_direction / np.linalg.norm(os_direction) + direct_weights = pls.direct_weights_.to_numpy() # W* maps scaled X -> scores + os_points = np.array( + [(x_design_scaled + s * os_direction) @ direct_weights for s in np.linspace(-3.0, 3.0, 25)] + ) + + scores = pls.scores_.to_numpy() + + fig, ax = plt.subplots(figsize=(6.4, 5.2)) + cal = ax.scatter(scores[:, 0], scores[:, 1], color=DARK_BLUE, s=30, label="Calibration cheeses") + (line3,) = ax.plot( + ns_line3[:, 0], ns_line3[:, 1], color=PURPLE, lw=2, linestyle=":", label=f"Null space (taste {THIRD_TASTE})" + ) + (line1,) = ax.plot( + ns_line[:, 0], ns_line[:, 1], color=ORANGE, lw=2, label=f"Null space (taste {TARGET_TASTE})" + ) + (line2,) = ax.plot( + ns_line2[:, 0], ns_line2[:, 1], color=RED, lw=2, linestyle="--", label=f"Null space (taste {SECOND_TASTE})" + ) + os_circles = ax.scatter( + os_points[:, 0], + os_points[:, 1], + facecolors="none", + edgecolors=GREEN, + s=70, + linewidths=1.8, + label="Orthogonal space (O-PLS), projected", + ) + di = ax.scatter( + tau[0], tau[1], color=BLACK, marker="s", s=70, zorder=5, label="Direct-inversion solution" + ) + ax.axhline(0, color="0.7", lw=0.8) + ax.axvline(0, color="0.7", lw=0.8) + ax.set_xlabel("$t_1$") + ax.set_ylabel("$t_2$") + ax.set_title(f"Cheddar cheese: designing toward a taste of {TARGET_TASTE}") + + # Two legends: the null-space lines (lowest to highest taste) at top left, and + # the markers at bottom right. + lines_legend = ax.legend( + handles=[line3, line1, line2], loc="upper left", framealpha=0.9, fontsize=9, title="Null spaces" + ) + ax.add_artist(lines_legend) + ax.legend(handles=[cal, os_circles, di], loc="lower right", framealpha=0.9, fontsize=9) + fig.tight_layout() + fig.savefig(out_dir / "pls-model-inversion-null-space.png") + plt.close(fig) + + +def main() -> None: + """Entry point: render into the given directory (default: this script's own).""" + out_dir = Path(sys.argv[1]) if len(sys.argv) > 1 else Path(__file__).resolve().parent + build_figure(out_dir) + + +if __name__ == "__main__": + main()