Multistep K-Forecast, changed description and input parameters for L …

src/py_ssa_lib/MSSA.py

@@ -338,51 +338,71 @@ def estimate_LRR(self, idx_components):
         else :
             print(nu_sq)
 
-    def L_Forecast(self, ts, M, idx_components, mode='forward'):
+    def L_Forecast(self, ts, M, idx_components, return_full=True):
         """
-                    Forecasts  M values for a given time series using LRR 
+                    Forecasts  M values for a given time series using LRR and U space 
                     Parameters
                     ----------
                     idx_components:list or numpy.arange of positive integer numbers, denotes the indices of elementary components 
                     ts: numpy.array, input time series 
                     M:int, number of  values to forecast or estimate
-                    mode:str, forecasts M future values for S time series if mode is "forward", 
+                    return_full:bool, if True returns  original time series + M forecasted values,
 
 
                     Returns
                     -------
-                    y_pred: numpy.array, original time series + M forecasted values, or original time series, where the last M values are estimated
+                    y_pred: numpy.array, (original time series +) M forecasted values,
         """
         R = self.estimate_LRR(idx_components)
         R = R.reshape(-1,1)
-        if M<=0:
-            y_pred = ts[:,:self.N]
-            return y_pred
-
-        if mode == 'forward':
-            y_pred = np.zeros((ts.shape[0],self.N+M))
-            y_pred[:,:self.N] = ts[:,:self.N]
-
-            for m in range(0,M):
-                    y_pred[:,self.N+m] = (y_pred[:,self.N-self.L+m+1:y_pred.shape[1]-M+m]@R).flatten()
+
+        y_pred = np.zeros((min(ts.shape),self.N+M))
+        y_pred[:,:self.N] = ts[:,:self.N]
 
-        #elif mode == 'retrospective':
-            # This method is no longer supported
-            #y_pred = np.zeros((ts.shape[0],self.N))
-            #y_pred[:,:self.N-M] = ts[:,:self.N-M]
-
-            #for m in range(0,M):
-
-                    #y_pred[:,self.N+m-M] = (y_pred[:,self.N-self.L+m+1-M:y_pred.shape[1]-M+m] @ R).flatten()
-
-
+        for m in range(0,M):
+                y_pred[:,self.N+m] = (y_pred[:,self.N-self.L+m+1:y_pred.shape[1]-M+m]@R).flatten()
+
+        if return_full:
+            return y_pred
         else:
-
-            raise  ValueError('Wrong type of the forecasting mode')
+            return y_pred[:,-M:]
 
+    def K_forecast(self, ts, M, idx_components,  return_full=True):
+          """"
+                    Forecasts  M values for a given time series using V eigenspace
+                    Parameters
+                    ----------
+                    idx_components:list or numpy.arange of positive integer numbers, denotes the indices of elementary components 
+                    ts: numpy.array, input time series 
+                    M:int, number of  values to forecast or estimate
+                    return_full:bool, if True returns  original time series + M forecasted values,
+
+
+                    Returns
+                    -------
+                    y_pred: numpy.array, (original time series +) M forecasted values,
+        """
 
+        W = self.V[:,idx_components]
+        W = W[[(self.K) * s for s in range(0, min(ts.shape))],: ]
+
+        Q = self.V[:,idx_components]
+        Q = np.delete(Q,[(self.K) * s for s in range(0, min(ts.shape))],0 )
+
+        inverse_m = np.linalg.inv(np.eye(min(ts.shape)) - W @ W.T )
+
+        y_pred = np.zeros((min(ts.shape),self.N+M))
+        y_pred[:,:self.N] = ts 
+        for m in range(0,M): 
 
-        return y_pred
+            Z = y_pred[:, self.N-self.K+m+1: self.N+m].flatten()
+            y_pred[:,self.N + m] = inverse_m @ W @ Q.T @ Z
+
+        if return_full:
+            return y_pred
+        else:
+            return y_pred[:,-M:]
+
 
     def estimate_ESPRIT(self, idx_components=[0], decompose_rho_omega=False):
 

src/py_ssa_lib/SSA.py

@@ -373,75 +373,73 @@ def estimate_LRR(self, idx_components):
             return R
         else :
             print(nu_sq)
-    
-    def L_Forecast(self, ts, M, idx_components, mode='forward'):
-        """
-                    Forecasts M values for a given time series using LRR (L-Forecasting)
+
+    def K_forecast(self, ts, idx_components, M, return_full=True):
+        """"
+                    Forecasts  M values for a given time series using V eigenspace
                     Parameters
                     ----------
                     idx_components:list or numpy.arange of positive integer numbers, denotes the indices of elementary components 
-                    ts: numpy.array, input time series 
+                    ts: numpy.array, single input time series . Check that its dimension is 1xN, or apply np.reshape(1,-1) on the input array
                     M:int, number of  values to forecast or estimate
-                    mode:str, forecasts M future values for S time series if mode is "forward"
+                    return_full:bool, if True returns  original time series + M forecasted values,
+
+
                     Returns
                     -------
-                    y_pred: numpy.array, original time series + M forecasted values, or original time series, where the last M values are estimated
+                    y_pred: numpy.array, (original time series +) M forecasted values,
         """
-        R = self.estimate_LRR(idx_components)
-        R = R.reshape(-1,1)
-        if M<=0:
-            y_pred = ts[:self.N]
-            return y_pred
-
-        if mode == 'forward':
-            y_pred = np.zeros((1,self.N+M))
-            y_pred[0,:self.N] = ts[:self.N]
-            for m in range(0,M):
-
-                    y_pred[:,self.N+m] = (y_pred[:,self.N-self.L+m+1:y_pred.shape[0]-M+m-1]) @ R
-        #elif mode == 'retrospective':
-            #this is an deprecated method
-            #y_pred = np.zeros((1,self.N))
-            #y_pred[0,:self.N-M] = ts[:self.N-M]
-            #for m in range(0,M):
-
-                    #y_pred[:,self.N+m-M] = (y_pred[:,self.N-self.L+m+1-M:y_pred.shape[0]-M+m-1]) @ R
+        W = self.V[:,idx_components]
+        W = W[[(self.K) * s for s in range(0, min(ts.shape))]]
 
-
-        else:
-
-            raise  ValueError('Wrong type of the forecasting mode')
-
+        Q = self.V[:,idx_components]
+        Q = np.delete(Q,[(self.K) * s for s in range(0, min(ts.shape))],0 )
+
+        inverse_m = np.linalg.inv(np.eye(min(ts.shape)) - W @ W.T )
 
+        y_pred = np.zeros((min(ts.shape),self.N+M))
+        y_pred[:,:self.N] = ts 
+        for m in range(0,M): 
 
-        return y_pred
+            Z = y_pred[:, self.N-self.K+m+1: self.N+m].flatten()
+            y_pred[:,self.N + m] = inverse_m @ W @ Q.T @ Z
 
-    def K_Forecast(model, ts, M, idx_components):
+        if return_full:
+            return y_pred
+        else:
+            return y_pred[:,-M:]
+
+    def L_Forecast(self, ts, M, idx_components, return_full=True):
         """
-                    Forecasts M values for a given time series using K-Forecasting 
+                    Forecasts  M values for a given time series using LRR and U space 
                     Parameters
                     ----------
                     idx_components:list or numpy.arange of positive integer numbers, denotes the indices of elementary components 
                     ts: numpy.array, input time series 
                     M:int, number of  values to forecast or estimate
+                    return_full:bool, if True returns  original time series + M forecasted values,
 
+
                     Returns
                     -------
-                    Z: numpy.array, original time series + M forecasted values, or original time series, where the last M values are estimated
+                    y_pred: numpy.array, (original time series +) M forecasted values,
         """
-
+        R = self.estimate_LRR(idx_components)
+        R = R.reshape(-1,1)
 
-        print(Warning("This function is under development! \nIt might contain some errors due to indexing of the signal space!"))
-        Q = np.delete(model.V[:, idx_components], model.K-1, axis=0)
-        W = model.V[model.K-1, idx_components].reshape(1,-1)
-        Inverse = np.linalg.inv(np.eye(1) - W @ W.T) 
-        Z = ts[-model.K+1:]
-        Z = np.append(Z, np.zeros((M))).reshape(-1,1)
-
+        y_pred = np.zeros((1,self.N+M))
+        y_pred[0,:self.N] = ts[:self.N]
         for m in range(0,M):
-            Z[-M+m,:] = Inverse @ W @ Q.T @ Z[-M+m - model.K+1:-M + m,:]
-
-        return Z.T
+
+                y_pred[:,self.N+m] = (y_pred[:,self.N-self.L+m+1:y_pred.shape[0]-M+m-1]) @ R
+
+
+        if return_full:
+            return y_pred
+        else:
+            return y_pred[:,-M:]
+
+
 
     def estimate_ESPRIT(self, idx_components=[0], decompose_rho_omega=False):
         """

test/.ipynb_checkpoints/test_MSSA-checkpoint.py

@@ -0,0 +1,102 @@
+import sys
+import py_ssa_lib
+import pytest
+import pandas as pd 
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy as sp
+from sklearn.utils.extmath import randomized_svd
+
+from py_ssa_lib.MSSA import MSSA
+from py_ssa_lib.datasets.load_energy_consumption_df import load_energy_consumption_df
+
+
+def test_mssa_init():
+
+
+            mssa_test = MSSA(Verbose=True)
+
+            assert mssa_test.Verbose == True
+
+
+def test_mssa_fit():
+
+            i_start_ts = 2
+            decomposition = "svd"
+            #decomposition = "rand_svd"
+            window_s = 7
+            demo_df = load_energy_consumption_df()
+
+            mssa_test = MSSA(Verbose=True)
+            mssa_test.fit(df=demo_df,decomposition=decomposition, idx_start_ts=i_start_ts, L=window_s)
+
+            assert (int(mssa_test.L) == window_s) 
+            assert  (mssa_test.decomposition == decomposition )
+            assert  (int(mssa_test.N) == demo_df.iloc[:,i_start_ts:].shape[1])
+
+def test_mssa_lrr():
+
+            i_start_ts = 2
+            decomposition = "svd"
+
+            window_s = 7
+            demo_df = load_energy_consumption_df()
+
+            mssa_test = MSSA()
+            mssa_test.fit(df=demo_df,decomposition=decomposition, idx_start_ts=i_start_ts, L=window_s)
+            ts_array = mssa_test.reconstruct_ts(idx_chosen_components=[0],return_as_df=False)
+            M_ = 2
+            fact_lrr = np.array([0.1756382, 0.1845678, 0.1981957, 0.2129944, 0.2277441, 0.2433179])
+
+            assert np.allclose(mssa_test.estimate_LRR( idx_components=[0]) , fact_lrr)==True
+
+
+
+def test_mssa_L_Forecast():
+            ts_index = 0
+            i_start_ts = 2
+            decomposition = "svd"
+
+            window_s = 7
+            demo_df = load_energy_consumption_df()
+
+            mssa_test = MSSA()
+            mssa_test.fit(df=demo_df,decomposition=decomposition, idx_start_ts=i_start_ts,L=window_s)
+            ts_array = mssa_test.reconstruct_ts(idx_chosen_components=[0],return_as_df=False)
+            M_ = 2
+
+            fact_lrr = np.array([0.1756382, 0.1845678, 0.1981957, 0.2129944, 0.2277441, 0.2433179])
+            fact_L_forecast = np.array([
+                                        [15.38531, 16.306364],
+                                        [18.77296, 19.915842],
+                                        [19.89430, 21.064436],
+                                        [83.41526, 89.702892],
+                                        [ 5.63838, 5.986697]
+                                        ])
+
+            assert np.allclose(mssa_test.estimate_LRR( idx_components=[0]) , fact_lrr )==True
+            assert np.allclose( 
+                                mssa_test.L_Forecast( ts_array, M=M_, idx_components=[0])[:,-M_:], fact_L_forecast ) == True
+
+def test_mssa_estimate_ESPRIT():
+
+            i_start_ts = 2
+            decomposition = "svd"
+
+            window_s = 7
+            demo_df = load_energy_consumption_df()
+
+            mssa_test = MSSA()
+            mssa_test.fit(df=demo_df,decomposition=decomposition, idx_start_ts=i_start_ts, L=window_s)
+            ts_array = mssa_test.reconstruct_ts(idx_chosen_components=[0],return_as_df=False)
+            M_ = 2
+            fact_esprit =np.array([ 
+                                   1.5373144+0.j,  1.0207511+0.j,  
+                                   0.91938+0.8687121j,  0.91938-0.8687121j,
+                                   -0.7997768+0.84053236j, -0.7997768-0.84053236j
+                                ])
+            assert np.allclose(mssa_test.estimate_ESPRIT(idx_components=np.arange(mssa_test.d-1), decompose_rho_omega=False),  fact_esprit) ==True
+
+
+
+