Deprecate np.martix usage (#81)

* Initial cleaning of all the np.matrix usage

* add tests for matrix tools

* update README

* add more tests in matrix tools

* fix typo in README

README.md

@@ -8,6 +8,7 @@ Welcome to [pydlm](https://pydlm.github.io/), a flexible time series modeling li
 Updates
 -------------------------------------------
 * Updates in the current Github version:
+    * Deprecated the usage of `np.matrix` to suppress the deprecation warning.
     * Fixed coveralls so all PR merges will report coverage change.
     * Updated the documentations in [pydlm.github.io](https://pydlm.github.io/) using `sphinx`. Exposed more complete APIs in the class reference.
     * Simplified the implementation of `longSeason` component and turns `longSeason` and `autoReg` to be stateless.

pydlm/access/_dlmGet.py

@@ -39,7 +39,7 @@ def _getLatentState(self, name, filterType, start, end):
         """
         end += 1
         indx = self.builder.componentIndex[name]
-        patten = lambda x: x if x is None else x[indx[0]:(indx[1] + 1), 0]
+        patten = lambda x: x if x is None else x[indx[0]:(indx[1] + 1), 0:1]
 
         if filterType == 'forwardFilter':
             return list(map(patten, self.result.filteredState[start:end]))

pydlm/access/dlmAccessMod.py

@@ -229,27 +229,21 @@ def getLatentState(self, filterType='forwardFilter', name='all'):
         start, end = self._checkAndGetWorkingDates(filterType=filterType)
         end += 1
         # to return the full latent states
+        to_array = lambda x: None if x is None else x.flatten().tolist()
         if name == 'all':
             if filterType == 'forwardFilter':
-                return list(map(lambda x: x if x is None
-                                else self._1DmatrixToArray(x),
-                                self.result.filteredState[start:end]))
+                return list(map(to_array, self.result.filteredState[start:end]))
             elif filterType == 'backwardSmoother':
-                return list(map(lambda x: x if x is None
-                                else self._1DmatrixToArray(x),
-                                self.result.smoothedState[start:end]))
+                return list(map(to_array, self.result.smoothedState[start:end]))
             elif filterType == 'predict':
-                return list(map(lambda x: x if x is None
-                                else self._1DmatrixToArray(x),
-                                self.result.predictedState[start:end]))
+                return list(map(to_array, self.result.predictedState[start:end]))
             else:
                 raise NameError('Incorrect filter type.')
 
         # to return the latent state for a given component
         self._checkComponent(name)
-        return list(map(lambda x: x if x is None else self._1DmatrixToArray(x),
-                        self._getLatentState(name=name, filterType=filterType,
-                                             start=start, end=(end - 1))))
+        return list(map(to_array, self._getLatentState(name=name, filterType=filterType,
+                                                       start=start, end=(end - 1))))
 
 
     def getLatentCov(self, filterType='forwardFilter', name='all'):

pydlm/base/kalmanFilter.py

@@ -51,7 +51,7 @@ def __init__(self, discount=[0.99], \
         """
 
         self.__checkDiscount__(discount)
-        self.discount = np.matrix(np.diag(1 / np.sqrt(np.array(discount))))
+        self.discount = np.diag(1 / np.sqrt(np.array(discount)))
         self.updateInnovation = updateInnovation
         self.index = index
 
@@ -74,6 +74,7 @@ def predict(self, model, dealWithMissingEvaluation = False):
 
         # if the step number == 0, we use result from the model state
         if model.prediction.step == 0:
+
             model.prediction.state = np.dot(model.transition, model.state)
             model.prediction.obs = np.dot(model.evaluation, model.prediction.state)
             model.prediction.sysVar = np.dot(np.dot(model.transition, model.sysVar),
@@ -88,20 +89,20 @@ def predict(self, model, dealWithMissingEvaluation = False):
             # add the innovation to the system variance
             model.prediction.sysVar += model.innovation
 
-            model.prediction.obsVar = np.dot(np.dot(model.evaluation, \
-                                                    model.prediction.sysVar), \
+            model.prediction.obsVar = np.dot(np.dot(model.evaluation,
+                                                    model.prediction.sysVar),
                                              model.evaluation.T) + model.noiseVar
             model.prediction.step = 1
 
         # otherwise, we use previous result to predict next time stamp
         else:
             model.prediction.state = np.dot(model.transition, model.prediction.state)
             model.prediction.obs = np.dot(model.evaluation, model.prediction.state)
-            model.prediction.sysVar = np.dot(np.dot(model.transition, \
-                                                    model.prediction.sysVar),\
+            model.prediction.sysVar = np.dot(np.dot(model.transition,
+                                                    model.prediction.sysVar),
                                              model.transition.T)
-            model.prediction.obsVar = np.dot(np.dot(model.evaluation, \
-                                                    model.prediction.sysVar), \
+            model.prediction.obsVar = np.dot(np.dot(model.evaluation,
+                                                    model.prediction.sysVar),
                                              model.evaluation.T) + model.noiseVar
             model.prediction.step += 1
 
@@ -139,25 +140,26 @@ def forwardFilter(self, model, y, dealWithMissingEvaluation = False):
 
             # the prediction error and the correction matrix
             err = y - model.prediction.obs
-            correction = np.dot(model.prediction.sysVar, model.evaluation.T) \
-                         / model.prediction.obsVar
+            correction = (np.dot(model.prediction.sysVar, model.evaluation.T)
+                          / model.prediction.obsVar)
 
             # update new states
             model.df += 1
             lastNoiseVar = model.noiseVar # for updating model.sysVar
-            model.noiseVar = model.noiseVar * \
-                             (1.0 - 1.0 / model.df + \
-                              err * err / model.df / model.prediction.obsVar)
+            model.noiseVar = (model.noiseVar *
+                              (1.0 - 1.0 / model.df +
+                               err * err / model.df / model.prediction.obsVar))
 
             model.state = model.prediction.state + correction * err
 
-            model.sysVar = model.noiseVar[0, 0] / lastNoiseVar[0, 0] * \
-                           (model.prediction.sysVar - np.dot(correction, correction.T) * \
-                            model.prediction.obsVar[0, 0])
+            model.sysVar = (model.noiseVar[0, 0] / lastNoiseVar[0, 0] *
+                            (model.prediction.sysVar - np.dot(correction, correction.T) *
+                             model.prediction.obsVar[0, 0]))
 
             model.obs = np.dot(model.evaluation, model.state)
-            model.obsVar = np.dot(np.dot(model.evaluation, model.sysVar), \
+            model.obsVar = np.dot(np.dot(model.evaluation, model.sysVar),
                                   model.evaluation.T) + model.noiseVar
+
             # update the innovation using discount
             # model.innovation = model.sysVar * (1 / self.discount - 1)
 
@@ -226,11 +228,11 @@ def backwardSmoother(self, model, rawState, rawSysVar):
 
         backward = np.dot(np.dot(rawSysVar, model.transition.T), predSysVarInv)
         model.state = rawState + np.dot(backward, (model.state - model.prediction.state))
-        model.sysVar = rawSysVar + \
-                       np.dot(np.dot(backward, \
-                                     (model.sysVar - model.prediction.sysVar)), backward.T)
+        model.sysVar = (rawSysVar + 
+                        np.dot(np.dot(backward,
+                                      (model.sysVar - model.prediction.sysVar)), backward.T))
         model.obs = np.dot(model.evaluation, model.state)
-        model.obsVar = np.dot(np.dot(model.evaluation, model.sysVar), \
+        model.obsVar = np.dot(np.dot(model.evaluation, model.sysVar),
                               model.evaluation.T) + model.noiseVar
 
         # recover the evaluation and the transition matrix
@@ -266,16 +268,16 @@ def backwardSampler(self, model, rawState, rawSysVar):
 
         backward = np.dot(np.dot(rawSysVar, model.transition.T), predSysVarInv)
         model.state = rawState + np.dot(backward, (model.state - model.prediction.state))
-        model.sysVar = rawSysVar + \
-                       np.dot(np.dot(backward, \
-                                     (model.sysVar - model.prediction.sysVar)), backward.T)
-        model.state = np.matrix(np.random.multivariate_normal(model.state.A1, \
-                                                              model.sysVar)).T
+        model.sysVar = (rawSysVar + 
+                        np.dot(np.dot(backward,
+                                      (model.sysVar - model.prediction.sysVar)), backward.T))
+        model.state = np.random.multivariate_normal(model.state.A1,
+                                                    model.sysVar).T
         model.obs = np.dot(model.evaluation, model.state)
-        model.obsVar = np.dot(np.dot(model.evaluation, model.sysVar), \
+        model.obsVar = np.dot(np.dot(model.evaluation, model.sysVar),
                               model.evaluation.T) + model.noiseVar
-        model.obs = np.matrix(np.random.multivariate_normal(model.obs.A1, \
-                                                              model.obsVar)).T
+        model.obs = np.random.multivariate_normal(model.obs.A1,
+                                                  model.obsVar).T
 
 
     # for updating the discounting factor
@@ -287,7 +289,7 @@ def updateDiscount(self, newDiscount):
         """
 
         self.__checkDiscount__(newDiscount)
-        self.discount = np.matrix(np.diag(1 / np.sqrt(newDiscount)))
+        self.discount = np.diag(1 / np.sqrt(newDiscount))
 
 
     def __checkDiscount__(self, discount):
@@ -306,8 +308,8 @@ def __updateInnovation__(self, model):
 
         """
 
-        model.innovation = np.dot(np.dot(self.discount, model.prediction.sysVar), \
-                                      self.discount) - model.prediction.sysVar
+        model.innovation = np.dot(np.dot(self.discount, model.prediction.sysVar),
+                                  self.discount) - model.prediction.sysVar
 
 
     # update the innovation
@@ -318,13 +320,13 @@ def __updateInnovation2__(self, model):
 
         """
 
-        innovation = np.dot(np.dot(self.discount, model.prediction.sysVar), \
-                                      self.discount) - model.prediction.sysVar
-        model.innovation = np.matrix(np.zeros(innovation.shape))
+        innovation = np.dot(np.dot(self.discount, model.prediction.sysVar),
+                            self.discount) - model.prediction.sysVar
+        model.innovation = np.zeros(innovation.shape)
         for name in self.index:
             indx = self.index[name]
-            model.innovation[indx[0]: (indx[1] + 1), indx[0]: (indx[1] + 1)] = \
-                    innovation[indx[0]: (indx[1] + 1), indx[0]: (indx[1] + 1)]
+            model.innovation[indx[0]: (indx[1] + 1), indx[0]: (indx[1] + 1)] = (
+                innovation[indx[0]: (indx[1] + 1), indx[0]: (indx[1] + 1)])
 
 
     # a generalized inverse of matrix A

pydlm/core/_dlm.py

@@ -485,7 +485,8 @@ def _1DmatrixToArray(self, arrayOf1dMatrix):
         """ Change an array of 1 x 1 matrix to normal array.
 
         """
-        return [item.tolist()[0][0] for item in arrayOf1dMatrix]
+        to_array = lambda x : None if x is None else x.tolist()[0][0]
+        return [to_array(item) for item in arrayOf1dMatrix]
 
 
     # function to judge whether a component is in the model

pydlm/dlm.py

@@ -24,7 +24,7 @@
 >>> myDlm.fitForwardFilter()
 
 >>> # extract the filtered result
->>> myDlm.getFilteredObs()
+>>> myDlm.getMean()
 
 """
 # This is the major class for fitting time series data using the

pydlm/modeler/autoReg.py

@@ -95,7 +95,7 @@ def createEvaluation(self, step, data):
         if step > len(data):
             raise NameError("There is no sufficient data for creating autoregressor.")
         # We pad numbers if the step is too early
-        self.evaluation = np.matrix([[self.padding] * (self.d - step) +
+        self.evaluation = np.array([[self.padding] * (self.d - step) +
                                     list(data[max(0, (step - self.d)) : step])])
 
 

pydlm/modeler/builder.py

@@ -293,10 +293,10 @@ def initialize(self, data=[], noise=1):
 
         self.statePrior = state
         self.sysVarPrior = sysVar
-        self.noiseVar = np.matrix(noise)
+        self.noiseVar = np.array([[noise]])
         self.model = baseModel(transition=transition,
                                evaluation=evaluation,
-                               noiseVar=np.matrix(noise),
+                               noiseVar=np.array([[noise]]),
                                sysVar=sysVar,
                                state=state,
                                df=self.initialDegreeFreedom)

pydlm/modeler/dynamic.py

@@ -79,7 +79,7 @@ def __init__(self,
                             "in the features.")
 
         self.features = deepcopy(features)
-        if isinstance(features, np.matrix):
+        if isinstance(features, np.matrix) or isinstance(features, np.ndarray):
             self.features = self.features.tolist()
         self.componentType = 'dynamic'
         self.name = name
@@ -104,7 +104,7 @@ def createEvaluation(self, step):
         given date
 
         """
-        self.evaluation = np.matrix([self.features[step]])
+        self.evaluation = np.array([self.features[step]])
 
 
     def createTransition(self):

pydlm/modeler/longSeason.py

@@ -81,7 +81,7 @@ def __init__(self,
         self.componentType = 'longSeason'
 
         # Initialize the evaluation vector
-        self.evaluation = np.matrix([0] * self.period)
+        self.evaluation = np.array([[0] * self.period])
 
 
     def updateEvaluation(self, step, data=None):

pydlm/modeler/matrixTools.py

@@ -6,15 +6,15 @@ class matrixTools:
     @staticmethod
     def matrixAddInDiag(A, B):
         if A is None:
-            return np.matrix(B)
+            return B
         elif B is None:
-            return np.matrix(A)
+            return A
         else:
             (An, Ap) = A.shape
             (Bn, Bp) = B.shape
 
-            newMatrixA = np.concatenate((A, np.matrix(np.zeros((An, Bp)))), 1)
-            newMatrixB = np.concatenate((np.matrix(np.zeros((Bn, Ap))), B), 1)
+            newMatrixA = np.concatenate((A, np.zeros((An, Bp))), 1)
+            newMatrixB = np.concatenate((np.zeros((Bn, Ap)), B), 1)
             return np.concatenate((newMatrixA, newMatrixB), 0)
 
 
@@ -33,18 +33,18 @@ def matrixAddByRow(A, B):
     @staticmethod
     def matrixAddByCol(A, B):
         if A is None:
-            return np.matrix(B)
+            return B
         elif B is None:
-            return np.matrix(A)
+            return A
         else:
             return np.concatenate((A, B), 1)
 
 
     @staticmethod
     def AddTwoVectors(a, b):
         if a is None:
-            return np.array(b)
+            return b
         elif b is None:
-            return np.array(a)
+            return a
         else:
             return np.concatenate((a, b))

pydlm/modeler/seasonality.py

@@ -94,7 +94,7 @@ def createEvaluation(self):
         """ Create the evaluation matrix
 
         """
-        self.evaluation = np.matrix(np.zeros((1, self.d)))
+        self.evaluation = np.zeros((1, self.d))
         self.evaluation[0, 0] = 1
 
 
@@ -116,22 +116,22 @@ def createTransition(self):
         [1 0 0 0]]
 
         """
-        self.transition = np.matrix(np.diag(np.ones(self.d - 1), 1))
+        self.transition = np.diag(np.ones(self.d - 1), 1)
         self.transition[self.d - 1, 0] = 1
 
 
     def createCovPrior(self, cov = 1e7):
         """Create the prior covariance matrix for the latent states.
 
         """
-        self.covPrior = np.matrix(np.eye(self.d)) * cov
+        self.covPrior = np.eye(self.d) * cov
 
 
     def createMeanPrior(self, mean = 0):
         """ Create the prior latent state
 
         """
-        self.meanPrior = np.matrix(np.ones((self.d, 1))) * mean
+        self.meanPrior = np.ones((self.d, 1)) * mean
 
 
     # Form free seasonality component ensures that sum(mean) = 0
@@ -147,9 +147,9 @@ def freeForm(self):
         if self.covPrior is None or self.meanPrior is None:
             raise NameError('freeForm can only be called after prior created.')
         else:
-            u = np.sum(np.sum(self.covPrior, 0), 1)[0, 0]
-            A = np.sum(self.covPrior, 1) / u
-            self.meanPrior = self.meanPrior - A * np.sum(self.meanPrior, 0)[0, 0]
+            u = np.sum(self.covPrior)
+            A = np.sum(self.covPrior, axis=1, keepdims=True) / u
+            self.meanPrior = self.meanPrior - A * np.sum(self.meanPrior)
             self.covPrior = self.covPrior - np.dot(A, A.T) * u