draft scikitlearn interface

requirements.txt

@@ -2,12 +2,10 @@ numpy>=1.16.4
 scipy>=1.2.1
 scikit-learn>=0.22.2
 pandas>=0.24.2
-matplotlib>=3.0.3
 statsmodels>=0.11.1
 
-# scripts
-ttictoc>=0.4.1
-seaborn>=0.10.0
+# examples
+matplotlib>=3.0.3
 
 # documentation
 Sphinx>=2.1.2

setup.py

@@ -44,13 +44,11 @@ setup(
                       "scipy",
                       "pandas",
                       "scikit-learn",
-                      "seaborn",
-                      "matplotlib",
                       "statsmodels"],
     author=__author__,
     author_email=__email__,
     url="https://gitlab.com/hsr-iet/smartwaterproject/swconsumption",
-    description="Utilities to run water cnsumption models",
+    description="Utilities to run water consumption models",
     long_description=long_description,
     long_description_content_type='text/x-rst',
     scripts=[],

swconsumption/__init__.py

@@ -22,3 +22,4 @@ __author__ = "Juan Pablo Carbajal"
 __maintainer__ = ["Juan Pablo Carbajal"]
 __email__ = "ajuanpi+dev@gmail.com"
 __status__ = "Development"
+

swconsumption/models.py

+"""
+Models for consupmtion time series
+"""
+
+# Copyright (C) 2020 HSR Hochschule für Technik Rapperswil
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+# Author: Juan Pablo Carbajal <ajuanpi+dev@gmail.com>
+
+from collections import namedtuple
+from functools import partial
+
+import numpy as np
+import pandas as pd
+
+import statsmodels.api as sm
+from statsmodels.tsa.ar_model import AutoReg
+from statsmodels.tsa.deterministic import (
+    DeterministicProcess,
+    CalendarFourier,
+    Fourier,
+    TimeTrend
+)
+
+from sklearn.base import BaseEstimator
+from sklearn.utils.validation import check_X_y, check_array, check_is_fitted
+
+class ARX(BaseEstimator):
+    """ ARX estimator
+    
+    Parameters
+    ----------
+    demo_param : str, default='demo_param'
+        A parameter used for demonstation of how to pass and store paramters.
+
+    Examples
+    --------
+    >>> from skltemplate import TemplateEstimator
+    >>> import numpy as np
+    >>> X = np.arange(100).reshape(100, 1)
+    >>> y = np.zeros((100, ))
+    >>> estimator = TemplateEstimator()
+    >>> estimator.fit(X, y)
+    TemplateEstimator(demo_param='demo_param')
+    """
+
+    def __init__(self, lags=1, subperiods=1, datetime=False, period='Y'):
+        self.lags = lags
+        self.subperiods = subperiods
+        self.datetime = datetime
+        self.period = period
+
+    def fit(self, X, y):
+        """A reference implementation of a fitting function.
+
+        Parameters
+        ----------
+        X : {array-like, sparse matrix}, shape (n_samples, n_features)
+            The training input samples.
+        y : array-like, shape (n_samples,) or (n_samples, n_outputs)
+            The target values (class labels in classification, real numbers in
+            regression).
+
+        Returns
+        -------
+        self : object
+            Returns self.
+        """
+#        X, y = check_X_y(X, y, ensure_min_features=2,
+#                              ensure_min_samples=max(self.lags)+1,
+#                              dtype=('object','numeric'),
+#                              accept_sparse=True)
+
+        self._fdata = self._data2pandas(X, y)
+
+        tt = TimeTrend(constant=True)
+        if self.datetime:
+            seasonal = CalendarFourier(freq=self.period, order=self.subperiods)
+        else:
+            if isinstance(self.period, str):
+                period = pd.Timedelta(f'1{self.period}').total_seconds() / 3600 / 24
+            else:
+                period = self.period
+            seasonal = Fourier(period=period, order=self.subperiods)
+
+        dp = DeterministicProcess(self._fdata.index, 
+                                  additional_terms=[tt, seasonal])
+        self._model = partial(AutoReg, 
+                              lags=self.lags, trend="n", deterministic=dp)
+
+        Exog = self._fdata[self._exognam]
+        self._fmodel = self._model(self._fdata.Q, exog=Exog).fit()
+
+        self.is_fitted_ = True
+        # `fit` should always return `self`
+        return self
+
+    def predict(self, X):
+        """ A reference implementation of a predicting function.
+
+        Parameters
+        ----------
+        X : {array-like, sparse matrix}, shape (n_samples, n_features)
+            The training input samples.
+
+        Returns
+        -------
+        y : ndarray, shape (n_samples,)
+            Returns an array of ones.
+        """
+        #X = check_array(X, ensure_min_features=2, dtype='numeric', accept_sparse=True)
+        check_is_fitted(self, 'is_fitted_')
+        
+        X_ = self._data2pandas(X)
+        oos = (X_.index > self._fdata.index[-1])
+
+        #import pdb; pdb.set_trace()
+        if self.datetime:
+            ini = X_.index[0]
+            end = X_.index[-1]
+        else:
+            ini = 0
+            end = len(X_.index)-1
+        Exog = X_[self._exognam][oos]
+        y_ = self._fmodel.predict(
+                        start=ini, 
+                        end=end,
+                        #exog= X_.Temp[~oos],
+                        exog_oos=Exog)
+        return y_
+
+    def _data2pandas(self, X, y=None):
+        """ organize inputs into dataframe """
+        self._exognam = [f'Ex{i}' for i in range(1,X.shape[1])]
+
+        df = pd.DataFrame(
+                        data=
+                        {k:X[:,i].astype(float) for i,k in enumerate(self._exognam,1)}, 
+                        index=X[:,0])
+        if self.datetime:
+            df.index = pd.DatetimeIndex(df.index, freq='infer')
+        df.sort_index(inplace=True)
+
+        if y is not None:
+            df['Q'] = y
+
+        return df
+
+def arx_tTQ(time, T, Q, *, lags=None, yearly_periods : int =24):
+    """ Autoregressive model with exogenous variables
+
+        Parameters
+        ----------
+        time:
+            Time vector
+        T:
+            Temperature data
+        Q:
+            Consumption data
+        lags: array_like, optional
+            lags for the AR model. Default [1, 5]
+        yearly_periods: integer, optional
+            Subdivision of the annual seasonal component
+
+        Returns
+        -------
+        :obj:`statsmodels.tsa.ar_model.AutoReg`
+            The model structure
+    """
+    lags = [1,5] if lags is None else lags
+
+    tt = TimeTrend(constant=True)
+    seasonal = CalendarFourier(freq='Y', order=yearly_periods)
+    dp = DeterministicProcess(time, additional_terms=[tt, seasonal])
+    mod = AutoReg(Q, lags, exog=T, trend="n", deterministic=dp)
+    return mod
+
+def mean_positive_error(y_true, y_pred):
+    """ Mean positive error
+    """
+    try:
+        resid = y_true - y_pred
+    except ValueError:
+        import pdb; pdb.set_trace()
+    return np.mean(resid[resid>0])
+
+def max_positive_error(y_true, y_pred):
+    """ Maximum positive error
+    """
+    resid = y_true - y_pred
+    return np.max(resid)