""" ****************************************************************************** SA_BioML_Pyrolysis.py ******************************************************************************* Supplementary analysis script for: "Mechanistic Machine Learning for Biomass Pyrolysis: Integrating Lignocellulosic Structure, Thermal Severity, and Multi-Objective Optimization for Organic Bio-Oil Yield" Hnin Yee Aye, Fabrizio Scala, Libor Stepanec, Ohn Zin Lin, Dagmar Juchelkova Energy Conversion and Management: X Purpose ------- Script for all supplementary analyses. Analyses -------- S1. Descriptive statistics for all input and output variables S2. Spearman rank correlation matrix (process variables + TSI) S3. Variance Inflation Factor (VIF) for all features S4. Held-out test set evaluation (80/20 split, LightGBM) S5. Five-model cross-validation: base features vs. base + TSI S6. LightGBM hyperparameter configuration S7. MLP hyperparameter configuration S8. Model training and single-sample prediction times S9. Residual diagnostics (predicted vs. residual; normal Q-Q plot) S10. ALE curves with 95% bootstrap confidence intervals Outputs (written to OUTDIR) --------------------------- descriptive_statistics.csv spearman_correlation.csv spearman_correlation_heatmap.png vif_results.csv holdout_test_results.csv tsi_contribution_all_models.csv tsi_contribution_lgbm_bar.png lgbm_hyperparameters.csv mlp_hyperparameters.csv model_timing.csv residual_plots_bio_liquid.png residual_plots_oil_org.png ale_bootstrap_ci_bio_liquid.png ale_bootstrap_ci_oil_org.png Requirements ------------ Python >= 3.9 numpy, pandas, matplotlib, scipy, scikit-learn >= 1.3, xgboost >= 1.7, lightgbm >= 4.0, statsmodels >= 0.14 alibi >= 0.9 (optional; built-in manual ALE used if not installed) Install all: pip install numpy pandas matplotlib scipy scikit-learn xgboost lightgbm statsmodels alibi Dataset ------- Publicly available: https://doi.org/10.5281/zenodo.19205977 Expected filename : pyrolysis.xlsx Usage ----- 1. Set EXCEL_PATH below to the location of pyrolysis.xlsx. 2. Optionally adjust OUTDIR, RANDOM_STATE, N_BOOT_ALE. 3. Run: python SA_BioML_Pyrolysis.py ****************************************************************************** """ # Standard library import os import time import warnings # Third-party import numpy as np import pandas as pd import scipy.stats as stats import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from sklearn.model_selection import KFold, train_test_split from sklearn.metrics import r2_score, mean_absolute_error, mean_squared_error from sklearn.linear_model import LinearRegression from sklearn.ensemble import RandomForestRegressor from sklearn.neural_network import MLPRegressor from xgboost import XGBRegressor from lightgbm import LGBMRegressor from statsmodels.stats.outliers_influence import variance_inflation_factor # Optional: alibi for ALE try: from alibi.explainers import ALE as AlibiALE HAS_ALIBI = True except ImportError: HAS_ALIBI = False print("[INFO] alibi not installed — built-in manual ALE will be used.") print(" Install with: pip install alibi\n") warnings.filterwarnings("ignore") # ******************************************** # CONFIGURATION — only edit this block # ******************************************** EXCEL_PATH = r"C:\Users\admin\Desktop\Bio-ML\V14122025\pyrolysis.xlsx" OUTDIR = os.path.join(os.path.dirname(EXCEL_PATH), "output_SA") RANDOM_STATE = 42 N_SPLITS = 5 N_BOOT_ALE = 100 # reduce to 50 if runtime is a concern # ****************************************************************** plt.rcParams.update({ "font.family": "serif", "font.serif": ["Times New Roman", "DejaVu Serif"], "mathtext.fontset": "dejavuserif", "font.size": 11, "axes.labelsize": 12, "axes.titlesize": 13, "xtick.labelsize": 10, "ytick.labelsize": 10, "legend.fontsize": 10, "lines.linewidth": 2.0, "axes.linewidth": 1.2, "grid.linewidth": 0.5, "axes.grid": False, "savefig.dpi": 300, "savefig.bbox": "tight", "savefig.pad_inches": 0.05, }) os.makedirs(OUTDIR, exist_ok=True) # ****************************************************************************** # SECTION 0 — DATA LOADING AND PREPROCESSING # ****************************************************************************** print("=" * 70) print("SA_BioML_Pyrolysis.py — Supplementary Analysis") print("=" * 70) print(f"\nLoading: {EXCEL_PATH}") df_raw = pd.read_excel(EXCEL_PATH) col_map = { "Ash(wt%)": "Ash", "FixedCarbon(wt%)": "fc", "Volatiles(wt%)": "vol", "C(wt%)": "C", "H/Ceff.(wt%)": "hceff", "H(wt%)": "H", "O(wt%)": "O", "N(wt%)": "N", "Cellulose(wt%)": "cellulose", "Hemicellulose(wt%)": "hemicellulose", "Lignin(wt%)": "lignin", "T(°C)": "temp", "HeatingRate(°C/min)": "HR", "ParticleSize(mm)": "PS", "FlowRate-Nitrogen(ml/min)": "Flow-N2", "bio-liquid yield(wt%)": "bio_liquid", "Water content": "water", } df = df_raw.rename(columns=col_map) required_cols = [ "Ash", "fc", "vol", "C", "H", "O", "N", "hceff", "cellulose", "hemicellulose", "lignin", "temp", "HR", "PS", "Flow-N2", "bio_liquid", "water", ] df = df.dropna(subset=required_cols).reset_index(drop=True) eps = 1e-6 df["oil_org"] = df["bio_liquid"] - df["water"] df["TSI"] = df["temp"] * np.log(df["HR"] + 1.0) / (df["PS"] + eps) BASE_FEATURES = [ "Ash", "fc", "vol", "C", "H", "O", "N", "hceff", "cellulose", "hemicellulose", "lignin", "temp", "HR", "PS", "Flow-N2", ] FEATURES_WITH_TSI = BASE_FEATURES + ["TSI"] PROCESS_FEATURES = ["temp", "HR", "PS", "Flow-N2"] TARGETS = { "bio_liquid": "Total bio-liquid yield (wt%)", "oil_org": "Organic oil yield (wt%)", } print(f"Cleaned dataset : {df.shape[0]} samples, {len(FEATURES_WITH_TSI)} features") print(f"Outputs written : {OUTDIR}\n") # ****************************************************************************** # HELPER FUNCTIONS # ****************************************************************************** def make_lgbm(): return LGBMRegressor( n_estimators=500, learning_rate=0.05, num_leaves=31, subsample=0.9, colsample_bytree=0.9, random_state=RANDOM_STATE, verbosity=-1, ) def make_all_models(): return { "Linear Regression": LinearRegression(), "Random Forest": RandomForestRegressor( n_estimators=500, random_state=RANDOM_STATE, n_jobs=-1), "XGBoost": XGBRegressor( n_estimators=500, learning_rate=0.05, max_depth=6, subsample=0.9, colsample_bytree=0.9, objective="reg:squarederror", random_state=RANDOM_STATE, n_jobs=-1, verbosity=0), "MLP": MLPRegressor( hidden_layer_sizes=(64, 32), activation="relu", solver="adam", max_iter=500, random_state=RANDOM_STATE), "LightGBM": make_lgbm(), } def cv_evaluate(X, y, model, n_splits=N_SPLITS): """K-fold CV; returns mean/std of R², MAE, RMSE.""" kf = KFold(n_splits=n_splits, shuffle=True, random_state=RANDOM_STATE) r2s, maes, rmses = [], [], [] for tr, te in kf.split(X): model.fit(X[tr], y[tr]) yp = model.predict(X[te]) r2s.append(r2_score(y[te], yp)) maes.append(mean_absolute_error(y[te], yp)) rmses.append(np.sqrt(mean_squared_error(y[te], yp))) return (np.mean(r2s), np.std(r2s), np.mean(maes), np.std(maes), np.mean(rmses), np.std(rmses)) def compute_ale_manual(model, X, feat_idx, n_bins=20): """ 1-D ALE curve following Apley & Zhu (2020). No external dependencies. Used as primary or fallback ALE engine. Returns ------- centres : ndarray, bin midpoints ale_cum : ndarray, centred cumulative ALE (wt%) Both are empty arrays if computation is not possible. """ x_col = X[:, feat_idx] quantiles = np.unique( np.percentile(x_col, np.linspace(0, 100, n_bins + 1)) ) if len(quantiles) < 2: return np.array([]), np.array([]) ale_vals, centres = [], [] for lo, hi in zip(quantiles[:-1], quantiles[1:]): mask = (x_col >= lo) & (x_col <= hi) if mask.sum() == 0: continue X_lo = X[mask].copy(); X_lo[:, feat_idx] = lo X_hi = X[mask].copy(); X_hi[:, feat_idx] = hi ale_vals.append((model.predict(X_hi) - model.predict(X_lo)).mean()) centres.append((lo + hi) / 2.0) if len(ale_vals) == 0: return np.array([]), np.array([]) ale_cum = np.cumsum(ale_vals) ale_cum -= ale_cum.mean() return np.array(centres), ale_cum def run_ale_with_bootstrap(model, X_full, feat, feat_idx_map, n_boot=N_BOOT_ALE, use_alibi=False, ale_explainer=None): """ Central ALE curve + 95% bootstrap CI for one feature. Bootstrap strategy ------------------ Each iteration tries alibi first. On any alibi failure, the same bootstrap sample is immediately retried with compute_ale_manual. This guarantees boot_ales is populated even when alibi is unreliable, preventing np.percentile from receiving an empty array. Returns ------- (xs, ys, ci_lo, ci_hi) — ci_lo/ci_hi are None if < 5 samples succeeded. None — if the central ALE itself fails. """ f_idx = feat_idx_map[feat] # Central ALE xs = ys = None if use_alibi and ale_explainer is not None: try: exp = ale_explainer.explain(X_full, features=[f_idx]) xs = np.asarray(exp.feature_values[0]).ravel() ys = np.asarray(exp.ale_values[0]).ravel() except Exception as e: print(f" [alibi central ALE failed for {feat}: {e}]" " — switching to manual") use_alibi = False if xs is None: xs, ys = compute_ale_manual(model, X_full, f_idx) if len(xs) == 0: print(f" [ALE empty for {feat}] — panel skipped") return None # Bootstrap rng = np.random.default_rng(RANDOM_STATE) boot_ales = [] n_alibi_fail = 0 n_total_fail = 0 for b in range(n_boot): idx = rng.choice(len(X_full), size=len(X_full), replace=True) Xb = X_full[idx] xb = yb = None # Step 1: try alibi if use_alibi and ale_explainer is not None: try: exp_b = ale_explainer.explain(Xb, features=[f_idx]) xb = np.asarray(exp_b.feature_values[0]).ravel() yb = np.asarray(exp_b.ale_values[0]).ravel() except Exception: n_alibi_fail += 1 xb = yb = None # Step 2: manual fallback if alibi failed or returned empty if xb is None or len(xb) < 2: try: xb, yb = compute_ale_manual(model, Xb, f_idx) xb = np.asarray(xb).ravel() yb = np.asarray(yb).ravel() except Exception as e: n_total_fail += 1 if n_total_fail <= 3: print(f" [iter {b} failed for {feat}: {e}]") continue if len(xb) < 2: n_total_fail += 1 continue boot_ales.append(np.interp(xs.ravel(), xb, yb)) # Diagnostics if n_alibi_fail > 0: print(f" [{feat}] alibi failed on {n_alibi_fail}/{n_boot} iters; " "manual used for those.") if n_total_fail > 0: print(f" [{feat}] {n_total_fail}/{n_boot} iters failed completely.") if len(boot_ales) < 5: print(f" [WARNING] Only {len(boot_ales)} successful bootstrap " f"samples for {feat}. CI band omitted.") return xs, ys, None, None mat = np.array(boot_ales) ci_lo = np.percentile(mat, 2.5, axis=0) ci_hi = np.percentile(mat, 97.5, axis=0) return xs, ys, ci_lo, ci_hi # ****************************************************************************** # S1 — DESCRIPTIVE STATISTICS # ****************************************************************************** print("─" * 70) print("S1. Descriptive statistics") print("─" * 70) all_vars = FEATURES_WITH_TSI + ["bio_liquid", "water", "oil_org"] desc = df[all_vars].agg( ["count", "mean", "median", "std", "min", "max", lambda x: x.max() - x.min(), lambda x: x.skew()] ).T desc.columns = ["N", "Mean", "Median", "Std", "Min", "Max", "Range", "Skewness"] desc = desc.round(3) desc.index.name = "Variable" p = os.path.join(OUTDIR, "descriptive_statistics.csv") desc.to_csv(p) print(desc.to_string()) print(f"\nSaved: {p}\n") # ****************************************************************************** # S2 — SPEARMAN RANK CORRELATION # ****************************************************************************** print("─" * 70) print("S2. Spearman rank correlation (process variables + TSI)") print("─" * 70) proc_tsi = PROCESS_FEATURES + ["TSI"] sp_mat, _ = stats.spearmanr(df[proc_tsi]) sp_df = pd.DataFrame(sp_mat, index=proc_tsi, columns=proc_tsi).round(3) p_csv = os.path.join(OUTDIR, "spearman_correlation.csv") sp_df.to_csv(p_csv) print(sp_df.to_string()) fig, ax = plt.subplots(figsize=(6, 5)) im = ax.imshow(sp_mat, vmin=-1, vmax=1, cmap="RdBu_r") ax.set_xticks(range(len(proc_tsi))); ax.set_xticklabels(proc_tsi, rotation=45, ha="right") ax.set_yticks(range(len(proc_tsi))); ax.set_yticklabels(proc_tsi) for i in range(len(proc_tsi)): for j in range(len(proc_tsi)): ax.text(j, i, f"{sp_mat[i,j]:.2f}", ha="center", va="center", fontsize=9) plt.colorbar(im, ax=ax, label="Spearman r") ax.set_title("Spearman Correlation — Process Variables + TSI") plt.tight_layout() p_png = os.path.join(OUTDIR, "spearman_correlation_heatmap.png") plt.savefig(p_png); plt.close() print(f"\nSaved: {p_csv}\nSaved: {p_png}\n") # ****************************************************************************** # S3 — VARIANCE INFLATION FACTOR # ****************************************************************************** print("─" * 70) print("S3. Variance Inflation Factor (all features including TSI)") print("─" * 70) X_vif = df[FEATURES_WITH_TSI].astype(float) vif_data = pd.DataFrame({ "Feature": FEATURES_WITH_TSI, "VIF": [variance_inflation_factor(X_vif.values, i) for i in range(X_vif.shape[1])], }).round(2) p = os.path.join(OUTDIR, "vif_results.csv") vif_data.to_csv(p, index=False) print(vif_data.to_string(index=False)) print("\nNote: VIF > 10 indicates potentially problematic multicollinearity.") print(f"Saved: {p}\n") # ****************************************************************************** # S4 — HELD-OUT TEST SET (80/20) # ****************************************************************************** print("─" * 70) print("S4. Held-out test set — LightGBM, 80/20 split") print("─" * 70) holdout_rows = [] for tgt_key in TARGETS: y = df[tgt_key].values for flabel, fcols in [("Base features", BASE_FEATURES), ("Base + TSI", FEATURES_WITH_TSI)]: X = df[fcols].values.astype(float) Xtr, Xte, ytr, yte = train_test_split( X, y, test_size=0.20, random_state=RANDOM_STATE) mdl = make_lgbm(); mdl.fit(Xtr, ytr) yp = mdl.predict(Xte) r2 = r2_score(yte, yp) mae = mean_absolute_error(yte, yp) rmse= np.sqrt(mean_squared_error(yte, yp)) holdout_rows.append({ "Target": tgt_key, "Feature set": flabel, "n_train": len(ytr), "n_test": len(yte), "R2": round(r2, 4), "MAE (wt%)": round(mae, 4), "RMSE (wt%)": round(rmse, 4), }) print(f" [{tgt_key} | {flabel}] R²={r2:.4f} MAE={mae:.4f} RMSE={rmse:.4f}") p = os.path.join(OUTDIR, "holdout_test_results.csv") pd.DataFrame(holdout_rows).to_csv(p, index=False) print(f"\nSaved: {p}\n") # ****************************************************************************** # S5 — FIVE-MODEL CV: BASE vs. BASE + TSI # ****************************************************************************** print("─" * 70) print("S5. Five-model CV: base features vs. base + TSI") print("─" * 70) tsi_rows = [] for tgt_key in TARGETS: y = df[tgt_key].values for flabel, fcols in [("base", BASE_FEATURES), ("with_tsi", FEATURES_WITH_TSI)]: X = df[fcols].values.astype(float) for mname, mdl in make_all_models().items(): print(f" CV: {tgt_key} | {flabel} | {mname}") r2m,r2s,maem,maes,rmsem,rmses = cv_evaluate(X, y, mdl) tsi_rows.append({ "Target": tgt_key, "Feature set": flabel, "Model": mname, "R2_mean": round(r2m, 4), "R2_std": round(r2s, 4), "MAE_mean": round(maem, 4), "MAE_std": round(maes, 4), "RMSE_mean": round(rmsem, 4), "RMSE_std": round(rmses, 4), }) tsi_df = pd.DataFrame(tsi_rows) p = os.path.join(OUTDIR, "tsi_contribution_all_models.csv") tsi_df.to_csv(p, index=False) print(f"\nSaved: {p}") # Bar chart: LightGBM base vs. with_tsi lgbm_r = tsi_df[tsi_df["Model"] == "LightGBM"] fig, axes = plt.subplots(1, 2, figsize=(11, 5)) mets = ["R2_mean", "MAE_mean", "RMSE_mean"] xp = np.arange(len(mets)); w = 0.35 for ax, tk in zip(axes, ["bio_liquid", "oil_org"]): sub = lgbm_r[lgbm_r["Target"] == tk] base = sub[sub["Feature set"] == "base"].iloc[0] wtsi = sub[sub["Feature set"] == "with_tsi"].iloc[0] ax.bar(xp - w/2, [base[m] for m in mets], w, label="Base features", color="#3B82F6") ax.bar(xp + w/2, [wtsi[m] for m in mets], w, label="Base + TSI", color="#F97316") ax.set_xticks(xp); ax.set_xticklabels(["R²", "MAE (wt%)", "RMSE (wt%)"]) ax.set_title(f"LightGBM — {tk.replace('_',' ').title()}") ax.legend(); ax.set_ylabel("Metric value") plt.suptitle("Effect of TSI on LightGBM performance (5-fold CV)", fontsize=13) plt.tight_layout() p = os.path.join(OUTDIR, "tsi_contribution_lgbm_bar.png") plt.savefig(p); plt.close() print(f"Saved: {p}\n") # ****************************************************************************** # S6 — LIGHTGBM HYPERPARAMETERS # ****************************************************************************** print("─" * 70) print("S6. LightGBM hyperparameter configuration") print("─" * 70) lgbm_cfg = { "n_estimators": 500, "learning_rate": 0.05, "num_leaves": 31, "subsample": 0.9, "colsample_bytree": 0.9, "min_child_samples": "20 (default)", "reg_alpha": "0.0 (default, L1)", "reg_lambda": "0.0 (default, L2)", "boosting_type": "gbdt (default)", "objective": "regression (default)", "n_jobs": -1, "random_state": RANDOM_STATE, "Note": ("Final surrogate trained with n_estimators=800, learning_rate=0.03; " "all other parameters identical. No automated tuning applied."), } lgbm_df = pd.DataFrame(list(lgbm_cfg.items()), columns=["Hyperparameter", "Value"]) p = os.path.join(OUTDIR, "lgbm_hyperparameters.csv") lgbm_df.to_csv(p, index=False) print(lgbm_df.to_string(index=False)) print(f"\nSaved: {p}\n") # ****************************************************************************** # S7 — MLP HYPERPARAMETERS # ****************************************************************************** print("─" * 70) print("S7. MLP hyperparameter configuration") print("─" * 70) mlp_cfg = { "hidden_layer_sizes": "(64, 32) — two hidden layers", "activation": "relu", "solver": "adam", "learning_rate": "adaptive", "learning_rate_init": "0.001 (default)", "max_iter": 500, "alpha (L2)": "0.0001 (default)", "batch_size": "auto [min(200, n_samples)]", "early_stopping": "False (default)", "random_state": RANDOM_STATE, "Note": ("No hyperparameter tuning applied. Results serve as a baseline. " "High R² std across folds reflects weight-initialisation sensitivity; " "tuning would be required for MLP to match tree-based methods."), } mlp_df = pd.DataFrame(list(mlp_cfg.items()), columns=["Hyperparameter", "Value"]) p = os.path.join(OUTDIR, "mlp_hyperparameters.csv") mlp_df.to_csv(p, index=False) print(mlp_df.to_string(index=False)) print(f"\nSaved: {p}\n") # ****************************************************************************** # S8 — MODEL TIMING # ****************************************************************************** print("─" * 70) print("S8. Model training and single-sample prediction times") print("─" * 70) timing_rows = [] for tgt_key in TARGETS: y = df[tgt_key].values X = df[FEATURES_WITH_TSI].values.astype(float) for mname, mdl in make_all_models().items(): t0 = time.perf_counter(); mdl.fit(X, y) tr_t = time.perf_counter() - t0 X1 = X[[0]] t0 = time.perf_counter() for _ in range(100): mdl.predict(X1) pr_t = (time.perf_counter() - t0) / 100 timing_rows.append({ "Target": tgt_key, "Model": mname, "Train time (s)": round(tr_t, 4), "Predict time (ms)": round(pr_t * 1000, 4), }) print(f" [{tgt_key} | {mname}] train={tr_t:.3f}s " f"predict={pr_t*1000:.3f}ms") p = os.path.join(OUTDIR, "model_timing.csv") pd.DataFrame(timing_rows).to_csv(p, index=False) print(f"\nSaved: {p}\n") # ****************************************************************************** # S9 — RESIDUAL DIAGNOSTICS # ****************************************************************************** print("─" * 70) print("S9. Residual diagnostics (LightGBM, held-out 20%)") print("─" * 70) for tgt_key in TARGETS: y = df[tgt_key].values X = df[FEATURES_WITH_TSI].values.astype(float) Xtr, Xte, ytr, yte = train_test_split( X, y, test_size=0.20, random_state=RANDOM_STATE) mdl = make_lgbm(); mdl.fit(Xtr, ytr) yp = mdl.predict(Xte) res = yte - yp r2 = r2_score(yte, yp) rmse= np.sqrt(mean_squared_error(yte, yp)) fig, axes = plt.subplots(1, 2, figsize=(12, 5)) ax = axes[0] ax.scatter(yp, res, alpha=0.6, s=30, color="#3B82F6", edgecolors="none") ax.axhline(0, color="red", linewidth=1.5, linestyle="--") ax.set_xlabel(f"Predicted {tgt_key.replace('_',' ')} (wt%)") ax.set_ylabel("Residual [actual − predicted] (wt%)") ax.set_title("(a) Predicted vs. Residual") ax.text(0.05, 0.93, f"R² = {r2:.3f}\nRMSE = {rmse:.3f} wt%", transform=ax.transAxes, fontsize=10, verticalalignment="top", bbox=dict(boxstyle="round", facecolor="white", alpha=0.75)) ax = axes[1] (osm, osr), (slope, intercept, _) = stats.probplot(res, dist="norm") ax.scatter(osm, osr, alpha=0.6, s=30, color="#F97316", edgecolors="none") xl = np.array([osm[0], osm[-1]]) ax.plot(xl, slope*xl + intercept, color="red", linewidth=1.5, linestyle="--", label="Normal reference") ax.set_xlabel("Theoretical quantiles") ax.set_ylabel("Sample quantiles (residuals, wt%)") ax.set_title("(b) Normal Q-Q Plot of Residuals") ax.legend() plt.suptitle( f"Residual diagnostics — LightGBM | {tgt_key} | held-out 20% (n={len(yte)})", fontsize=12) plt.tight_layout() p = os.path.join(OUTDIR, f"residual_plots_{tgt_key}.png") plt.savefig(p); plt.close() print(f" Saved: {p}") print() # ****************************************************************************** # S10 — ALE WITH 95% BOOTSTRAP CONFIDENCE INTERVALS # ****************************************************************************** method_label = "alibi" if HAS_ALIBI else "manual (built-in)" print("─" * 70) print(f"S10. ALE + 95% bootstrap CI (n_boot={N_BOOT_ALE}, method={method_label})") print("─" * 70) # Train final LightGBM models on full dataset final_models = {} for tgt_key in TARGETS: m = LGBMRegressor( n_estimators=800, learning_rate=0.03, num_leaves=31, subsample=0.9, colsample_bytree=0.9, random_state=RANDOM_STATE, verbosity=-1, ) m.fit(df[FEATURES_WITH_TSI].values.astype(float), df[tgt_key].values) final_models[tgt_key] = m X_full = df[FEATURES_WITH_TSI].values.astype(float) feat_idx_map = {f: i for i, f in enumerate(FEATURES_WITH_TSI)} for tgt_key in TARGETS: model = final_models[tgt_key] use_alibi = HAS_ALIBI ale_expl = None if use_alibi: try: ale_expl = AlibiALE(predictor=model.predict, feature_names=FEATURES_WITH_TSI) except Exception as e: print(f" [alibi init failed: {e}] — using manual fallback") use_alibi = False fig, axes = plt.subplots(2, 2, figsize=(12, 10)) axes = axes.ravel() for ax, feat in zip(axes, PROCESS_FEATURES): print(f" [{tgt_key} | {feat}] computing ALE + CI ...") result = run_ale_with_bootstrap( model, X_full, feat, feat_idx_map, n_boot=N_BOOT_ALE, use_alibi=use_alibi, ale_explainer=ale_expl, ) if result is None: ax.set_visible(False); continue xs, ys, ci_lo, ci_hi = result ax.plot(xs, ys, color="#1D4ED8", linewidth=2, label="ALE") if ci_lo is not None: ax.fill_between(xs, ci_lo, ci_hi, alpha=0.25, color="#3B82F6", label="95% CI (bootstrap)") ax.axhline(0, color="gray", linewidth=0.8, linestyle="--") ax.set_xlabel(feat, fontsize=13) ax.set_ylabel("ALE (wt%)", fontsize=13) ax.set_title(feat, fontsize=14) ax.legend(fontsize=10) fig.suptitle( f"ALE with 95% bootstrap CI — {tgt_key.replace('_',' ').title()} " f"(n_boot={N_BOOT_ALE}, method: {method_label})", fontsize=13, fontweight="bold", ) fig.tight_layout() p = os.path.join(OUTDIR, f"ale_bootstrap_ci_{tgt_key}.png") fig.savefig(p); plt.close(fig) print(f" Saved: {p}") print() # ****************************************************************************** # COMPLETION SUMMARY # ****************************************************************************** print("=" * 70) print("COMPLETE") print("=" * 70) print(f"\nAll outputs in: {OUTDIR}\n") manifest = [ ("descriptive_statistics.csv", "S1 — Descriptive statistics"), ("spearman_correlation.csv", "S2 — Spearman correlation matrix"), ("spearman_correlation_heatmap.png", "S2 — Spearman heatmap"), ("vif_results.csv", "S3 — Variance Inflation Factors"), ("holdout_test_results.csv", "S4 — Held-out 20% test results"), ("tsi_contribution_all_models.csv", "S5 — 5-model CV, base vs. base+TSI"), ("tsi_contribution_lgbm_bar.png", "S5 — TSI contribution bar chart"), ("lgbm_hyperparameters.csv", "S6 — LightGBM hyperparameters"), ("mlp_hyperparameters.csv", "S7 — MLP hyperparameters"), ("model_timing.csv", "S8 — Training and prediction times"), ("residual_plots_bio_liquid.png", "S9 — Residual diagnostics (bio_liquid)"), ("residual_plots_oil_org.png", "S9 — Residual diagnostics (oil_org)"), ("ale_bootstrap_ci_bio_liquid.png", "S10 — ALE + CI (bio_liquid)"), ("ale_bootstrap_ci_oil_org.png", "S10 — ALE + CI (oil_org)"), ] w = max(len(f) for f, _ in manifest) for fname, desc in manifest: status = "OK" if os.path.exists(os.path.join(OUTDIR, fname)) else "MISSING" print(f" [{status}] {fname:<{w+2}}{desc}") print()