diff --git a/sysidentpy/basis_function/_basis_function.py b/sysidentpy/basis_function/_basis_function.py index 05541fed..4d850f89 100644 --- a/sysidentpy/basis_function/_basis_function.py +++ b/sysidentpy/basis_function/_basis_function.py @@ -45,6 +45,9 @@ def fit( self, data: np.ndarray, max_lag: int = 1, + ylag: int = 1, + xlag: int = 1, + model_type: str = "NARMAX", predefined_regressors: Optional[np.ndarray] = None, ): """Build the Polynomial information matrix. @@ -59,6 +62,12 @@ def fit( The lagged matrix built with respect to each lag and column. max_lag : int Target data used on training phase. + ylag : ndarray of int + The range of lags according to user definition. + xlag : ndarray of int + The range of lags according to user definition. + model_type : str + The type of the model (NARMAX, NAR or NFIR). predefined_regressors : ndarray of int The index of the selected regressors by the Model Structure Selection algorithm. @@ -70,7 +79,7 @@ def fit( """ # Create combinations of all columns based on its index - iterable_list = range(data.shape[1]) + iterable_list = self.get_iterable_list(ylag, xlag, model_type) combinations = list(combinations_with_replacement(iterable_list, self.degree)) if predefined_regressors is not None: combinations = [combinations[index] for index in predefined_regressors] @@ -88,6 +97,9 @@ def transform( self, data: np.ndarray, max_lag: int = 1, + ylag: int=1, + xlag: int=1, + model_type: str = "NARMAX", predefined_regressors: Optional[np.ndarray] = None, ): """Build Polynomial Basis Functions. @@ -98,6 +110,12 @@ def transform( The lagged matrix built with respect to each lag and column. max_lag : int Maximum lag of list of regressors. + ylag : ndarray of int + The range of lags according to user definition. + xlag : ndarray of int + The range of lags according to user definition. + model_type : str + The type of the model (NARMAX, NAR or NFIR). predefined_regressors: ndarray Regressors to be filtered in the transformation. @@ -107,7 +125,7 @@ def transform( Transformed array. """ - return self.fit(data, max_lag, predefined_regressors) + return self.fit(data, max_lag, ylag, xlag, model_type, predefined_regressors) class Fourier(BaseBasisFunction): @@ -149,6 +167,9 @@ def fit( self, data: np.ndarray, max_lag: int = 1, + ylag: int = 1, + xlag: int = 1, + model_type: str = "NARMAX", predefined_regressors: Optional[np.ndarray] = None, ): """Build the Polynomial information matrix. @@ -163,6 +184,12 @@ def fit( The lagged matrix built with respect to each lag and column. max_lag : int Target data used on training phase. + ylag : ndarray of int + The range of lags according to user definition. + xlag : ndarray of int + The range of lags according to user definition. + model_type : str + The type of the model (NARMAX, NAR or NFIR). predefined_regressors : ndarray of int The index of the selected regressors by the Model Structure Selection algorithm. @@ -175,7 +202,7 @@ def fit( """ # remove intercept (because the data always have the intercept) if self.degree > 1: - data = Polynomial().fit(data, max_lag, predefined_regressors=None) + data = Polynomial().fit(data, max_lag, ylag, xlag, model_type, predefined_regressors=None) data = data[:, 1:] else: data = data[max_lag:, 1:] @@ -205,6 +232,9 @@ def transform( self, data: np.ndarray, max_lag: int = 1, + ylag: int = 1, + xlag: int = 1, + model_type: str = "NARMAX", predefined_regressors: Optional[np.ndarray] = None, ): """Build Fourier Basis Functions. @@ -215,6 +245,12 @@ def transform( The lagged matrix built with respect to each lag and column. max_lag : int Maximum lag of list of regressors. + ylag : ndarray of int + The range of lags according to user definition. + xlag : ndarray of int + The range of lags according to user definition. + model_type : str + The type of the model (NARMAX, NAR or NFIR). predefined_regressors: ndarray Regressors to be filtered in the transformation. @@ -224,7 +260,7 @@ def transform( Transformed array. """ - return self.fit(data, max_lag, predefined_regressors) + return self.fit(data, max_lag, ylag, xlag, model_type, predefined_regressors) class Bersntein(BaseBasisFunction): @@ -294,11 +330,14 @@ def fit( self, data: np.ndarray, max_lag: int = 1, + ylag: int = 1, + xlag: int = 1, + model_type: str = "NARMAX", predefined_regressors: Optional[np.ndarray] = None, ): # remove intercept (because the data always have the intercept) if self.degree > 1: - data = Polynomial().fit(data, max_lag, predefined_regressors=None) + data = Polynomial().fit(data, max_lag, ylag, xlag, model_type, predefined_regressors=None) data = data[:, 1:] else: data = data[max_lag:, 1:] @@ -325,6 +364,169 @@ def transform( self, data: np.ndarray, max_lag: int = 1, + ylag: int = 1, + xlag: int = 1, + model_type: str = "NARMAX", predefined_regressors: Optional[np.ndarray] = None, ): - return self.fit(data, max_lag, predefined_regressors) + """Build Bersntein Basis Functions. + + Parameters + ---------- + data : ndarray of floats + The lagged matrix built with respect to each lag and column. + max_lag : int + Maximum lag of list of regressors. + ylag : ndarray of int + The range of lags according to user definition. + xlag : ndarray of int + The range of lags according to user definition. + predefined_regressors: ndarray + Regressors to be filtered in the transformation. + + Returns + ------- + X_tr : {ndarray, sparse matrix} of shape (n_samples, n_features) + Transformed array. + + """ + return self.fit(data, max_lag, ylag, xlag, model_type, predefined_regressors) + +class Bilinear(BaseBasisFunction): + r"""Build Bilinear basis function. + + A general bilinear input-output model takes the form + + .. math:: + y(k) = a_0 + \sum_{i=1}^{n_y} a_i y(k-i) + \sum_{i=1}^{n_u} b_i u(k-i) + + \sum_{i=1}^{n_y} \sum_{j=1}^{n_u} c_{ij} y(k-i) u(k-j) + + This is a special case of the Polynomial NARMAX model. + + Bilinear system theory has been widely studied and it plays an important role in the context of continuous-time + systems. This is because, roughly speaking, the set of bilinear + systems is dense in the space of continuous-time systems and any continuous causal + functional can be arbitrarily well approximated by bilinear systems within any + bounded time interval (see for example Fliess and Normand-Cyrot 1982). Moreover, + many real continuous-time processes are naturally in bilinear form. A few examples + are distillation columns (España and Landau 1978), nuclear and thermal control + processes (Mohler 1973). + + Sampling the continuous-time bilinear system, however, produces a NARMAX model + which is more complex than a discrete-time bilinear model. + + Parameters + ---------- + degree : int (max_degree), default=2 + The maximum degree of the polynomial features. + + Notes + ----- + Be aware that the number of features in the output array scales + significantly as the number of inputs, the max lag of the input and output, and + degree increases. High degrees can cause overfitting. + """ + + + def __init__( + self, + degree: int = 2, + ): + self.degree = degree + + def fit( + self, + data: np.ndarray, + max_lag: int = 1, + ylag: int = 1, + xlag: int = 1, + model_type: str = "NARMAX", + predefined_regressors: Optional[np.ndarray] = None, + ): + """Build the Bilinear information matrix. + + Each columns of the information matrix represents a candidate + regressor. The set of candidate regressors are based on xlag, + ylag, and degree defined by the user. + + Parameters + ---------- + data : ndarray of floats + The lagged matrix built with respect to each lag and column. + max_lag : int + Target data used on training phase. + ylag : ndarray of int + The range of lags according to user definition. + xlag : ndarray of int + The range of lags according to user definition. + model_type : str + The type of the model (NARMAX, NAR or NFIR). + predefined_regressors : ndarray of int + The index of the selected regressors by the Model Structure + Selection algorithm. + + Returns + ------- + psi = ndarray of floats + The lagged matrix built in respect with each lag and column. + + """ + # Create combinations of all columns based on its index + iterable_list = self.get_iterable_list(ylag, xlag, model_type) + combinations = list(combinations_with_replacement(iterable_list, self.degree)) + if self.degree == 1: + Warning('You choose a bilinear basis function and nonlinear degree = 1. In this case, you have a linear polynomial model.') + else: + ny = self.get_max_ylag(ylag) + nx = self.get_max_xlag(xlag) + combination_ylag = list(combinations_with_replacement(list(range(1, ny + 1)), self.degree)) + combination_xlag = list(combinations_with_replacement(list(range(ny + 1, nx + ny + 1)), self.degree)) + combinations_xy = combination_xlag + combination_ylag + combinations = list(set(combinations)-set(combinations_xy)) + + if predefined_regressors is not None: + combinations = [combinations[index] for index in predefined_regressors] + + + psi = np.column_stack( + [ + np.prod(data[:, combinations[i]], axis=1) + for i in range(len(combinations)) + ] + ) + psi = psi[max_lag:, :] + return psi + + def transform( + self, + data: np.ndarray, + max_lag: int = 1, + ylag: int = 1, + xlag: int = 1, + model_type: str = "NARMAX", + predefined_regressors: Optional[np.ndarray] = None, + ): + """Build Polynomial Basis Functions. + + Parameters + ---------- + data : ndarray of floats + The lagged matrix built with respect to each lag and column. + max_lag : int + Maximum lag of list of regressors. + ylag : ndarray of int + The range of lags according to user definition. + xlag : ndarray of int + The range of lags according to user definition. + model_type : str + The type of the model (NARMAX, NAR or NFIR). + predefined_regressors: ndarray + Regressors to be filtered in the transformation. + + Returns + ------- + X_tr : {ndarray, sparse matrix} of shape (n_samples, n_features) + Transformed array. + + """ + return self.fit(data, max_lag, ylag, xlag, model_type, predefined_regressors) diff --git a/sysidentpy/basis_function/basis_function_base.py b/sysidentpy/basis_function/basis_function_base.py index e5414299..19a37cf4 100644 --- a/sysidentpy/basis_function/basis_function_base.py +++ b/sysidentpy/basis_function/basis_function_base.py @@ -1,5 +1,6 @@ """Base class for Basis Function.""" +from itertools import combinations_with_replacement, chain from abc import ABCMeta, abstractmethod from typing import Optional @@ -18,6 +19,8 @@ def fit( self, data: np.ndarray, max_lag: int = 1, + ylag: int = 1, + xlag: int = 1, predefined_regressors: Optional[np.ndarray] = None, ): """Abstract method.""" @@ -27,6 +30,79 @@ def transform( self, data: np.ndarray, max_lag: int = 1, + ylag: int = 1, + xlag: int = 1, predefined_regressors: Optional[np.ndarray] = None, ) -> np.ndarray: """Abstract methods.""" + + def get_max_ylag(self, ylag: int = 1): + """Get maximum ylag. + + Parameters + ---------- + ylag : ndarray of int + The range of lags according to user definition. + + Returns + ------- + ny : list + Maximum value of ylag. + + """ + ny = np.max(list(chain.from_iterable([[ylag]]))) + return ny + + def get_max_xlag(self, xlag: int = 1): + """Get maximum xlag. + + Parameters + ---------- + xlag : ndarray of int + The range of lags according to user definition. + + Returns + ------- + nx : list + Maximum value of xlag. + + """ + nx = np.max(list(chain.from_iterable([[np.array(xlag, dtype=object)]]))) + return nx + + def get_iterable_list( + self, + ylag: int = 1, + xlag: int = 1, + model_type: str = "NARMAX" + ): + """Get iterable list. + + Parameters + ---------- + ylag : ndarray of int + The range of lags according to user definition. + xlag : ndarray of int + The range of lags according to user definition. + model_type : str + The type of the model (NARMAX, NAR or NFIR). + + Returns + ------- + iterable_list : list + List of tuples of the regressor combinations. + + """ + if model_type == "NARMAX": + ny = self.get_max_ylag(ylag) + nx = self.get_max_xlag(xlag) + iterable_list = list(range(ny + nx + 1)) + combinations = list(combinations_with_replacement(iterable_list, self.degree)) + elif model_type == "NAR": + ny = self.get_max_ylag(ylag) + iterable_list = list(range(ny + 1)) + else: + nx = self.get_max_xlag(xlag) + iterable_list = list(range(nx + 1)) + return iterable_list + diff --git a/sysidentpy/general_estimators/narx.py b/sysidentpy/general_estimators/narx.py index df81147e..576721ad 100644 --- a/sysidentpy/general_estimators/narx.py +++ b/sysidentpy/general_estimators/narx.py @@ -160,7 +160,7 @@ def fit(self, *, X=None, y=None): self.max_lag = self._get_max_lag() lagged_data = self.build_matrix(X, y) reg_matrix = self.basis_function.fit( - lagged_data, self.max_lag, predefined_regressors=None + lagged_data, self.max_lag, self.ylag, self.xlag, self.model_type, predefined_regressors=None ) if X is not None: @@ -265,7 +265,9 @@ def _one_step_ahead_prediction(self, X, y): X_base = self.basis_function.transform( lagged_data, self.max_lag, - # predefined_regressors=self.pivv[: len(self.final_model)], + self.ylag, + self.xlag, + self.model_type ) yhat = self.base_estimator.predict(X_base) @@ -485,7 +487,9 @@ def _basis_function_predict(self, X, y_initial, forecast_horizon=None): X_tmp = self.basis_function.transform( lagged_data, self.max_lag, - # predefined_regressors=self.pivv[: len(self.final_model)], + self.ylag, + self.xlag, + self.model_type ) a = self.base_estimator.predict(X_tmp) diff --git a/sysidentpy/model_structure_selection/accelerated_orthogonal_least_squares.py b/sysidentpy/model_structure_selection/accelerated_orthogonal_least_squares.py index 7e54f0b6..77ec0c94 100644 --- a/sysidentpy/model_structure_selection/accelerated_orthogonal_least_squares.py +++ b/sysidentpy/model_structure_selection/accelerated_orthogonal_least_squares.py @@ -308,7 +308,7 @@ def fit(self, *, X: Optional[np.ndarray] = None, y: Optional[np.ndarray] = None) self.max_lag = self._get_max_lag() lagged_data = self.build_matrix(X, y) reg_matrix = self.basis_function.fit( - lagged_data, self.max_lag, predefined_regressors=None + lagged_data, self.max_lag, self.ylag, self.xlag, self.model_type, predefined_regressors=None ) if X is not None: @@ -424,6 +424,9 @@ def _one_step_ahead_prediction( X_base = self.basis_function.transform( lagged_data, self.max_lag, + self.ylag, + self.xlag, + self.model_type, predefined_regressors=self.pivv[: len(self.final_model)], ) diff --git a/sysidentpy/model_structure_selection/entropic_regression.py b/sysidentpy/model_structure_selection/entropic_regression.py index a6077968..4c6db987 100644 --- a/sysidentpy/model_structure_selection/entropic_regression.py +++ b/sysidentpy/model_structure_selection/entropic_regression.py @@ -584,7 +584,7 @@ def fit(self, *, X=None, y=None): lagged_data = self.build_matrix(X, y) reg_matrix = self.basis_function.fit( - lagged_data, self.max_lag, predefined_regressors=None + lagged_data, self.max_lag, self.ylag, self.xlag, self.model_type, predefined_regressors=None ) if X is not None: @@ -739,6 +739,9 @@ def _one_step_ahead_prediction(self, X, y): X_base = self.basis_function.transform( lagged_data, self.max_lag, + self.ylag, + self.xlag, + self.model_type, predefined_regressors=self.pivv[: len(self.final_model)], ) diff --git a/sysidentpy/model_structure_selection/forward_regression_orthogonal_least_squares.py b/sysidentpy/model_structure_selection/forward_regression_orthogonal_least_squares.py index 571a8f3c..5ee4eb84 100644 --- a/sysidentpy/model_structure_selection/forward_regression_orthogonal_least_squares.py +++ b/sysidentpy/model_structure_selection/forward_regression_orthogonal_least_squares.py @@ -602,7 +602,7 @@ def fit(self, *, X: Optional[np.ndarray] = None, y: np.ndarray): lagged_data = self.build_matrix(X, y) reg_matrix = self.basis_function.fit( - lagged_data, self.max_lag, predefined_regressors=None + lagged_data, self.max_lag, self.ylag, self.xlag, self.model_type, predefined_regressors=None ) if X is not None: @@ -743,6 +743,9 @@ def _one_step_ahead_prediction( X_base = self.basis_function.transform( lagged_data, self.max_lag, + self.ylag, + self.xlag, + self.model_type, predefined_regressors=self.pivv[: len(self.final_model)], ) diff --git a/sysidentpy/model_structure_selection/meta_model_structure_selection.py b/sysidentpy/model_structure_selection/meta_model_structure_selection.py index 405fbbf9..ec5b4174 100644 --- a/sysidentpy/model_structure_selection/meta_model_structure_selection.py +++ b/sysidentpy/model_structure_selection/meta_model_structure_selection.py @@ -385,7 +385,7 @@ def evaluate_objective_function( lagged_data = self.build_matrix(X_train, y_train) psi = self.basis_function.fit( - lagged_data, self.max_lag, predefined_regressors=self.pivv + lagged_data, self.max_lag, self.xlag, self.ylag, self.model_type, predefined_regressors=self.pivv ) pos_insignificant_terms, _, _ = self.perform_t_test( diff --git a/sysidentpy/multiobjective_parameter_estimation/estimators.py b/sysidentpy/multiobjective_parameter_estimation/estimators.py index 1d11d045..f76710bb 100644 --- a/sysidentpy/multiobjective_parameter_estimation/estimators.py +++ b/sysidentpy/multiobjective_parameter_estimation/estimators.py @@ -384,7 +384,7 @@ def build_psi(self, X: np.ndarray, y: np.ndarray) -> np.ndarray: ) psi = Polynomial(degree=self.degree).fit( - lagged_data, max_lag=self.max_lag, predefined_regressors=pivv + lagged_data, max_lag=self.max_lag, ylag=ylag, xlag=xlag, predefined_regressors=pivv ) return psi diff --git a/sysidentpy/narmax_base.py b/sysidentpy/narmax_base.py index ddfbce7e..7738f3f0 100644 --- a/sysidentpy/narmax_base.py +++ b/sysidentpy/narmax_base.py @@ -902,6 +902,9 @@ def _basis_function_predict( X_tmp = self.basis_function.transform( lagged_data, self.max_lag, + self.ylag, + self.xlag, + self.model_type, predefined_regressors=self.pivv[: len(self.final_model)], ) diff --git a/sysidentpy/neural_network/narx_nn.py b/sysidentpy/neural_network/narx_nn.py index d1a4374e..ee0dc3eb 100644 --- a/sysidentpy/neural_network/narx_nn.py +++ b/sysidentpy/neural_network/narx_nn.py @@ -274,12 +274,12 @@ def split_data(self, X, y): basis_name = self.basis_function.__class__.__name__ if basis_name == "Polynomial": reg_matrix = self.basis_function.fit( - lagged_data, self.max_lag, predefined_regressors=None + lagged_data, self.max_lag, self.ylag, self.xlag, self.model_type, predefined_regressors=None ) reg_matrix = reg_matrix[:, 1:] else: reg_matrix = self.basis_function.fit( - lagged_data, self.max_lag, predefined_regressors=None + lagged_data, self.max_lag, self.ylag, self.xlag, self.model_type, predefined_regressors=None ) if X is not None: @@ -518,12 +518,18 @@ def _one_step_ahead_prediction(self, X, y): X_base = self.basis_function.transform( lagged_data, self.max_lag, + self.ylag, + self.xlag, + self.model_type ) X_base = X_base[:, 1:] else: X_base = self.basis_function.transform( lagged_data, self.max_lag, + self.ylag, + self.xlag, + self.model_type ) yhat = np.zeros(X.shape[0], dtype=float) @@ -712,6 +718,9 @@ def _basis_function_predict(self, X, y_initial, forecast_horizon=None): X_tmp = self.basis_function.transform( lagged_data, self.max_lag, + self.ylag, + self.xlag, + self.model_type ) X_tmp = np.atleast_1d(X_tmp).astype(np.float32) yhat = yhat.astype(np.float32) diff --git a/sysidentpy/parameter_estimation/estimators_base.py b/sysidentpy/parameter_estimation/estimators_base.py index 90a45a4c..90ae8df9 100644 --- a/sysidentpy/parameter_estimation/estimators_base.py +++ b/sysidentpy/parameter_estimation/estimators_base.py @@ -88,7 +88,7 @@ def unbiased_estimator( lagged_data = im.build_output_matrix(None, e) e_regressors = basis_function.fit( - lagged_data, max_lag, predefined_regressors=None + lagged_data, max_lag, ylag=elag, predefined_regressors=None ) psi_extended = np.concatenate([psi, e_regressors], axis=1) diff --git a/sysidentpy/simulation/_simulation.py b/sysidentpy/simulation/_simulation.py index e0c797b8..f8ce542b 100644 --- a/sysidentpy/simulation/_simulation.py +++ b/sysidentpy/simulation/_simulation.py @@ -277,7 +277,7 @@ def simulate( self.max_lag = self._get_max_lag() lagged_data = self.build_matrix(X_train, y_train) psi = self.basis_function.fit( - lagged_data, self.max_lag, predefined_regressors=self.pivv + lagged_data, self.max_lag, self.ylag, self.xlag, self.model_type, predefined_regressors=self.pivv ) self.theta = self.estimator.optimize( @@ -303,7 +303,7 @@ def simulate( self.max_lag = self._get_max_lag() lagged_data = self.build_matrix(X_train, y_train) psi = self.basis_function.fit( - lagged_data, self.max_lag, predefined_regressors=self.pivv + lagged_data, self.max_lag, self.ylag, self.xlag, self.model_type, predefined_regressors=self.pivv ) _, self.err, _, _ = self.error_reduction_ratio( @@ -485,6 +485,9 @@ def _one_step_ahead_prediction(self, X, y): X_base = self.basis_function.transform( lagged_data, self.max_lag, + self.ylag, + self.xlag, + self.model_type, predefined_regressors=self.pivv[: len(self.final_model)], )