from __future__ import annotations import json from dataclasses import dataclass from pathlib import Path import joblib import numpy as np import pandas as pd from gfis.config import FEATURE_COLUMNS, MODEL_DIR, REPORT_DIR from gfis.data import theoretical_methane_upper_bound from gfis.db import log_prediction, log_simulation from gfis.preprocessing import add_engineered_features, clean_frame DEFAULT_INPUT = { "temperature": 37.0, "pH": 7.1, "OLR": 3.2, "HRT": 25.0, "TS": 9.0, "VS": 6.8, "C_N_ratio": 25.0, "ambient_temperature": 28.0, "moisture": 82.0, "methane_yield_lag1": 200.0, "methane_yield_roll3": 200.0, } @dataclass class GFISResult: methane_yield: float physics_upper_bound: float physics_violation: bool vfa_alk_ratio: float stability_label: str def classify_stability(vfa_alk_ratio: float) -> str: if vfa_alk_ratio < 0.30: return "Stable" if vfa_alk_ratio < 0.45: return "Warning" return "Critical" def _fallback_prediction(values: dict) -> tuple[float, float]: vs = float(values["VS"]) temp_score = np.exp(-((float(values["temperature"]) - 37.0) ** 2) / 95.0) ph_score = np.exp(-((float(values["pH"]) - 7.15) ** 2) / 0.20) hrt_effect = 1.0 - np.exp(-float(values["HRT"]) / 18.0) methane = 85.0 + (36.0 * vs + 3.0 * float(values["TS"])) * temp_score * ph_score * hrt_effect methane += 7.5 * float(values["OLR"]) - max(0.0, float(values["OLR"]) - 4.2) * 18.0 vfa_alk = ( 0.18 + 0.055 * max(0.0, float(values["OLR"]) - 3.0) + 0.018 * abs(float(values["pH"]) - 7.1) + 0.010 * max(0.0, 32.0 - float(values["temperature"])) - 0.003 * max(0.0, float(values["HRT"]) - 22.0) ) return float(methane), float(np.clip(vfa_alk, 0.08, 0.75)) class GFISService: def __init__(self, model_dir: Path = MODEL_DIR): self.model_dir = model_dir self.bundle = None bundle_path = model_dir / "gfis_tabular_bundle.joblib" if bundle_path.exists(): try: self.bundle = joblib.load(bundle_path) except Exception: self.bundle = None def _frame_from_values(self, values: dict) -> pd.DataFrame: row = {**DEFAULT_INPUT, **values} frame = pd.DataFrame([row]) for col in FEATURE_COLUMNS: frame[f"{col}_lag1"] = frame[col] frame[f"{col}_roll3"] = frame[col] frame["methane_yield_lag1"] = float(row.get("methane_yield_lag1", 200.0)) frame["methane_yield_roll3"] = float(row.get("methane_yield_roll3", frame["methane_yield_lag1"].iloc[0])) frame["OLR_HRT_ratio"] = frame["OLR"] / frame["HRT"].replace(0, np.nan) frame["VS_TS_ratio"] = frame["VS"] / frame["TS"].replace(0, np.nan) frame["temperature_pH_interaction"] = frame["temperature"] * frame["pH"] return clean_frame(frame) def predict(self, values: dict) -> GFISResult: merged = {**DEFAULT_INPUT, **values} upper = float(theoretical_methane_upper_bound(np.asarray([merged["VS"]]))[0]) if self.bundle is not None: frame = self._frame_from_values(merged) feature_names = self.bundle["feature_names"] x = self.bundle["scaler"].transform(frame[feature_names]) methane = float(self.bundle["tabular_model"].predict(x)[0]) vfa_alk = float(self.bundle["soft_sensor"].predict(x)[0]) else: methane, vfa_alk = _fallback_prediction(merged) physics_violation = methane > upper methane = min(methane, upper) result = GFISResult( methane_yield=round(methane, 3), physics_upper_bound=round(upper, 3), physics_violation=bool(physics_violation), vfa_alk_ratio=round(vfa_alk, 3), stability_label=classify_stability(vfa_alk), ) log_prediction(merged, result.__dict__) return result def simulate(self, base_values: dict, scenarios: list[dict]) -> list[dict]: outputs = [] for scenario in scenarios: merged = {**base_values, **scenario} result = self.predict(merged) outputs.append({**merged, **result.__dict__}) log_simulation(base_values, scenarios, outputs) return outputs def optimize(self, base_values: dict) -> dict: best = None for temp in np.linspace(33, 40, 8): for ph in np.linspace(6.8, 7.5, 8): for olr in np.linspace(2.0, 4.6, 7): result = self.predict({**base_values, "temperature": temp, "pH": ph, "OLR": olr}) if result.stability_label == "Critical": continue candidate = {**result.__dict__, "temperature": round(float(temp), 2), "pH": round(float(ph), 2), "OLR": round(float(olr), 2)} if best is None or candidate["methane_yield"] > best["methane_yield"]: best = candidate return best or self.predict(base_values).__dict__ def plant_run(self, base_values: dict, hours: int = 48) -> list[dict]: """Run a deterministic multi-hour digital-twin trace. This simulates plant operation with mild operational drift and feeds the previous methane output back as methane history. It is intentionally lightweight but useful for viva/demo evidence. """ outputs: list[dict] = [] state = {**DEFAULT_INPUT, **base_values} methane_history = float(state.get("methane_yield_lag1", 200.0)) rolling_history = methane_history for hour in range(hours): drift = np.sin(hour / 7.5) scenario = { **state, "temperature": float(state["temperature"]) + 0.7 * drift, "pH": float(state["pH"]) + 0.04 * np.sin(hour / 11.0), "OLR": max(1.0, float(state["OLR"]) + 0.18 * np.sin(hour / 5.0)), "methane_yield_lag1": methane_history, "methane_yield_roll3": rolling_history, } result = self.predict(scenario) methane_history = result.methane_yield rolling_history = 0.65 * rolling_history + 0.35 * methane_history outputs.append( { "hour": hour, "temperature": round(scenario["temperature"], 3), "pH": round(scenario["pH"], 3), "OLR": round(scenario["OLR"], 3), **result.__dict__, } ) log_simulation(base_values, [{"plant_run_hours": hours}], outputs) return outputs def evaluation(self) -> dict: metrics_path = REPORT_DIR / "evaluation_metrics.json" if metrics_path.exists(): return json.loads(metrics_path.read_text()) return {"status": "not_trained", "message": "Run python -m gfis.train first."}