"""Type-1 Takagi–Sugeno–Kang fuzzy inference system (see SPEC.md).""" from __future__ import annotations import numpy as np from sklearn.cluster import KMeans # Minimum denominator for normalized firing strengths (see predict). _FIRE_EPS = 1e-12 class T1TSK: """Type-1 first-order TSK fuzzy system with Gaussian antecedents. Inference (SPEC.md): f_k = prod_j mu_kj(x_j) g_k(x) = a_k0 + sum_j a_kj x_j y = sum_k f_k g_k(x) / sum_k f_k Training follows SPEC.md §4 steps 1–4 on the provided (X, y): standardize features, k-means centers, cluster-based spreads, then least squares on the same normalized weighted consequent form used at prediction time. """ def __init__(self, n_rules: int, random_state: int = 42) -> None: """Create an unfitted T1 TSK model. Args: n_rules: Number of fuzzy rules (K). random_state: Seed for k-means (reproducible initialization). """ if n_rules < 1: raise ValueError(f"n_rules must be >= 1, got {n_rules}") self.n_rules = n_rules self.random_state = random_state self._fitted = False self._x_mean: np.ndarray | None = None self._x_scale: np.ndarray | None = None # Centers m[k, j] and spreads sigma[k, j] (one Gaussian per rule/feature). self.m_: np.ndarray | None = None self.sigma_: np.ndarray | None = None # Consequents: row k is [a_k0, a_k1, ..., a_kp]. self.coef_: np.ndarray | None = None def _check_fitted(self) -> None: if not self._fitted: raise RuntimeError("Call fit(X, y) before predict(X).") def _standardize(self, X: np.ndarray) -> np.ndarray: assert self._x_mean is not None and self._x_scale is not None return (X - self._x_mean) / self._x_scale def _memberships(self, X_scaled: np.ndarray) -> np.ndarray: """Gaussian memberships mu[n, k, j] matching GaussianMF in membership.py.""" assert self.m_ is not None and self.sigma_ is not None # X_scaled: (n, p), m_/sigma_: (K, p) -> (n, K, p) diff = X_scaled[:, None, :] - self.m_[None, :, :] exponent = -(diff**2) / (2.0 * self.sigma_[None, :, :] ** 2) return np.exp(exponent) def _firing_strengths(self, X_scaled: np.ndarray) -> np.ndarray: """Rule firing strengths f[n, k] = prod_j mu_kj(x_j).""" mu = self._memberships(X_scaled) return np.prod(mu, axis=2) def _rule_outputs(self, X_scaled: np.ndarray) -> np.ndarray: """Linear consequents g[n, k] for each rule.""" assert self.coef_ is not None n_samples, n_features = X_scaled.shape X_aug = np.hstack( [np.ones((n_samples, 1), dtype=float), X_scaled] ) return X_aug @ self.coef_.T def predict(self, X: np.ndarray) -> np.ndarray: """Predict target values for feature matrix X (n_samples, n_features).""" self._check_fitted() X = np.asarray(X, dtype=float) if X.ndim != 2: raise ValueError(f"X must be 2D, got shape {X.shape}") X_scaled = self._standardize(X) firing = self._firing_strengths(X_scaled) rule_out = self._rule_outputs(X_scaled) numer = np.sum(firing * rule_out, axis=1) denom = np.sum(firing, axis=1) # The specification does not define an output when all rules underflow. # Use the neutral mean of the rule consequents instead of returning 0. zero_firing = denom <= _FIRE_EPS safe_denom = np.maximum(denom, _FIRE_EPS) predictions = numer / safe_denom if np.any(zero_firing): predictions[zero_firing] = np.mean(rule_out[zero_firing], axis=1) return predictions def fit(self, X: np.ndarray, y: np.ndarray) -> "T1TSK": """Fit antecedents and consequents from training data. Steps (SPEC.md §4): 1. Standardize features using training X only. 2. k-means cluster centers -> Gaussian means m_kj. 3. Per-cluster feature std -> sigma_kj (with a small floor). 4. Weighted least squares for all a_kj using normalized firing weights. """ X = np.asarray(X, dtype=float) y = np.asarray(y, dtype=float).reshape(-1) if X.ndim != 2: raise ValueError(f"X must be 2D, got shape {X.shape}") n_samples, n_features = X.shape if y.shape[0] != n_samples: raise ValueError("X and y must have the same number of samples") if n_samples < self.n_rules: raise ValueError( "Need at least as many samples as rules for k-means." ) # 1. Standardize features (training statistics only). self._x_mean = np.mean(X, axis=0) self._x_scale = np.std(X, axis=0, ddof=1) self._x_scale = np.where( self._x_scale < 1e-12, 1.0, self._x_scale ) X_scaled = self._standardize(X) # 2. Rule centers from k-means on standardized features. kmeans = KMeans( n_clusters=self.n_rules, random_state=self.random_state, n_init=10, ) labels = kmeans.fit_predict(X_scaled) self.m_ = kmeans.cluster_centers_.astype(float) # 3. Gaussian spreads from cluster dispersion (assumption: sample std # per feature within each cluster; empty/tiny clusters use global std). global_std = np.std(X_scaled, axis=0, ddof=1) global_std = np.maximum(global_std, 1e-8) self.sigma_ = np.zeros((self.n_rules, n_features), dtype=float) for k in range(self.n_rules): mask = labels == k if np.count_nonzero(mask) < 2: self.sigma_[k] = global_std else: self.sigma_[k] = np.std(X_scaled[mask], axis=0, ddof=1) # Keep every sigma positive (same idea as rejecting invalid sigma in GaussianMF). self.sigma_ = np.maximum(self.sigma_, 0.1 * global_std) self.sigma_ = np.maximum(self.sigma_, 1e-8) # 4. Least squares on the inference formula with fixed antecedents. firing = self._firing_strengths(X_scaled) weight_sum = np.sum(firing, axis=1, keepdims=True) weight_sum = np.maximum(weight_sum, _FIRE_EPS) w = firing / weight_sum X_aug = np.hstack( [np.ones((n_samples, 1), dtype=float), X_scaled] ) n_params = self.n_rules * (n_features + 1) design = np.zeros((n_samples, n_params), dtype=float) for k in range(self.n_rules): base = k * (n_features + 1) design[:, base] = w[:, k] design[:, base + 1 : base + 1 + n_features] = w[:, k : k + 1] * X_scaled coef_flat, _, _, _ = np.linalg.lstsq(design, y, rcond=None) self.coef_ = coef_flat.reshape(self.n_rules, n_features + 1) self._fitted = True return self