Fix all open issues: #69, #70, #71

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:

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]

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":

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)

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",
         )