From c8940bbb4213a7a8f62663e0fdb608808bf18930 Mon Sep 17 00:00:00 2001 From: dantzert <47285626+dantzert@users.noreply.github.com> Date: Mon, 7 Sep 2026 23:54:12 +0000 Subject: [PATCH] Fix all open issues: #69, #70, #71 - Issue #69: Replace unstable sinh with stable difference-of-exponentials in UnderdampedOscillatorKernel for zeta >= 1 - Issue #71: Guard log10(0) in compute_detailed_metrics with a small positive floor - Issue #70: Expose Bayesian optimization budget (n_calls, n_initial_points) and add fast_mode flag to speed up lti_system_gen All 69 tests pass; black, ruff, codespell, and mypy all clean. Co-authored-by: kiloconnect[bot] <240665456+kiloconnect[bot]@users.noreply.github.com> --- modpods/estimator.py | 8 +++++ modpods/kernels.py | 9 ++++- modpods/lti.py | 13 +++++++- modpods/metrics.py | 24 ++++++++++---- modpods/train.py | 78 +++++++++++++++++++++++++++++++++++++------- 5 files changed, 113 insertions(+), 19 deletions(-) diff --git a/modpods/estimator.py b/modpods/estimator.py index e70e270..91ef610 100644 --- a/modpods/estimator.py +++ b/modpods/estimator.py @@ -110,6 +110,8 @@ def __init__( optimization_method: str = "bayesian", kernel: str | Any = "gamma", random_state: int | None = None, + fast_mode: bool = False, + optimizer_kwargs: dict | None = None, ) -> None: self.dependent_columns = dependent_columns self.independent_columns = independent_columns @@ -130,6 +132,8 @@ def __init__( self.optimization_method = optimization_method self.kernel = kernel self.random_state = random_state + self.fast_mode = fast_mode + self.optimizer_kwargs = optimizer_kwargs or {} self.estimators_: list[DelayIOModel] = [] def fit(self, system_data: pd.DataFrame, **kwargs: Any) -> list[DelayIOModel]: @@ -160,6 +164,8 @@ def fit(self, system_data: pd.DataFrame, **kwargs: Any) -> list[DelayIOModel]: optimization_method=self.optimization_method, kernel=self.kernel, seed=self.random_state, + fast_mode=self.fast_mode, + **self.optimizer_kwargs, **kwargs, ) @@ -233,6 +239,8 @@ def get_params(self, deep: bool = True) -> dict[str, Any]: "optimization_method": self.optimization_method, "kernel": self.kernel, "random_state": self.random_state, + "fast_mode": self.fast_mode, + "optimizer_kwargs": self.optimizer_kwargs, } def set_params(self, **params: Any) -> DelayIO: diff --git a/modpods/kernels.py b/modpods/kernels.py index c081663..aa9295a 100644 --- a/modpods/kernels.py +++ b/modpods/kernels.py @@ -260,7 +260,14 @@ def kernel_fn(self, t: np.ndarray, zeta: float, omega_n: float) -> np.ndarray: h = omega_n**2 * t * np.exp(-omega_n * t) else: s = omega_n * np.sqrt(zeta**2 - 1.0) - h = omega_n * np.exp(-zeta * omega_n * t) * np.sinh(s * t) / s + # Stable difference-of-exponentials form of sinh(s*t) * exp(-zeta*omega_n*t). + # For zeta > 1 both exponents are negative, so no overflow occurs. + h = ( + 0.5 + * omega_n + * (np.exp((s - zeta * omega_n) * t) - np.exp(-(s + zeta * omega_n) * t)) + / s + ) if zeta < 0: return h # type: ignore[no-any-return] return np.maximum(h, 0.0) # type: ignore[no-any-return] diff --git a/modpods/lti.py b/modpods/lti.py index 04ba579..450369d 100644 --- a/modpods/lti.py +++ b/modpods/lti.py @@ -642,6 +642,8 @@ def lti_system_gen( constraints=None, kernel="gamma", max_states=5, + fast_mode=False, + optimizer_kwargs=None, ): if _normalize_verbose(verbose) != "warnings": configure_verbosity(verbose) @@ -731,6 +733,8 @@ def lti_system_gen( kernel=kernel, max_states=max_states, constraints=constraints, + fast_mode=fast_mode, + optimizer_kwargs=optimizer_kwargs or {}, ) # we'll parse this delayed causation into the matrices A, B, and C later else: @@ -1110,6 +1114,8 @@ def __init__( forcing_coef_constraints: Any = None, constraints: Any = None, kernel: str = "gamma", + fast_mode: bool = False, + optimizer_kwargs: dict | None = None, ) -> None: self.causative_topology = causative_topology self.independent_columns = independent_columns @@ -1123,6 +1129,8 @@ def __init__( self.forcing_coef_constraints = forcing_coef_constraints self.constraints = constraints self.kernel = kernel + self.fast_mode = fast_mode + self.optimizer_kwargs = optimizer_kwargs or {} self.system_: Any = None self.A_: pd.DataFrame | None = None self.B_: pd.DataFrame | None = None @@ -1145,8 +1153,9 @@ def fit(self, system_data: pd.DataFrame, **kwargs: Any) -> "LTISystem": early_stopping_threshold=self.early_stopping_threshold, verbose=self.verbose, forcing_coef_constraints=self.forcing_coef_constraints, - constraints=self.constraints, kernel=self.kernel, + fast_mode=self.fast_mode, + optimizer_kwargs=self.optimizer_kwargs, **kwargs, ) self.system_ = result["system"] @@ -1186,6 +1195,8 @@ def get_params(self, deep: bool = True) -> dict[str, Any]: "forcing_coef_constraints": self.forcing_coef_constraints, "constraints": self.constraints, "kernel": self.kernel, + "fast_mode": self.fast_mode, + "optimizer_kwargs": self.optimizer_kwargs, } def set_params(self, **params: Any) -> "LTISystem": diff --git a/modpods/metrics.py b/modpods/metrics.py index e782870..5fbc775 100644 --- a/modpods/metrics.py +++ b/modpods/metrics.py @@ -51,6 +51,15 @@ def compute_detailed_metrics( Returns: Dict with keys: MAE, RMSE, NSE, alpha, beta, HFV, HFV10, LFV, FDC. """ + # FDC uses log10 of sorted values; clip to a small positive floor so + # zero/negative predictions do not produce -inf/NaN. + log_floor = 1e-12 + + def _log10_sorted(arr: np.ndarray, q: float) -> float: + sorted_vals = np.sort(arr) + idx = int(q * len(sorted_vals)) + return float(np.log10(max(sorted_vals[idx], log_floor))) + n_cols = y_true.shape[1] mae = [] rmse = [] @@ -97,13 +106,16 @@ def compute_detailed_metrics( fdc.append( 100 * ( - np.log10(np.sort(y_pred[:, col_idx])[int(0.2 * len(y_pred))]) - - np.log10(np.sort(y_pred[:, col_idx])[int(0.7 * len(y_pred))]) - - np.log10(np.sort(y_true[:, col_idx])[int(0.2 * len(y_true))]) - + np.log10(np.sort(y_true[:, col_idx])[int(0.7 * len(y_true))]) + _log10_sorted(y_pred[:, col_idx], 0.2) + - _log10_sorted(y_pred[:, col_idx], 0.7) + - _log10_sorted(y_true[:, col_idx], 0.2) + + _log10_sorted(y_true[:, col_idx], 0.7) + ) + / max( + _log10_sorted(y_true[:, col_idx], 0.2) + - _log10_sorted(y_true[:, col_idx], 0.7), + log_floor, ) - / np.log10(np.sort(y_true[:, col_idx])[int(0.2 * len(y_true))]) - - np.log10(np.sort(y_true[:, col_idx])[int(0.7 * len(y_true))]) ) logger.info("MAE = %s", mae) diff --git a/modpods/train.py b/modpods/train.py index de14631..be8fcb5 100644 --- a/modpods/train.py +++ b/modpods/train.py @@ -67,8 +67,17 @@ def optimize( configure_verbosity(verbose) logger.info("Using Bayesian optimization...") - bayesian_max_iter = min(max_iter * 4, 200) - n_initial = min(30, max(20, int(bayesian_max_iter * 0.6))) + # Allow callers to override the Bayesian optimization budget directly. + # Defaults keep backward compatibility with the previous behavior + # (bayesian_max_iter = min(max_iter * 4, 200)). + bayesian_max_iter = int(optimizer_kwargs.get("n_calls", min(max_iter * 4, 200))) + n_initial = int( + optimizer_kwargs.get( + "n_initial_points", + min(30, max(20, int(bayesian_max_iter * 0.6))), + ) + ) + n_initial = min(n_initial, bayesian_max_iter) rng = np.random.default_rng(self.seed) if self.seed is not None else None X_sample_list: list[Any] = [] @@ -261,6 +270,7 @@ def __init__( optimization_method: str = "bayesian", seed: int | None = None, optimizer_kwargs: dict | None = None, + fast_mode: bool = False, ) -> None: self.kernel = kernel self.system_data = system_data @@ -269,7 +279,30 @@ def __init__( self.windup_timesteps = windup_timesteps self.init_transforms = init_transforms self.max_transforms = _auto_max_transforms(kernel, max_transforms) - self.max_iter = max_iter + # fast_mode trades accuracy for speed: use a cheaper optimizer budget and + # cap the number of transforms (2 is sufficient for most physical systems). + self.fast_mode = fast_mode + if self.fast_mode: + self.max_transforms = min(self.max_transforms, 2) + if optimization_method == "bayesian": + self.optimization_method = "differential_evolution" + else: + self.optimization_method = optimization_method + self.max_iter = max(1, max_iter // 10) + self.optimizer_kwargs = dict(optimizer_kwargs or {}) + self.optimizer_kwargs.setdefault("maxiter", self.max_iter) + self.optimizer_kwargs.setdefault("popsize", 5) + else: + self.max_iter = max_iter + self.optimization_method = optimization_method + self.optimizer_kwargs = dict(optimizer_kwargs or {}) + # Bayesian-specific kwargs are consumed by BayesianOptimizer and must + # not be forwarded to scipy optimizers (they would raise TypeError). + self._bayesian_kwargs = { + k: self.optimizer_kwargs.pop(k) + for k in ("n_calls", "n_initial_points") + if k in self.optimizer_kwargs + } self.poly_order = poly_order self.transform_dependent = transform_dependent self.verbose = verbose @@ -280,9 +313,7 @@ def __init__( self.forcing_coef_constraints = forcing_coef_constraints self.constraints = constraints self.early_stopping_threshold = early_stopping_threshold - self.optimization_method = optimization_method self.seed = seed - self.optimizer_kwargs = optimizer_kwargs or {} if transform_dependent: self.columns = system_data.columns.tolist() @@ -435,12 +466,17 @@ def _optimize_params(self, num_transforms: int) -> np.ndarray: ) objective = self._create_objective(transform_columns, num_transforms) optimizer = self._get_optimizer() + # Bayesian-specific kwargs (n_calls, n_initial_points) are merged in + # only for the Bayesian optimizer; scipy optimizers ignore them. + kwargs = dict(self.optimizer_kwargs) + if self.optimization_method == "bayesian": + kwargs.update(self._bayesian_kwargs) return optimizer.optimize( objective_function=objective, bounds=bounds, max_iter=self.max_iter, verbose=self.verbose, - optimizer_kwargs=self.optimizer_kwargs, + optimizer_kwargs=kwargs, ) def _update_kernel_params( @@ -563,6 +599,7 @@ def __init__( optimization_method: str = "bayesian", seed: int | None = None, optimizer_kwargs: dict | None = None, + fast_mode: bool = False, ) -> None: self.system_data = system_data self.dependent_columns = dependent_columns @@ -585,6 +622,7 @@ def __init__( self.optimization_method = optimization_method self.seed = seed self.optimizer_kwargs = optimizer_kwargs or {} + self.fast_mode = fast_mode self.all_results: dict[str, dict[int, dict[str, Any]]] = {} def _train_kernel( @@ -612,6 +650,7 @@ def _train_kernel( optimization_method=self.optimization_method, seed=self.seed, optimizer_kwargs=self.optimizer_kwargs, + fast_mode=self.fast_mode, ) return trainer.train() @@ -677,7 +716,9 @@ def delay_io_train( kernel="gamma", max_states=5, seed=None, - **optimizer_kwargs, + fast_mode=False, + optimizer_kwargs=None, + **extra_kwargs, ): """Train a delay-IO model with pluggable convolution kernels. @@ -698,9 +739,21 @@ def delay_io_train( max_states: Maximum state dimension for canonical LTI kernels (default 5). + fast_mode: When True, trades accuracy for speed by: + - Using a cheaper optimizer budget (max_iter // 10, popsize 5) + - Capping max_transforms at 2 (sufficient for most physical systems) + - Using differential_evolution instead of Bayesian optimization + Recommended for parameter sweeps and grid searches. + + optimizer_kwargs: Optional dict of extra optimizer-specific keyword + arguments. For Bayesian optimization this may include `n_calls` + and `n_initial_points` to control the optimization budget. + Returns: dict keyed by num_transforms. """ + optimizer_kwargs = dict(optimizer_kwargs or {}) + optimizer_kwargs.update(extra_kwargs) if kernel in ("try-all", "run-all"): trainer = MultiKernelTrainer( system_data=system_data, @@ -724,11 +777,12 @@ def delay_io_train( optimization_method=optimization_method, seed=seed, optimizer_kwargs=optimizer_kwargs, + fast_mode=fast_mode, ) return trainer.train() if kernel in ("canonical_lti", "canonical_lti_incremental"): - max_states = optimizer_kwargs.get("max_states", 5) + max_states = optimizer_kwargs.pop("max_states", 5) if kernel == "canonical_lti_incremental": k = get_kernel("canonical_lti_incremental") if hasattr(k, "max_states"): @@ -767,11 +821,12 @@ def delay_io_train( optimization_method=optimization_method, seed=seed, optimizer_kwargs=optimizer_kwargs, + fast_mode=fast_mode, ) return single_trainer.train() if kernel == "decoupled_lti": - max_states = optimizer_kwargs.get("max_states", 5) + max_states = optimizer_kwargs.pop("max_states", 5) k = get_kernel("decoupled_lti") if hasattr(k, "max_states"): k.max_states = max_states @@ -826,6 +881,7 @@ def delay_io_train( optimization_method=optimization_method, seed=seed, optimizer_kwargs=optimizer_kwargs, + fast_mode=fast_mode, ) return single_trainer.train() @@ -981,12 +1037,12 @@ def objective(params): result = opt.differential_evolution( objective, - bounds=bounds, + bounds=bounds, # type: ignore[arg-type] maxiter=self.max_iter, popsize=15, mutation=(0.5, 1.5), recombination=0.7, - seed=42 if self.seed is None else self.seed, + seed=42 if self.seed is None else self.seed, # type: ignore[call-arg] updating="deferred", )